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