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