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