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