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 for 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 "utils.h" 21 #include "utils_eigensolver_jacobi.h" 22 23 typedef struct SGS_DD_ModelContext_ *SGS_DDModelContext; 24 struct SGS_DD_ModelContext_ { 25 CeedInt num_inputs, num_outputs; 26 CeedInt num_layers; 27 CeedInt num_neurons; 28 CeedScalar alpha; 29 30 struct NewtonianIdealGasContext_ gas; 31 struct { 32 size_t bias1, bias2; 33 size_t weight1, weight2; 34 size_t out_scaling; 35 } offsets; 36 size_t total_bytes; 37 CeedScalar data[1]; 38 }; 39 40 // @brief Calculate the inverse of the multiplicity, reducing to a single component 41 CEED_QFUNCTION(InverseMultiplicity)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 42 const CeedScalar(*multiplicity)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 43 CeedScalar(*inv_multiplicity) = (CeedScalar(*))out[0]; 44 45 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) inv_multiplicity[i] = 1.0 / multiplicity[0][i]; 46 return 0; 47 } 48 49 // @brief Calculate Frobenius norm of velocity gradient from eigenframe quantities 50 CEED_QFUNCTION_HELPER CeedScalar VelocityGradientMagnitude(const CeedScalar strain_sframe[3], const CeedScalar vorticity_sframe[3]) { 51 return sqrt(Dot3(strain_sframe, strain_sframe) + 0.5 * Dot3(vorticity_sframe, vorticity_sframe)); 52 }; 53 54 // @brief Denormalize outputs using min-max (de-)normalization 55 CEED_QFUNCTION_HELPER void DenormalizeDDOutputs(CeedScalar output[6], const CeedScalar (*new_bounds)[2], const CeedScalar old_bounds[6][2]) { 56 CeedScalar bounds_ratio; 57 for (int i = 0; i < 6; i++) { 58 bounds_ratio = (new_bounds[i][1] - new_bounds[i][0]) / (old_bounds[i][1] - old_bounds[i][0]); 59 output[i] = bounds_ratio * (output[i] - old_bounds[i][1]) + new_bounds[i][1]; 60 } 61 } 62 63 // @brief Change the order of basis vectors so that they align with vector and obey right-hand rule 64 // @details The e_1 and e_3 basis vectors are the closest aligned to the vector. The e_2 is set via e_3 x e_1 65 // The basis vectors are assumed to form the rows of the basis matrix. 66 CEED_QFUNCTION_HELPER void OrientBasisWithVector(CeedScalar basis[3][3], const CeedScalar vector[3]) { 67 CeedScalar alignment[3] = {0.}, cross[3]; 68 69 MatVec3(basis, vector, CEED_NOTRANSPOSE, alignment); 70 71 if (alignment[0] < 0) ScaleN(basis[0], -1, 3); 72 if (alignment[2] < 0) ScaleN(basis[2], -1, 3); 73 74 Cross3(basis[2], basis[0], cross); 75 CeedScalar basis_1_orientation = Dot3(cross, basis[1]); 76 if (basis_1_orientation < 0) ScaleN(basis[1], -1, 3); 77 } 78 79 CEED_QFUNCTION_HELPER void LeakyReLU(CeedScalar *x, const CeedScalar alpha, const CeedInt N) { 80 for (CeedInt i = 0; i < N; i++) x[i] *= (x[i] < 0 ? alpha : 1.); 81 } 82 83 CEED_QFUNCTION_HELPER void DataDrivenInference(const CeedScalar *inputs, CeedScalar *outputs, SGS_DDModelContext sgsdd_ctx) { 84 const CeedInt num_neurons = sgsdd_ctx->num_neurons; 85 const CeedInt num_inputs = sgsdd_ctx->num_inputs; 86 const CeedInt num_outputs = sgsdd_ctx->num_outputs; 87 const CeedScalar alpha = sgsdd_ctx->alpha; 88 const CeedScalar *bias1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias1]; 89 const CeedScalar *bias2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias2]; 90 const CeedScalar *weight1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight1]; 91 const CeedScalar *weight2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight2]; 92 CeedScalar V[20] = {0.}; 93 94 CopyN(bias1, V, num_neurons); 95 MatVecNM(weight1, inputs, num_neurons, num_inputs, CEED_NOTRANSPOSE, V); 96 LeakyReLU(V, alpha, num_neurons); 97 CopyN(bias2, outputs, num_outputs); 98 MatVecNM(weight2, V, num_outputs, num_neurons, CEED_NOTRANSPOSE, outputs); 99 } 100 101 /** 102 * @brief Compute model inputs for anisotropic data-driven model 103 * 104 * @param[in] grad_velo_aniso Gradient of velocity in physical (anisotropic) coordinates 105 * @param[in] km_A_ij Anisotropy tensor, in Kelvin-Mandel notation 106 * @param[in] delta Length used to create anisotropy tensor 107 * @param[in] viscosity Kinematic viscosity 108 * @param[out] eigenvectors Eigenvectors of the (anisotropic) velocity gradient 109 * @param[out] inputs Data-driven model inputs 110 * @param[out] grad_velo_magnitude Frobenius norm of the velocity gradient 111 */ 112 CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropicInputs(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta, 113 const CeedScalar viscosity, CeedScalar eigenvectors[3][3], CeedScalar inputs[6], CeedScalar *grad_velo_magnitude) { 114 CeedScalar strain_sframe[3] = {0.}, vorticity_sframe[3] = {0.}; 115 CeedScalar A_ij[3][3] = {{0.}}, grad_velo_iso[3][3] = {{0.}}; 116 117 // -- Transform physical, anisotropic velocity gradient to isotropic 118 KMUnpack(km_A_ij, A_ij); 119 MatMat3(grad_velo_aniso, A_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, grad_velo_iso); 120 121 { // -- Get Eigenframe 122 CeedScalar kmstrain_iso[6], strain_iso[3][3]; 123 CeedInt work_vector[3] = {0}; 124 KMStrainRate(grad_velo_iso, kmstrain_iso); 125 KMUnpack(kmstrain_iso, strain_iso); 126 Diagonalize3(strain_iso, strain_sframe, eigenvectors, work_vector, SORT_DECREASING_EVALS, true, 5); 127 } 128 129 { // -- Get vorticity in S-frame 130 CeedScalar rotation_iso[3][3]; 131 RotationRate(grad_velo_iso, rotation_iso); 132 CeedScalar vorticity_iso[3] = {-2 * rotation_iso[1][2], 2 * rotation_iso[0][2], -2 * rotation_iso[0][1]}; 133 OrientBasisWithVector(eigenvectors, vorticity_iso); 134 MatVec3(eigenvectors, vorticity_iso, CEED_NOTRANSPOSE, vorticity_sframe); 135 } 136 137 // -- Calculate DD model inputs 138 *grad_velo_magnitude = VelocityGradientMagnitude(strain_sframe, vorticity_sframe); 139 inputs[0] = strain_sframe[0]; 140 inputs[1] = strain_sframe[1]; 141 inputs[2] = strain_sframe[2]; 142 inputs[3] = vorticity_sframe[0]; 143 inputs[4] = vorticity_sframe[1]; 144 inputs[5] = viscosity / Square(delta); 145 ScaleN(inputs, 1 / (*grad_velo_magnitude + CEED_EPSILON), 6); 146 } 147 148 CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropic(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta, 149 const CeedScalar viscosity, CeedScalar kmsgs_stress[6], SGS_DDModelContext sgsdd_ctx) { 150 CeedScalar inputs[6], grad_velo_magnitude, eigenvectors[3][3], sgs_sframe_sym[6] = {0.}; 151 152 ComputeSGS_DDAnisotropicInputs(grad_velo_aniso, km_A_ij, delta, viscosity, eigenvectors, inputs, &grad_velo_magnitude); 153 154 DataDrivenInference(inputs, sgs_sframe_sym, sgsdd_ctx); 155 156 CeedScalar old_bounds[6][2] = {{0}}; 157 for (int j = 0; j < 6; j++) old_bounds[j][1] = 1; 158 const CeedScalar(*new_bounds)[2] = (const CeedScalar(*)[2]) & sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling]; 159 DenormalizeDDOutputs(sgs_sframe_sym, new_bounds, old_bounds); 160 161 // Re-dimensionalize sgs_stress 162 ScaleN(sgs_sframe_sym, Square(delta) * Square(grad_velo_magnitude), 6); 163 164 CeedScalar sgs_stress[3][3] = {{0.}}; 165 { // Rotate SGS Stress back to physical frame, SGS_physical = E^T SGS_sframe E 166 CeedScalar Evec_sgs[3][3] = {{0.}}; 167 const CeedScalar sgs_sframe[3][3] = { 168 {sgs_sframe_sym[0], sgs_sframe_sym[3], sgs_sframe_sym[4]}, 169 {sgs_sframe_sym[3], sgs_sframe_sym[1], sgs_sframe_sym[5]}, 170 {sgs_sframe_sym[4], sgs_sframe_sym[5], sgs_sframe_sym[2]}, 171 }; 172 MatMat3(eigenvectors, sgs_sframe, CEED_TRANSPOSE, CEED_NOTRANSPOSE, Evec_sgs); 173 MatMat3(Evec_sgs, eigenvectors, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, sgs_stress); 174 } 175 176 KMPack(sgs_stress, kmsgs_stress); 177 } 178 179 // @brief Calculate subgrid stress at nodes using anisotropic data-driven model 180 CEED_QFUNCTION_HELPER int ComputeSGS_DDAnisotropicNodal(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 181 StateVariable state_var) { 182 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 183 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 184 const CeedScalar(*grad_velo)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[2]; 185 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 186 const CeedScalar(*inv_multiplicity) = (const CeedScalar(*))in[4]; 187 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 188 189 const SGS_DDModelContext sgsdd_ctx = (SGS_DDModelContext)ctx; 190 const NewtonianIdealGasContext gas = &sgsdd_ctx->gas; 191 192 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 193 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 194 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 195 const CeedScalar grad_velo_aniso[3][3] = { 196 {grad_velo[0][0][i], grad_velo[0][1][i], grad_velo[0][2][i]}, 197 {grad_velo[1][0][i], grad_velo[1][1][i], grad_velo[1][2][i]}, 198 {grad_velo[2][0][i], grad_velo[2][1][i], grad_velo[2][2][i]} 199 }; 200 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]}; 201 const CeedScalar delta = A_ij_delta[6][i]; 202 const State s = StateFromQ(gas, qi, x_i, state_var); 203 CeedScalar km_sgs[6]; 204 205 ComputeSGS_DDAnisotropic(grad_velo_aniso, km_A_ij, delta, gas->mu / s.U.density, km_sgs, sgsdd_ctx); 206 207 for (int j = 0; j < 6; j++) v[j][i] = inv_multiplicity[i] * km_sgs[j]; 208 } 209 return 0; 210 } 211 212 CEED_QFUNCTION(ComputeSGS_DDAnisotropicNodal_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 213 return ComputeSGS_DDAnisotropicNodal(ctx, Q, in, out, STATEVAR_PRIMITIVE); 214 } 215 216 CEED_QFUNCTION(ComputeSGS_DDAnisotropicNodal_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 217 return ComputeSGS_DDAnisotropicNodal(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 218 } 219 220 // @brief Adds subgrid stress to residual (during IFunction evaluation) 221 CEED_QFUNCTION_HELPER int FluxSubgridStress(const StatePrimitive Y, const CeedScalar km_sgs[6], CeedScalar Flux[5][3]) { 222 CeedScalar sgs[3][3]; 223 224 KMUnpack(km_sgs, sgs); 225 for (CeedInt j = 0; j < 3; j++) { 226 Flux[0][j] = 0.; 227 for (CeedInt k = 0; k < 3; k++) Flux[k + 1][j] = sgs[k][j]; 228 Flux[4][j] = Y.velocity[0] * sgs[0][j] + Y.velocity[1] * sgs[1][j] + Y.velocity[2] * sgs[2][j]; 229 } 230 return 0; 231 } 232 233 CEED_QFUNCTION_HELPER int IFunction_NodalSubgridStress(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 234 StateVariable state_var) { 235 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 236 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 237 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 238 const CeedScalar(*km_sgs)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 239 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[0]; 240 241 SGS_DDModelContext sgsdd_ctx = (SGS_DDModelContext)ctx; 242 NewtonianIdealGasContext gas = &sgsdd_ctx->gas; 243 244 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 245 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 246 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 247 const State s = StateFromQ(gas, qi, x_i, state_var); 248 249 const CeedScalar wdetJ = q_data[0][i]; 250 const CeedScalar dXdx[3][3] = { 251 {q_data[1][i], q_data[2][i], q_data[3][i]}, 252 {q_data[4][i], q_data[5][i], q_data[6][i]}, 253 {q_data[7][i], q_data[8][i], q_data[9][i]} 254 }; 255 256 CeedScalar Flux[5][3]; 257 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]}; 258 FluxSubgridStress(s.Y, km_sgs_i, Flux); 259 260 for (CeedInt k = 0; k < 3; k++) { 261 for (CeedInt j = 0; j < 5; j++) { 262 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 263 } 264 } 265 } 266 return 0; 267 } 268 269 CEED_QFUNCTION(IFunction_NodalSubgridStress_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 270 return IFunction_NodalSubgridStress(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 271 } 272 273 CEED_QFUNCTION(IFunction_NodalSubgridStress_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 274 return IFunction_NodalSubgridStress(ctx, Q, in, out, STATEVAR_PRIMITIVE); 275 } 276 277 #endif // sgs_dd_model_h 278