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 /// Structs and helper functions regarding the state of a newtonian simulation 10 11 12 #ifndef newtonian_state_h 13 #define newtonian_state_h 14 15 #include <math.h> 16 #include <ceed.h> 17 #include "newtonian_types.h" 18 19 typedef struct { 20 CeedScalar pressure; 21 CeedScalar velocity[3]; 22 CeedScalar temperature; 23 } StatePrimitive; 24 25 typedef struct { 26 CeedScalar density; 27 CeedScalar momentum[3]; 28 CeedScalar E_total; 29 } StateConservative; 30 31 typedef struct { 32 StateConservative U; 33 StatePrimitive Y; 34 } State; 35 36 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative( 37 NewtonianIdealGasContext gas, StateConservative U, const CeedScalar x[3]) { 38 StatePrimitive Y; 39 for (CeedInt i=0; i<3; i++) Y.velocity[i] = U.momentum[i] / U.density; 40 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 41 CeedScalar e_potential = -Dot3(gas->g, x); 42 CeedScalar e_total = U.E_total / U.density; 43 CeedScalar e_internal = e_total - e_kinetic - e_potential; 44 Y.temperature = e_internal / gas->cv; 45 Y.pressure = (gas->cp / gas->cv - 1) * U.density * e_internal; 46 return Y; 47 } 48 49 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd( 50 NewtonianIdealGasContext gas, State s, StateConservative dU, 51 const CeedScalar x[3], const CeedScalar dx[3]) { 52 StatePrimitive dY; 53 for (CeedInt i=0; i<3; i++) { 54 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 55 } 56 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 57 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 58 CeedScalar e_potential = -Dot3(gas->g, x); 59 CeedScalar de_potential = -Dot3(gas->g, dx); 60 CeedScalar e_total = s.U.E_total / s.U.density; 61 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 62 CeedScalar e_internal = e_total - e_kinetic - e_potential; 63 CeedScalar de_internal = de_total - de_kinetic - de_potential; 64 dY.temperature = de_internal / gas->cv; 65 dY.pressure = (gas->cp / gas->cv - 1) 66 * (dU.density * e_internal + s.U.density * de_internal); 67 return dY; 68 } 69 70 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, 71 const CeedScalar U[5], const CeedScalar x[3]) { 72 State s; 73 s.U.density = U[0]; 74 s.U.momentum[0] = U[1]; 75 s.U.momentum[1] = U[2]; 76 s.U.momentum[2] = U[3]; 77 s.U.E_total = U[4]; 78 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 79 return s; 80 } 81 82 CEED_QFUNCTION_HELPER State StateFromU_fwd(NewtonianIdealGasContext gas, 83 State s, const CeedScalar dU[5], 84 const CeedScalar x[3], const CeedScalar dx[3]) { 85 State ds; 86 ds.U.density = dU[0]; 87 ds.U.momentum[0] = dU[1]; 88 ds.U.momentum[1] = dU[2]; 89 ds.U.momentum[2] = dU[3]; 90 ds.U.E_total = dU[4]; 91 ds.Y = StatePrimitiveFromConservative_fwd(gas, s, ds.U, x, dx); 92 return ds; 93 } 94 95 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, 96 StateConservative Flux[3]) { 97 for (CeedInt i=0; i<3; i++) { 98 Flux[i].density = s.U.momentum[i]; 99 for (CeedInt j=0; j<3; j++) 100 Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] 101 + s.Y.pressure * (i == j); 102 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 103 } 104 } 105 106 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, 107 State s, State ds, StateConservative dFlux[3]) { 108 for (CeedInt i=0; i<3; i++) { 109 dFlux[i].density = ds.U.momentum[i]; 110 for (CeedInt j=0; j<3; j++) 111 dFlux[i].momentum[j] = ds.U.momentum[i] * s.Y.velocity[j] + 112 s.U.momentum[i] * ds.Y.velocity[j] + ds.Y.pressure * (i == j); 113 dFlux[i].E_total = (ds.U.E_total + ds.Y.pressure) * s.Y.velocity[i] + 114 (s.U.E_total + s.Y.pressure) * ds.Y.velocity[i]; 115 } 116 } 117 118 // Kelvin-Mandel notation 119 CEED_QFUNCTION_HELPER void KMStrainRate(const State grad_s[3], 120 CeedScalar strain_rate[6]) { 121 const CeedScalar weight = 1 / sqrt(2.); 122 strain_rate[0] = grad_s[0].Y.velocity[0]; 123 strain_rate[1] = grad_s[1].Y.velocity[1]; 124 strain_rate[2] = grad_s[2].Y.velocity[2]; 125 strain_rate[3] = weight * (grad_s[2].Y.velocity[1] + grad_s[1].Y.velocity[2]); 126 strain_rate[4] = weight * (grad_s[2].Y.velocity[0] + grad_s[0].Y.velocity[2]); 127 strain_rate[5] = weight * (grad_s[1].Y.velocity[0] + grad_s[0].Y.velocity[1]); 128 } 129 130 CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) { 131 const CeedScalar weight = 1 / sqrt(2.); 132 A[0][0] = v[0]; 133 A[1][1] = v[1]; 134 A[2][2] = v[2]; 135 A[2][1] = A[1][2] = weight * v[3]; 136 A[2][0] = A[0][2] = weight * v[4]; 137 A[1][0] = A[0][1] = weight * v[5]; 138 } 139 140 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, 141 const CeedScalar strain_rate[6], CeedScalar stress[6]) { 142 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 143 for (CeedInt i=0; i<6; i++) { 144 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 145 } 146 } 147 148 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, 149 StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 150 CeedScalar Fe[3]) { 151 for (CeedInt i=0; i<3; i++) { 152 Fe[i] = - Y.velocity[0] * stress[0][i] 153 - Y.velocity[1] * stress[1][i] 154 - Y.velocity[2] * stress[2][i] 155 - gas->k * grad_s[i].Y.temperature; 156 } 157 } 158 159 CEED_QFUNCTION_HELPER void ViscousEnergyFlux_fwd(NewtonianIdealGasContext gas, 160 StatePrimitive Y, StatePrimitive dY, const State grad_ds[3], 161 const CeedScalar stress[3][3], 162 const CeedScalar dstress[3][3], 163 CeedScalar dFe[3]) { 164 for (CeedInt i=0; i<3; i++) { 165 dFe[i] = - Y.velocity[0] * dstress[0][i] - dY.velocity[0] * stress[0][i] 166 - Y.velocity[1] * dstress[1][i] - dY.velocity[1] * stress[1][i] 167 - Y.velocity[2] * dstress[2][i] - dY.velocity[2] * stress[2][i] 168 - gas->k * grad_ds[i].Y.temperature; 169 } 170 } 171 172 #endif // newtonian_state_h 173