1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Structs and helper functions regarding the state of a newtonian simulation 10 11 12 #ifndef newtonian_state_h 13 #define newtonian_state_h 14 15 #include <ceed.h> 16 #include <math.h> 17 #include "newtonian_types.h" 18 #include "utils.h" 19 20 typedef struct { 21 CeedScalar pressure; 22 CeedScalar velocity[3]; 23 CeedScalar temperature; 24 } StatePrimitive; 25 26 typedef struct { 27 CeedScalar density; 28 CeedScalar momentum[3]; 29 CeedScalar E_total; 30 } StateConservative; 31 32 typedef struct { 33 StateConservative U; 34 StatePrimitive Y; 35 } State; 36 37 CEED_QFUNCTION_HELPER void UnpackState_U(StateConservative s, CeedScalar U[5]) { 38 U[0] = s.density; 39 for (int i=0; i<3; i++) U[i+1] = s.momentum[i]; 40 U[4] = s.E_total; 41 } 42 43 CEED_QFUNCTION_HELPER void UnpackState_Y(StatePrimitive s, CeedScalar Y[5]) { 44 Y[0] = s.pressure; 45 for (int i=0; i<3; i++) Y[i+1] = s.velocity[i]; 46 Y[4] = s.temperature; 47 } 48 49 CEED_QFUNCTION_HELPER CeedScalar HeatCapacityRatio( 50 NewtonianIdealGasContext gas) { 51 return gas->cp / gas->cv; 52 } 53 54 CEED_QFUNCTION_HELPER CeedScalar GasConstant( 55 NewtonianIdealGasContext gas) { 56 return gas->cp - gas->cv; 57 } 58 59 CEED_QFUNCTION_HELPER CeedScalar Prandtl(NewtonianIdealGasContext gas) { 60 return gas->cp * gas->mu / gas->k; 61 } 62 63 CEED_QFUNCTION_HELPER CeedScalar SoundSpeed(NewtonianIdealGasContext gas, 64 CeedScalar T) { 65 return sqrt(gas->cp * (HeatCapacityRatio(gas) - 1.) * T); 66 } 67 68 CEED_QFUNCTION_HELPER CeedScalar Mach(NewtonianIdealGasContext gas, 69 CeedScalar T, CeedScalar u) { 70 return u / SoundSpeed(gas, T); 71 } 72 73 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy( 74 NewtonianIdealGasContext gas, const State s) { 75 // Ignoring potential energy 76 CeedScalar e_internal = gas->cv*s.Y.temperature; 77 CeedScalar e_kinetic = 0.5*Dot3(s.Y.velocity, s.Y.velocity); 78 return e_internal + e_kinetic + s.Y.pressure/s.U.density; 79 } 80 81 CEED_QFUNCTION_HELPER CeedScalar TotalSpecificEnthalpy_fwd( 82 NewtonianIdealGasContext gas, const State s, const State ds) { 83 // Ignoring potential energy 84 CeedScalar de_kinetic = Dot3(ds.Y.velocity, s.Y.velocity); 85 CeedScalar de_internal = gas->cv * ds.Y.temperature; 86 return de_internal + de_kinetic + ds.Y.pressure/s.U.density 87 - s.Y.pressure/Square(s.U.density)*ds.U.density; 88 } 89 90 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative( 91 NewtonianIdealGasContext gas, StateConservative U, const CeedScalar x[3]) { 92 StatePrimitive Y; 93 for (CeedInt i=0; i<3; i++) Y.velocity[i] = U.momentum[i] / U.density; 94 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 95 CeedScalar e_potential = -Dot3(gas->g, x); 96 CeedScalar e_total = U.E_total / U.density; 97 CeedScalar e_internal = e_total - e_kinetic - e_potential; 98 Y.temperature = e_internal / gas->cv; 99 Y.pressure = (HeatCapacityRatio(gas) - 1) * U.density * e_internal; 100 return Y; 101 } 102 103 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd( 104 NewtonianIdealGasContext gas, State s, StateConservative dU, 105 const CeedScalar x[3], const CeedScalar dx[3]) { 106 StatePrimitive dY; 107 for (CeedInt i=0; i<3; i++) { 108 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 109 } 110 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 111 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 112 CeedScalar e_potential = -Dot3(gas->g, x); 113 CeedScalar de_potential = -Dot3(gas->g, dx); 114 CeedScalar e_total = s.U.E_total / s.U.density; 115 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 116 CeedScalar e_internal = e_total - e_kinetic - e_potential; 117 CeedScalar de_internal = de_total - de_kinetic - de_potential; 118 dY.temperature = de_internal / gas->cv; 119 dY.pressure = (HeatCapacityRatio(gas) - 1) 120 * (dU.density * e_internal + s.U.density * de_internal); 121 return dY; 122 } 123 124 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive( 125 NewtonianIdealGasContext gas, StatePrimitive Y, const CeedScalar x[3]) { 126 StateConservative U; 127 U.density = Y.pressure / (GasConstant(gas) * Y.temperature); 128 for (int i=0; i<3; i++) U.momentum[i] = U.density*Y.velocity[i]; 129 CeedScalar e_internal = gas->cv * Y.temperature; 130 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 131 CeedScalar e_potential = -Dot3(gas->g, x); 132 CeedScalar e_total = e_internal + e_kinetic + e_potential; 133 U.E_total = U.density*e_total; 134 return U; 135 } 136 137 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive_fwd( 138 NewtonianIdealGasContext gas, State s, StatePrimitive dY, 139 const CeedScalar x[3], const CeedScalar dx[3]) { 140 StateConservative dU; 141 dU.density = (dY.pressure * s.Y.temperature - s.Y.pressure * dY.temperature) / 142 (GasConstant(gas) * s.Y.temperature * s.Y.temperature); 143 for (int i=0; i<3; i++) { 144 dU.momentum[i] = dU.density * s.Y.velocity[i] + s.U.density * dY.velocity[i]; 145 } 146 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 147 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 148 CeedScalar e_potential = -Dot3(gas->g, x); 149 CeedScalar de_potential = -Dot3(gas->g, dx); 150 CeedScalar e_internal = gas->cv * s.Y.temperature; 151 CeedScalar de_internal = gas->cv * dY.temperature; 152 CeedScalar e_total = e_internal + e_kinetic + e_potential; 153 CeedScalar de_total = de_internal + de_kinetic + de_potential; 154 dU.E_total = dU.density*e_total + s.U.density*de_total; 155 return dU; 156 } 157 158 // Function pointer types for generic state array -> State struct functions 159 typedef State (*StateFromQi_t)(NewtonianIdealGasContext gas, 160 const CeedScalar qi[5], const CeedScalar x[3]); 161 typedef State (*StateFromQi_fwd_t)(NewtonianIdealGasContext gas, 162 State s, const CeedScalar dqi[5], 163 const CeedScalar x[3], const CeedScalar dx[3]); 164 165 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, 166 const CeedScalar U[5], const CeedScalar x[3]) { 167 State s; 168 s.U.density = U[0]; 169 s.U.momentum[0] = U[1]; 170 s.U.momentum[1] = U[2]; 171 s.U.momentum[2] = U[3]; 172 s.U.E_total = U[4]; 173 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 174 return s; 175 } 176 177 CEED_QFUNCTION_HELPER State StateFromU_fwd(NewtonianIdealGasContext gas, 178 State s, const CeedScalar dU[5], 179 const CeedScalar x[3], const CeedScalar dx[3]) { 180 State ds; 181 ds.U.density = dU[0]; 182 ds.U.momentum[0] = dU[1]; 183 ds.U.momentum[1] = dU[2]; 184 ds.U.momentum[2] = dU[3]; 185 ds.U.E_total = dU[4]; 186 ds.Y = StatePrimitiveFromConservative_fwd(gas, s, ds.U, x, dx); 187 return ds; 188 } 189 190 CEED_QFUNCTION_HELPER State StateFromY(NewtonianIdealGasContext gas, 191 const CeedScalar Y[5], const CeedScalar x[3]) { 192 State s; 193 s.Y.pressure = Y[0]; 194 s.Y.velocity[0] = Y[1]; 195 s.Y.velocity[1] = Y[2]; 196 s.Y.velocity[2] = Y[3]; 197 s.Y.temperature = Y[4]; 198 s.U = StateConservativeFromPrimitive(gas, s.Y, x); 199 return s; 200 } 201 202 CEED_QFUNCTION_HELPER State StateFromY_fwd(NewtonianIdealGasContext gas, 203 State s, const CeedScalar dY[5], 204 const CeedScalar x[3], const CeedScalar dx[3]) { 205 State ds; 206 ds.Y.pressure = dY[0]; 207 ds.Y.velocity[0] = dY[1]; 208 ds.Y.velocity[1] = dY[2]; 209 ds.Y.velocity[2] = dY[3]; 210 ds.Y.temperature = dY[4]; 211 ds.U = StateConservativeFromPrimitive_fwd(gas, s, ds.Y, x, dx); 212 return ds; 213 } 214 215 // Function pointer types for State struct -> generic state array 216 typedef void (*StateToQi_t)(NewtonianIdealGasContext gas, 217 const State input, CeedScalar qi[5]); 218 219 CEED_QFUNCTION_HELPER void StateToU(NewtonianIdealGasContext gas, 220 const State input, CeedScalar U[5]) { 221 UnpackState_U(input.U, U); 222 } 223 224 CEED_QFUNCTION_HELPER void StateToY(NewtonianIdealGasContext gas, 225 const State input, CeedScalar Y[5]) { 226 UnpackState_Y(input.Y, Y); 227 } 228 229 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, 230 StateConservative Flux[3]) { 231 for (CeedInt i=0; i<3; i++) { 232 Flux[i].density = s.U.momentum[i]; 233 for (CeedInt j=0; j<3; j++) 234 Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] 235 + s.Y.pressure * (i == j); 236 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 237 } 238 } 239 240 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, 241 State s, State ds, StateConservative dFlux[3]) { 242 for (CeedInt i=0; i<3; i++) { 243 dFlux[i].density = ds.U.momentum[i]; 244 for (CeedInt j=0; j<3; j++) 245 dFlux[i].momentum[j] = ds.U.momentum[i] * s.Y.velocity[j] + 246 s.U.momentum[i] * ds.Y.velocity[j] + ds.Y.pressure * (i == j); 247 dFlux[i].E_total = (ds.U.E_total + ds.Y.pressure) * s.Y.velocity[i] + 248 (s.U.E_total + s.Y.pressure) * ds.Y.velocity[i]; 249 } 250 } 251 252 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal( 253 NewtonianIdealGasContext gas, State s, const CeedScalar normal[3]) { 254 StateConservative Flux[3], Flux_dot_n = {0}; 255 FluxInviscid(gas, s, Flux); 256 for (CeedInt i=0; i<3; i++) { 257 Flux_dot_n.density += Flux[i].density * normal[i]; 258 for (CeedInt j=0; j<3; j++) 259 Flux_dot_n.momentum[j] += Flux[i].momentum[j] * normal[i]; 260 Flux_dot_n.E_total += Flux[i].E_total * normal[i]; 261 } 262 return Flux_dot_n; 263 } 264 265 CEED_QFUNCTION_HELPER StateConservative FluxInviscidDotNormal_fwd( 266 NewtonianIdealGasContext gas, State s, State ds, const CeedScalar normal[3]) { 267 StateConservative dFlux[3], Flux_dot_n = {0}; 268 FluxInviscid_fwd(gas, s, ds, dFlux); 269 for (CeedInt i=0; i<3; i++) { 270 Flux_dot_n.density += dFlux[i].density * normal[i]; 271 for (CeedInt j=0; j<3; j++) 272 Flux_dot_n.momentum[j] += dFlux[i].momentum[j] * normal[i]; 273 Flux_dot_n.E_total += dFlux[i].E_total * normal[i]; 274 } 275 return Flux_dot_n; 276 } 277 278 CEED_QFUNCTION_HELPER void FluxInviscidStrong(NewtonianIdealGasContext gas, 279 State s, State ds[3], CeedScalar strong_conv[5]) { 280 for (CeedInt i=0; i<5; i++) strong_conv[i] = 0; 281 for (CeedInt i=0; i<3; i++) { 282 StateConservative dF[3]; 283 FluxInviscid_fwd(gas, s, ds[i], dF); 284 CeedScalar dF_i[5]; 285 UnpackState_U(dF[i], dF_i); 286 for (CeedInt j=0; j<5; j++) 287 strong_conv[j] += dF_i[j]; 288 } 289 } 290 291 CEED_QFUNCTION_HELPER void FluxTotal(const StateConservative F_inviscid[3], 292 CeedScalar stress[3][3], CeedScalar Fe[3], CeedScalar Flux[5][3]) { 293 for (CeedInt j=0; j<3; j++) { 294 Flux[0][j] = F_inviscid[j].density; 295 for (CeedInt k=0; k<3; k++) 296 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 297 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 298 } 299 } 300 301 CEED_QFUNCTION_HELPER void FluxTotal_Boundary( 302 const StateConservative F_inviscid[3], const CeedScalar stress[3][3], 303 const CeedScalar Fe[3], const CeedScalar normal[3], CeedScalar Flux[5]) { 304 305 for (CeedInt j=0; j<5; j++) Flux[j] = 0.; 306 for (CeedInt j=0; j<3; j++) { 307 Flux[0] += F_inviscid[j].density * normal[j]; 308 for (CeedInt k=0; k<3; k++) { 309 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * normal[j]; 310 } 311 Flux[4] += (F_inviscid[j].E_total + Fe[j]) * normal[j]; 312 } 313 } 314 315 // Kelvin-Mandel notation 316 CEED_QFUNCTION_HELPER void KMStrainRate(const State grad_s[3], 317 CeedScalar strain_rate[6]) { 318 const CeedScalar weight = 1 / sqrt(2.); 319 strain_rate[0] = grad_s[0].Y.velocity[0]; 320 strain_rate[1] = grad_s[1].Y.velocity[1]; 321 strain_rate[2] = grad_s[2].Y.velocity[2]; 322 strain_rate[3] = weight * (grad_s[2].Y.velocity[1] + grad_s[1].Y.velocity[2]); 323 strain_rate[4] = weight * (grad_s[2].Y.velocity[0] + grad_s[0].Y.velocity[2]); 324 strain_rate[5] = weight * (grad_s[1].Y.velocity[0] + grad_s[0].Y.velocity[1]); 325 } 326 327 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, 328 const CeedScalar strain_rate[6], CeedScalar stress[6]) { 329 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 330 for (CeedInt i=0; i<6; i++) { 331 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 332 } 333 } 334 335 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, 336 StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 337 CeedScalar Fe[3]) { 338 for (CeedInt i=0; i<3; i++) { 339 Fe[i] = - Y.velocity[0] * stress[0][i] 340 - Y.velocity[1] * stress[1][i] 341 - Y.velocity[2] * stress[2][i] 342 - gas->k * grad_s[i].Y.temperature; 343 } 344 } 345 346 CEED_QFUNCTION_HELPER void ViscousEnergyFlux_fwd(NewtonianIdealGasContext gas, 347 StatePrimitive Y, StatePrimitive dY, const State grad_ds[3], 348 const CeedScalar stress[3][3], const CeedScalar dstress[3][3], 349 CeedScalar dFe[3]) { 350 for (CeedInt i=0; i<3; i++) { 351 dFe[i] = - Y.velocity[0] * dstress[0][i] - dY.velocity[0] * stress[0][i] 352 - Y.velocity[1] * dstress[1][i] - dY.velocity[1] * stress[1][i] 353 - Y.velocity[2] * dstress[2][i] - dY.velocity[2] * stress[2][i] 354 - gas->k * grad_ds[i].Y.temperature; 355 } 356 } 357 358 #endif // newtonian_state_h 359