1 // Copyright (c) 2017-2023, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Structs and helper functions to evaluate data-driven subgrid-stress modeling 10 /// See 'Invariant data-driven subgrid stress modeling in the strain-rate eigenframe for large eddy simulation' 2022 and 'S-frame discrepancy 11 /// correction models for data-informed Reynolds stress closure' 2022 12 13 #ifndef sgs_dd_model_h 14 #define sgs_dd_model_h 15 16 #include <ceed.h> 17 18 #include "newtonian_state.h" 19 #include "newtonian_types.h" 20 #include "sgs_dd_utils.h" 21 #include "utils.h" 22 #include "utils_eigensolver_jacobi.h" 23 24 typedef struct SGS_DD_ModelContext_ *SGS_DDModelContext; 25 struct SGS_DD_ModelContext_ { 26 CeedInt num_inputs, num_outputs; 27 CeedInt num_layers; 28 CeedInt num_neurons; 29 CeedScalar alpha; 30 31 struct NewtonianIdealGasContext_ gas; 32 struct { 33 size_t bias1, bias2; 34 size_t weight1, weight2; 35 size_t out_scaling; 36 } offsets; 37 size_t total_bytes; 38 CeedScalar data[1]; 39 }; 40 41 // @brief Denormalize outputs using min-max (de-)normalization 42 CEED_QFUNCTION_HELPER void DenormalizeDDOutputs(CeedScalar output[6], const CeedScalar (*new_bounds)[2], const CeedScalar old_bounds[6][2]) { 43 CeedScalar bounds_ratio; 44 for (int i = 0; i < 6; i++) { 45 bounds_ratio = (new_bounds[i][1] - new_bounds[i][0]) / (old_bounds[i][1] - old_bounds[i][0]); 46 output[i] = bounds_ratio * (output[i] - old_bounds[i][1]) + new_bounds[i][1]; 47 } 48 } 49 50 CEED_QFUNCTION_HELPER void LeakyReLU(CeedScalar *x, const CeedScalar alpha, const CeedInt N) { 51 for (CeedInt i = 0; i < N; i++) x[i] *= (x[i] < 0 ? alpha : 1.); 52 } 53 54 CEED_QFUNCTION_HELPER void DataDrivenInference(const CeedScalar *inputs, CeedScalar *outputs, SGS_DDModelContext sgsdd_ctx) { 55 const CeedInt num_neurons = sgsdd_ctx->num_neurons; 56 const CeedInt num_inputs = sgsdd_ctx->num_inputs; 57 const CeedInt num_outputs = sgsdd_ctx->num_outputs; 58 const CeedScalar alpha = sgsdd_ctx->alpha; 59 const CeedScalar *bias1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias1]; 60 const CeedScalar *bias2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias2]; 61 const CeedScalar *weight1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight1]; 62 const CeedScalar *weight2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight2]; 63 CeedScalar V[20] = {0.}; 64 65 CopyN(bias1, V, num_neurons); 66 MatVecNM(weight1, inputs, num_neurons, num_inputs, CEED_NOTRANSPOSE, V); 67 LeakyReLU(V, alpha, num_neurons); 68 CopyN(bias2, outputs, num_outputs); 69 MatVecNM(weight2, V, num_outputs, num_neurons, CEED_NOTRANSPOSE, outputs); 70 } 71 72 CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropic(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta, 73 const CeedScalar viscosity, CeedScalar kmsgs_stress[6], SGS_DDModelContext sgsdd_ctx) { 74 CeedScalar inputs[6], grad_velo_magnitude, eigenvectors[3][3], sgs_sframe_sym[6] = {0.}; 75 76 ComputeSGS_DDAnisotropicInputs(grad_velo_aniso, km_A_ij, delta, viscosity, eigenvectors, inputs, &grad_velo_magnitude); 77 78 DataDrivenInference(inputs, sgs_sframe_sym, sgsdd_ctx); 79 80 CeedScalar old_bounds[6][2] = {{0}}; 81 for (int j = 0; j < 6; j++) old_bounds[j][1] = 1; 82 const CeedScalar(*new_bounds)[2] = (const CeedScalar(*)[2]) & sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling]; 83 DenormalizeDDOutputs(sgs_sframe_sym, new_bounds, old_bounds); 84 85 // Re-dimensionalize sgs_stress 86 ScaleN(sgs_sframe_sym, Square(delta) * Square(grad_velo_magnitude), 6); 87 88 CeedScalar sgs_stress[3][3] = {{0.}}; 89 { // Rotate SGS Stress back to physical frame, SGS_physical = E^T SGS_sframe E 90 CeedScalar Evec_sgs[3][3] = {{0.}}; 91 const CeedScalar sgs_sframe[3][3] = { 92 {sgs_sframe_sym[0], sgs_sframe_sym[3], sgs_sframe_sym[4]}, 93 {sgs_sframe_sym[3], sgs_sframe_sym[1], sgs_sframe_sym[5]}, 94 {sgs_sframe_sym[4], sgs_sframe_sym[5], sgs_sframe_sym[2]}, 95 }; 96 MatMat3(eigenvectors, sgs_sframe, CEED_TRANSPOSE, CEED_NOTRANSPOSE, Evec_sgs); 97 MatMat3(Evec_sgs, eigenvectors, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, sgs_stress); 98 } 99 100 KMPack(sgs_stress, kmsgs_stress); 101 } 102 103 // @brief Calculate subgrid stress at nodes using anisotropic data-driven model 104 CEED_QFUNCTION_HELPER int ComputeSGS_DDAnisotropicNodal(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 105 StateVariable state_var) { 106 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 107 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 108 const CeedScalar(*grad_velo)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[2]; 109 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 110 const CeedScalar(*inv_multiplicity) = (const CeedScalar(*))in[4]; 111 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 112 113 const SGS_DDModelContext sgsdd_ctx = (SGS_DDModelContext)ctx; 114 const NewtonianIdealGasContext gas = &sgsdd_ctx->gas; 115 116 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 117 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 118 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 119 const CeedScalar grad_velo_aniso[3][3] = { 120 {grad_velo[0][0][i], grad_velo[0][1][i], grad_velo[0][2][i]}, 121 {grad_velo[1][0][i], grad_velo[1][1][i], grad_velo[1][2][i]}, 122 {grad_velo[2][0][i], grad_velo[2][1][i], grad_velo[2][2][i]} 123 }; 124 const CeedScalar km_A_ij[6] = {A_ij_delta[0][i], A_ij_delta[1][i], A_ij_delta[2][i], A_ij_delta[3][i], A_ij_delta[4][i], A_ij_delta[5][i]}; 125 const CeedScalar delta = A_ij_delta[6][i]; 126 const State s = StateFromQ(gas, qi, x_i, state_var); 127 CeedScalar km_sgs[6]; 128 129 ComputeSGS_DDAnisotropic(grad_velo_aniso, km_A_ij, delta, gas->mu / s.U.density, km_sgs, sgsdd_ctx); 130 131 for (int j = 0; j < 6; j++) v[j][i] = inv_multiplicity[i] * km_sgs[j]; 132 } 133 return 0; 134 } 135 136 CEED_QFUNCTION(ComputeSGS_DDAnisotropicNodal_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 137 return ComputeSGS_DDAnisotropicNodal(ctx, Q, in, out, STATEVAR_PRIMITIVE); 138 } 139 140 CEED_QFUNCTION(ComputeSGS_DDAnisotropicNodal_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 141 return ComputeSGS_DDAnisotropicNodal(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 142 } 143 144 // @brief Adds subgrid stress to residual (during IFunction evaluation) 145 CEED_QFUNCTION_HELPER int FluxSubgridStress(const StatePrimitive Y, const CeedScalar km_sgs[6], CeedScalar Flux[5][3]) { 146 CeedScalar sgs[3][3]; 147 148 KMUnpack(km_sgs, sgs); 149 for (CeedInt j = 0; j < 3; j++) { 150 Flux[0][j] = 0.; 151 for (CeedInt k = 0; k < 3; k++) Flux[k + 1][j] = sgs[k][j]; 152 Flux[4][j] = Y.velocity[0] * sgs[0][j] + Y.velocity[1] * sgs[1][j] + Y.velocity[2] * sgs[2][j]; 153 } 154 return 0; 155 } 156 157 CEED_QFUNCTION_HELPER int IFunction_NodalSubgridStress(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 158 StateVariable state_var) { 159 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 160 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 161 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 162 const CeedScalar(*km_sgs)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 163 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[0]; 164 165 SGS_DDModelContext sgsdd_ctx = (SGS_DDModelContext)ctx; 166 NewtonianIdealGasContext gas = &sgsdd_ctx->gas; 167 168 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 169 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 170 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 171 const State s = StateFromQ(gas, qi, x_i, state_var); 172 173 const CeedScalar wdetJ = q_data[0][i]; 174 const CeedScalar dXdx[3][3] = { 175 {q_data[1][i], q_data[2][i], q_data[3][i]}, 176 {q_data[4][i], q_data[5][i], q_data[6][i]}, 177 {q_data[7][i], q_data[8][i], q_data[9][i]} 178 }; 179 180 CeedScalar Flux[5][3]; 181 const CeedScalar km_sgs_i[6] = {km_sgs[0][i], km_sgs[1][i], km_sgs[2][i], km_sgs[3][i], km_sgs[4][i], km_sgs[5][i]}; 182 FluxSubgridStress(s.Y, km_sgs_i, Flux); 183 184 for (CeedInt k = 0; k < 3; k++) { 185 for (CeedInt j = 0; j < 5; j++) { 186 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 187 } 188 } 189 } 190 return 0; 191 } 192 193 CEED_QFUNCTION(IFunction_NodalSubgridStress_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 194 return IFunction_NodalSubgridStress(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 195 } 196 197 CEED_QFUNCTION(IFunction_NodalSubgridStress_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 198 return IFunction_NodalSubgridStress(ctx, Q, in, out, STATEVAR_PRIMITIVE); 199 } 200 201 #endif // sgs_dd_model_h 202