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