1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Structs and helper functions to evaluate data-driven subgrid-stress modeling 6 /// See 'Invariant data-driven subgrid stress modeling in the strain-rate eigenframe for large eddy simulation' 2022 and 'S-frame discrepancy 7 /// correction models for data-informed Reynolds stress closure' 2022 8 #include <ceed/types.h> 9 10 #include "newtonian_state.h" 11 #include "newtonian_types.h" 12 #include "sgs_dd_utils.h" 13 #include "utils.h" 14 #include "utils_eigensolver_jacobi.h" 15 16 typedef struct SgsDDContext_ *SgsDDContext; 17 struct SgsDDContext_ { 18 CeedInt num_inputs, num_outputs; 19 CeedInt num_layers; 20 CeedInt num_neurons; 21 CeedScalar alpha; 22 23 struct NewtonianIdealGasContext_ newt_ctx; 24 struct { 25 size_t bias1, bias2; 26 size_t weight1, weight2; 27 size_t out_scaling; 28 } offsets; 29 size_t total_bytes; 30 CeedScalar data[1]; 31 }; 32 33 CEED_QFUNCTION_HELPER void LeakyReLU(CeedScalar *x, const CeedScalar alpha, const CeedInt N) { 34 for (CeedInt i = 0; i < N; i++) x[i] *= (x[i] < 0 ? alpha : 1.); 35 } 36 37 CEED_QFUNCTION_HELPER void DataDrivenInference(const CeedScalar *inputs, CeedScalar *outputs, SgsDDContext sgsdd_ctx) { 38 const CeedInt num_neurons = sgsdd_ctx->num_neurons; 39 const CeedInt num_inputs = sgsdd_ctx->num_inputs; 40 const CeedInt num_outputs = sgsdd_ctx->num_outputs; 41 const CeedScalar alpha = sgsdd_ctx->alpha; 42 const CeedScalar *bias1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias1]; 43 const CeedScalar *bias2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.bias2]; 44 const CeedScalar *weight1 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight1]; 45 const CeedScalar *weight2 = &sgsdd_ctx->data[sgsdd_ctx->offsets.weight2]; 46 CeedScalar V[20] = {0.}; 47 48 CopyN(bias1, V, num_neurons); 49 MatVecNM(weight1, inputs, num_neurons, num_inputs, CEED_NOTRANSPOSE, V); 50 LeakyReLU(V, alpha, num_neurons); 51 CopyN(bias2, outputs, num_outputs); 52 MatVecNM(weight2, V, num_outputs, num_neurons, CEED_NOTRANSPOSE, outputs); 53 } 54 55 CEED_QFUNCTION_HELPER void ComputeSgsDD_Fused(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta, 56 const CeedScalar viscosity, CeedScalar kmsgs_stress[6], SgsDDContext sgsdd_ctx) { 57 CeedScalar inputs[6], grad_velo_magnitude, eigenvectors[3][3], sgs_sframe_sym[6] = {0.}, new_bounds[6][2]; 58 // Copying new_bounds because Sycl online compiler doesn't like direct casting the pointer 59 CopyN(&sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling], (CeedScalar *)new_bounds, 12); 60 61 ComputeSgsDDInputs(grad_velo_aniso, km_A_ij, delta, viscosity, eigenvectors, inputs, &grad_velo_magnitude); 62 DataDrivenInference(inputs, sgs_sframe_sym, sgsdd_ctx); 63 ComputeSgsDDOutputs(sgs_sframe_sym, delta, eigenvectors, new_bounds, grad_velo_magnitude, kmsgs_stress); 64 } 65 66 // @brief Calculate subgrid stress at nodes using anisotropic data-driven model 67 CEED_QFUNCTION_HELPER int ComputeSgsDDNodal_Fused(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 68 StateVariable state_var) { 69 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 70 const CeedScalar(*grad_velo)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[2]; 71 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 72 const CeedScalar(*inv_multiplicity) = (const CeedScalar(*))in[4]; 73 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 74 75 const SgsDDContext sgsdd_ctx = (SgsDDContext)ctx; 76 const NewtonianIGProperties gas = sgsdd_ctx->newt_ctx.gas; 77 78 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 79 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 80 const CeedScalar grad_velo_aniso[3][3] = { 81 {grad_velo[0][0][i], grad_velo[0][1][i], grad_velo[0][2][i]}, 82 {grad_velo[1][0][i], grad_velo[1][1][i], grad_velo[1][2][i]}, 83 {grad_velo[2][0][i], grad_velo[2][1][i], grad_velo[2][2][i]} 84 }; 85 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]}; 86 const CeedScalar delta = A_ij_delta[6][i]; 87 const State s = StateFromQ(gas, qi, state_var); 88 CeedScalar km_sgs[6]; 89 90 ComputeSgsDD_Fused(grad_velo_aniso, km_A_ij, delta, gas.mu / s.U.density, km_sgs, sgsdd_ctx); 91 92 for (int j = 0; j < 6; j++) v[j][i] = inv_multiplicity[i] * km_sgs[j]; 93 } 94 return 0; 95 } 96 97 CEED_QFUNCTION(ComputeSgsDDNodal_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 98 return ComputeSgsDDNodal_Fused(ctx, Q, in, out, STATEVAR_PRIMITIVE); 99 } 100 101 CEED_QFUNCTION(ComputeSgsDDNodal_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 102 return ComputeSgsDDNodal_Fused(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 103 } 104 105 CEED_QFUNCTION(ComputeSgsDDNodal_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 106 return ComputeSgsDDNodal_Fused(ctx, Q, in, out, STATEVAR_ENTROPY); 107 } 108 109 // @brief Calculate inputs to anisotropic data-driven model 110 CEED_QFUNCTION_HELPER int ComputeSgsDDNodal_Sequential_Inputs(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 111 StateVariable state_var) { 112 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 113 const CeedScalar(*grad_velo)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1]; 114 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 115 const CeedScalar(*inv_multiplicity) = (const CeedScalar(*))in[3]; 116 CeedScalar(*eigenvectors_stored) = out[0]; 117 CeedScalar(*model_inputs)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 118 119 const SgsDDContext sgsdd_ctx = (SgsDDContext)ctx; 120 const NewtonianIGProperties gas = sgsdd_ctx->newt_ctx.gas; 121 122 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 123 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 124 const CeedScalar grad_velo_aniso[3][3] = { 125 {grad_velo[0][0][i], grad_velo[0][1][i], grad_velo[0][2][i]}, 126 {grad_velo[1][0][i], grad_velo[1][1][i], grad_velo[1][2][i]}, 127 {grad_velo[2][0][i], grad_velo[2][1][i], grad_velo[2][2][i]} 128 }; 129 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]}; 130 const CeedScalar delta = A_ij_delta[6][i]; 131 const State s = StateFromQ(gas, qi, state_var); 132 133 CeedScalar model_inputs_i[6], grad_velo_magnitude, eigenvectors[3][3]; 134 ComputeSgsDDInputs(grad_velo_aniso, km_A_ij, delta, gas.mu / s.U.density, eigenvectors, model_inputs_i, &grad_velo_magnitude); 135 136 ScaleN(model_inputs_i, inv_multiplicity[i], 6); 137 StoredValuesPack(Q, i, 0, 6, model_inputs_i, (CeedScalar *)model_inputs); 138 StoredValuesPack(Q, i, 0, 9, (const CeedScalar *)eigenvectors, eigenvectors_stored); 139 StoredValuesPack(Q, i, 9, 1, &grad_velo_magnitude, eigenvectors_stored); 140 } 141 return CEED_ERROR_SUCCESS; 142 } 143 144 CEED_QFUNCTION(ComputeSgsDDNodal_Sequential_Inputs_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 145 return ComputeSgsDDNodal_Sequential_Inputs(ctx, Q, in, out, STATEVAR_PRIMITIVE); 146 } 147 148 CEED_QFUNCTION(ComputeSgsDDNodal_Sequential_Inputs_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 149 return ComputeSgsDDNodal_Sequential_Inputs(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 150 } 151 152 CEED_QFUNCTION(ComputeSgsDDNodal_Sequential_Inputs_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 153 return ComputeSgsDDNodal_Sequential_Inputs(ctx, Q, in, out, STATEVAR_ENTROPY); 154 } 155 156 // @brief Runs inference on the data-driven model, used predominantsly for testing and validation 157 CEED_QFUNCTION(ComputeSgsDDNodal_Sequential_Inference)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 158 const CeedScalar(*model_inputs) = in[0]; 159 const CeedScalar(*inv_multiplicity) = in[1]; 160 CeedScalar(*model_outputs) = out[0]; 161 162 const SgsDDContext sgsdd_ctx = (SgsDDContext)ctx; 163 164 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 165 CeedScalar model_inputs_i[6], model_outputs_i[6]; 166 167 StoredValuesUnpack(Q, i, 0, 6, (const CeedScalar *)model_inputs, model_inputs_i); 168 DataDrivenInference(model_inputs_i, model_outputs_i, sgsdd_ctx); 169 ScaleN(model_outputs_i, inv_multiplicity[i], 6); 170 StoredValuesPack(Q, i, 0, 6, model_outputs_i, model_outputs); 171 } 172 return CEED_ERROR_SUCCESS; 173 } 174 175 // @brief Calculates SGS from outputs of anisotropic data-driven model 176 CEED_QFUNCTION(ComputeSgsDDNodal_Sequential_Outputs)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 177 const CeedScalar(*model_outputs) = in[0]; 178 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 179 const CeedScalar(*inv_multiplicity) = (const CeedScalar(*))in[2]; 180 const CeedScalar(*eigenvectors_stored) = in[3]; 181 CeedScalar(*kmsgs_stress)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 182 183 const SgsDDContext sgsdd_ctx = (SgsDDContext)ctx; 184 CeedScalar new_bounds[6][2]; 185 CopyN(&sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling], (CeedScalar *)new_bounds, 12); 186 187 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 188 CeedScalar model_outputs_i[6]; 189 const CeedScalar delta = A_ij_delta[6][i]; 190 191 StoredValuesUnpack(Q, i, 0, 6, model_outputs, model_outputs_i); 192 CeedScalar grad_velo_magnitude, eigenvectors[3][3], kmsgs_stress_i[6]; 193 StoredValuesUnpack(Q, i, 0, 9, eigenvectors_stored, (CeedScalar *)eigenvectors); 194 StoredValuesUnpack(Q, i, 9, 1, eigenvectors_stored, &grad_velo_magnitude); 195 ComputeSgsDDOutputs(model_outputs_i, delta, eigenvectors, new_bounds, grad_velo_magnitude, kmsgs_stress_i); 196 197 for (int j = 0; j < 6; j++) kmsgs_stress[j][i] = inv_multiplicity[i] * kmsgs_stress_i[j]; 198 } 199 return CEED_ERROR_SUCCESS; 200 } 201 202 // @brief Adds subgrid stress to residual (during IFunction evaluation) 203 CEED_QFUNCTION_HELPER int FluxSubgridStress(const StatePrimitive Y, const CeedScalar km_sgs[6], CeedScalar Flux[5][3]) { 204 CeedScalar sgs[3][3]; 205 206 KMUnpack(km_sgs, sgs); 207 for (CeedInt j = 0; j < 3; j++) { 208 Flux[0][j] = 0.; 209 for (CeedInt k = 0; k < 3; k++) Flux[k + 1][j] = sgs[k][j]; 210 Flux[4][j] = Y.velocity[0] * sgs[0][j] + Y.velocity[1] * sgs[1][j] + Y.velocity[2] * sgs[2][j]; 211 } 212 return 0; 213 } 214 215 CEED_QFUNCTION_HELPER int IFunction_NodalSgs(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 216 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 217 const CeedScalar(*q_data) = in[1]; 218 const CeedScalar(*km_sgs)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 219 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[0]; 220 221 NewtonianIdealGasContext newt_ctx = (NewtonianIdealGasContext)ctx; 222 223 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 224 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 225 const State s = StateFromQ(newt_ctx->gas, qi, state_var); 226 227 CeedScalar wdetJ, dXdx[3][3]; 228 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 229 230 CeedScalar Flux[5][3]; 231 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]}; 232 FluxSubgridStress(s.Y, km_sgs_i, Flux); 233 234 for (CeedInt k = 0; k < 3; k++) { 235 for (CeedInt j = 0; j < 5; j++) { 236 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 237 } 238 } 239 } 240 return 0; 241 } 242 243 CEED_QFUNCTION(IFunction_NodalSgs_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 244 return IFunction_NodalSgs(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 245 } 246 247 CEED_QFUNCTION(IFunction_NodalSgs_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 248 return IFunction_NodalSgs(ctx, Q, in, out, STATEVAR_PRIMITIVE); 249 } 250 251 CEED_QFUNCTION(IFunction_NodalSgs_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 252 return IFunction_NodalSgs(ctx, Q, in, out, STATEVAR_ENTROPY); 253 } 254