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 "utils.h" 20 #include "utils_eigensolver_jacobi.h" 21 22 typedef struct SGS_DD_ModelContext_ *SGS_DDModelContext; 23 struct SGS_DD_ModelContext_ { 24 CeedInt num_inputs, num_outputs; 25 CeedInt num_layers; 26 CeedInt num_neurons; 27 CeedScalar alpha; 28 29 struct { 30 size_t bias1, bias2; 31 size_t weight1, weight2; 32 size_t out_scaling; 33 } offsets; 34 size_t total_bytes; 35 CeedScalar data[1]; 36 }; 37 38 // @brief Calculate the inverse of the multiplicity, reducing to a single component 39 CEED_QFUNCTION(InverseMultiplicity)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 40 const CeedScalar(*multiplicity)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 41 CeedScalar(*inv_multiplicity) = (CeedScalar(*))out[0]; 42 43 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) inv_multiplicity[i] = 1.0 / multiplicity[0][i]; 44 return 0; 45 } 46 47 // @brief Calculate Frobenius norm of velocity gradient from eigenframe quantities 48 CEED_QFUNCTION_HELPER CeedScalar VelocityGradientMagnitude(const CeedScalar strain_sframe[3], const CeedScalar vorticity_sframe[3]) { 49 return sqrt(Dot3(strain_sframe, strain_sframe) + 0.5 * Dot3(vorticity_sframe, vorticity_sframe)); 50 }; 51 52 // @brief Denormalize outputs using min-max (de-)normalization 53 CEED_QFUNCTION_HELPER void DenormalizeDDOutputs(CeedScalar output[6], const CeedScalar new_bounds[6][2], const CeedScalar old_bounds[6][2]) { 54 CeedScalar bounds_ratio; 55 for (int i = 0; i < 6; i++) { 56 bounds_ratio = (new_bounds[i][1] - new_bounds[i][0]) / (old_bounds[i][1] - old_bounds[i][0]); 57 output[i] = bounds_ratio * (output[i] - old_bounds[i][1]) + new_bounds[i][1]; 58 } 59 } 60 61 // @brief Change the order of basis vectors so that they align with vector and obey right-hand rule 62 // @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 63 // The basis vectors are assumed to form the rows of the basis matrix. 64 CEED_QFUNCTION_HELPER void OrientBasisWithVector(CeedScalar basis[3][3], const CeedScalar vector[3]) { 65 CeedScalar alignment[3] = {0.}, cross[3]; 66 67 MatVec3(basis, vector, CEED_NOTRANSPOSE, alignment); 68 69 if (alignment[0] < 0) ScaleN(basis[0], -1, 3); 70 if (alignment[2] < 0) ScaleN(basis[2], -1, 3); 71 72 Cross3(basis[2], basis[0], cross); 73 CeedScalar basis_1_orientation = Dot3(cross, basis[1]); 74 if (basis_1_orientation < 0) ScaleN(basis[1], -1, 3); 75 } 76 77 CEED_QFUNCTION_HELPER void LeakyReLU(CeedScalar *x, const CeedScalar alpha, const CeedInt N) { 78 for (CeedInt i = 0; i < N; i++) x[i] *= (x[i] < 0 ? alpha : 1.); 79 } 80 81 CEED_QFUNCTION_HELPER void DataDrivenInference(const CeedScalar *inputs, CeedScalar *outputs, SGS_DDModelContext sgsdd_ctx) { 82 const CeedInt num_neurons = sgsdd_ctx->num_neurons; 83 const CeedInt num_inputs = sgsdd_ctx->num_inputs; 84 const CeedInt num_outputs = sgsdd_ctx->num_outputs; 85 const CeedScalar alpha = sgsdd_ctx->alpha; 86 const CeedScalar *bias1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias1]; 87 const CeedScalar *bias2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias2]; 88 const CeedScalar *weight1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight1]; 89 const CeedScalar *weight2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight2]; 90 CeedScalar V[20] = {0.}; 91 92 CopyN(bias1, V, num_neurons); 93 MatVecNM(weight1, inputs, num_neurons, num_inputs, CEED_NOTRANSPOSE, V); 94 LeakyReLU(V, alpha, num_neurons); 95 CopyN(bias2, outputs, num_outputs); 96 MatVecNM(weight2, V, num_outputs, num_neurons, CEED_NOTRANSPOSE, outputs); 97 } 98 99 CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropic(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta, 100 const CeedScalar viscosity, CeedScalar kmsgs_stress[6], SGS_DDModelContext sgsdd_ctx) { 101 CeedScalar strain_sframe[3] = {0.}, vorticity_sframe[3] = {0.}, eigenvectors[3][3]; 102 CeedScalar A_ij[3][3] = {{0.}}, grad_velo_iso[3][3] = {{0.}}; 103 104 // -- Unpack anisotropy tensor 105 KMUnpack(km_A_ij, A_ij); 106 107 // -- Transform physical, anisotropic velocity gradient to isotropic 108 MatMat3(grad_velo_aniso, A_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, grad_velo_iso); 109 110 { // -- Get Eigenframe 111 CeedScalar kmstrain_iso[6], strain_iso[3][3]; 112 CeedInt work_vector[3] = {0}; 113 KMStrainRate(grad_velo_iso, kmstrain_iso); 114 KMUnpack(kmstrain_iso, strain_iso); 115 Diagonalize3(strain_iso, strain_sframe, eigenvectors, work_vector, SORT_DECREASING_EVALS, true, 5); 116 } 117 118 { // -- Get vorticity in S-frame 119 CeedScalar rotation_iso[3][3]; 120 RotationRate(grad_velo_iso, rotation_iso); 121 CeedScalar vorticity_iso[3] = {-2 * rotation_iso[1][2], 2 * rotation_iso[0][2], -2 * rotation_iso[0][1]}; 122 OrientBasisWithVector(eigenvectors, vorticity_iso); 123 MatVec3(eigenvectors, vorticity_iso, CEED_NOTRANSPOSE, vorticity_sframe); 124 } 125 126 // -- Setup DD model inputs 127 const CeedScalar grad_velo_magnitude = VelocityGradientMagnitude(strain_sframe, vorticity_sframe); 128 CeedScalar inputs[6] = {strain_sframe[0], strain_sframe[1], strain_sframe[2], vorticity_sframe[0], vorticity_sframe[1], viscosity / Square(delta)}; 129 ScaleN(inputs, 1 / (grad_velo_magnitude + CEED_EPSILON), 6); 130 131 CeedScalar sgs_sframe_sym[6] = {0.}; 132 DataDrivenInference(inputs, sgs_sframe_sym, sgsdd_ctx); 133 134 CeedScalar old_bounds[6][2] = {{0}}; 135 for (int j = 0; j < 6; j++) old_bounds[j][1] = 1; 136 const CeedScalar(*new_bounds)[2] = (const CeedScalar(*)[2]) & sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling]; 137 DenormalizeDDOutputs(sgs_sframe_sym, new_bounds, old_bounds); 138 139 // Re-dimensionalize sgs_stress 140 ScaleN(sgs_sframe_sym, Square(delta) * Square(grad_velo_magnitude), 6); 141 142 CeedScalar sgs_stress[3][3] = {{0.}}; 143 { // Rotate SGS Stress back to physical frame, SGS_physical = E^T SGS_sframe E 144 CeedScalar Evec_sgs[3][3] = {{0.}}; 145 const CeedScalar sgs_sframe[3][3] = { 146 {sgs_sframe_sym[0], sgs_sframe_sym[3], sgs_sframe_sym[4]}, 147 {sgs_sframe_sym[3], sgs_sframe_sym[1], sgs_sframe_sym[5]}, 148 {sgs_sframe_sym[4], sgs_sframe_sym[5], sgs_sframe_sym[2]}, 149 }; 150 MatMat3(eigenvectors, sgs_sframe, CEED_TRANSPOSE, CEED_NOTRANSPOSE, Evec_sgs); 151 MatMat3(Evec_sgs, eigenvectors, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, sgs_stress); 152 } 153 154 KMPack(sgs_stress, kmsgs_stress); 155 } 156 157 #endif // sgs_dd_model_h 158