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