1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Operator for Navier-Stokes example using PETSc 10 #include <ceed/types.h> 11 #ifndef CEED_RUNNING_JIT_PASS 12 #include <math.h> 13 #include <stdlib.h> 14 #endif 15 16 #include "newtonian_state.h" 17 #include "newtonian_types.h" 18 #include "stabilization.h" 19 #include "utils.h" 20 21 CEED_QFUNCTION_HELPER void InternalDampingLayer(const NewtonianIdealGasContext context, const State s, const CeedScalar sigma, CeedScalar damp_Y[5], 22 CeedScalar damp_residual[5]) { 23 ScaleN(damp_Y, sigma, 5); 24 State damp_s = StateFromY_fwd(context, s, damp_Y); 25 26 CeedScalar U[5]; 27 UnpackState_U(damp_s.U, U); 28 for (int i = 0; i < 5; i++) damp_residual[i] += U[i]; 29 } 30 31 // ***************************************************************************** 32 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 33 // ***************************************************************************** 34 CEED_QFUNCTION_HELPER int ICsNewtonianIG(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 35 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 36 37 const SetupContext context = (SetupContext)ctx; 38 39 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 40 CeedScalar q[5] = {0.}; 41 State s = StateFromPrimitive(&context->gas, context->reference); 42 StateToQ(&context->gas, s, q, state_var); 43 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 44 } 45 return 0; 46 } 47 48 CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 49 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 50 } 51 52 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 53 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_PRIMITIVE); 54 } 55 56 CEED_QFUNCTION(ICsNewtonianIG_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 57 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_ENTROPY); 58 } 59 60 CEED_QFUNCTION_HELPER void MassFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 61 StateVariable state_var) { 62 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 63 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 64 const CeedScalar(*q_data) = in[2]; 65 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 66 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 67 68 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 69 70 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 71 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 72 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]}; 73 const State s = StateFromQ(context, qi, state_var); 74 const State s_dot = StateFromQ(context, qi_dot, state_var); 75 CeedScalar wdetJ, dXdx[3][3]; 76 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 77 78 // Standard mass matrix term 79 for (CeedInt f = 0; f < 5; f++) { 80 v[f][i] = wdetJ * qi_dot[f]; 81 } 82 83 // Stabilization method: none (Galerkin), SU, or SUPG 84 State grad_s[3] = {{{0.}}}; 85 CeedScalar Tau_d[3], stab[5][3], body_force[5] = {0.}, U_dot[5]; 86 UnpackState_U(s_dot.U, U_dot); 87 Tau_diagPrim(context, s, dXdx, context->dt, Tau_d); 88 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab); 89 90 // Stabilized mass term 91 for (CeedInt j = 0; j < 5; j++) { 92 for (CeedInt k = 0; k < 3; k++) { 93 Grad_v[k][j][i] = wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 94 } 95 } 96 } 97 } 98 99 CEED_QFUNCTION(MassFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 100 MassFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 101 return 0; 102 } 103 104 // ***************************************************************************** 105 // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method 106 // 107 // This is 3D compressible Navier-Stokes in conservation form with state variables of density, momentum density, and total energy density. 108 // 109 // State Variables: q = ( rho, U1, U2, U3, E ) 110 // rho - Mass Density 111 // Ui - Momentum Density, Ui = rho ui 112 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 113 // 114 // Navier-Stokes Equations: 115 // drho/dt + div( U ) = 0 116 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 117 // dE/dt + div( (E + P) u ) = div( Fe ) 118 // 119 // Viscous Stress: 120 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 121 // 122 // Thermal Stress: 123 // Fe = u Fu + k grad( T ) 124 // Equation of State 125 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 126 // 127 // Stabilization: 128 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 129 // f1 = rho sqrt(ui uj gij) 130 // gij = dXi/dX * dXi/dX 131 // TauC = Cc f1 / (8 gii) 132 // TauM = min( 1 , 1 / f1 ) 133 // TauE = TauM / (Ce cv) 134 // 135 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 136 // 137 // Constants: 138 // lambda = - 2 / 3, From Stokes hypothesis 139 // mu , Dynamic viscosity 140 // k , Thermal conductivity 141 // cv , Specific heat, constant volume 142 // cp , Specific heat, constant pressure 143 // g , Gravity 144 // gamma = cp / cv, Specific heat ratio 145 // 146 // 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 147 // gradu ) 148 // ***************************************************************************** 149 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 150 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 151 const CeedScalar(*Grad_q) = in[1]; 152 const CeedScalar(*q_data) = in[2]; 153 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 154 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 155 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 156 157 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 158 const CeedScalar *g = context->g; 159 const CeedScalar dt = context->dt; 160 const CeedScalar P0 = context->idl_pressure; 161 162 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 163 CeedScalar U[5], wdetJ, dXdx[3][3]; 164 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 165 for (int j = 0; j < 5; j++) U[j] = q[j][i]; 166 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 167 State s = StateFromU(context, U); 168 169 State grad_s[3]; 170 StatePhysicalGradientFromReference(Q, i, context, s, STATEVAR_CONSERVATIVE, Grad_q, dXdx, grad_s); 171 172 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 173 KMStrainRate_State(grad_s, strain_rate); 174 NewtonianStress(context, strain_rate, kmstress); 175 KMUnpack(kmstress, stress); 176 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 177 178 StateConservative F_inviscid[3]; 179 FluxInviscid(context, s, F_inviscid); 180 181 // Total flux 182 CeedScalar Flux[5][3]; 183 FluxTotal(F_inviscid, stress, Fe, Flux); 184 185 for (CeedInt j = 0; j < 5; j++) { 186 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]); 187 } 188 189 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)}; 190 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 191 192 if (context->idl_enable) { 193 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 194 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 195 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 196 for (int j = 0; j < 5; j++) v[j][i] -= wdetJ * idl_residual[j]; 197 } 198 199 // -- Stabilization method: none (Galerkin), SU, or SUPG 200 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 201 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 202 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab); 203 204 for (CeedInt j = 0; j < 5; j++) { 205 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]); 206 } 207 } 208 return 0; 209 } 210 211 // ***************************************************************************** 212 // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepping method 213 // 214 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 215 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 216 // (diffusive terms will be added later) 217 // ***************************************************************************** 218 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 219 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 220 const CeedScalar(*Grad_q) = in[1]; 221 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 222 const CeedScalar(*q_data) = in[3]; 223 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 224 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 225 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 226 CeedScalar(*jac_data) = out[2]; 227 228 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 229 const CeedScalar *g = context->g; 230 const CeedScalar dt = context->dt; 231 const CeedScalar P0 = context->idl_pressure; 232 233 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 234 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 235 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 236 const State s = StateFromQ(context, qi, state_var); 237 238 CeedScalar wdetJ, dXdx[3][3]; 239 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 240 State grad_s[3]; 241 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 242 243 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 244 KMStrainRate_State(grad_s, strain_rate); 245 NewtonianStress(context, strain_rate, kmstress); 246 KMUnpack(kmstress, stress); 247 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 248 249 StateConservative F_inviscid[3]; 250 FluxInviscid(context, s, F_inviscid); 251 252 // Total flux 253 CeedScalar Flux[5][3]; 254 FluxTotal(F_inviscid, stress, Fe, Flux); 255 256 for (CeedInt j = 0; j < 5; j++) { 257 for (CeedInt k = 0; k < 3; k++) { 258 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 259 } 260 } 261 262 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)}; 263 264 // -- Stabilization method: none (Galerkin), SU, or SUPG 265 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5]; 266 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 267 State s_dot = StateFromQ_fwd(context, s, qi_dot, state_var); 268 UnpackState_U(s_dot.U, U_dot); 269 270 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 271 if (context->idl_enable) { 272 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 273 StoredValuesPack(Q, i, 14, 1, &sigma, jac_data); 274 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 275 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 276 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 277 } 278 279 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 280 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab); 281 282 for (CeedInt j = 0; j < 5; j++) { 283 for (CeedInt k = 0; k < 3; k++) { 284 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 285 } 286 } 287 StoredValuesPack(Q, i, 0, 5, qi, jac_data); 288 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data); 289 StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data); 290 } 291 return 0; 292 } 293 294 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 295 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 296 } 297 298 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 299 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 300 } 301 302 CEED_QFUNCTION(IFunction_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 303 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 304 } 305 306 // ***************************************************************************** 307 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 308 // ***************************************************************************** 309 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 310 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 311 const CeedScalar(*Grad_dq) = in[1]; 312 const CeedScalar(*q_data) = in[2]; 313 const CeedScalar(*jac_data) = in[3]; 314 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 315 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 316 317 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 318 const CeedScalar *g = context->g; 319 320 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 321 CeedScalar wdetJ, dXdx[3][3]; 322 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 323 324 CeedScalar qi[5], kmstress[6], Tau_d[3]; 325 StoredValuesUnpack(Q, i, 0, 5, jac_data, qi); 326 StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress); 327 StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d); 328 State s = StateFromQ(context, qi, state_var); 329 330 CeedScalar dqi[5]; 331 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 332 State ds = StateFromQ_fwd(context, s, dqi, state_var); 333 334 State grad_ds[3]; 335 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 336 337 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 338 KMStrainRate_State(grad_ds, dstrain_rate); 339 NewtonianStress(context, dstrain_rate, dkmstress); 340 KMUnpack(dkmstress, dstress); 341 KMUnpack(kmstress, stress); 342 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 343 344 StateConservative dF_inviscid[3]; 345 FluxInviscid_fwd(context, s, ds, dF_inviscid); 346 347 // Total flux 348 CeedScalar dFlux[5][3]; 349 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 350 351 for (int j = 0; j < 5; j++) { 352 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]); 353 } 354 355 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)}; 356 CeedScalar dU[5] = {0.}; 357 UnpackState_U(ds.U, dU); 358 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 359 360 if (context->idl_enable) { 361 const CeedScalar sigma = jac_data[14 * Q + i]; 362 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 363 // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 364 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 365 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 366 } 367 368 // -- Stabilization method: none (Galerkin), SU, or SUPG 369 CeedScalar dstab[5][3], U_dot[5] = {0}; 370 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 371 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, dstab); 372 373 for (int j = 0; j < 5; j++) { 374 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]); 375 } 376 } 377 return 0; 378 } 379 380 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 381 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 382 } 383 384 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 385 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 386 } 387 388 CEED_QFUNCTION(IJacobian_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 389 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 390 } 391 392 // ***************************************************************************** 393 // Compute boundary integral (ie. for strongly set inflows) 394 // ***************************************************************************** 395 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 396 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 397 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 398 const CeedScalar(*Grad_q) = in[1]; 399 const CeedScalar(*q_data_sur) = in[2]; 400 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 401 CeedScalar(*jac_data_sur) = context->is_implicit ? out[1] : NULL; 402 403 const bool is_implicit = context->is_implicit; 404 405 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 406 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 407 State s = StateFromQ(context, qi, state_var); 408 409 CeedScalar wdetJb, dXdx[2][3], norm[3]; 410 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 411 wdetJb *= is_implicit ? -1. : 1.; 412 413 State grad_s[3]; 414 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 415 416 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 417 KMStrainRate_State(grad_s, strain_rate); 418 NewtonianStress(context, strain_rate, kmstress); 419 KMUnpack(kmstress, stress); 420 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 421 422 StateConservative F_inviscid[3]; 423 FluxInviscid(context, s, F_inviscid); 424 425 CeedScalar Flux[5]; 426 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 427 428 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 429 430 if (is_implicit) { 431 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 432 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 433 } 434 } 435 return 0; 436 } 437 438 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 439 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 440 } 441 442 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 443 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE); 444 } 445 446 CEED_QFUNCTION(BoundaryIntegral_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 447 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_ENTROPY); 448 } 449 450 // ***************************************************************************** 451 // Jacobian for "set nothing" boundary integral 452 // ***************************************************************************** 453 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 454 StateVariable state_var) { 455 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 456 const CeedScalar(*Grad_dq) = in[1]; 457 const CeedScalar(*q_data_sur) = in[2]; 458 const CeedScalar(*jac_data_sur) = in[4]; 459 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 460 461 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 462 const bool is_implicit = context->is_implicit; 463 464 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 465 CeedScalar wdetJb, dXdx[2][3], norm[3]; 466 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 467 wdetJb *= is_implicit ? -1. : 1.; 468 469 CeedScalar qi[5], kmstress[6], dqi[5]; 470 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 471 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 472 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 473 474 State s = StateFromQ(context, qi, state_var); 475 State ds = StateFromQ_fwd(context, s, dqi, state_var); 476 477 State grad_ds[3]; 478 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 479 480 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 481 KMStrainRate_State(grad_ds, dstrain_rate); 482 NewtonianStress(context, dstrain_rate, dkmstress); 483 KMUnpack(dkmstress, dstress); 484 KMUnpack(kmstress, stress); 485 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 486 487 StateConservative dF_inviscid[3]; 488 FluxInviscid_fwd(context, s, ds, dF_inviscid); 489 490 CeedScalar dFlux[5]; 491 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 492 493 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 494 } 495 return 0; 496 } 497 498 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 499 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 500 } 501 502 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 503 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 504 } 505 506 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 507 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY); 508 } 509