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 density; 22 CeedScalar momentum[3]; 23 CeedScalar E_total; 24 } StateConservative; 25 26 typedef struct { 27 StateConservative U; 28 StatePrimitive Y; 29 } State; 30 31 CEED_QFUNCTION_HELPER void UnpackState_U(StateConservative s, CeedScalar U[5]) { 32 U[0] = s.density; 33 for (int i = 0; i < 3; i++) U[i + 1] = s.momentum[i]; 34 U[4] = s.E_total; 35 } 36 37 CEED_QFUNCTION_HELPER void UnpackState_Y(StatePrimitive s, CeedScalar Y[5]) { 38 Y[0] = s.pressure; 39 for (int i = 0; i < 3; i++) Y[i + 1] = s.velocity[i]; 40 Y[4] = s.temperature; 41 } 42 43 CEED_QFUNCTION_HELPER CeedScalar HeatCapacityRatio(NewtonianIdealGasContext gas) { return gas->cp / gas->cv; } 44 45 CEED_QFUNCTION_HELPER CeedScalar GasConstant(NewtonianIdealGasContext gas) { return gas->cp - gas->cv; } 46 47 CEED_QFUNCTION_HELPER CeedScalar Prandtl(NewtonianIdealGasContext gas) { return gas->cp * gas->mu / gas->k; } 48 49 CEED_QFUNCTION_HELPER CeedScalar SoundSpeed(NewtonianIdealGasContext gas, CeedScalar T) { return sqrt(gas->cp * (HeatCapacityRatio(gas) - 1.) * T); } 50 51 CEED_QFUNCTION_HELPER CeedScalar Mach(NewtonianIdealGasContext gas, CeedScalar T, CeedScalar u) { return u / SoundSpeed(gas, T); } 52 53 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy(NewtonianIdealGasContext gas, const State s) { 54 // Ignoring potential energy 55 CeedScalar e_internal = gas->cv * s.Y.temperature; 56 CeedScalar e_kinetic = 0.5 * Dot3(s.Y.velocity, s.Y.velocity); 57 return e_internal + e_kinetic + s.Y.pressure / s.U.density; 58 } 59 60 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy_fwd(NewtonianIdealGasContext gas, const State s, const State ds) { 61 // Ignoring potential energy 62 CeedScalar de_kinetic = Dot3(ds.Y.velocity, s.Y.velocity); 63 CeedScalar de_internal = gas->cv * ds.Y.temperature; 64 return de_internal + de_kinetic + ds.Y.pressure / s.U.density - s.Y.pressure / Square(s.U.density) * ds.U.density; 65 } 66 67 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative(NewtonianIdealGasContext gas, StateConservative U, const CeedScalar x[3]) { 68 StatePrimitive Y; 69 for (CeedInt i = 0; i < 3; i++) Y.velocity[i] = U.momentum[i] / U.density; 70 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 71 CeedScalar e_potential = -Dot3(gas->g, x); 72 CeedScalar e_total = U.E_total / U.density; 73 CeedScalar e_internal = e_total - e_kinetic - e_potential; 74 Y.temperature = e_internal / gas->cv; 75 Y.pressure = (HeatCapacityRatio(gas) - 1) * U.density * e_internal; 76 return Y; 77 } 78 79 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd(NewtonianIdealGasContext gas, State s, StateConservative dU, 80 const CeedScalar x[3], const CeedScalar dx[3]) { 81 StatePrimitive dY; 82 for (CeedInt i = 0; i < 3; i++) { 83 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 84 } 85 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 86 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 87 CeedScalar e_potential = -Dot3(gas->g, x); 88 CeedScalar de_potential = -Dot3(gas->g, dx); 89 CeedScalar e_total = s.U.E_total / s.U.density; 90 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 91 CeedScalar e_internal = e_total - e_kinetic - e_potential; 92 CeedScalar de_internal = de_total - de_kinetic - de_potential; 93 dY.temperature = de_internal / gas->cv; 94 dY.pressure = (HeatCapacityRatio(gas) - 1) * (dU.density * e_internal + s.U.density * de_internal); 95 return dY; 96 } 97 98 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive(NewtonianIdealGasContext gas, StatePrimitive Y, const CeedScalar x[3]) { 99 StateConservative U; 100 U.density = Y.pressure / (GasConstant(gas) * Y.temperature); 101 for (int i = 0; i < 3; i++) U.momentum[i] = U.density * Y.velocity[i]; 102 CeedScalar e_internal = gas->cv * Y.temperature; 103 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 104 CeedScalar e_potential = -Dot3(gas->g, x); 105 CeedScalar e_total = e_internal + e_kinetic + e_potential; 106 U.E_total = U.density * e_total; 107 return U; 108 } 109 110 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive_fwd(NewtonianIdealGasContext gas, State s, StatePrimitive dY, 111 const CeedScalar x[3], const CeedScalar dx[3]) { 112 StateConservative dU; 113 dU.density = (dY.pressure * s.Y.temperature - s.Y.pressure * dY.temperature) / (GasConstant(gas) * s.Y.temperature * s.Y.temperature); 114 for (int i = 0; i < 3; i++) { 115 dU.momentum[i] = dU.density * s.Y.velocity[i] + s.U.density * dY.velocity[i]; 116 } 117 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 118 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 119 CeedScalar e_potential = -Dot3(gas->g, x); 120 CeedScalar de_potential = -Dot3(gas->g, dx); 121 CeedScalar e_internal = gas->cv * s.Y.temperature; 122 CeedScalar de_internal = gas->cv * dY.temperature; 123 CeedScalar e_total = e_internal + e_kinetic + e_potential; 124 CeedScalar de_total = de_internal + de_kinetic + de_potential; 125 dU.E_total = dU.density * e_total + s.U.density * de_total; 126 return dU; 127 } 128 129 CEED_QFUNCTION_HELPER State StateFromPrimitive(NewtonianIdealGasContext gas, StatePrimitive Y, const CeedScalar x[3]) { 130 StateConservative U = StateConservativeFromPrimitive(gas, Y, x); 131 State s; 132 s.U = U; 133 s.Y = Y; 134 return s; 135 } 136 137 CEED_QFUNCTION_HELPER State StateFromPrimitive_fwd(NewtonianIdealGasContext gas, State s, StatePrimitive dY, const CeedScalar x[3], 138 const CeedScalar dx[3]) { 139 StateConservative dU = StateConservativeFromPrimitive_fwd(gas, s, dY, x, dx); 140 State ds; 141 ds.U = dU; 142 ds.Y = dY; 143 return ds; 144 } 145 146 // linear combination of n states 147 CEED_QFUNCTION_HELPER StateConservative StateConservativeMult(CeedInt n, const CeedScalar a[], const StateConservative X[]) { 148 StateConservative R = {0}; 149 for (CeedInt i = 0; i < n; i++) { 150 R.density += a[i] * X[i].density; 151 for (int j = 0; j < 3; j++) R.momentum[j] += a[i] * X[i].momentum[j]; 152 R.E_total += a[i] * X[i].E_total; 153 } 154 return R; 155 } 156 157 CEED_QFUNCTION_HELPER StateConservative StateConservativeAXPBYPCZ(CeedScalar a, StateConservative X, CeedScalar b, StateConservative Y, CeedScalar c, 158 StateConservative Z) { 159 StateConservative R; 160 R.density = a * X.density + b * Y.density + c * Z.density; 161 for (int i = 0; i < 3; i++) R.momentum[i] = a * X.momentum[i] + b * Y.momentum[i] + c * Z.momentum[i]; 162 R.E_total = a * X.E_total + b * Y.E_total + c * Z.E_total; 163 return R; 164 } 165 166 CEED_QFUNCTION_HELPER void StateToU(NewtonianIdealGasContext gas, const State input, CeedScalar U[5]) { UnpackState_U(input.U, U); } 167 168 CEED_QFUNCTION_HELPER void StateToY(NewtonianIdealGasContext gas, const State input, CeedScalar Y[5]) { UnpackState_Y(input.Y, Y); } 169 170 CEED_QFUNCTION_HELPER void StateToQ(NewtonianIdealGasContext gas, const State input, CeedScalar Q[5], StateVariable state_var) { 171 switch (state_var) { 172 case STATEVAR_CONSERVATIVE: 173 StateToU(gas, input, Q); 174 break; 175 case STATEVAR_PRIMITIVE: 176 StateToY(gas, input, Q); 177 break; 178 } 179 } 180 181 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, const CeedScalar U[5], const CeedScalar x[3]) { 182 State s; 183 s.U.density = U[0]; 184 s.U.momentum[0] = U[1]; 185 s.U.momentum[1] = U[2]; 186 s.U.momentum[2] = U[3]; 187 s.U.E_total = U[4]; 188 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 189 return s; 190 } 191 192 CEED_QFUNCTION_HELPER State StateFromU_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dU[5], const CeedScalar x[3], 193 const CeedScalar dx[3]) { 194 State ds; 195 ds.U.density = dU[0]; 196 ds.U.momentum[0] = dU[1]; 197 ds.U.momentum[1] = dU[2]; 198 ds.U.momentum[2] = dU[3]; 199 ds.U.E_total = dU[4]; 200 ds.Y = StatePrimitiveFromConservative_fwd(gas, s, ds.U, x, dx); 201 return ds; 202 } 203 204 CEED_QFUNCTION_HELPER State StateFromY(NewtonianIdealGasContext gas, const CeedScalar Y[5], const CeedScalar x[3]) { 205 State s; 206 s.Y.pressure = Y[0]; 207 s.Y.velocity[0] = Y[1]; 208 s.Y.velocity[1] = Y[2]; 209 s.Y.velocity[2] = Y[3]; 210 s.Y.temperature = Y[4]; 211 s.U = StateConservativeFromPrimitive(gas, s.Y, x); 212 return s; 213 } 214 215 CEED_QFUNCTION_HELPER State StateFromY_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dY[5], const CeedScalar x[3], 216 const CeedScalar dx[3]) { 217 State ds; 218 ds.Y.pressure = dY[0]; 219 ds.Y.velocity[0] = dY[1]; 220 ds.Y.velocity[1] = dY[2]; 221 ds.Y.velocity[2] = dY[3]; 222 ds.Y.temperature = dY[4]; 223 ds.U = StateConservativeFromPrimitive_fwd(gas, s, ds.Y, x, dx); 224 return ds; 225 } 226 227 CEED_QFUNCTION_HELPER State StateFromQ(NewtonianIdealGasContext gas, const CeedScalar Q[5], const CeedScalar x[3], StateVariable state_var) { 228 State s; 229 switch (state_var) { 230 case STATEVAR_CONSERVATIVE: 231 s = StateFromU(gas, Q, x); 232 break; 233 case STATEVAR_PRIMITIVE: 234 s = StateFromY(gas, Q, x); 235 break; 236 } 237 return s; 238 } 239 240 CEED_QFUNCTION_HELPER State StateFromQ_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dQ[5], const CeedScalar x[3], 241 const CeedScalar dx[3], StateVariable state_var) { 242 State ds; 243 switch (state_var) { 244 case STATEVAR_CONSERVATIVE: 245 ds = StateFromU_fwd(gas, s, dQ, x, dx); 246 break; 247 case STATEVAR_PRIMITIVE: 248 ds = StateFromY_fwd(gas, s, dQ, x, dx); 249 break; 250 } 251 return ds; 252 } 253 254 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, StateConservative Flux[3]) { 255 for (CeedInt i = 0; i < 3; i++) { 256 Flux[i].density = s.U.momentum[i]; 257 for (CeedInt j = 0; j < 3; j++) Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] + s.Y.pressure * (i == j); 258 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 259 } 260 } 261 262 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, State s, State ds, StateConservative dFlux[3]) { 263 for (CeedInt i = 0; i < 3; i++) { 264 dFlux[i].density = ds.U.momentum[i]; 265 for (CeedInt j = 0; j < 3; j++) { 266 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); 267 } 268 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]; 269 } 270 } 271 272 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal(NewtonianIdealGasContext gas, State s, const CeedScalar normal[3]) { 273 StateConservative Flux[3], Flux_dot_n = {0}; 274 FluxInviscid(gas, s, Flux); 275 for (CeedInt i = 0; i < 3; i++) { 276 Flux_dot_n.density += Flux[i].density * normal[i]; 277 for (CeedInt j = 0; j < 3; j++) Flux_dot_n.momentum[j] += Flux[i].momentum[j] * normal[i]; 278 Flux_dot_n.E_total += Flux[i].E_total * normal[i]; 279 } 280 return Flux_dot_n; 281 } 282 283 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal_fwd(NewtonianIdealGasContext gas, State s, State ds, const CeedScalar normal[3]) { 284 StateConservative dFlux[3], Flux_dot_n = {0}; 285 FluxInviscid_fwd(gas, s, ds, dFlux); 286 for (CeedInt i = 0; i < 3; i++) { 287 Flux_dot_n.density += dFlux[i].density * normal[i]; 288 for (CeedInt j = 0; j < 3; j++) Flux_dot_n.momentum[j] += dFlux[i].momentum[j] * normal[i]; 289 Flux_dot_n.E_total += dFlux[i].E_total * normal[i]; 290 } 291 return Flux_dot_n; 292 } 293 294 CEED_QFUNCTION_HELPER void FluxInviscidStrong(NewtonianIdealGasContext gas, State s, State ds[3], CeedScalar strong_conv[5]) { 295 for (CeedInt i = 0; i < 5; i++) strong_conv[i] = 0; 296 for (CeedInt i = 0; i < 3; i++) { 297 StateConservative dF[3]; 298 FluxInviscid_fwd(gas, s, ds[i], dF); 299 CeedScalar dF_i[5]; 300 UnpackState_U(dF[i], dF_i); 301 for (CeedInt j = 0; j < 5; j++) strong_conv[j] += dF_i[j]; 302 } 303 } 304 305 CEED_QFUNCTION_HELPER void FluxTotal(const StateConservative F_inviscid[3], CeedScalar stress[3][3], CeedScalar Fe[3], CeedScalar Flux[5][3]) { 306 for (CeedInt j = 0; j < 3; j++) { 307 Flux[0][j] = F_inviscid[j].density; 308 for (CeedInt k = 0; k < 3; k++) Flux[k + 1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 309 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 310 } 311 } 312 313 CEED_QFUNCTION_HELPER void FluxTotal_Boundary(const StateConservative F_inviscid[3], const CeedScalar stress[3][3], const CeedScalar Fe[3], 314 const CeedScalar normal[3], CeedScalar Flux[5]) { 315 for (CeedInt j = 0; j < 5; j++) Flux[j] = 0.; 316 for (CeedInt j = 0; j < 3; j++) { 317 Flux[0] += F_inviscid[j].density * normal[j]; 318 for (CeedInt k = 0; k < 3; k++) { 319 Flux[k + 1] += (F_inviscid[j].momentum[k] - stress[k][j]) * normal[j]; 320 } 321 Flux[4] += (F_inviscid[j].E_total + Fe[j]) * normal[j]; 322 } 323 } 324 325 CEED_QFUNCTION_HELPER void FluxTotal_RiemannBoundary(const StateConservative F_inviscid_normal, const CeedScalar stress[3][3], const CeedScalar Fe[3], 326 const CeedScalar normal[3], CeedScalar Flux[5]) { 327 Flux[0] = F_inviscid_normal.density; 328 for (CeedInt k = 0; k < 3; k++) Flux[k + 1] = F_inviscid_normal.momentum[k]; 329 Flux[4] = F_inviscid_normal.E_total; 330 for (CeedInt j = 0; j < 3; j++) { 331 for (CeedInt k = 0; k < 3; k++) { 332 Flux[k + 1] -= stress[k][j] * normal[j]; 333 } 334 Flux[4] += Fe[j] * normal[j]; 335 } 336 } 337 338 CEED_QFUNCTION_HELPER void VelocityGradient(const State grad_s[3], CeedScalar grad_velocity[3][3]) { 339 grad_velocity[0][0] = grad_s[0].Y.velocity[0]; 340 grad_velocity[0][1] = grad_s[1].Y.velocity[0]; 341 grad_velocity[0][2] = grad_s[2].Y.velocity[0]; 342 grad_velocity[1][0] = grad_s[0].Y.velocity[1]; 343 grad_velocity[1][1] = grad_s[1].Y.velocity[1]; 344 grad_velocity[1][2] = grad_s[2].Y.velocity[1]; 345 grad_velocity[2][0] = grad_s[0].Y.velocity[2]; 346 grad_velocity[2][1] = grad_s[1].Y.velocity[2]; 347 grad_velocity[2][2] = grad_s[2].Y.velocity[2]; 348 } 349 350 CEED_QFUNCTION_HELPER void KMStrainRate(const CeedScalar grad_velocity[3][3], CeedScalar strain_rate[6]) { 351 const CeedScalar weight = 1 / sqrt(2.); // Really sqrt(2.) / 2 352 strain_rate[0] = grad_velocity[0][0]; 353 strain_rate[1] = grad_velocity[1][1]; 354 strain_rate[2] = grad_velocity[2][2]; 355 strain_rate[3] = weight * (grad_velocity[1][2] + grad_velocity[2][1]); 356 strain_rate[4] = weight * (grad_velocity[0][2] + grad_velocity[2][0]); 357 strain_rate[5] = weight * (grad_velocity[0][1] + grad_velocity[1][0]); 358 } 359 360 // Kelvin-Mandel notation 361 CEED_QFUNCTION_HELPER void KMStrainRate_State(const State grad_s[3], CeedScalar strain_rate[6]) { 362 CeedScalar grad_velocity[3][3]; 363 VelocityGradient(grad_s, grad_velocity); 364 KMStrainRate(grad_velocity, strain_rate); 365 } 366 367 CEED_QFUNCTION_HELPER void RotationRate(const CeedScalar grad_velocity[3][3], CeedScalar rotation_rate[3][3]) { 368 rotation_rate[0][0] = 0; 369 rotation_rate[1][1] = 0; 370 rotation_rate[2][2] = 0; 371 rotation_rate[1][2] = 0.5 * (grad_velocity[1][2] - grad_velocity[2][1]); 372 rotation_rate[0][2] = 0.5 * (grad_velocity[0][2] - grad_velocity[2][0]); 373 rotation_rate[0][1] = 0.5 * (grad_velocity[0][1] - grad_velocity[1][0]); 374 rotation_rate[2][1] = -rotation_rate[1][2]; 375 rotation_rate[2][0] = -rotation_rate[0][2]; 376 rotation_rate[1][0] = -rotation_rate[0][1]; 377 } 378 379 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, const CeedScalar strain_rate[6], CeedScalar stress[6]) { 380 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 381 for (CeedInt i = 0; i < 6; i++) { 382 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 383 } 384 } 385 386 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 387 CeedScalar Fe[3]) { 388 for (CeedInt i = 0; i < 3; i++) { 389 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; 390 } 391 } 392 393 CEED_QFUNCTION_HELPER void ViscousEnergyFlux_fwd(NewtonianIdealGasContext gas, StatePrimitive Y, StatePrimitive dY, const State grad_ds[3], 394 const CeedScalar stress[3][3], const CeedScalar dstress[3][3], CeedScalar dFe[3]) { 395 for (CeedInt i = 0; i < 3; i++) { 396 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] - 397 Y.velocity[2] * dstress[2][i] - dY.velocity[2] * stress[2][i] - gas->k * grad_ds[i].Y.temperature; 398 } 399 } 400 401 CEED_QFUNCTION_HELPER void Vorticity(const State grad_s[3], CeedScalar vorticity[3]) { 402 CeedScalar grad_velocity[3][3]; 403 VelocityGradient(grad_s, grad_velocity); 404 Curl3(grad_velocity, vorticity); 405 } 406 407 #endif // newtonian_state_h 408