xref: /honee/qfunctions/sgs_dd_training.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 for training data-driven subgrid-stress models
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 "differential_filter_enums.h"
15 #include "newtonian_state.h"
16 #include "newtonian_types.h"
17 #include "sgs_dd_utils.h"
18 #include "utils.h"
19 #include "utils_eigensolver_jacobi.h"
20 
21 typedef struct SGS_DD_TrainingContext_ *SGS_DDTrainingContext;
22 struct SGS_DD_TrainingContext_ {
23   struct NewtonianIdealGasContext_ gas;
24 };
25 
26 // @brief Calculate Data-Driven SGS model training data at nodes
27 CEED_QFUNCTION_HELPER int ComputeSGS_DDAnisotropicTrainingDataNodal(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
28                                                                     StateVariable state_var) {
29   const CeedScalar(*q)[CEED_Q_VLA]            = (const CeedScalar(*)[CEED_Q_VLA])in[0];
30   const CeedScalar(*velo_prod)[CEED_Q_VLA]    = (const CeedScalar(*)[CEED_Q_VLA])in[1];
31   const CeedScalar(*grad_velo)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[2];
32   const CeedScalar(*A_ij_delta)[CEED_Q_VLA]   = (const CeedScalar(*)[CEED_Q_VLA])in[3];
33   const CeedScalar(*inv_multiplicity)         = (const CeedScalar(*))in[4];
34   CeedScalar(*v)[CEED_Q_VLA]                  = (CeedScalar(*)[CEED_Q_VLA])out[0];
35 
36   const SGS_DDTrainingContext    sgsdd_ctx = (SGS_DDTrainingContext)ctx;
37   const NewtonianIdealGasContext gas       = &sgsdd_ctx->gas;
38 
39   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
40     const CeedScalar qi[5]                 = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
41     const CeedScalar grad_velo_aniso[3][3] = {
42         {grad_velo[0][0][i], grad_velo[0][1][i], grad_velo[0][2][i]},
43         {grad_velo[1][0][i], grad_velo[1][1][i], grad_velo[1][2][i]},
44         {grad_velo[2][0][i], grad_velo[2][1][i], grad_velo[2][2][i]}
45     };
46     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]};
47     const CeedScalar delta      = A_ij_delta[6][i];
48     const State      s          = StateFromQ(gas, qi, state_var);
49     CeedScalar       inputs[6];
50     CeedScalar       eigenvectors[3][3], grad_velo_magnitude;  // dummy variables, don't actually use them
51 
52     ComputeSgsDDInputs(grad_velo_aniso, km_A_ij, delta, gas->mu / s.U.density, eigenvectors, inputs, &grad_velo_magnitude);
53 
54     for (int j = 0; j < 6; j++) v[j][i] = inv_multiplicity[i] * inputs[j];
55 
56     v[0 + 6][i] = (velo_prod[DIFF_FILTER_VELOCITY_SQUARED_XX][i] - Square(s.Y.velocity[0])) * inv_multiplicity[i];
57     v[1 + 6][i] = (velo_prod[DIFF_FILTER_VELOCITY_SQUARED_YY][i] - Square(s.Y.velocity[1])) * inv_multiplicity[i];
58     v[2 + 6][i] = (velo_prod[DIFF_FILTER_VELOCITY_SQUARED_ZZ][i] - Square(s.Y.velocity[2])) * inv_multiplicity[i];
59     v[3 + 6][i] = (velo_prod[DIFF_FILTER_VELOCITY_SQUARED_YZ][i] - s.Y.velocity[1] * s.Y.velocity[2]) * inv_multiplicity[i];
60     v[4 + 6][i] = (velo_prod[DIFF_FILTER_VELOCITY_SQUARED_XZ][i] - s.Y.velocity[0] * s.Y.velocity[2]) * inv_multiplicity[i];
61     v[5 + 6][i] = (velo_prod[DIFF_FILTER_VELOCITY_SQUARED_XY][i] - s.Y.velocity[0] * s.Y.velocity[1]) * inv_multiplicity[i];
62   }
63   return 0;
64 }
65 
66 CEED_QFUNCTION(ComputeSGS_DDAnisotropicTrainingDataNodal_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
67   return ComputeSGS_DDAnisotropicTrainingDataNodal(ctx, Q, in, out, STATEVAR_PRIMITIVE);
68 }
69