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