1 // Copyright (c) 2017-2026, 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 /// Helper functions for computing stabilization terms of a newtonian simulation
10 #include <ceed/types.h>
11
12 #include "newtonian_state.h"
13
14 // *****************************************************************************
15 // Helper function for computing the variation in primitive variables, given Tau_d
16 // *****************************************************************************
dYFromTau(CeedScalar Y[5],CeedScalar Tau_d[3],CeedScalar dY[5])17 CEED_QFUNCTION_HELPER void dYFromTau(CeedScalar Y[5], CeedScalar Tau_d[3], CeedScalar dY[5]) {
18 dY[0] = Tau_d[0] * Y[0];
19 dY[1] = Tau_d[1] * Y[1];
20 dY[2] = Tau_d[1] * Y[2];
21 dY[3] = Tau_d[1] * Y[3];
22 dY[4] = Tau_d[2] * Y[4];
23 }
24
25 // *****************************************************************************
26 // Helper functions for computing the stabilization terms
27 // *****************************************************************************
StabilizationMatrix(NewtonianIdealGasContext gas,State s,CeedScalar Tau_d[3],CeedScalar strong_residual[5],CeedScalar stab[5][3])28 CEED_QFUNCTION_HELPER void StabilizationMatrix(NewtonianIdealGasContext gas, State s, CeedScalar Tau_d[3], CeedScalar strong_residual[5],
29 CeedScalar stab[5][3]) {
30 CeedScalar dY[5];
31 StateConservative dF[3];
32 // Zero stab so all future terms can safely sum into it
33 for (CeedInt i = 0; i < 5; i++) {
34 for (CeedInt j = 0; j < 3; j++) stab[i][j] = 0;
35 }
36 dYFromTau(strong_residual, Tau_d, dY);
37 State ds = StateFromY_fwd(gas, s, dY);
38 FluxInviscid_fwd(gas, s, ds, dF);
39 for (CeedInt i = 0; i < 3; i++) {
40 CeedScalar dF_i[5];
41 UnpackState_U(dF[i], dF_i);
42 for (CeedInt j = 0; j < 5; j++) stab[j][i] += dF_i[j];
43 }
44 }
45
Stabilization(NewtonianIdealGasContext gas,State s,CeedScalar Tau_d[3],State ds[3],CeedScalar U_dot[5],const CeedScalar body_force[5],CeedScalar stab[5][3])46 CEED_QFUNCTION_HELPER void Stabilization(NewtonianIdealGasContext gas, State s, CeedScalar Tau_d[3], State ds[3], CeedScalar U_dot[5],
47 const CeedScalar body_force[5], CeedScalar stab[5][3]) {
48 // -- Stabilization method: none (Galerkin), SU, or SUPG
49 CeedScalar strong_residual[5] = {0};
50 switch (gas->stabilization) {
51 case STAB_NONE:
52 break;
53 case STAB_SU:
54 FluxInviscidStrong(gas, s, ds, strong_residual);
55 break;
56 case STAB_SUPG:
57 FluxInviscidStrong(gas, s, ds, strong_residual);
58 for (CeedInt j = 0; j < 5; j++) strong_residual[j] += U_dot[j] - body_force[j];
59 break;
60 }
61 StabilizationMatrix(gas, s, Tau_d, strong_residual, stab);
62 }
63
64 // *****************************************************************************
65 // Helper function for computing Tau elements (stabilization constant)
66 // Model from:
67 // PHASTA
68 //
69 // Tau[i] = itau=0 which is diagonal-Shakib (3 values still but not spatial)
70 // *****************************************************************************
Tau_diagPrim(NewtonianIdealGasContext gas,State s,const CeedScalar dXdx[3][3],const CeedScalar dt,CeedScalar Tau_d[3])71 CEED_QFUNCTION_HELPER void Tau_diagPrim(NewtonianIdealGasContext gas, State s, const CeedScalar dXdx[3][3], const CeedScalar dt,
72 CeedScalar Tau_d[3]) {
73 // Context
74 const CeedScalar Ctau_t = gas->Ctau_t;
75 const CeedScalar Ctau_v = gas->Ctau_v;
76 const CeedScalar Ctau_C = gas->Ctau_C;
77 const CeedScalar Ctau_M = gas->Ctau_M;
78 const CeedScalar Ctau_E = gas->Ctau_E;
79 const CeedScalar cv = gas->cv;
80 const CeedScalar mu = gas->mu;
81 const CeedScalar rho = s.U.density;
82
83 CeedScalar tau;
84 CeedScalar dts;
85 CeedScalar fact;
86
87 CeedScalar gijd_mat[3][3] = {{0.}}, velocity_term;
88 MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, gijd_mat);
89
90 dts = Ctau_t / dt;
91
92 { // u_i g_ij u_j
93 CeedScalar gij_uj[3] = {0.};
94 MatVec3(gijd_mat, s.Y.velocity, CEED_NOTRANSPOSE, gij_uj);
95 velocity_term = Dot3(s.Y.velocity, gij_uj);
96 }
97
98 tau = Square(rho) * (4. * Square(dts) + velocity_term) + Ctau_v * Square(mu) * DotN((CeedScalar *)gijd_mat, (CeedScalar *)gijd_mat, 9);
99
100 fact = sqrt(tau);
101
102 Tau_d[0] = Ctau_C * fact / (rho * (gijd_mat[0][0] + gijd_mat[1][1] + gijd_mat[2][2])) * 0.125;
103 Tau_d[1] = Ctau_M / fact;
104 Tau_d[2] = Ctau_E / (fact * cv);
105
106 // consider putting back the way I initially had it
107 // Ctau_E * Tau_d[1] /cv to avoid a division if the compiler is smart enough to see that cv IS a constant that it could invert once for all elements
108 // but in that case energy tau is scaled by the product of Ctau_E * Ctau_M
109 // OR we could absorb cv into Ctau_E but this puts more burden on user to know how to change constants with a change of fluid or units. Same for
110 // Ctau_v * mu * mu IF AND ONLY IF we don't add viscosity law =f(T)
111 }
112