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