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 #include "utils.h" 19 20 typedef struct { 21 CeedScalar pressure; 22 CeedScalar velocity[3]; 23 CeedScalar temperature; 24 } StatePrimitive; 25 26 typedef struct { 27 CeedScalar density; 28 CeedScalar momentum[3]; 29 CeedScalar E_total; 30 } StateConservative; 31 32 typedef struct { 33 StateConservative U; 34 StatePrimitive Y; 35 } State; 36 37 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative( 38 NewtonianIdealGasContext gas, StateConservative U, const CeedScalar x[3]) { 39 StatePrimitive Y; 40 for (CeedInt i=0; i<3; i++) Y.velocity[i] = U.momentum[i] / U.density; 41 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 42 CeedScalar e_potential = -Dot3(gas->g, x); 43 CeedScalar e_total = U.E_total / U.density; 44 CeedScalar e_internal = e_total - e_kinetic - e_potential; 45 Y.temperature = e_internal / gas->cv; 46 Y.pressure = (gas->cp / gas->cv - 1) * U.density * e_internal; 47 return Y; 48 } 49 50 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd( 51 NewtonianIdealGasContext gas, State s, StateConservative dU, 52 const CeedScalar x[3], const CeedScalar dx[3]) { 53 StatePrimitive dY; 54 for (CeedInt i=0; i<3; i++) { 55 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 56 } 57 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 58 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 59 CeedScalar e_potential = -Dot3(gas->g, x); 60 CeedScalar de_potential = -Dot3(gas->g, dx); 61 CeedScalar e_total = s.U.E_total / s.U.density; 62 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 63 CeedScalar e_internal = e_total - e_kinetic - e_potential; 64 CeedScalar de_internal = de_total - de_kinetic - de_potential; 65 dY.temperature = de_internal / gas->cv; 66 dY.pressure = (gas->cp / gas->cv - 1) 67 * (dU.density * e_internal + s.U.density * de_internal); 68 return dY; 69 } 70 71 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive( 72 NewtonianIdealGasContext gas, StatePrimitive Y, const CeedScalar x[3]) { 73 StateConservative U; 74 CeedScalar R = gas->cp - gas->cv; 75 U.density = Y.pressure / (R * Y.temperature); 76 for (int i=0; i<3; i++) U.momentum[i] = U.density*Y.velocity[i]; 77 CeedScalar e_internal = gas->cv * Y.temperature; 78 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 79 CeedScalar e_potential = -Dot3(gas->g, x); 80 CeedScalar e_total = e_internal + e_kinetic + e_potential; 81 U.E_total = U.density*e_total; 82 return U; 83 } 84 85 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive_fwd( 86 NewtonianIdealGasContext gas, State s, StatePrimitive dY, 87 const CeedScalar x[3], const CeedScalar dx[3]) { 88 StateConservative dU; 89 CeedScalar R = gas->cp - gas->cv; 90 dU.density = (dY.pressure * s.Y.temperature - s.Y.pressure * dY.temperature) / 91 (R * s.Y.temperature * s.Y.temperature); 92 for (int i=0; i<3; i++) { 93 dU.momentum[i] = dU.density * s.Y.velocity[i] + s.U.density * dY.velocity[i]; 94 } 95 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 96 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 97 CeedScalar e_potential = -Dot3(gas->g, x); 98 CeedScalar de_potential = -Dot3(gas->g, dx); 99 CeedScalar e_internal = gas->cv * s.Y.temperature; 100 CeedScalar de_internal = gas->cv * dY.temperature; 101 CeedScalar e_total = e_internal + e_kinetic + e_potential; 102 CeedScalar de_total = de_internal + de_kinetic + de_potential; 103 dU.E_total = dU.density*e_total + s.U.density*de_total; 104 return dU; 105 } 106 107 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, 108 const CeedScalar U[5], const CeedScalar x[3]) { 109 State s; 110 s.U.density = U[0]; 111 s.U.momentum[0] = U[1]; 112 s.U.momentum[1] = U[2]; 113 s.U.momentum[2] = U[3]; 114 s.U.E_total = U[4]; 115 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 116 return s; 117 } 118 119 CEED_QFUNCTION_HELPER State StateFromU_fwd(NewtonianIdealGasContext gas, 120 State s, const CeedScalar dU[5], 121 const CeedScalar x[3], const CeedScalar dx[3]) { 122 State ds; 123 ds.U.density = dU[0]; 124 ds.U.momentum[0] = dU[1]; 125 ds.U.momentum[1] = dU[2]; 126 ds.U.momentum[2] = dU[3]; 127 ds.U.E_total = dU[4]; 128 ds.Y = StatePrimitiveFromConservative_fwd(gas, s, ds.U, x, dx); 129 return ds; 130 } 131 132 CEED_QFUNCTION_HELPER State StateFromY(NewtonianIdealGasContext gas, 133 const CeedScalar Y[5], const CeedScalar x[3]) { 134 State s; 135 s.Y.pressure = Y[0]; 136 s.Y.velocity[0] = Y[1]; 137 s.Y.velocity[1] = Y[2]; 138 s.Y.velocity[2] = Y[3]; 139 s.Y.temperature = Y[4]; 140 s.U = StateConservativeFromPrimitive(gas, s.Y, x); 141 return s; 142 } 143 144 CEED_QFUNCTION_HELPER State StateFromY_fwd(NewtonianIdealGasContext gas, 145 State s, const CeedScalar dY[5], 146 const CeedScalar x[3], const CeedScalar dx[3]) { 147 State ds; 148 ds.Y.pressure = dY[0]; 149 ds.Y.velocity[0] = dY[1]; 150 ds.Y.velocity[1] = dY[2]; 151 ds.Y.velocity[2] = dY[3]; 152 ds.Y.temperature = dY[4]; 153 ds.U = StateConservativeFromPrimitive_fwd(gas, s, ds.Y, x, dx); 154 return ds; 155 } 156 157 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, 158 StateConservative Flux[3]) { 159 for (CeedInt i=0; i<3; i++) { 160 Flux[i].density = s.U.momentum[i]; 161 for (CeedInt j=0; j<3; j++) 162 Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] 163 + s.Y.pressure * (i == j); 164 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 165 } 166 } 167 168 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, 169 State s, State ds, StateConservative dFlux[3]) { 170 for (CeedInt i=0; i<3; i++) { 171 dFlux[i].density = ds.U.momentum[i]; 172 for (CeedInt j=0; j<3; j++) 173 dFlux[i].momentum[j] = ds.U.momentum[i] * s.Y.velocity[j] + 174 s.U.momentum[i] * ds.Y.velocity[j] + ds.Y.pressure * (i == j); 175 dFlux[i].E_total = (ds.U.E_total + ds.Y.pressure) * s.Y.velocity[i] + 176 (s.U.E_total + s.Y.pressure) * ds.Y.velocity[i]; 177 } 178 } 179 180 // Kelvin-Mandel notation 181 CEED_QFUNCTION_HELPER void KMStrainRate(const State grad_s[3], 182 CeedScalar strain_rate[6]) { 183 const CeedScalar weight = 1 / sqrt(2.); 184 strain_rate[0] = grad_s[0].Y.velocity[0]; 185 strain_rate[1] = grad_s[1].Y.velocity[1]; 186 strain_rate[2] = grad_s[2].Y.velocity[2]; 187 strain_rate[3] = weight * (grad_s[2].Y.velocity[1] + grad_s[1].Y.velocity[2]); 188 strain_rate[4] = weight * (grad_s[2].Y.velocity[0] + grad_s[0].Y.velocity[2]); 189 strain_rate[5] = weight * (grad_s[1].Y.velocity[0] + grad_s[0].Y.velocity[1]); 190 } 191 192 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, 193 const CeedScalar strain_rate[6], CeedScalar stress[6]) { 194 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 195 for (CeedInt i=0; i<6; i++) { 196 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 197 } 198 } 199 200 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, 201 StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 202 CeedScalar Fe[3]) { 203 for (CeedInt i=0; i<3; i++) { 204 Fe[i] = - Y.velocity[0] * stress[0][i] 205 - Y.velocity[1] * stress[1][i] 206 - Y.velocity[2] * stress[2][i] 207 - gas->k * grad_s[i].Y.temperature; 208 } 209 } 210 211 CEED_QFUNCTION_HELPER void ViscousEnergyFlux_fwd(NewtonianIdealGasContext gas, 212 StatePrimitive Y, StatePrimitive dY, const State grad_ds[3], 213 const CeedScalar stress[3][3], 214 const CeedScalar dstress[3][3], 215 CeedScalar dFe[3]) { 216 for (CeedInt i=0; i<3; i++) { 217 dFe[i] = - Y.velocity[0] * dstress[0][i] - dY.velocity[0] * stress[0][i] 218 - Y.velocity[1] * dstress[1][i] - dY.velocity[1] * stress[1][i] 219 - Y.velocity[2] * dstress[2][i] - dY.velocity[2] * stress[2][i] 220 - gas->k * grad_ds[i].Y.temperature; 221 } 222 } 223 224 #endif // newtonian_state_h 225