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