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