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 // ***************************************************************************** 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 // ***************************************************************************** 24 CEED_QFUNCTION_HELPER void StabilizationMatrix(const NewtonianIdealGasContext 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 for (CeedInt i = 0; i < 5; i++) { 30 for (CeedInt j = 0; j < 3; j++) stab[i][j] = 0; 31 } 32 dYFromTau(strong_residual, Tau_d, dY); 33 State ds = StateFromY_fwd(gas, s, dY); 34 FluxInviscid_fwd(gas, s, ds, dF); 35 for (CeedInt i = 0; i < 3; i++) { 36 CeedScalar dF_i[5]; 37 UnpackState_U(dF[i], dF_i); 38 for (CeedInt j = 0; j < 5; j++) stab[j][i] += dF_i[j]; 39 } 40 } 41 42 CEED_QFUNCTION_HELPER void Stabilization(const NewtonianIdealGasContext gas, const State s, const CeedScalar Tau_d[3], const State ds[3], 43 const CeedScalar U_dot[5], const CeedScalar body_force[5], const CeedScalar divFdiff[5], 44 CeedScalar stab[5][3]) { 45 // -- Stabilization method: none (Galerkin), SU, or SUPG 46 CeedScalar strong_residual[5] = {0}; 47 switch (gas->stabilization) { 48 case STAB_NONE: 49 break; 50 case STAB_SU: 51 FluxInviscidStrong(gas, s, ds, strong_residual); 52 break; 53 case STAB_SUPG: 54 FluxInviscidStrong(gas, s, ds, strong_residual); 55 for (CeedInt j = 0; j < 5; j++) strong_residual[j] += U_dot[j] - body_force[j] + divFdiff[j]; 56 break; 57 } 58 StabilizationMatrix(gas, s, Tau_d, strong_residual, stab); 59 } 60 61 // ***************************************************************************** 62 // Helper function for computing Tau elements (stabilization constant) 63 // Model from: 64 // PHASTA 65 // 66 // Tau[i] = itau=0 which is diagonal-Shakib (3 values still but not spatial) 67 // ***************************************************************************** 68 CEED_QFUNCTION_HELPER void Tau_diagPrim(NewtonianIdealGasContext gas, State s, const CeedScalar dXdx[3][3], const CeedScalar dt, 69 CeedScalar Tau_d[3]) { 70 // Context 71 const CeedScalar Ctau_t = gas->Ctau_t; 72 const CeedScalar Ctau_v = gas->Ctau_v; 73 const CeedScalar Ctau_C = gas->Ctau_C; 74 const CeedScalar Ctau_M = gas->Ctau_M; 75 const CeedScalar Ctau_E = gas->Ctau_E; 76 const CeedScalar cv = gas->cv; 77 const CeedScalar mu = gas->mu; 78 const CeedScalar rho = s.U.density; 79 80 CeedScalar tau; 81 CeedScalar dts; 82 CeedScalar fact; 83 84 CeedScalar gijd_mat[3][3] = {{0.}}, velocity_term; 85 MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, gijd_mat); 86 87 dts = Ctau_t / dt; 88 89 { // u_i g_ij u_j 90 CeedScalar gij_uj[3] = {0.}; 91 MatVec3(gijd_mat, s.Y.velocity, CEED_NOTRANSPOSE, gij_uj); 92 velocity_term = Dot3(s.Y.velocity, gij_uj); 93 } 94 95 tau = Square(rho) * (4. * Square(dts) + velocity_term) + Ctau_v * Square(mu) * DotN((CeedScalar *)gijd_mat, (CeedScalar *)gijd_mat, 9); 96 97 fact = sqrt(tau); 98 99 Tau_d[0] = Ctau_C * fact / (rho * (gijd_mat[0][0] + gijd_mat[1][1] + gijd_mat[2][2])) * 0.125; 100 Tau_d[1] = Ctau_M / fact; 101 Tau_d[2] = Ctau_E / (fact * cv); 102 103 // consider putting back the way I initially had it 104 // 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 105 // but in that case energy tau is scaled by the product of Ctau_E * Ctau_M 106 // 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 107 // Ctau_v * mu * mu IF AND ONLY IF we don't add viscosity law =f(T) 108 } 109