1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Newtonian fluids operator for HONEE 6 #include <ceed/types.h> 7 8 #include "newtonian_state.h" 9 #include "newtonian_types.h" 10 #include "stabilization.h" 11 #include "utils.h" 12 13 CEED_QFUNCTION_HELPER void InternalDampingLayer(const NewtonianIdealGasContext context, const State s, const CeedScalar sigma, CeedScalar damp_Y[5], 14 CeedScalar damp_residual[5]) { 15 ScaleN(damp_Y, sigma, 5); 16 State damp_s = StateFromY_fwd(context, s, damp_Y); 17 18 CeedScalar U[5]; 19 UnpackState_U(damp_s.U, U); 20 for (int i = 0; i < 5; i++) damp_residual[i] += U[i]; 21 } 22 23 // ***************************************************************************** 24 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 25 // ***************************************************************************** 26 CEED_QFUNCTION_HELPER int ICsNewtonianIG(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 27 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 28 29 const SetupContext context = (SetupContext)ctx; 30 31 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 32 CeedScalar q[5]; 33 State s = StateFromPrimitive(&context->gas, context->reference); 34 StateToQ(&context->gas, s, q, state_var); 35 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 36 } 37 return 0; 38 } 39 40 CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 41 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 42 } 43 44 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 45 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_PRIMITIVE); 46 } 47 48 CEED_QFUNCTION(ICsNewtonianIG_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 49 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_ENTROPY); 50 } 51 52 CEED_QFUNCTION_HELPER int MassFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 53 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 54 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 55 const CeedScalar(*q_data) = in[2]; 56 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 57 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 58 59 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 60 61 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 62 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 63 const CeedScalar qi_dot[5] = {q_dot[0][i], q_dot[1][i], q_dot[2][i], q_dot[3][i], q_dot[4][i]}; 64 const State s = StateFromQ(context, qi, state_var); 65 const State s_dot = StateFromQ(context, qi_dot, state_var); 66 CeedScalar wdetJ, dXdx[3][3]; 67 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 68 69 // Standard mass matrix term 70 for (CeedInt f = 0; f < 5; f++) { 71 v[f][i] = wdetJ * qi_dot[f]; 72 } 73 74 // Stabilization method: none (Galerkin), SU, or SUPG 75 State grad_s[3] = {{{0.}}}; 76 CeedScalar Tau_d[3], stab[5][3], body_force[5] = {0.}, divFdiff[5] = {0.}, U_dot[5]; 77 UnpackState_U(s_dot.U, U_dot); 78 Tau_diagPrim(context, s, dXdx, context->dt, Tau_d); 79 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, divFdiff, stab); 80 81 // Stabilized mass term 82 for (CeedInt j = 0; j < 5; j++) { 83 for (CeedInt k = 0; k < 3; k++) { 84 Grad_v[k][j][i] = wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 85 } 86 } 87 } 88 return 0; 89 } 90 91 CEED_QFUNCTION(MassFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 92 return MassFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 93 } 94 95 // ***************************************************************************** 96 // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method 97 // 98 // This is 3D compressible Navier-Stokes in conservation form with state variables of density, momentum density, and total energy density. 99 // 100 // State Variables: q = ( rho, U1, U2, U3, E ) 101 // rho - Mass Density 102 // Ui - Momentum Density, Ui = rho ui 103 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 104 // 105 // Navier-Stokes Equations: 106 // drho/dt + div( U ) = 0 107 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 108 // dE/dt + div( (E + P) u ) = div( Fe ) 109 // 110 // Viscous Stress: 111 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 112 // 113 // Thermal Stress: 114 // Fe = u Fu + k grad( T ) 115 // Equation of State 116 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 117 // 118 // Stabilization: 119 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 120 // f1 = rho sqrt(ui uj gij) 121 // gij = dXi/dX * dXi/dX 122 // TauC = Cc f1 / (8 gii) 123 // TauM = min( 1 , 1 / f1 ) 124 // TauE = TauM / (Ce cv) 125 // 126 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 127 // 128 // Constants: 129 // lambda = - 2 / 3, From Stokes hypothesis 130 // mu , Dynamic viscosity 131 // k , Thermal conductivity 132 // cv , Specific heat, constant volume 133 // cp , Specific heat, constant pressure 134 // g , Gravity 135 // gamma = cp / cv, Specific heat ratio 136 // 137 // We require the product of the inverse of the Jacobian (dXdx_j,k) and its transpose (dXdx_k,j) to properly compute integrals of the form: int( gradv 138 // gradu ) 139 // ***************************************************************************** 140 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 141 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 142 const bool use_divFdiff = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE; 143 144 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 145 const CeedScalar(*Grad_q) = in[1]; 146 const CeedScalar(*q_data) = in[2]; 147 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 148 const CeedScalar(*divFdiff)[CEED_Q_VLA] = use_divFdiff ? (const CeedScalar(*)[CEED_Q_VLA])in[4] : NULL; 149 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 150 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 151 152 const CeedScalar *g = context->g; 153 const CeedScalar dt = context->dt; 154 const CeedScalar idl_pressure = context->idl_pressure; 155 156 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 157 CeedScalar U[5], wdetJ, dXdx[3][3]; 158 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 159 for (int j = 0; j < 5; j++) U[j] = q[j][i]; 160 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 161 State s = StateFromU(context, U); 162 163 State grad_s[3]; 164 StatePhysicalGradientFromReference(Q, i, context, s, STATEVAR_CONSERVATIVE, Grad_q, dXdx, grad_s); 165 166 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 167 KMStrainRate_State(grad_s, strain_rate); 168 NewtonianStress(context, strain_rate, kmstress); 169 KMUnpack(kmstress, stress); 170 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 171 172 StateConservative F_inviscid[3]; 173 FluxInviscid(context, s, F_inviscid); 174 175 // Total flux 176 CeedScalar Flux[5][3]; 177 FluxTotal(F_inviscid, stress, Fe, Flux); 178 179 for (CeedInt j = 0; j < 5; j++) { 180 for (CeedInt k = 0; k < 3; k++) Grad_v[k][j][i] = wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 181 } 182 183 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], Dot3(s.U.momentum, g)}; 184 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 185 186 if (context->idl_enable) { 187 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 188 CeedScalar damp_state[5] = {s.Y.pressure - idl_pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 189 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 190 for (int j = 0; j < 5; j++) v[j][i] -= wdetJ * idl_residual[j]; 191 } 192 193 CeedScalar divFdiff_i[5] = {0.}; 194 if (use_divFdiff) 195 for (int j = 1; j < 5; j++) divFdiff_i[j] = divFdiff[j - 1][i]; 196 197 // -- Stabilization method: none (Galerkin), SU, or SUPG 198 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 199 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 200 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, divFdiff_i, stab); 201 202 for (CeedInt j = 0; j < 5; j++) { 203 for (CeedInt k = 0; k < 3; k++) Grad_v[k][j][i] -= wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 204 } 205 } 206 return 0; 207 } 208 209 // ***************************************************************************** 210 // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepping method 211 // 212 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 213 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 214 // (diffusive terms will be added later) 215 // ***************************************************************************** 216 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 217 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 218 const bool use_divFdiff = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE; 219 220 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 221 const CeedScalar(*Grad_q) = in[1]; 222 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 223 const CeedScalar(*q_data) = in[3]; 224 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 225 const CeedScalar(*divFdiff)[CEED_Q_VLA] = use_divFdiff ? (const CeedScalar(*)[CEED_Q_VLA])in[5] : NULL; 226 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 227 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 228 CeedScalar(*jac_data) = out[2]; 229 230 const CeedScalar *g = context->g; 231 const CeedScalar dt = context->dt; 232 const CeedScalar idl_pressure = context->idl_pressure; 233 234 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 235 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 236 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 237 const State s = StateFromQ(context, qi, state_var); 238 239 CeedScalar wdetJ, dXdx[3][3]; 240 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 241 State grad_s[3]; 242 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 243 244 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 245 KMStrainRate_State(grad_s, strain_rate); 246 NewtonianStress(context, strain_rate, kmstress); 247 KMUnpack(kmstress, stress); 248 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 249 250 StateConservative F_inviscid[3]; 251 FluxInviscid(context, s, F_inviscid); 252 253 // Total flux 254 CeedScalar Flux[5][3]; 255 FluxTotal(F_inviscid, stress, Fe, Flux); 256 257 for (CeedInt j = 0; j < 5; j++) { 258 for (CeedInt k = 0; k < 3; k++) { 259 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 260 } 261 } 262 263 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], Dot3(s.U.momentum, g)}; 264 265 // -- Stabilization method: none (Galerkin), SU, or SUPG 266 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5]; 267 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 268 State s_dot = StateFromQ_fwd(context, s, qi_dot, state_var); 269 UnpackState_U(s_dot.U, U_dot); 270 271 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 272 if (context->idl_enable) { 273 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 274 StoredValuesPack(Q, i, 14, 1, &sigma, jac_data); 275 CeedScalar damp_state[5] = {s.Y.pressure - idl_pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 276 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 277 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 278 } 279 280 CeedScalar divFdiff_i[5] = {0.}; 281 if (use_divFdiff) { 282 for (int j = 1; j < 5; j++) divFdiff_i[j] = divFdiff[j - 1][i]; 283 } 284 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 285 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, divFdiff_i, stab); 286 287 for (CeedInt j = 0; j < 5; j++) { 288 for (CeedInt k = 0; k < 3; k++) { 289 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 290 } 291 } 292 StoredValuesPack(Q, i, 0, 5, qi, jac_data); 293 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data); 294 StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data); 295 } 296 return 0; 297 } 298 299 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 300 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 301 } 302 303 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 304 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 305 } 306 307 CEED_QFUNCTION(IFunction_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 308 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 309 } 310 311 // ***************************************************************************** 312 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 313 // ***************************************************************************** 314 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 315 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 316 const CeedScalar(*Grad_dq) = in[1]; 317 const CeedScalar(*q_data) = in[2]; 318 const CeedScalar(*jac_data) = in[3]; 319 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 320 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 321 322 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 323 const CeedScalar *g = context->g; 324 325 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 326 CeedScalar wdetJ, dXdx[3][3]; 327 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 328 329 CeedScalar qi[5], kmstress[6], Tau_d[3]; 330 StoredValuesUnpack(Q, i, 0, 5, jac_data, qi); 331 StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress); 332 StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d); 333 State s = StateFromQ(context, qi, state_var); 334 335 CeedScalar dqi[5]; 336 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 337 State ds = StateFromQ_fwd(context, s, dqi, state_var); 338 339 State grad_ds[3]; 340 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 341 342 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 343 KMStrainRate_State(grad_ds, dstrain_rate); 344 NewtonianStress(context, dstrain_rate, dkmstress); 345 KMUnpack(dkmstress, dstress); 346 KMUnpack(kmstress, stress); 347 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 348 349 StateConservative dF_inviscid[3]; 350 FluxInviscid_fwd(context, s, ds, dF_inviscid); 351 352 // Total flux 353 CeedScalar dFlux[5][3]; 354 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 355 356 for (int j = 0; j < 5; j++) { 357 for (int k = 0; k < 3; k++) Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * dFlux[j][0] + dXdx[k][1] * dFlux[j][1] + dXdx[k][2] * dFlux[j][2]); 358 } 359 360 const CeedScalar dbody_force[5] = {0, ds.U.density * g[0], ds.U.density * g[1], ds.U.density * g[2], Dot3(ds.U.momentum, g)}; 361 CeedScalar dU[5] = {0.}; 362 UnpackState_U(ds.U, dU); 363 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 364 365 if (context->idl_enable) { 366 const CeedScalar sigma = jac_data[14 * Q + i]; 367 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 368 // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 369 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 370 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 371 } 372 373 // -- Stabilization method: none (Galerkin), SU, or SUPG 374 CeedScalar dstab[5][3], U_dot[5] = {0}; 375 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 376 const CeedScalar zeroFlux[5] = {0.}; 377 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, zeroFlux, dstab); 378 379 for (int j = 0; j < 5; j++) { 380 for (int k = 0; k < 3; k++) Grad_v[k][j][i] += wdetJ * (dstab[j][0] * dXdx[k][0] + dstab[j][1] * dXdx[k][1] + dstab[j][2] * dXdx[k][2]); 381 } 382 } 383 return 0; 384 } 385 386 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 387 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 388 } 389 390 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 391 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 392 } 393 394 CEED_QFUNCTION(IJacobian_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 395 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 396 } 397 398 // ***************************************************************************** 399 // Compute boundary integral (ie. for strongly set inflows) 400 // ***************************************************************************** 401 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 402 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 403 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 404 const CeedScalar(*Grad_q) = in[1]; 405 const CeedScalar(*q_data_sur) = in[2]; 406 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 407 CeedScalar(*jac_data_sur) = context->is_implicit ? out[1] : NULL; 408 409 const bool is_implicit = context->is_implicit; 410 411 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 412 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 413 State s = StateFromQ(context, qi, state_var); 414 415 CeedScalar wdetJb, dXdx[2][3], normal[3]; 416 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 417 wdetJb *= is_implicit ? -1. : 1.; 418 419 State grad_s[3]; 420 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 421 422 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 423 KMStrainRate_State(grad_s, strain_rate); 424 NewtonianStress(context, strain_rate, kmstress); 425 KMUnpack(kmstress, stress); 426 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 427 428 StateConservative F_inviscid[3]; 429 FluxInviscid(context, s, F_inviscid); 430 431 CeedScalar Flux[5]; 432 FluxTotal_Boundary(F_inviscid, stress, Fe, normal, Flux); 433 434 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 435 436 if (is_implicit) { 437 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 438 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 439 } 440 } 441 return 0; 442 } 443 444 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 445 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 446 } 447 448 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 449 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE); 450 } 451 452 CEED_QFUNCTION(BoundaryIntegral_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 453 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_ENTROPY); 454 } 455 456 // ***************************************************************************** 457 // Jacobian for "set nothing" boundary integral 458 // ***************************************************************************** 459 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 460 StateVariable state_var) { 461 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 462 const CeedScalar(*Grad_dq) = in[1]; 463 const CeedScalar(*q_data_sur) = in[2]; 464 const CeedScalar(*jac_data_sur) = in[4]; 465 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 466 467 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 468 const bool is_implicit = context->is_implicit; 469 470 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 471 CeedScalar wdetJb, dXdx[2][3], normal[3]; 472 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 473 wdetJb *= is_implicit ? -1. : 1.; 474 475 CeedScalar qi[5], kmstress[6], dqi[5]; 476 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 477 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 478 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 479 480 State s = StateFromQ(context, qi, state_var); 481 State ds = StateFromQ_fwd(context, s, dqi, state_var); 482 483 State grad_ds[3]; 484 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 485 486 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 487 KMStrainRate_State(grad_ds, dstrain_rate); 488 NewtonianStress(context, dstrain_rate, dkmstress); 489 KMUnpack(dkmstress, dstress); 490 KMUnpack(kmstress, stress); 491 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 492 493 StateConservative dF_inviscid[3]; 494 FluxInviscid_fwd(context, s, ds, dF_inviscid); 495 496 CeedScalar dFlux[5]; 497 FluxTotal_Boundary(dF_inviscid, dstress, dFe, normal, dFlux); 498 499 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 500 } 501 return 0; 502 } 503 504 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 505 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 506 } 507 508 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 509 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 510 } 511 512 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 513 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY); 514 } 515 516 // @brief Volume integral for RHS of divergence of diffusive flux direct projection 517 CEED_QFUNCTION_HELPER int DivDiffusiveFluxVolumeRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 518 StateVariable state_var) { 519 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 520 const CeedScalar(*Grad_q) = in[1]; 521 const CeedScalar(*q_data) = in[2]; 522 CeedScalar(*Grad_v)[4][CEED_Q_VLA] = (CeedScalar(*)[4][CEED_Q_VLA])out[0]; 523 524 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 525 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 526 527 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 528 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 529 const State s = StateFromQ(context, qi, state_var); 530 CeedScalar wdetJ, dXdx[3][3]; 531 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 532 533 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 534 { // Get stress and Fe 535 State grad_s[3]; 536 CeedScalar strain_rate[6], kmstress[6]; 537 538 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 539 KMStrainRate_State(grad_s, strain_rate); 540 NewtonianStress(context, strain_rate, kmstress); 541 KMUnpack(kmstress, stress); 542 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 543 } 544 545 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 546 547 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 548 for (CeedInt k = 0; k < 3; k++) { 549 Grad_v[k][j - 1][i] = -wdetJ * Dot3(dXdx[k], Fdiff[j]); 550 } 551 } 552 } 553 return 0; 554 } 555 556 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 557 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 558 } 559 560 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 561 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 562 } 563 564 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 565 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 566 } 567 568 // @brief Boundary integral for RHS of divergence of diffusive flux direct projection 569 CEED_QFUNCTION_HELPER int DivDiffusiveFluxBoundaryRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 570 StateVariable state_var) { 571 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 572 const CeedScalar(*Grad_q) = in[1]; 573 const CeedScalar(*q_data) = in[2]; 574 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 575 576 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 577 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 578 579 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 580 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 581 const State s = StateFromQ(context, qi, state_var); 582 CeedScalar wdetJ, dXdx[3][3], normal[3]; 583 CeedScalar stress[3][3], Fe[3], Fdiff[5]; 584 585 QdataBoundaryGradientUnpack_3D(Q, i, q_data, &wdetJ, dXdx, normal); 586 { // Get stress and Fe 587 State grad_s[3]; 588 CeedScalar strain_rate[6], kmstress[6]; 589 590 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 591 KMStrainRate_State(grad_s, strain_rate); 592 NewtonianStress(context, strain_rate, kmstress); 593 KMUnpack(kmstress, stress); 594 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 595 } 596 597 FluxTotal_Boundary(ZeroInviscidFluxes, stress, Fe, normal, Fdiff); 598 599 // Continuity has no diffusive flux, therefore skip 600 for (CeedInt j = 1; j < 5; j++) v[j - 1][i] = wdetJ * Fdiff[j]; 601 } 602 return 0; 603 } 604 605 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 606 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 607 } 608 609 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 610 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 611 } 612 613 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 614 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 615 } 616 617 // @brief Integral for RHS of diffusive flux indirect projection 618 CEED_QFUNCTION_HELPER int DiffusiveFluxRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 619 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 620 const CeedScalar(*Grad_q) = in[1]; 621 const CeedScalar(*q_data) = in[2]; 622 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 623 624 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 625 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 626 627 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 628 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 629 const State s = StateFromQ(context, qi, state_var); 630 CeedScalar wdetJ, dXdx[3][3]; 631 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 632 633 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 634 { // Get stress and Fe 635 State grad_s[3]; 636 CeedScalar strain_rate[6], kmstress[6]; 637 638 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 639 KMStrainRate_State(grad_s, strain_rate); 640 NewtonianStress(context, strain_rate, kmstress); 641 KMUnpack(kmstress, stress); 642 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 643 } 644 645 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 646 647 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 648 for (CeedInt k = 0; k < 3; k++) { 649 v[(j - 1) * 3 + k][i] = wdetJ * Fdiff[j][k]; 650 } 651 } 652 } 653 return 0; 654 } 655 656 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 657 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 658 } 659 660 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 661 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 662 } 663 664 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 665 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 666 } 667