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 StatePrimitive StatePrimitiveFromConservative( 50 NewtonianIdealGasContext gas, StateConservative U, const CeedScalar x[3]) { 51 StatePrimitive Y; 52 for (CeedInt i=0; i<3; i++) Y.velocity[i] = U.momentum[i] / U.density; 53 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 54 CeedScalar e_potential = -Dot3(gas->g, x); 55 CeedScalar e_total = U.E_total / U.density; 56 CeedScalar e_internal = e_total - e_kinetic - e_potential; 57 Y.temperature = e_internal / gas->cv; 58 Y.pressure = (gas->cp / gas->cv - 1) * U.density * e_internal; 59 return Y; 60 } 61 62 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd( 63 NewtonianIdealGasContext gas, State s, StateConservative dU, 64 const CeedScalar x[3], const CeedScalar dx[3]) { 65 StatePrimitive dY; 66 for (CeedInt i=0; i<3; i++) { 67 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 68 } 69 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 70 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 71 CeedScalar e_potential = -Dot3(gas->g, x); 72 CeedScalar de_potential = -Dot3(gas->g, dx); 73 CeedScalar e_total = s.U.E_total / s.U.density; 74 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 75 CeedScalar e_internal = e_total - e_kinetic - e_potential; 76 CeedScalar de_internal = de_total - de_kinetic - de_potential; 77 dY.temperature = de_internal / gas->cv; 78 dY.pressure = (gas->cp / gas->cv - 1) 79 * (dU.density * e_internal + s.U.density * de_internal); 80 return dY; 81 } 82 83 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive( 84 NewtonianIdealGasContext gas, StatePrimitive Y, const CeedScalar x[3]) { 85 StateConservative U; 86 CeedScalar R = gas->cp - gas->cv; 87 U.density = Y.pressure / (R * Y.temperature); 88 for (int i=0; i<3; i++) U.momentum[i] = U.density*Y.velocity[i]; 89 CeedScalar e_internal = gas->cv * Y.temperature; 90 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 91 CeedScalar e_potential = -Dot3(gas->g, x); 92 CeedScalar e_total = e_internal + e_kinetic + e_potential; 93 U.E_total = U.density*e_total; 94 return U; 95 } 96 97 CEED_QFUNCTION_HELPER StateConservative StateConservativeFromPrimitive_fwd( 98 NewtonianIdealGasContext gas, State s, StatePrimitive dY, 99 const CeedScalar x[3], const CeedScalar dx[3]) { 100 StateConservative dU; 101 CeedScalar R = gas->cp - gas->cv; 102 dU.density = (dY.pressure * s.Y.temperature - s.Y.pressure * dY.temperature) / 103 (R * s.Y.temperature * s.Y.temperature); 104 for (int i=0; i<3; i++) { 105 dU.momentum[i] = dU.density * s.Y.velocity[i] + s.U.density * dY.velocity[i]; 106 } 107 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 108 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 109 CeedScalar e_potential = -Dot3(gas->g, x); 110 CeedScalar de_potential = -Dot3(gas->g, dx); 111 CeedScalar e_internal = gas->cv * s.Y.temperature; 112 CeedScalar de_internal = gas->cv * dY.temperature; 113 CeedScalar e_total = e_internal + e_kinetic + e_potential; 114 CeedScalar de_total = de_internal + de_kinetic + de_potential; 115 dU.E_total = dU.density*e_total + s.U.density*de_total; 116 return dU; 117 } 118 119 // Function pointer types for generic state array -> State struct functions 120 typedef State (*StateFromQi_t)(NewtonianIdealGasContext gas, 121 const CeedScalar qi[5], const CeedScalar x[3]); 122 typedef State (*StateFromQi_fwd_t)(NewtonianIdealGasContext gas, 123 State s, const CeedScalar dqi[5], 124 const CeedScalar x[3], const CeedScalar dx[3]); 125 126 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, 127 const CeedScalar U[5], const CeedScalar x[3]) { 128 State s; 129 s.U.density = U[0]; 130 s.U.momentum[0] = U[1]; 131 s.U.momentum[1] = U[2]; 132 s.U.momentum[2] = U[3]; 133 s.U.E_total = U[4]; 134 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 135 return s; 136 } 137 138 CEED_QFUNCTION_HELPER State StateFromU_fwd(NewtonianIdealGasContext gas, 139 State s, const CeedScalar dU[5], 140 const CeedScalar x[3], const CeedScalar dx[3]) { 141 State ds; 142 ds.U.density = dU[0]; 143 ds.U.momentum[0] = dU[1]; 144 ds.U.momentum[1] = dU[2]; 145 ds.U.momentum[2] = dU[3]; 146 ds.U.E_total = dU[4]; 147 ds.Y = StatePrimitiveFromConservative_fwd(gas, s, ds.U, x, dx); 148 return ds; 149 } 150 151 CEED_QFUNCTION_HELPER State StateFromY(NewtonianIdealGasContext gas, 152 const CeedScalar Y[5], const CeedScalar x[3]) { 153 State s; 154 s.Y.pressure = Y[0]; 155 s.Y.velocity[0] = Y[1]; 156 s.Y.velocity[1] = Y[2]; 157 s.Y.velocity[2] = Y[3]; 158 s.Y.temperature = Y[4]; 159 s.U = StateConservativeFromPrimitive(gas, s.Y, x); 160 return s; 161 } 162 163 CEED_QFUNCTION_HELPER State StateFromY_fwd(NewtonianIdealGasContext gas, 164 State s, const CeedScalar dY[5], 165 const CeedScalar x[3], const CeedScalar dx[3]) { 166 State ds; 167 ds.Y.pressure = dY[0]; 168 ds.Y.velocity[0] = dY[1]; 169 ds.Y.velocity[1] = dY[2]; 170 ds.Y.velocity[2] = dY[3]; 171 ds.Y.temperature = dY[4]; 172 ds.U = StateConservativeFromPrimitive_fwd(gas, s, ds.Y, x, dx); 173 return ds; 174 } 175 176 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, 177 StateConservative Flux[3]) { 178 for (CeedInt i=0; i<3; i++) { 179 Flux[i].density = s.U.momentum[i]; 180 for (CeedInt j=0; j<3; j++) 181 Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] 182 + s.Y.pressure * (i == j); 183 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 184 } 185 } 186 187 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, 188 State s, State ds, StateConservative dFlux[3]) { 189 for (CeedInt i=0; i<3; i++) { 190 dFlux[i].density = ds.U.momentum[i]; 191 for (CeedInt j=0; j<3; j++) 192 dFlux[i].momentum[j] = ds.U.momentum[i] * s.Y.velocity[j] + 193 s.U.momentum[i] * ds.Y.velocity[j] + ds.Y.pressure * (i == j); 194 dFlux[i].E_total = (ds.U.E_total + ds.Y.pressure) * s.Y.velocity[i] + 195 (s.U.E_total + s.Y.pressure) * ds.Y.velocity[i]; 196 } 197 } 198 199 CEED_QFUNCTION_HELPER void FluxInviscidStrong(NewtonianIdealGasContext gas, 200 State s, State ds[3], CeedScalar strong_conv[5]) { 201 for (CeedInt i=0; i<5; i++) strong_conv[i] = 0; 202 for (CeedInt i=0; i<3; i++) { 203 StateConservative dF[3]; 204 FluxInviscid_fwd(gas, s, ds[i], dF); 205 CeedScalar dF_i[5]; 206 UnpackState_U(dF[i], dF_i); 207 for (CeedInt j=0; j<5; j++) 208 strong_conv[j] += dF_i[j]; 209 } 210 } 211 212 CEED_QFUNCTION_HELPER void FluxTotal(const StateConservative F_inviscid[3], 213 CeedScalar stress[3][3], CeedScalar Fe[3], CeedScalar Flux[5][3]) { 214 for (CeedInt j=0; j<3; j++) { 215 Flux[0][j] = F_inviscid[j].density; 216 for (CeedInt k=0; k<3; k++) 217 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 218 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 219 } 220 } 221 222 CEED_QFUNCTION_HELPER void FluxTotal_Boundary( 223 const StateConservative F_inviscid[3], const CeedScalar stress[3][3], 224 const CeedScalar Fe[3], const CeedScalar normal[3], CeedScalar Flux[5]) { 225 226 for(CeedInt j=0; j<5; j++) Flux[j] = 0.; 227 for (CeedInt j=0; j<3; j++) { 228 Flux[0] += F_inviscid[j].density * normal[j]; 229 for (CeedInt k=0; k<3; k++) { 230 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * normal[j]; 231 } 232 Flux[4] += (F_inviscid[j].E_total + Fe[j]) * normal[j]; 233 } 234 } 235 236 // Kelvin-Mandel notation 237 CEED_QFUNCTION_HELPER void KMStrainRate(const State grad_s[3], 238 CeedScalar strain_rate[6]) { 239 const CeedScalar weight = 1 / sqrt(2.); 240 strain_rate[0] = grad_s[0].Y.velocity[0]; 241 strain_rate[1] = grad_s[1].Y.velocity[1]; 242 strain_rate[2] = grad_s[2].Y.velocity[2]; 243 strain_rate[3] = weight * (grad_s[2].Y.velocity[1] + grad_s[1].Y.velocity[2]); 244 strain_rate[4] = weight * (grad_s[2].Y.velocity[0] + grad_s[0].Y.velocity[2]); 245 strain_rate[5] = weight * (grad_s[1].Y.velocity[0] + grad_s[0].Y.velocity[1]); 246 } 247 248 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, 249 const CeedScalar strain_rate[6], CeedScalar stress[6]) { 250 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 251 for (CeedInt i=0; i<6; i++) { 252 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 253 } 254 } 255 256 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, 257 StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 258 CeedScalar Fe[3]) { 259 for (CeedInt i=0; i<3; i++) { 260 Fe[i] = - Y.velocity[0] * stress[0][i] 261 - Y.velocity[1] * stress[1][i] 262 - Y.velocity[2] * stress[2][i] 263 - gas->k * grad_s[i].Y.temperature; 264 } 265 } 266 267 CEED_QFUNCTION_HELPER void ViscousEnergyFlux_fwd(NewtonianIdealGasContext gas, 268 StatePrimitive Y, StatePrimitive dY, const State grad_ds[3], 269 const CeedScalar stress[3][3], 270 const CeedScalar dstress[3][3], 271 CeedScalar dFe[3]) { 272 for (CeedInt i=0; i<3; i++) { 273 dFe[i] = - Y.velocity[0] * dstress[0][i] - dY.velocity[0] * stress[0][i] 274 - Y.velocity[1] * dstress[1][i] - dY.velocity[1] * stress[1][i] 275 - Y.velocity[2] * dstress[2][i] - dY.velocity[2] * stress[2][i] 276 - gas->k * grad_ds[i].Y.temperature; 277 } 278 } 279 280 #endif // newtonian_state_h 281