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