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