1 // Copyright (c) 2017-2022, 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 11 12 #ifndef stabilization_h 13 #define stabilization_h 14 15 #include <ceed.h> 16 #include "newtonian_state.h" 17 18 // ***************************************************************************** 19 // Helper function for computing the variation in primitive variables, 20 // given Tau_d 21 // ***************************************************************************** 22 CEED_QFUNCTION_HELPER void dYFromTau(CeedScalar Y[5], CeedScalar Tau_d[3], 23 CeedScalar dY[5]) { 24 dY[0] = Tau_d[0] * Y[0]; 25 dY[1] = Tau_d[1] * Y[1]; 26 dY[2] = Tau_d[1] * Y[2]; 27 dY[3] = Tau_d[1] * Y[3]; 28 dY[4] = Tau_d[2] * Y[4]; 29 } 30 31 // ***************************************************************************** 32 // Helper functions for computing the stabilization terms 33 // ***************************************************************************** 34 CEED_QFUNCTION_HELPER void StabilizationMatrix(NewtonianIdealGasContext gas, 35 State s, CeedScalar Tau_d[3], CeedScalar R[5], const CeedScalar x[3], 36 CeedScalar stab[5][3]) { 37 CeedScalar dY[5]; 38 const CeedScalar dx_i[3] = {0}; 39 StateConservative dF[3]; 40 // Zero stab so all future terms can safely sum into it 41 for (CeedInt i=0; i<5; i++) 42 for (CeedInt j=0; j<3; j++) 43 stab[i][j] = 0; 44 dYFromTau(R, Tau_d, dY); 45 State ds = StateFromY_fwd(gas, s, dY, x, dx_i); 46 FluxInviscid_fwd(gas, s, ds, dF); 47 for (CeedInt i=0; i<3; i++) { 48 CeedScalar dF_i[5]; 49 UnpackState_U(dF[i], dF_i); 50 for (CeedInt j=0; j<5; j++) 51 stab[j][i] += dF_i[j]; 52 } 53 } 54 55 CEED_QFUNCTION_HELPER void Stabilization(NewtonianIdealGasContext gas, State s, 56 CeedScalar Tau_d[3], State ds[3], CeedScalar U_dot[5], 57 const CeedScalar body_force[5], const CeedScalar x[3], CeedScalar stab[5][3]) { 58 // -- Stabilization method: none (Galerkin), SU, or SUPG 59 CeedScalar R[5] = {0}; 60 switch (gas->stabilization) { 61 case STAB_NONE: 62 break; 63 case STAB_SU: 64 FluxInviscidStrong(gas, s, ds, R); 65 break; 66 case STAB_SUPG: 67 FluxInviscidStrong(gas, s, ds, R); 68 for (CeedInt j=0; j<5; j++) R[j] += U_dot[j] - body_force[j]; 69 break; 70 } 71 StabilizationMatrix(gas, s, Tau_d, R, x, stab); 72 } 73 74 // ***************************************************************************** 75 // Helper function for computing Tau elements (stabilization constant) 76 // Model from: 77 // PHASTA 78 // 79 // Tau[i] = itau=0 which is diagonal-Shakib (3 values still but not spatial) 80 // 81 // ***************************************************************************** 82 CEED_QFUNCTION_HELPER void Tau_diagPrim(NewtonianIdealGasContext gas, State s, 83 const CeedScalar dXdx[3][3], 84 const CeedScalar dt, CeedScalar Tau_d[3]) { 85 // Context 86 const CeedScalar Ctau_t = gas->Ctau_t; 87 const CeedScalar Ctau_v = gas->Ctau_v; 88 const CeedScalar Ctau_C = gas->Ctau_C; 89 const CeedScalar Ctau_M = gas->Ctau_M; 90 const CeedScalar Ctau_E = gas->Ctau_E; 91 const CeedScalar cv = gas->cv; 92 const CeedScalar mu = gas->mu; 93 const CeedScalar u[3] = {s.Y.velocity[0], s.Y.velocity[1], s.Y.velocity[2]}; 94 const CeedScalar rho = s.U.density; 95 96 CeedScalar gijd[6]; 97 CeedScalar tau; 98 CeedScalar dts; 99 CeedScalar fact; 100 101 //*INDENT-OFF* 102 gijd[0] = dXdx[0][0] * dXdx[0][0] 103 + dXdx[1][0] * dXdx[1][0] 104 + dXdx[2][0] * dXdx[2][0]; 105 106 gijd[1] = dXdx[0][0] * dXdx[0][1] 107 + dXdx[1][0] * dXdx[1][1] 108 + dXdx[2][0] * dXdx[2][1]; 109 110 gijd[2] = dXdx[0][1] * dXdx[0][1] 111 + dXdx[1][1] * dXdx[1][1] 112 + dXdx[2][1] * dXdx[2][1]; 113 114 gijd[3] = dXdx[0][0] * dXdx[0][2] 115 + dXdx[1][0] * dXdx[1][2] 116 + dXdx[2][0] * dXdx[2][2]; 117 118 gijd[4] = dXdx[0][1] * dXdx[0][2] 119 + dXdx[1][1] * dXdx[1][2] 120 + dXdx[2][1] * dXdx[2][2]; 121 122 gijd[5] = dXdx[0][2] * dXdx[0][2] 123 + dXdx[1][2] * dXdx[1][2] 124 + dXdx[2][2] * dXdx[2][2]; 125 //*INDENT-ON* 126 127 dts = Ctau_t / dt ; 128 129 tau = rho*rho*((4. * dts * dts) 130 + u[0] * ( u[0] * gijd[0] + 2. * ( u[1] * gijd[1] + u[2] * gijd[3])) 131 + u[1] * ( u[1] * gijd[2] + 2. * u[2] * gijd[4]) 132 + u[2] * u[2] * gijd[5]) 133 + Ctau_v* mu * mu * 134 (gijd[0]*gijd[0] + gijd[2]*gijd[2] + gijd[5]*gijd[5] + 135 + 2. * (gijd[1]*gijd[1] + gijd[3]*gijd[3] + gijd[4]*gijd[4])); 136 137 fact = sqrt(tau); 138 139 Tau_d[0] = Ctau_C * fact / (rho*(gijd[0] + gijd[2] + gijd[5]))*0.125; 140 141 Tau_d[1] = Ctau_M / fact; 142 Tau_d[2] = Ctau_E / ( fact * cv ); 143 144 // consider putting back the way I initially had it Ctau_E * Tau_d[1] /cv 145 // to avoid a division if the compiler is smart enough to see that cv IS 146 // a constant that it could invert once for all elements 147 // but in that case energy tau is scaled by the product of Ctau_E * Ctau_M 148 // OR we could absorb cv into Ctau_E but this puts more burden on user to 149 // know how to change constants with a change of fluid or units. Same for 150 // Ctau_v * mu * mu IF AND ONLY IF we don't add viscosity law =f(T) 151 } 152 153 // ***************************************************************************** 154 155 #endif // stabilization_h 156