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