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