xref: /honee/qfunctions/newtonian.h (revision 4b96a86bba11b72cac7c28056075834ab8e47fc9)
1dc936754SJeremy L Thompson // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors.
2727da7e7SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
33a8779fbSJames Wright //
4727da7e7SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause
53a8779fbSJames Wright //
6727da7e7SJeremy L Thompson // This file is part of CEED:  http://github.com/ceed
73a8779fbSJames Wright 
83a8779fbSJames Wright /// @file
93a8779fbSJames Wright /// Operator for Navier-Stokes example using PETSc
103a8779fbSJames Wright #include <ceed.h>
11d0cce58aSJeremy L Thompson #include <math.h>
127b530f2aSAdelekeBankole #include <stdlib.h>
132b916ea7SJeremy L Thompson 
14475b2820SJames Wright #include "newtonian_state.h"
15d0cce58aSJeremy L Thompson #include "newtonian_types.h"
16d1b9ef12SLeila Ghaffari #include "stabilization.h"
17d0cce58aSJeremy L Thompson #include "utils.h"
18bb8a0c61SJames Wright 
1994a7b3d2SKenneth E. Jansen CEED_QFUNCTION_HELPER void InternalDampingLayer(const NewtonianIdealGasContext context, const State s, const CeedScalar sigma, CeedScalar damp_Y[5],
20e7754af5SKenneth E. Jansen                                                 CeedScalar damp_residual[5]) {
21e7754af5SKenneth E. Jansen   ScaleN(damp_Y, sigma, 5);
22edcfef1bSKenneth E. Jansen   State damp_s = StateFromY_fwd(context, s, damp_Y);
23e7754af5SKenneth E. Jansen 
24e7754af5SKenneth E. Jansen   CeedScalar U[5];
25e7754af5SKenneth E. Jansen   UnpackState_U(damp_s.U, U);
26e7754af5SKenneth E. Jansen   for (int i = 0; i < 5; i++) damp_residual[i] += U[i];
27e7754af5SKenneth E. Jansen }
28e7754af5SKenneth E. Jansen 
29bb8a0c61SJames Wright // *****************************************************************************
303a8779fbSJames Wright // This QFunction sets a "still" initial condition for generic Newtonian IG problems
313a8779fbSJames Wright // *****************************************************************************
328fff8293SJames Wright CEED_QFUNCTION_HELPER int ICsNewtonianIG(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
333a8779fbSJames Wright   // Inputs
343a8779fbSJames Wright 
353a8779fbSJames Wright   // Outputs
363a8779fbSJames Wright   CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
373a8779fbSJames Wright 
38bb8a0c61SJames Wright   // Context
39bb8a0c61SJames Wright   const SetupContext context = (SetupContext)ctx;
40bb8a0c61SJames Wright 
413a8779fbSJames Wright   // Quadrature Point Loop
422b916ea7SJeremy L Thompson   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
433a8779fbSJames Wright     CeedScalar q[5] = {0.};
44edcfef1bSKenneth E. Jansen     State      s    = StateFromPrimitive(&context->gas, context->reference);
458fff8293SJames Wright     StateToQ(&context->gas, s, q, state_var);
462b916ea7SJeremy L Thompson     for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
473a8779fbSJames Wright   }  // End of Quadrature Point Loop
483a8779fbSJames Wright   return 0;
493a8779fbSJames Wright }
503a8779fbSJames Wright 
512b916ea7SJeremy L Thompson CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
528fff8293SJames Wright   return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_PRIMITIVE);
53b8fb7609SAdeleke O. Bankole }
54b8fb7609SAdeleke O. Bankole CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
558fff8293SJames Wright   return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
56cbe60e31SLeila Ghaffari }
57cbe60e31SLeila Ghaffari 
58cbe60e31SLeila Ghaffari // *****************************************************************************
5904e40bb6SJeremy L Thompson // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method
603a8779fbSJames Wright //
6104e40bb6SJeremy L Thompson // This is 3D compressible Navier-Stokes in conservation form with state variables of density, momentum density, and total energy density.
623a8779fbSJames Wright //
633a8779fbSJames Wright // State Variables: q = ( rho, U1, U2, U3, E )
643a8779fbSJames Wright //   rho - Mass Density
653a8779fbSJames Wright //   Ui  - Momentum Density,      Ui = rho ui
663a8779fbSJames Wright //   E   - Total Energy Density,  E  = rho (cv T + (u u)/2 + g z)
673a8779fbSJames Wright //
683a8779fbSJames Wright // Navier-Stokes Equations:
693a8779fbSJames Wright //   drho/dt + div( U )                               = 0
703a8779fbSJames Wright //   dU/dt   + div( rho (u x u) + P I3 ) + rho g khat = div( Fu )
713a8779fbSJames Wright //   dE/dt   + div( (E + P) u )                       = div( Fe )
723a8779fbSJames Wright //
733a8779fbSJames Wright // Viscous Stress:
743a8779fbSJames Wright //   Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3)
753a8779fbSJames Wright //
763a8779fbSJames Wright // Thermal Stress:
773a8779fbSJames Wright //   Fe = u Fu + k grad( T )
78bb8a0c61SJames Wright // Equation of State
793a8779fbSJames Wright //   P = (gamma - 1) (E - rho (u u) / 2 - rho g z)
803a8779fbSJames Wright //
813a8779fbSJames Wright // Stabilization:
823a8779fbSJames Wright //   Tau = diag(TauC, TauM, TauM, TauM, TauE)
833a8779fbSJames Wright //     f1 = rho  sqrt(ui uj gij)
843a8779fbSJames Wright //     gij = dXi/dX * dXi/dX
853a8779fbSJames Wright //     TauC = Cc f1 / (8 gii)
863a8779fbSJames Wright //     TauM = min( 1 , 1 / f1 )
873a8779fbSJames Wright //     TauE = TauM / (Ce cv)
883a8779fbSJames Wright //
893a8779fbSJames Wright //  SU   = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) )
903a8779fbSJames Wright //
913a8779fbSJames Wright // Constants:
923a8779fbSJames Wright //   lambda = - 2 / 3,  From Stokes hypothesis
933a8779fbSJames Wright //   mu              ,  Dynamic viscosity
943a8779fbSJames Wright //   k               ,  Thermal conductivity
953a8779fbSJames Wright //   cv              ,  Specific heat, constant volume
963a8779fbSJames Wright //   cp              ,  Specific heat, constant pressure
973a8779fbSJames Wright //   g               ,  Gravity
983a8779fbSJames Wright //   gamma  = cp / cv,  Specific heat ratio
993a8779fbSJames Wright //
10004e40bb6SJeremy L Thompson // 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
10104e40bb6SJeremy L Thompson // gradu )
1023a8779fbSJames Wright // *****************************************************************************
1032b916ea7SJeremy L Thompson CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
1043a8779fbSJames Wright   // Inputs
1053d65b166SJames Wright   const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
10687bd45e7SJames Wright   const CeedScalar(*Grad_q)        = in[1];
107ade49511SJames Wright   const CeedScalar(*q_data)        = in[2];
1083d65b166SJames Wright 
1093a8779fbSJames Wright   // Outputs
1103d65b166SJames Wright   CeedScalar(*v)[CEED_Q_VLA]         = (CeedScalar(*)[CEED_Q_VLA])out[0];
1113d65b166SJames Wright   CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
1123a8779fbSJames Wright 
1133a8779fbSJames Wright   // Context
1143a8779fbSJames Wright   NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
115bb8a0c61SJames Wright   const CeedScalar        *g       = context->g;
116bb8a0c61SJames Wright   const CeedScalar         dt      = context->dt;
1173a8779fbSJames Wright 
1183a8779fbSJames Wright   // Quadrature Point Loop
1193d65b166SJames Wright   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
120ade49511SJames Wright     CeedScalar U[5], wdetJ, dXdx[3][3];
121c1a52365SJed Brown     for (int j = 0; j < 5; j++) U[j] = q[j][i];
122ade49511SJames Wright     StoredValuesUnpack(Q, i, 0, 1, q_data, &wdetJ);
123ade49511SJames Wright     StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx);
124edcfef1bSKenneth E. Jansen     State s = StateFromU(context, U);
125c1a52365SJed Brown 
126c1a52365SJed Brown     State grad_s[3];
127edcfef1bSKenneth E. Jansen     StatePhysicalGradientFromReference(Q, i, context, s, STATEVAR_CONSERVATIVE, Grad_q, dXdx, grad_s);
128c1a52365SJed Brown 
129c1a52365SJed Brown     CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
13040a33f2dSJames Wright     KMStrainRate_State(grad_s, strain_rate);
131c1a52365SJed Brown     NewtonianStress(context, strain_rate, kmstress);
132c1a52365SJed Brown     KMUnpack(kmstress, stress);
133c1a52365SJed Brown     ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe);
134c1a52365SJed Brown 
135c1a52365SJed Brown     StateConservative F_inviscid[3];
136c1a52365SJed Brown     FluxInviscid(context, s, F_inviscid);
137c1a52365SJed Brown 
138c1a52365SJed Brown     // Total flux
139c1a52365SJed Brown     CeedScalar Flux[5][3];
140d1b9ef12SLeila Ghaffari     FluxTotal(F_inviscid, stress, Fe, Flux);
141c1a52365SJed Brown 
1427523f6aaSJames Wright     for (CeedInt j = 0; j < 5; j++) {
1437523f6aaSJames Wright       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]);
1442b916ea7SJeremy L Thompson     }
145c1a52365SJed Brown 
14660dbb574SKenneth E. Jansen     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)};
1472b916ea7SJeremy L Thompson     for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j];
1483a8779fbSJames Wright 
149d1b9ef12SLeila Ghaffari     // -- Stabilization method: none (Galerkin), SU, or SUPG
150d1b9ef12SLeila Ghaffari     CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0};
151d1b9ef12SLeila Ghaffari     Tau_diagPrim(context, s, dXdx, dt, Tau_d);
152edcfef1bSKenneth E. Jansen     Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab);
1533a8779fbSJames Wright 
1542b916ea7SJeremy L Thompson     for (CeedInt j = 0; j < 5; j++) {
1552b916ea7SJeremy L Thompson       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]);
1562b916ea7SJeremy L Thompson     }
1573a8779fbSJames Wright   }  // End Quadrature Point Loop
1583a8779fbSJames Wright 
1593a8779fbSJames Wright   // Return
1603a8779fbSJames Wright   return 0;
1613a8779fbSJames Wright }
1623a8779fbSJames Wright 
1633a8779fbSJames Wright // *****************************************************************************
16404e40bb6SJeremy L Thompson // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepping method
1653a8779fbSJames Wright //
1663a8779fbSJames Wright //  SU   = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) )
1673a8779fbSJames Wright //  SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) )
16804e40bb6SJeremy L Thompson //                                       (diffusive terms will be added later)
1693a8779fbSJames Wright // *****************************************************************************
1708fff8293SJames Wright CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
1713a8779fbSJames Wright   // Inputs
1723d65b166SJames Wright   const CeedScalar(*q)[CEED_Q_VLA]     = (const CeedScalar(*)[CEED_Q_VLA])in[0];
17387bd45e7SJames Wright   const CeedScalar(*Grad_q)            = in[1];
1743d65b166SJames Wright   const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
175ade49511SJames Wright   const CeedScalar(*q_data)            = in[3];
1763d65b166SJames Wright   const CeedScalar(*x)[CEED_Q_VLA]     = (const CeedScalar(*)[CEED_Q_VLA])in[4];
1773d65b166SJames Wright 
1783a8779fbSJames Wright   // Outputs
1793d65b166SJames Wright   CeedScalar(*v)[CEED_Q_VLA]         = (CeedScalar(*)[CEED_Q_VLA])out[0];
1803d65b166SJames Wright   CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
181ade49511SJames Wright   CeedScalar(*jac_data)              = out[2];
1823d65b166SJames Wright 
1833a8779fbSJames Wright   // Context
1843a8779fbSJames Wright   NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
185bb8a0c61SJames Wright   const CeedScalar        *g       = context->g;
186bb8a0c61SJames Wright   const CeedScalar         dt      = context->dt;
187e7754af5SKenneth E. Jansen   const CeedScalar         P0      = context->P0;
1883a8779fbSJames Wright 
1893a8779fbSJames Wright   // Quadrature Point Loop
1903d65b166SJames Wright   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
1913d65b166SJames Wright     const CeedScalar qi[5]  = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
192c1a52365SJed Brown     const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
193edcfef1bSKenneth E. Jansen     const State      s      = StateFromQ(context, qi, state_var);
194c1a52365SJed Brown 
195ade49511SJames Wright     CeedScalar wdetJ, dXdx[3][3];
196ade49511SJames Wright     QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
197c1a52365SJed Brown     State grad_s[3];
198edcfef1bSKenneth E. Jansen     StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s);
199c1a52365SJed Brown 
200c1a52365SJed Brown     CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
20140a33f2dSJames Wright     KMStrainRate_State(grad_s, strain_rate);
202c1a52365SJed Brown     NewtonianStress(context, strain_rate, kmstress);
203c1a52365SJed Brown     KMUnpack(kmstress, stress);
204c1a52365SJed Brown     ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe);
205c1a52365SJed Brown 
206c1a52365SJed Brown     StateConservative F_inviscid[3];
207c1a52365SJed Brown     FluxInviscid(context, s, F_inviscid);
208c1a52365SJed Brown 
209c1a52365SJed Brown     // Total flux
210c1a52365SJed Brown     CeedScalar Flux[5][3];
211d1b9ef12SLeila Ghaffari     FluxTotal(F_inviscid, stress, Fe, Flux);
212c1a52365SJed Brown 
2137523f6aaSJames Wright     for (CeedInt j = 0; j < 5; j++) {
2147523f6aaSJames Wright       for (CeedInt k = 0; k < 3; k++) {
2157523f6aaSJames Wright         Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]);
2163d65b166SJames Wright       }
2172b916ea7SJeremy L Thompson     }
218c1a52365SJed Brown 
21960dbb574SKenneth E. Jansen     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)};
2203a8779fbSJames Wright 
221d1b9ef12SLeila Ghaffari     // -- Stabilization method: none (Galerkin), SU, or SUPG
222edcfef1bSKenneth E. Jansen     CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5];
22376555becSJames Wright     for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i];
224edcfef1bSKenneth E. Jansen     State s_dot = StateFromQ_fwd(context, s, qi_dot, state_var);
22576555becSJames Wright     UnpackState_U(s_dot.U, U_dot);
22676555becSJames Wright 
2272b916ea7SJeremy L Thompson     for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]);
228e7754af5SKenneth E. Jansen     if (context->idl_enable) {
22994a7b3d2SKenneth E. Jansen       const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]);
23094a7b3d2SKenneth E. Jansen       StoredValuesPack(Q, i, 14, 1, &sigma, jac_data);
231e7754af5SKenneth E. Jansen       CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.};
23294a7b3d2SKenneth E. Jansen       InternalDampingLayer(context, s, sigma, damp_state, idl_residual);
233e7754af5SKenneth E. Jansen       for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j];
234e7754af5SKenneth E. Jansen     }
235e7754af5SKenneth E. Jansen 
236d1b9ef12SLeila Ghaffari     Tau_diagPrim(context, s, dXdx, dt, Tau_d);
237edcfef1bSKenneth E. Jansen     Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab);
2383a8779fbSJames Wright 
2392b916ea7SJeremy L Thompson     for (CeedInt j = 0; j < 5; j++) {
2403d65b166SJames Wright       for (CeedInt k = 0; k < 3; k++) {
2413d65b166SJames Wright         Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]);
2423d65b166SJames Wright       }
2432b916ea7SJeremy L Thompson     }
244ade49511SJames Wright     StoredValuesPack(Q, i, 0, 5, qi, jac_data);
245ade49511SJames Wright     StoredValuesPack(Q, i, 5, 6, kmstress, jac_data);
246ade49511SJames Wright     StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data);
2473a8779fbSJames Wright 
2483a8779fbSJames Wright   }  // End Quadrature Point Loop
2493a8779fbSJames Wright 
2503a8779fbSJames Wright   // Return
2513a8779fbSJames Wright   return 0;
2523a8779fbSJames Wright }
253f0b65372SJed Brown 
2542b916ea7SJeremy L Thompson CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
2558fff8293SJames Wright   return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
25676555becSJames Wright }
25776555becSJames Wright 
2582b916ea7SJeremy L Thompson CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
2598fff8293SJames Wright   return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
26076555becSJames Wright }
26176555becSJames Wright 
262cbe60e31SLeila Ghaffari // *****************************************************************************
26304e40bb6SJeremy L Thompson // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method.
264cbe60e31SLeila Ghaffari // *****************************************************************************
2658fff8293SJames Wright CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
266f0b65372SJed Brown   // Inputs
2673d65b166SJames Wright   const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
26887bd45e7SJames Wright   const CeedScalar(*Grad_dq)        = in[1];
269ade49511SJames Wright   const CeedScalar(*q_data)         = in[2];
27094a7b3d2SKenneth E. Jansen   const CeedScalar(*jac_data)       = in[3];
2713d65b166SJames Wright 
272f0b65372SJed Brown   // Outputs
2733d65b166SJames Wright   CeedScalar(*v)[CEED_Q_VLA]         = (CeedScalar(*)[CEED_Q_VLA])out[0];
2743d65b166SJames Wright   CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
2753d65b166SJames Wright 
276f0b65372SJed Brown   // Context
277f0b65372SJed Brown   NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
278f0b65372SJed Brown   const CeedScalar        *g       = context->g;
279f0b65372SJed Brown 
280f0b65372SJed Brown   // Quadrature Point Loop
2813d65b166SJames Wright   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
282ade49511SJames Wright     CeedScalar wdetJ, dXdx[3][3];
283ade49511SJames Wright     QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
284f0b65372SJed Brown 
2858789e95fSJames Wright     CeedScalar qi[5], kmstress[6], Tau_d[3];
286ade49511SJames Wright     StoredValuesUnpack(Q, i, 0, 5, jac_data, qi);
287ade49511SJames Wright     StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress);
288ade49511SJames Wright     StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d);
289edcfef1bSKenneth E. Jansen     State s = StateFromQ(context, qi, state_var);
290f0b65372SJed Brown 
291edcfef1bSKenneth E. Jansen     CeedScalar dqi[5];
29276555becSJames Wright     for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];
293edcfef1bSKenneth E. Jansen     State ds = StateFromQ_fwd(context, s, dqi, state_var);
294f0b65372SJed Brown 
295f0b65372SJed Brown     State grad_ds[3];
296edcfef1bSKenneth E. Jansen     StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds);
297f0b65372SJed Brown 
298f0b65372SJed Brown     CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3];
29940a33f2dSJames Wright     KMStrainRate_State(grad_ds, dstrain_rate);
300f0b65372SJed Brown     NewtonianStress(context, dstrain_rate, dkmstress);
301f0b65372SJed Brown     KMUnpack(dkmstress, dstress);
302f0b65372SJed Brown     KMUnpack(kmstress, stress);
303f0b65372SJed Brown     ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe);
304f0b65372SJed Brown 
305f0b65372SJed Brown     StateConservative dF_inviscid[3];
306f0b65372SJed Brown     FluxInviscid_fwd(context, s, ds, dF_inviscid);
307f0b65372SJed Brown 
308f0b65372SJed Brown     // Total flux
309f0b65372SJed Brown     CeedScalar dFlux[5][3];
310d1b9ef12SLeila Ghaffari     FluxTotal(dF_inviscid, dstress, dFe, dFlux);
311f0b65372SJed Brown 
31222387d3aSJames Wright     for (int j = 0; j < 5; j++) {
31322387d3aSJames Wright       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]);
3142b916ea7SJeremy L Thompson     }
315f0b65372SJed Brown 
31660dbb574SKenneth E. Jansen     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)};
31776555becSJames Wright     CeedScalar       dU[5]          = {0.};
31876555becSJames Wright     UnpackState_U(ds.U, dU);
3192b916ea7SJeremy L Thompson     for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]);
320f0b65372SJed Brown 
321e7754af5SKenneth E. Jansen     if (context->idl_enable) {
32294a7b3d2SKenneth E. Jansen       const CeedScalar sigma         = jac_data[14 * Q + i];
323e7754af5SKenneth E. Jansen       CeedScalar       damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.};
324e7754af5SKenneth E. Jansen       // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds.
32594a7b3d2SKenneth E. Jansen       InternalDampingLayer(context, s, sigma, damp_state, idl_residual);
326e7754af5SKenneth E. Jansen       for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j];
327e7754af5SKenneth E. Jansen     }
328e7754af5SKenneth E. Jansen 
329d1b9ef12SLeila Ghaffari     // -- Stabilization method: none (Galerkin), SU, or SUPG
330d1b9ef12SLeila Ghaffari     CeedScalar dstab[5][3], U_dot[5] = {0};
331d1b9ef12SLeila Ghaffari     for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j];
332edcfef1bSKenneth E. Jansen     Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, dstab);
333d1b9ef12SLeila Ghaffari 
3342b916ea7SJeremy L Thompson     for (int j = 0; j < 5; j++) {
3352b916ea7SJeremy L Thompson       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]);
3362b916ea7SJeremy L Thompson     }
337f0b65372SJed Brown   }  // End Quadrature Point Loop
338f0b65372SJed Brown   return 0;
339f0b65372SJed Brown }
3408085925cSJames Wright 
3412b916ea7SJeremy L Thompson CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
3428fff8293SJames Wright   return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
34376555becSJames Wright }
34476555becSJames Wright 
3452b916ea7SJeremy L Thompson CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
3468fff8293SJames Wright   return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
34776555becSJames Wright }
34876555becSJames Wright 
349d1b9ef12SLeila Ghaffari // *****************************************************************************
3508085925cSJames Wright // Compute boundary integral (ie. for strongly set inflows)
351d1b9ef12SLeila Ghaffari // *****************************************************************************
3528fff8293SJames Wright CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
353*4b96a86bSJames Wright   const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
3543d65b166SJames Wright   const CeedScalar(*q)[CEED_Q_VLA]       = (const CeedScalar(*)[CEED_Q_VLA])in[0];
35587bd45e7SJames Wright   const CeedScalar(*Grad_q)              = in[1];
356ade49511SJames Wright   const CeedScalar(*q_data_sur)          = in[2];
3573d65b166SJames Wright   CeedScalar(*v)[CEED_Q_VLA]             = (CeedScalar(*)[CEED_Q_VLA])out[0];
358*4b96a86bSJames Wright   CeedScalar(*jac_data_sur)              = context->is_implicit ? out[1] : NULL;
3598085925cSJames Wright 
360d3b25f3aSJames Wright   const bool is_implicit = context->is_implicit;
3618085925cSJames Wright 
3622b916ea7SJeremy L Thompson   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
36341e73928SJames Wright     const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
364edcfef1bSKenneth E. Jansen     State            s     = StateFromQ(context, qi, state_var);
3658085925cSJames Wright 
366ade49511SJames Wright     CeedScalar wdetJb, dXdx[2][3], norm[3];
367ade49511SJames Wright     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm);
368ade49511SJames Wright     wdetJb *= is_implicit ? -1. : 1.;
3698085925cSJames Wright 
370d3b25f3aSJames Wright     State grad_s[3];
371edcfef1bSKenneth E. Jansen     StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_q, dXdx, grad_s);
3728085925cSJames Wright 
373d3b25f3aSJames Wright     CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
37440a33f2dSJames Wright     KMStrainRate_State(grad_s, strain_rate);
375d3b25f3aSJames Wright     NewtonianStress(context, strain_rate, kmstress);
376d3b25f3aSJames Wright     KMUnpack(kmstress, stress);
377d3b25f3aSJames Wright     ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe);
378d3b25f3aSJames Wright 
379d3b25f3aSJames Wright     StateConservative F_inviscid[3];
380d3b25f3aSJames Wright     FluxInviscid(context, s, F_inviscid);
381d3b25f3aSJames Wright 
382c5740391SJames Wright     CeedScalar Flux[5];
383c5740391SJames Wright     FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux);
384d3b25f3aSJames Wright 
385c5740391SJames Wright     for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];
3868085925cSJames Wright 
387*4b96a86bSJames Wright     if (is_implicit) {
388ade49511SJames Wright       StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
389ade49511SJames Wright       StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
3908085925cSJames Wright     }
391*4b96a86bSJames Wright   }
3928085925cSJames Wright   return 0;
3938085925cSJames Wright }
3948085925cSJames Wright 
3952b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
3968fff8293SJames Wright   return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
397d4559bbeSJames Wright }
398d4559bbeSJames Wright 
3992b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
4008fff8293SJames Wright   return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE);
401d4559bbeSJames Wright }
402d4559bbeSJames Wright 
403d1b9ef12SLeila Ghaffari // *****************************************************************************
40468ae065aSJames Wright // Jacobian for "set nothing" boundary integral
405d1b9ef12SLeila Ghaffari // *****************************************************************************
4062b916ea7SJeremy L Thompson CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
4078fff8293SJames Wright                                                     StateVariable state_var) {
40868ae065aSJames Wright   // Inputs
4093d65b166SJames Wright   const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
41087bd45e7SJames Wright   const CeedScalar(*Grad_dq)        = in[1];
411ade49511SJames Wright   const CeedScalar(*q_data_sur)     = in[2];
412c1484fadSKenneth E. Jansen   const CeedScalar(*jac_data_sur)   = in[4];
4133d65b166SJames Wright 
41468ae065aSJames Wright   // Outputs
41568ae065aSJames Wright   CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
41668ae065aSJames Wright 
41768ae065aSJames Wright   const NewtonianIdealGasContext context     = (NewtonianIdealGasContext)ctx;
418ade49511SJames Wright   const bool                     is_implicit = context->is_implicit;
41968ae065aSJames Wright 
42068ae065aSJames Wright   // Quadrature Point Loop
4213d65b166SJames Wright   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
422ade49511SJames Wright     CeedScalar wdetJb, dXdx[2][3], norm[3];
423ade49511SJames Wright     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm);
424ade49511SJames Wright     wdetJb *= is_implicit ? -1. : 1.;
42568ae065aSJames Wright 
426edcfef1bSKenneth E. Jansen     CeedScalar qi[5], kmstress[6], dqi[5];
427ade49511SJames Wright     StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi);
428ade49511SJames Wright     StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress);
42941e73928SJames Wright     for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];
4303934e2b1SJames Wright 
431edcfef1bSKenneth E. Jansen     State s  = StateFromQ(context, qi, state_var);
432edcfef1bSKenneth E. Jansen     State ds = StateFromQ_fwd(context, s, dqi, state_var);
43368ae065aSJames Wright 
43468ae065aSJames Wright     State grad_ds[3];
435edcfef1bSKenneth E. Jansen     StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds);
43668ae065aSJames Wright 
43768ae065aSJames Wright     CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3];
43840a33f2dSJames Wright     KMStrainRate_State(grad_ds, dstrain_rate);
43968ae065aSJames Wright     NewtonianStress(context, dstrain_rate, dkmstress);
44068ae065aSJames Wright     KMUnpack(dkmstress, dstress);
44168ae065aSJames Wright     KMUnpack(kmstress, stress);
44268ae065aSJames Wright     ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe);
44368ae065aSJames Wright 
44468ae065aSJames Wright     StateConservative dF_inviscid[3];
44568ae065aSJames Wright     FluxInviscid_fwd(context, s, ds, dF_inviscid);
44668ae065aSJames Wright 
447c5740391SJames Wright     CeedScalar dFlux[5];
448c5740391SJames Wright     FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux);
44968ae065aSJames Wright 
450c5740391SJames Wright     for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j];
45168ae065aSJames Wright   }  // End Quadrature Point Loop
45268ae065aSJames Wright   return 0;
45368ae065aSJames Wright }
45468ae065aSJames Wright 
4552b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
4568fff8293SJames Wright   return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
457d4559bbeSJames Wright }
458d4559bbeSJames Wright 
4592b916ea7SJeremy L Thompson CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
4608fff8293SJames Wright   return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
461d4559bbeSJames Wright }
462