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 StateEntropy V; 28 } State; 29 30 CEED_QFUNCTION_HELPER void UnpackState_U(StateConservative s, CeedScalar U[5]) { 31 U[0] = s.density; 32 for (int i = 0; i < 3; i++) U[i + 1] = s.momentum[i]; 33 U[4] = s.E_total; 34 } 35 36 CEED_QFUNCTION_HELPER void UnpackState_Y(StatePrimitive s, CeedScalar Y[5]) { 37 Y[0] = s.pressure; 38 for (int i = 0; i < 3; i++) Y[i + 1] = s.velocity[i]; 39 Y[4] = s.temperature; 40 } 41 42 CEED_QFUNCTION_HELPER void UnpackState_V(StateEntropy s, CeedScalar V[5]) { 43 V[0] = s.S_density; 44 for (int i = 0; i < 3; i++) V[i + 1] = s.S_momentum[i]; 45 V[4] = s.S_energy; 46 } 47 48 CEED_QFUNCTION_HELPER CeedScalar HeatCapacityRatio(NewtonianIdealGasContext gas) { return gas->cp / gas->cv; } 49 50 CEED_QFUNCTION_HELPER CeedScalar GasConstant(NewtonianIdealGasContext gas) { return gas->cp - gas->cv; } 51 52 CEED_QFUNCTION_HELPER CeedScalar Prandtl(NewtonianIdealGasContext gas) { return gas->cp * gas->mu / gas->k; } 53 54 CEED_QFUNCTION_HELPER CeedScalar SoundSpeed(NewtonianIdealGasContext gas, CeedScalar T) { return sqrt(gas->cp * (HeatCapacityRatio(gas) - 1.) * T); } 55 56 CEED_QFUNCTION_HELPER CeedScalar Mach(NewtonianIdealGasContext gas, CeedScalar T, CeedScalar u) { return u / SoundSpeed(gas, T); } 57 58 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy(NewtonianIdealGasContext gas, const State s) { 59 // Ignoring potential energy 60 CeedScalar e_internal = gas->cv * s.Y.temperature; 61 CeedScalar e_kinetic = 0.5 * Dot3(s.Y.velocity, s.Y.velocity); 62 return e_internal + e_kinetic + s.Y.pressure / s.U.density; 63 } 64 65 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy_fwd(NewtonianIdealGasContext gas, const State s, const State ds) { 66 // Ignoring potential energy 67 CeedScalar de_kinetic = Dot3(ds.Y.velocity, s.Y.velocity); 68 CeedScalar de_internal = gas->cv * ds.Y.temperature; 69 return de_internal + de_kinetic + ds.Y.pressure / s.U.density - s.Y.pressure / Square(s.U.density) * ds.U.density; 70 } 71 72 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative(NewtonianIdealGasContext gas, StateConservative U) { 73 StatePrimitive Y; 74 for (CeedInt i = 0; i < 3; i++) Y.velocity[i] = U.momentum[i] / U.density; 75 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 76 CeedScalar e_total = U.E_total / U.density; 77 CeedScalar e_internal = e_total - e_kinetic; 78 Y.temperature = e_internal / gas->cv; 79 Y.pressure = (HeatCapacityRatio(gas) - 1) * U.density * e_internal; 80 return Y; 81 } 82 83 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd(NewtonianIdealGasContext gas, State s, StateConservative dU) { 84 StatePrimitive dY; 85 for (CeedInt i = 0; i < 3; i++) { 86 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 87 } 88 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 89 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 90 CeedScalar e_total = s.U.E_total / s.U.density; 91 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 92 CeedScalar e_internal = e_total - e_kinetic; 93 CeedScalar de_internal = de_total - de_kinetic; 94 dY.temperature = de_internal / gas->cv; 95 dY.pressure = (HeatCapacityRatio(gas) - 1) * (dU.density * e_internal + s.U.density * de_internal); 96 return dY; 97 } 98 99 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive(NewtonianIdealGasContext gas, StatePrimitive Y) { 100 StateConservative U; 101 U.density = Y.pressure / (GasConstant(gas) * Y.temperature); 102 for (int i = 0; i < 3; i++) U.momentum[i] = U.density * Y.velocity[i]; 103 CeedScalar e_internal = gas->cv * Y.temperature; 104 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 105 CeedScalar e_total = e_internal + e_kinetic; 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 StateConservative dU; 112 dU.density = (dY.pressure * s.Y.temperature - s.Y.pressure * dY.temperature) / (GasConstant(gas) * s.Y.temperature * s.Y.temperature); 113 for (int i = 0; i < 3; i++) { 114 dU.momentum[i] = dU.density * s.Y.velocity[i] + s.U.density * dY.velocity[i]; 115 } 116 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 117 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 118 CeedScalar e_internal = gas->cv * s.Y.temperature; 119 CeedScalar de_internal = gas->cv * dY.temperature; 120 CeedScalar e_total = e_internal + e_kinetic; 121 CeedScalar de_total = de_internal + de_kinetic; 122 dU.E_total = dU.density * e_total + s.U.density * de_total; 123 return dU; 124 } 125 126 CEED_QFUNCTION_HELPER StateEntropy StateEntropyFromConservative(NewtonianIdealGasContext gas, StateConservative U) { 127 StateEntropy V; 128 const CeedScalar gamma = HeatCapacityRatio(gas); 129 const CeedScalar e_kinetic = .5 * Dot3(U.momentum, U.momentum) / U.density; 130 const CeedScalar e_internal = U.E_total - e_kinetic; 131 const CeedScalar p = (gamma - 1) * e_internal; 132 const CeedScalar entropy = log(p) - gamma * log(U.density); 133 134 V.S_density = (gamma - entropy) / (gamma - 1) - e_kinetic / p; 135 for (int i = 0; i < 3; i++) V.S_momentum[i] = U.momentum[i] / p; 136 V.S_energy = -U.density / p; 137 return V; 138 } 139 140 CEED_QFUNCTION_HELPER StateEntropy StateEntropyFromConservative_fwd(NewtonianIdealGasContext gas, State s, StateConservative dU) { 141 StateEntropy dV; 142 const CeedScalar gamma = HeatCapacityRatio(gas); 143 const CeedScalar e_kinetic = .5 * Dot3(s.U.momentum, s.U.momentum) / s.U.density; 144 const CeedScalar de_kinetic = (Dot3(s.U.momentum, dU.momentum) - e_kinetic * dU.density) / s.U.density; 145 const CeedScalar de_internal = dU.E_total - de_kinetic; 146 const CeedScalar p = s.Y.pressure; 147 const CeedScalar dp = (gamma - 1) * de_internal; 148 149 CeedScalar dentropy = dp / p - gamma * dU.density / s.U.density; 150 151 dV.S_density = -dentropy / (gamma - 1) - de_kinetic / p + dp * e_kinetic / Square(p); 152 for (CeedInt i = 0; i < 3; i++) { 153 dV.S_momentum[i] = (dU.momentum[i] - s.U.momentum[i] * dp / p) / p; 154 } 155 dV.S_energy = -(dU.density - s.U.density * dp / p) / p; 156 return dV; 157 } 158 159 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromEntropy(NewtonianIdealGasContext gas, StateEntropy V) { 160 StateConservative U; 161 CeedScalar velocity[3]; 162 for (int i = 0; i < 3; i++) velocity[i] = -V.S_momentum[i] / V.S_energy; 163 const CeedScalar gamma = HeatCapacityRatio(gas); 164 const CeedScalar e_kinetic = 0.5 * Dot3(velocity, velocity); 165 const CeedScalar entropy = gamma - (gamma - 1) * (V.S_density - e_kinetic * V.S_energy); 166 const CeedScalar log_rho = -(entropy + log(-V.S_energy)) / (gamma - 1); 167 U.density = exp(log_rho); 168 for (int i = 0; i < 3; i++) U.momentum[i] = U.density * velocity[i]; 169 170 const CeedScalar e_internal = -gas->cv / (GasConstant(gas) * V.S_energy); 171 U.E_total = U.density * (e_internal + e_kinetic); 172 return U; 173 } 174 175 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromEntropy_fwd(NewtonianIdealGasContext gas, State s, StateEntropy dV) { 176 StateConservative dU; 177 CeedScalar dvelocity[3]; 178 for (int i = 0; i < 3; i++) dvelocity[i] = (-dV.S_momentum[i] - s.Y.velocity[i] * dV.S_energy) / s.V.S_energy; 179 const CeedScalar gamma = HeatCapacityRatio(gas); 180 const CeedScalar e_kinetic = 0.5 * Dot3(s.Y.velocity, s.Y.velocity); 181 const CeedScalar de_kinetic = Dot3(dvelocity, s.Y.velocity); 182 const CeedScalar entropy = gamma - (gamma - 1) * (s.V.S_density - e_kinetic * s.V.S_energy); 183 const CeedScalar dentropy = -(gamma - 1) * (dV.S_density - (de_kinetic * s.V.S_energy + e_kinetic * dV.S_energy)); 184 const CeedScalar log_rho = -(entropy + log(-s.V.S_energy)) / (gamma - 1); 185 const CeedScalar rho = exp(log_rho); 186 dU.density = -rho / (gamma - 1) * (dentropy + dV.S_energy / s.V.S_energy); 187 for (int i = 0; i < 3; i++) dU.momentum[i] = dU.density * s.Y.velocity[i] + s.U.density * dvelocity[i]; 188 189 const CeedScalar e_internal = -gas->cv / (GasConstant(gas) * s.V.S_energy); 190 const CeedScalar de_internal = gas->cv * dV.S_energy / (GasConstant(gas) * s.V.S_energy * s.V.S_energy); 191 const CeedScalar e_total = e_internal + e_kinetic; 192 dU.E_total = dU.density * e_total + s.U.density * (de_internal + de_kinetic); 193 return dU; 194 } 195 196 CEED_QFUNCTION_HELPER StateEntropy StateEntropyFromPrimitive(NewtonianIdealGasContext gas, StatePrimitive Y) { 197 StateEntropy V; 198 const CeedScalar gamma = HeatCapacityRatio(gas); 199 const CeedScalar rho = Y.pressure / (GasConstant(gas) * Y.temperature); 200 const CeedScalar entropy = log(Y.pressure) - gamma * log(rho); 201 const CeedScalar rho_div_p = rho / Y.pressure; 202 const CeedScalar e_kinetic = 0.5 * Dot3(Y.velocity, Y.velocity); 203 204 V.S_density = (gamma - entropy) / (gamma - 1) - rho_div_p * e_kinetic; 205 for (int i = 0; i < 3; i++) V.S_momentum[i] = rho_div_p * Y.velocity[i]; 206 V.S_energy = -rho_div_p; 207 return V; 208 } 209 210 CEED_QFUNCTION_HELPER StateEntropy StateEntropyFromPrimitive_fwd(NewtonianIdealGasContext gas, State s, StatePrimitive dY) { 211 StateEntropy dV; 212 const CeedScalar gamma = HeatCapacityRatio(gas); 213 CeedScalar drho = (dY.pressure * s.Y.temperature - s.Y.pressure * dY.temperature) / (GasConstant(gas) * s.Y.temperature * s.Y.temperature); 214 215 const CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 216 const CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 217 const CeedScalar rho_div_p = s.U.density / s.Y.pressure; 218 const CeedScalar drho_div_p = (drho * s.Y.pressure - s.U.density * dY.pressure) / Square(s.Y.pressure); 219 220 CeedScalar dentropy = dY.pressure / s.Y.pressure - gamma * drho / s.U.density; 221 222 dV.S_density = -dentropy / (gamma - 1) - de_kinetic * rho_div_p - e_kinetic * drho_div_p; 223 for (CeedInt i = 0; i < 3; i++) dV.S_momentum[i] = rho_div_p * dY.velocity[i] + drho_div_p * s.Y.velocity[i]; 224 dV.S_energy = -drho_div_p; 225 return dV; 226 } 227 228 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromEntropy(NewtonianIdealGasContext gas, StateEntropy V) { 229 StatePrimitive Y; 230 for (int i = 0; i < 3; i++) Y.velocity[i] = -V.S_momentum[i] / V.S_energy; 231 Y.temperature = -1 / (GasConstant(gas) * V.S_energy); 232 const CeedScalar gamma = HeatCapacityRatio(gas); 233 const CeedScalar e_kinetic = 0.5 * Dot3(Y.velocity, Y.velocity); 234 const CeedScalar entropy = gamma - (gamma - 1) * (V.S_density - e_kinetic * V.S_energy); 235 const CeedScalar log_P = -(entropy + gamma * log(-V.S_energy)) / (gamma - 1); 236 Y.pressure = exp(log_P); 237 return Y; 238 } 239 240 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromEntropy_fwd(NewtonianIdealGasContext gas, State s, StateEntropy dV) { 241 StatePrimitive dY; 242 for (int i = 0; i < 3; i++) dY.velocity[i] = -(dV.S_momentum[i] - s.V.S_momentum[i] * dV.S_energy / s.V.S_energy) / s.V.S_energy; 243 dY.temperature = dV.S_energy / (GasConstant(gas) * s.V.S_energy * s.V.S_energy); 244 const CeedScalar gamma = HeatCapacityRatio(gas); 245 const CeedScalar e_kinetic = 0.5 * Dot3(s.Y.velocity, s.Y.velocity); 246 const CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 247 const CeedScalar dentropy = (1 - gamma) * (dV.S_density - e_kinetic * dV.S_energy - de_kinetic * s.V.S_energy); 248 dY.pressure = s.Y.pressure * (-dentropy - gamma * dV.S_energy / s.V.S_energy) / (gamma - 1); 249 return dY; 250 } 251 252 CEED_QFUNCTION_HELPER State StateFromPrimitive(NewtonianIdealGasContext gas, StatePrimitive Y) { 253 StateConservative U = StateConservativeFromPrimitive(gas, Y); 254 StateEntropy V = StateEntropyFromPrimitive(gas, Y); 255 State s; 256 s.U = U; 257 s.Y = Y; 258 s.V = V; 259 return s; 260 } 261 262 CEED_QFUNCTION_HELPER State StateFromPrimitive_fwd(NewtonianIdealGasContext gas, State s, StatePrimitive dY) { 263 StateConservative dU = StateConservativeFromPrimitive_fwd(gas, s, dY); 264 StateEntropy dV = StateEntropyFromPrimitive_fwd(gas, s, dY); 265 State ds; 266 ds.U = dU; 267 ds.Y = dY; 268 ds.V = dV; 269 return ds; 270 } 271 272 // linear combination of n states 273 CEED_QFUNCTION_HELPER StateConservative StateConservativeMult(CeedInt n, const CeedScalar a[], const StateConservative X[]) { 274 StateConservative R = {0}; 275 for (CeedInt i = 0; i < n; i++) { 276 R.density += a[i] * X[i].density; 277 for (int j = 0; j < 3; j++) R.momentum[j] += a[i] * X[i].momentum[j]; 278 R.E_total += a[i] * X[i].E_total; 279 } 280 return R; 281 } 282 283 CEED_QFUNCTION_HELPER StateConservative StateConservativeAXPBYPCZ(CeedScalar a, StateConservative X, CeedScalar b, StateConservative Y, CeedScalar c, 284 StateConservative Z) { 285 StateConservative R; 286 R.density = a * X.density + b * Y.density + c * Z.density; 287 for (int i = 0; i < 3; i++) R.momentum[i] = a * X.momentum[i] + b * Y.momentum[i] + c * Z.momentum[i]; 288 R.E_total = a * X.E_total + b * Y.E_total + c * Z.E_total; 289 return R; 290 } 291 292 CEED_QFUNCTION_HELPER void StateToU(NewtonianIdealGasContext gas, const State input, CeedScalar U[5]) { UnpackState_U(input.U, U); } 293 294 CEED_QFUNCTION_HELPER void StateToY(NewtonianIdealGasContext gas, const State input, CeedScalar Y[5]) { UnpackState_Y(input.Y, Y); } 295 296 CEED_QFUNCTION_HELPER void StateToV(NewtonianIdealGasContext gas, const State input, CeedScalar V[5]) { UnpackState_V(input.V, V); } 297 298 CEED_QFUNCTION_HELPER void StateToQ(NewtonianIdealGasContext gas, const State input, CeedScalar Q[5], StateVariable state_var) { 299 switch (state_var) { 300 case STATEVAR_CONSERVATIVE: 301 StateToU(gas, input, Q); 302 break; 303 case STATEVAR_PRIMITIVE: 304 StateToY(gas, input, Q); 305 break; 306 case STATEVAR_ENTROPY: 307 StateToV(gas, input, Q); 308 break; 309 } 310 } 311 312 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, const CeedScalar U[5]) { 313 State s; 314 s.U.density = U[0]; 315 s.U.momentum[0] = U[1]; 316 s.U.momentum[1] = U[2]; 317 s.U.momentum[2] = U[3]; 318 s.U.E_total = U[4]; 319 s.Y = StatePrimitiveFromConservative(gas, s.U); 320 s.V = StateEntropyFromConservative(gas, s.U); 321 return s; 322 } 323 324 CEED_QFUNCTION_HELPER State StateFromU_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dU[5]) { 325 State ds; 326 ds.U.density = dU[0]; 327 ds.U.momentum[0] = dU[1]; 328 ds.U.momentum[1] = dU[2]; 329 ds.U.momentum[2] = dU[3]; 330 ds.U.E_total = dU[4]; 331 ds.Y = StatePrimitiveFromConservative_fwd(gas, s, ds.U); 332 ds.V = StateEntropyFromConservative_fwd(gas, s, ds.U); 333 return ds; 334 } 335 336 CEED_QFUNCTION_HELPER State StateFromY(NewtonianIdealGasContext gas, const CeedScalar Y[5]) { 337 State s; 338 s.Y.pressure = Y[0]; 339 s.Y.velocity[0] = Y[1]; 340 s.Y.velocity[1] = Y[2]; 341 s.Y.velocity[2] = Y[3]; 342 s.Y.temperature = Y[4]; 343 s.U = StateConservativeFromPrimitive(gas, s.Y); 344 s.V = StateEntropyFromPrimitive(gas, s.Y); 345 return s; 346 } 347 348 CEED_QFUNCTION_HELPER State StateFromY_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dY[5]) { 349 State ds; 350 ds.Y.pressure = dY[0]; 351 ds.Y.velocity[0] = dY[1]; 352 ds.Y.velocity[1] = dY[2]; 353 ds.Y.velocity[2] = dY[3]; 354 ds.Y.temperature = dY[4]; 355 ds.U = StateConservativeFromPrimitive_fwd(gas, s, ds.Y); 356 ds.V = StateEntropyFromPrimitive_fwd(gas, s, ds.Y); 357 return ds; 358 } 359 360 CEED_QFUNCTION_HELPER State StateFromV(NewtonianIdealGasContext gas, const CeedScalar V[5]) { 361 State s; 362 s.V.S_density = V[0]; 363 s.V.S_momentum[0] = V[1]; 364 s.V.S_momentum[1] = V[2]; 365 s.V.S_momentum[2] = V[3]; 366 s.V.S_energy = V[4]; 367 s.U = StateConservativeFromEntropy(gas, s.V); 368 s.Y = StatePrimitiveFromEntropy(gas, s.V); 369 return s; 370 } 371 372 CEED_QFUNCTION_HELPER State StateFromV_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dV[5]) { 373 State ds; 374 ds.V.S_density = dV[0]; 375 ds.V.S_momentum[0] = dV[1]; 376 ds.V.S_momentum[1] = dV[2]; 377 ds.V.S_momentum[2] = dV[3]; 378 ds.V.S_energy = dV[4]; 379 ds.U = StateConservativeFromEntropy_fwd(gas, s, ds.V); 380 ds.Y = StatePrimitiveFromEntropy_fwd(gas, s, ds.V); 381 return ds; 382 } 383 384 CEED_QFUNCTION_HELPER State StateFromQ(NewtonianIdealGasContext gas, const CeedScalar Q[5], StateVariable state_var) { 385 State s; 386 switch (state_var) { 387 case STATEVAR_CONSERVATIVE: 388 s = StateFromU(gas, Q); 389 break; 390 case STATEVAR_PRIMITIVE: 391 s = StateFromY(gas, Q); 392 break; 393 case STATEVAR_ENTROPY: 394 s = StateFromV(gas, Q); 395 break; 396 } 397 return s; 398 } 399 400 CEED_QFUNCTION_HELPER State StateFromQ_fwd(NewtonianIdealGasContext gas, State s, const CeedScalar dQ[5], StateVariable state_var) { 401 State ds; 402 switch (state_var) { 403 case STATEVAR_CONSERVATIVE: 404 ds = StateFromU_fwd(gas, s, dQ); 405 break; 406 case STATEVAR_PRIMITIVE: 407 ds = StateFromY_fwd(gas, s, dQ); 408 break; 409 case STATEVAR_ENTROPY: 410 ds = StateFromV_fwd(gas, s, dQ); 411 break; 412 } 413 return ds; 414 } 415 416 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, StateConservative Flux[3]) { 417 for (CeedInt i = 0; i < 3; i++) { 418 Flux[i].density = s.U.momentum[i]; 419 for (CeedInt j = 0; j < 3; j++) Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] + s.Y.pressure * (i == j); 420 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 421 } 422 } 423 424 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, State s, State ds, StateConservative dFlux[3]) { 425 for (CeedInt i = 0; i < 3; i++) { 426 dFlux[i].density = ds.U.momentum[i]; 427 for (CeedInt j = 0; j < 3; j++) { 428 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); 429 } 430 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]; 431 } 432 } 433 434 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal(NewtonianIdealGasContext gas, State s, const CeedScalar normal[3]) { 435 StateConservative Flux[3], Flux_dot_n = {0}; 436 FluxInviscid(gas, s, Flux); 437 for (CeedInt i = 0; i < 3; i++) { 438 Flux_dot_n.density += Flux[i].density * normal[i]; 439 for (CeedInt j = 0; j < 3; j++) Flux_dot_n.momentum[j] += Flux[i].momentum[j] * normal[i]; 440 Flux_dot_n.E_total += Flux[i].E_total * normal[i]; 441 } 442 return Flux_dot_n; 443 } 444 445 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal_fwd(NewtonianIdealGasContext gas, State s, State ds, const CeedScalar normal[3]) { 446 StateConservative dFlux[3], Flux_dot_n = {0}; 447 FluxInviscid_fwd(gas, s, ds, dFlux); 448 for (CeedInt i = 0; i < 3; i++) { 449 Flux_dot_n.density += dFlux[i].density * normal[i]; 450 for (CeedInt j = 0; j < 3; j++) Flux_dot_n.momentum[j] += dFlux[i].momentum[j] * normal[i]; 451 Flux_dot_n.E_total += dFlux[i].E_total * normal[i]; 452 } 453 return Flux_dot_n; 454 } 455 456 CEED_QFUNCTION_HELPER void FluxInviscidStrong(NewtonianIdealGasContext gas, State s, State ds[3], CeedScalar strong_conv[5]) { 457 for (CeedInt i = 0; i < 5; i++) strong_conv[i] = 0; 458 for (CeedInt i = 0; i < 3; i++) { 459 StateConservative dF[3]; 460 FluxInviscid_fwd(gas, s, ds[i], dF); 461 CeedScalar dF_i[5]; 462 UnpackState_U(dF[i], dF_i); 463 for (CeedInt j = 0; j < 5; j++) strong_conv[j] += dF_i[j]; 464 } 465 } 466 467 CEED_QFUNCTION_HELPER void FluxTotal(const StateConservative F_inviscid[3], CeedScalar stress[3][3], CeedScalar Fe[3], CeedScalar Flux[5][3]) { 468 for (CeedInt j = 0; j < 3; j++) { 469 Flux[0][j] = F_inviscid[j].density; 470 for (CeedInt k = 0; k < 3; k++) Flux[k + 1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 471 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 472 } 473 } 474 475 CEED_QFUNCTION_HELPER void FluxTotal_Boundary(const StateConservative F_inviscid[3], const CeedScalar stress[3][3], const CeedScalar Fe[3], 476 const CeedScalar normal[3], CeedScalar Flux[5]) { 477 for (CeedInt j = 0; j < 5; j++) Flux[j] = 0.; 478 for (CeedInt j = 0; j < 3; j++) { 479 Flux[0] += F_inviscid[j].density * normal[j]; 480 for (CeedInt k = 0; k < 3; k++) { 481 Flux[k + 1] += (F_inviscid[j].momentum[k] - stress[k][j]) * normal[j]; 482 } 483 Flux[4] += (F_inviscid[j].E_total + Fe[j]) * normal[j]; 484 } 485 } 486 487 CEED_QFUNCTION_HELPER void FluxTotal_RiemannBoundary(const StateConservative F_inviscid_normal, const CeedScalar stress[3][3], const CeedScalar Fe[3], 488 const CeedScalar normal[3], CeedScalar Flux[5]) { 489 Flux[0] = F_inviscid_normal.density; 490 for (CeedInt k = 0; k < 3; k++) Flux[k + 1] = F_inviscid_normal.momentum[k]; 491 Flux[4] = F_inviscid_normal.E_total; 492 for (CeedInt j = 0; j < 3; j++) { 493 for (CeedInt k = 0; k < 3; k++) { 494 Flux[k + 1] -= stress[k][j] * normal[j]; 495 } 496 Flux[4] += Fe[j] * normal[j]; 497 } 498 } 499 500 CEED_QFUNCTION_HELPER void VelocityGradient(const State grad_s[3], CeedScalar grad_velocity[3][3]) { 501 grad_velocity[0][0] = grad_s[0].Y.velocity[0]; 502 grad_velocity[0][1] = grad_s[1].Y.velocity[0]; 503 grad_velocity[0][2] = grad_s[2].Y.velocity[0]; 504 grad_velocity[1][0] = grad_s[0].Y.velocity[1]; 505 grad_velocity[1][1] = grad_s[1].Y.velocity[1]; 506 grad_velocity[1][2] = grad_s[2].Y.velocity[1]; 507 grad_velocity[2][0] = grad_s[0].Y.velocity[2]; 508 grad_velocity[2][1] = grad_s[1].Y.velocity[2]; 509 grad_velocity[2][2] = grad_s[2].Y.velocity[2]; 510 } 511 512 CEED_QFUNCTION_HELPER void KMStrainRate(const CeedScalar grad_velocity[3][3], CeedScalar strain_rate[6]) { 513 const CeedScalar weight = 1 / sqrt(2.); // Really sqrt(2.) / 2 514 strain_rate[0] = grad_velocity[0][0]; 515 strain_rate[1] = grad_velocity[1][1]; 516 strain_rate[2] = grad_velocity[2][2]; 517 strain_rate[3] = weight * (grad_velocity[1][2] + grad_velocity[2][1]); 518 strain_rate[4] = weight * (grad_velocity[0][2] + grad_velocity[2][0]); 519 strain_rate[5] = weight * (grad_velocity[0][1] + grad_velocity[1][0]); 520 } 521 522 // Kelvin-Mandel notation 523 CEED_QFUNCTION_HELPER void KMStrainRate_State(const State grad_s[3], CeedScalar strain_rate[6]) { 524 CeedScalar grad_velocity[3][3]; 525 VelocityGradient(grad_s, grad_velocity); 526 KMStrainRate(grad_velocity, strain_rate); 527 } 528 529 //@brief Given velocity gradient du_i/dx_j, return 0.5*(du_i/dx_j - du_j/dx_i) 530 CEED_QFUNCTION_HELPER void RotationRate(const CeedScalar grad_velocity[3][3], CeedScalar rotation_rate[3][3]) { 531 rotation_rate[0][0] = 0; 532 rotation_rate[1][1] = 0; 533 rotation_rate[2][2] = 0; 534 rotation_rate[1][2] = 0.5 * (grad_velocity[1][2] - grad_velocity[2][1]); 535 rotation_rate[0][2] = 0.5 * (grad_velocity[0][2] - grad_velocity[2][0]); 536 rotation_rate[0][1] = 0.5 * (grad_velocity[0][1] - grad_velocity[1][0]); 537 rotation_rate[2][1] = -rotation_rate[1][2]; 538 rotation_rate[2][0] = -rotation_rate[0][2]; 539 rotation_rate[1][0] = -rotation_rate[0][1]; 540 } 541 542 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, const CeedScalar strain_rate[6], CeedScalar stress[6]) { 543 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 544 for (CeedInt i = 0; i < 6; i++) { 545 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 546 } 547 } 548 549 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 550 CeedScalar Fe[3]) { 551 for (CeedInt i = 0; i < 3; i++) { 552 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; 553 } 554 } 555 556 CEED_QFUNCTION_HELPER void ViscousEnergyFlux_fwd(NewtonianIdealGasContext gas, StatePrimitive Y, StatePrimitive dY, const State grad_ds[3], 557 const CeedScalar stress[3][3], const CeedScalar dstress[3][3], CeedScalar dFe[3]) { 558 for (CeedInt i = 0; i < 3; i++) { 559 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] - 560 Y.velocity[2] * dstress[2][i] - dY.velocity[2] * stress[2][i] - gas->k * grad_ds[i].Y.temperature; 561 } 562 } 563 564 CEED_QFUNCTION_HELPER void Vorticity(const State grad_s[3], CeedScalar vorticity[3]) { 565 CeedScalar grad_velocity[3][3]; 566 VelocityGradient(grad_s, grad_velocity); 567 Curl3(grad_velocity, vorticity); 568 } 569 570 CEED_QFUNCTION_HELPER void StatePhysicalGradientFromReference(CeedInt Q, CeedInt i, NewtonianIdealGasContext gas, State s, StateVariable state_var, 571 const CeedScalar *grad_q, const CeedScalar dXdx[3][3], State grad_s[3]) { 572 for (CeedInt k = 0; k < 3; k++) { 573 CeedScalar dqi[5]; 574 for (CeedInt j = 0; j < 5; j++) { 575 dqi[j] = 576 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]; 577 } 578 grad_s[k] = StateFromQ_fwd(gas, s, dqi, state_var); 579 } 580 } 581 582 CEED_QFUNCTION_HELPER void StatePhysicalGradientFromReference_Boundary(CeedInt Q, CeedInt i, NewtonianIdealGasContext gas, State s, 583 StateVariable state_var, const CeedScalar *grad_q, const CeedScalar dXdx[2][3], 584 State grad_s[3]) { 585 for (CeedInt k = 0; k < 3; k++) { 586 CeedScalar dqi[5]; 587 for (CeedInt j = 0; j < 5; j++) { 588 dqi[j] = grad_q[(Q * 5) * 0 + Q * j + i] * dXdx[0][k] + grad_q[(Q * 5) * 1 + Q * j + i] * dXdx[1][k]; 589 } 590 grad_s[k] = StateFromQ_fwd(gas, s, dqi, state_var); 591 } 592 } 593