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 #ifndef newtonian_state_h 12 #define newtonian_state_h 13 14 #include <ceed.h> 15 #include <math.h> 16 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 void UnpackState_U(StateConservative s, CeedScalar U[5]) { 38 U[0] = s.density; 39 for (int i = 0; i < 3; i++) U[i + 1] = s.momentum[i]; 40 U[4] = s.E_total; 41 } 42 43 CEED_QFUNCTION_HELPER void UnpackState_Y(StatePrimitive s, CeedScalar Y[5]) { 44 Y[0] = s.pressure; 45 for (int i = 0; i < 3; i++) Y[i + 1] = s.velocity[i]; 46 Y[4] = s.temperature; 47 } 48 49 CEED_QFUNCTION_HELPER CeedScalar HeatCapacityRatio(NewtonianIdealGasContext gas) { return gas->cp / gas->cv; } 50 51 CEED_QFUNCTION_HELPER CeedScalar GasConstant(NewtonianIdealGasContext gas) { return gas->cp - gas->cv; } 52 53 CEED_QFUNCTION_HELPER CeedScalar Prandtl(NewtonianIdealGasContext gas) { return gas->cp * gas->mu / gas->k; } 54 55 CEED_QFUNCTION_HELPER CeedScalar SoundSpeed(NewtonianIdealGasContext gas, CeedScalar T) { return sqrt(gas->cp * (HeatCapacityRatio(gas) - 1.) * T); } 56 57 CEED_QFUNCTION_HELPER CeedScalar Mach(NewtonianIdealGasContext gas, CeedScalar T, CeedScalar u) { return u / SoundSpeed(gas, T); } 58 59 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy(NewtonianIdealGasContext gas, const State s) { 60 // Ignoring potential energy 61 CeedScalar e_internal = gas->cv * s.Y.temperature; 62 CeedScalar e_kinetic = 0.5 * Dot3(s.Y.velocity, s.Y.velocity); 63 return e_internal + e_kinetic + s.Y.pressure / s.U.density; 64 } 65 66 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy_fwd(NewtonianIdealGasContext gas, const State s, const State ds) { 67 // Ignoring potential energy 68 CeedScalar de_kinetic = Dot3(ds.Y.velocity, s.Y.velocity); 69 CeedScalar de_internal = gas->cv * ds.Y.temperature; 70 return de_internal + de_kinetic + ds.Y.pressure / s.U.density - s.Y.pressure / Square(s.U.density) * ds.U.density; 71 } 72 73 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative(NewtonianIdealGasContext gas, StateConservative U, const CeedScalar x[3]) { 74 StatePrimitive Y; 75 for (CeedInt i = 0; i < 3; i++) Y.velocity[i] = U.momentum[i] / U.density; 76 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 77 CeedScalar e_potential = -Dot3(gas->g, x); 78 CeedScalar e_total = U.E_total / U.density; 79 CeedScalar e_internal = e_total - e_kinetic - e_potential; 80 Y.temperature = e_internal / gas->cv; 81 Y.pressure = (HeatCapacityRatio(gas) - 1) * U.density * e_internal; 82 return Y; 83 } 84 85 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd(NewtonianIdealGasContext gas, State s, StateConservative dU, 86 const CeedScalar x[3], const CeedScalar dx[3]) { 87 StatePrimitive dY; 88 for (CeedInt i = 0; i < 3; i++) { 89 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 90 } 91 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 92 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 93 CeedScalar e_potential = -Dot3(gas->g, x); 94 CeedScalar de_potential = -Dot3(gas->g, dx); 95 CeedScalar e_total = s.U.E_total / s.U.density; 96 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 97 CeedScalar e_internal = e_total - e_kinetic - e_potential; 98 CeedScalar de_internal = de_total - de_kinetic - de_potential; 99 dY.temperature = de_internal / gas->cv; 100 dY.pressure = (HeatCapacityRatio(gas) - 1) * (dU.density * e_internal + s.U.density * de_internal); 101 return dY; 102 } 103 104 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive(NewtonianIdealGasContext gas, StatePrimitive Y, const CeedScalar x[3]) { 105 StateConservative U; 106 U.density = Y.pressure / (GasConstant(gas) * Y.temperature); 107 for (int i = 0; i < 3; i++) U.momentum[i] = U.density * Y.velocity[i]; 108 CeedScalar e_internal = gas->cv * Y.temperature; 109 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 110 CeedScalar e_potential = -Dot3(gas->g, x); 111 CeedScalar e_total = e_internal + e_kinetic + e_potential; 112 U.E_total = U.density * e_total; 113 return U; 114 } 115 116 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive_fwd(NewtonianIdealGasContext gas, State s, StatePrimitive dY, 117 const CeedScalar x[3], const CeedScalar dx[3]) { 118 StateConservative dU; 119 dU.density = (dY.pressure * s.Y.temperature - s.Y.pressure * dY.temperature) / (GasConstant(gas) * s.Y.temperature * s.Y.temperature); 120 for (int i = 0; i < 3; i++) { 121 dU.momentum[i] = dU.density * s.Y.velocity[i] + s.U.density * dY.velocity[i]; 122 } 123 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 124 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 125 CeedScalar e_potential = -Dot3(gas->g, x); 126 CeedScalar de_potential = -Dot3(gas->g, dx); 127 CeedScalar e_internal = gas->cv * s.Y.temperature; 128 CeedScalar de_internal = gas->cv * dY.temperature; 129 CeedScalar e_total = e_internal + e_kinetic + e_potential; 130 CeedScalar de_total = de_internal + de_kinetic + de_potential; 131 dU.E_total = dU.density * e_total + s.U.density * de_total; 132 return dU; 133 } 134 135 // linear combination of n states 136 CEED_QFUNCTION_HELPER StateConservative StateConservativeMult(CeedInt n, const CeedScalar a[], const StateConservative X[]) { 137 StateConservative R = {0}; 138 for (CeedInt i = 0; i < n; i++) { 139 R.density += a[i] * X[i].density; 140 for (int j = 0; j < 3; j++) R.momentum[j] += a[i] * X[i].momentum[j]; 141 R.E_total += a[i] * X[i].E_total; 142 } 143 return R; 144 } 145 146 CEED_QFUNCTION_HELPER StateConservative StateConservativeAXPBYPCZ(CeedScalar a, StateConservative X, CeedScalar b, StateConservative Y, CeedScalar c, 147 StateConservative Z) { 148 StateConservative R; 149 R.density = a * X.density + b * Y.density + c * Z.density; 150 for (int i = 0; i < 3; i++) R.momentum[i] = a * X.momentum[i] + b * Y.momentum[i] + c * Z.momentum[i]; 151 R.E_total = a * X.E_total + b * Y.E_total + c * Z.E_total; 152 return R; 153 } 154 155 // Function pointer types for generic state array -> State struct functions 156 typedef State (*StateFromQi_t)(NewtonianIdealGasContext gas, const CeedScalar qi[5], const CeedScalar x[3]); 157 typedef State (*StateFromQi_fwd_t)(NewtonianIdealGasContext gas, State s, const CeedScalar dqi[5], const CeedScalar x[3], const CeedScalar dx[3]); 158 159 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, const CeedScalar U[5], const CeedScalar x[3]) { 160 State s; 161 s.U.density = U[0]; 162 s.U.momentum[0] = U[1]; 163 s.U.momentum[1] = U[2]; 164 s.U.momentum[2] = U[3]; 165 s.U.E_total = U[4]; 166 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 167 return s; 168 } 169 170 CEED_QFUNCTION_HELPER State StateFromU_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dU[5], const CeedScalar x[3], 171 const CeedScalar dx[3]) { 172 State ds; 173 ds.U.density = dU[0]; 174 ds.U.momentum[0] = dU[1]; 175 ds.U.momentum[1] = dU[2]; 176 ds.U.momentum[2] = dU[3]; 177 ds.U.E_total = dU[4]; 178 ds.Y = StatePrimitiveFromConservative_fwd(gas, s, ds.U, x, dx); 179 return ds; 180 } 181 182 CEED_QFUNCTION_HELPER State StateFromY(NewtonianIdealGasContext gas, const CeedScalar Y[5], const CeedScalar x[3]) { 183 State s; 184 s.Y.pressure = Y[0]; 185 s.Y.velocity[0] = Y[1]; 186 s.Y.velocity[1] = Y[2]; 187 s.Y.velocity[2] = Y[3]; 188 s.Y.temperature = Y[4]; 189 s.U = StateConservativeFromPrimitive(gas, s.Y, x); 190 return s; 191 } 192 193 CEED_QFUNCTION_HELPER State StateFromY_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dY[5], const CeedScalar x[3], 194 const CeedScalar dx[3]) { 195 State ds; 196 ds.Y.pressure = dY[0]; 197 ds.Y.velocity[0] = dY[1]; 198 ds.Y.velocity[1] = dY[2]; 199 ds.Y.velocity[2] = dY[3]; 200 ds.Y.temperature = dY[4]; 201 ds.U = StateConservativeFromPrimitive_fwd(gas, s, ds.Y, x, dx); 202 return ds; 203 } 204 205 // Function pointer types for State struct -> generic state array 206 typedef void (*StateToQi_t)(NewtonianIdealGasContext gas, const State input, CeedScalar qi[5]); 207 208 CEED_QFUNCTION_HELPER void StateToU(NewtonianIdealGasContext gas, const State input, CeedScalar U[5]) { UnpackState_U(input.U, U); } 209 210 CEED_QFUNCTION_HELPER void StateToY(NewtonianIdealGasContext gas, const State input, CeedScalar Y[5]) { UnpackState_Y(input.Y, Y); } 211 212 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, StateConservative Flux[3]) { 213 for (CeedInt i = 0; i < 3; i++) { 214 Flux[i].density = s.U.momentum[i]; 215 for (CeedInt j = 0; j < 3; j++) Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] + s.Y.pressure * (i == j); 216 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 217 } 218 } 219 220 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, State s, State ds, StateConservative dFlux[3]) { 221 for (CeedInt i = 0; i < 3; i++) { 222 dFlux[i].density = ds.U.momentum[i]; 223 for (CeedInt j = 0; j < 3; j++) { 224 dFlux[i].momentum[j] = ds.U.momentum[i] * s.Y.velocity[j] + s.U.momentum[i] * ds.Y.velocity[j] + ds.Y.pressure * (i == j); 225 } 226 dFlux[i].E_total = (ds.U.E_total + ds.Y.pressure) * s.Y.velocity[i] + (s.U.E_total + s.Y.pressure) * ds.Y.velocity[i]; 227 } 228 } 229 230 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal(NewtonianIdealGasContext gas, State s, const CeedScalar normal[3]) { 231 StateConservative Flux[3], Flux_dot_n = {0}; 232 FluxInviscid(gas, s, Flux); 233 for (CeedInt i = 0; i < 3; i++) { 234 Flux_dot_n.density += Flux[i].density * normal[i]; 235 for (CeedInt j = 0; j < 3; j++) Flux_dot_n.momentum[j] += Flux[i].momentum[j] * normal[i]; 236 Flux_dot_n.E_total += Flux[i].E_total * normal[i]; 237 } 238 return Flux_dot_n; 239 } 240 241 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal_fwd(NewtonianIdealGasContext gas, State s, State ds, const CeedScalar normal[3]) { 242 StateConservative dFlux[3], Flux_dot_n = {0}; 243 FluxInviscid_fwd(gas, s, ds, dFlux); 244 for (CeedInt i = 0; i < 3; i++) { 245 Flux_dot_n.density += dFlux[i].density * normal[i]; 246 for (CeedInt j = 0; j < 3; j++) Flux_dot_n.momentum[j] += dFlux[i].momentum[j] * normal[i]; 247 Flux_dot_n.E_total += dFlux[i].E_total * normal[i]; 248 } 249 return Flux_dot_n; 250 } 251 252 CEED_QFUNCTION_HELPER void FluxInviscidStrong(NewtonianIdealGasContext gas, State s, State ds[3], CeedScalar strong_conv[5]) { 253 for (CeedInt i = 0; i < 5; i++) strong_conv[i] = 0; 254 for (CeedInt i = 0; i < 3; i++) { 255 StateConservative dF[3]; 256 FluxInviscid_fwd(gas, s, ds[i], dF); 257 CeedScalar dF_i[5]; 258 UnpackState_U(dF[i], dF_i); 259 for (CeedInt j = 0; j < 5; j++) strong_conv[j] += dF_i[j]; 260 } 261 } 262 263 CEED_QFUNCTION_HELPER void FluxTotal(const StateConservative F_inviscid[3], CeedScalar stress[3][3], CeedScalar Fe[3], CeedScalar Flux[5][3]) { 264 for (CeedInt j = 0; j < 3; j++) { 265 Flux[0][j] = F_inviscid[j].density; 266 for (CeedInt k = 0; k < 3; k++) Flux[k + 1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 267 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 268 } 269 } 270 271 CEED_QFUNCTION_HELPER void FluxTotal_Boundary(const StateConservative F_inviscid[3], const CeedScalar stress[3][3], const CeedScalar Fe[3], 272 const CeedScalar normal[3], CeedScalar Flux[5]) { 273 for (CeedInt j = 0; j < 5; j++) Flux[j] = 0.; 274 for (CeedInt j = 0; j < 3; j++) { 275 Flux[0] += F_inviscid[j].density * normal[j]; 276 for (CeedInt k = 0; k < 3; k++) { 277 Flux[k + 1] += (F_inviscid[j].momentum[k] - stress[k][j]) * normal[j]; 278 } 279 Flux[4] += (F_inviscid[j].E_total + Fe[j]) * normal[j]; 280 } 281 } 282 283 // Kelvin-Mandel notation 284 CEED_QFUNCTION_HELPER void KMStrainRate(const State grad_s[3], CeedScalar strain_rate[6]) { 285 const CeedScalar weight = 1 / sqrt(2.); 286 strain_rate[0] = grad_s[0].Y.velocity[0]; 287 strain_rate[1] = grad_s[1].Y.velocity[1]; 288 strain_rate[2] = grad_s[2].Y.velocity[2]; 289 strain_rate[3] = weight * (grad_s[2].Y.velocity[1] + grad_s[1].Y.velocity[2]); 290 strain_rate[4] = weight * (grad_s[2].Y.velocity[0] + grad_s[0].Y.velocity[2]); 291 strain_rate[5] = weight * (grad_s[1].Y.velocity[0] + grad_s[0].Y.velocity[1]); 292 } 293 294 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, const CeedScalar strain_rate[6], CeedScalar stress[6]) { 295 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 296 for (CeedInt i = 0; i < 6; i++) { 297 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 298 } 299 } 300 301 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 302 CeedScalar Fe[3]) { 303 for (CeedInt i = 0; i < 3; i++) { 304 Fe[i] = -Y.velocity[0] * stress[0][i] - Y.velocity[1] * stress[1][i] - Y.velocity[2] * stress[2][i] - gas->k * grad_s[i].Y.temperature; 305 } 306 } 307 308 CEED_QFUNCTION_HELPER void ViscousEnergyFlux_fwd(NewtonianIdealGasContext gas, StatePrimitive Y, StatePrimitive dY, const State grad_ds[3], 309 const CeedScalar stress[3][3], const CeedScalar dstress[3][3], CeedScalar dFe[3]) { 310 for (CeedInt i = 0; i < 3; i++) { 311 dFe[i] = -Y.velocity[0] * dstress[0][i] - dY.velocity[0] * stress[0][i] - Y.velocity[1] * dstress[1][i] - dY.velocity[1] * stress[1][i] - 312 Y.velocity[2] * dstress[2][i] - dY.velocity[2] * stress[2][i] - gas->k * grad_ds[i].Y.temperature; 313 } 314 } 315 316 #endif // newtonian_state_h 317