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