xref: /honee/qfunctions/bc_freestream.h (revision 8a8cb6e06ce4728cc6d80ca92f8de31da49852e5)
1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 
4 /// @file
5 /// QFunctions for the `bc_freestream` and `bc_outflow` boundary conditions
6 #include "bc_freestream_type.h"
7 #include "newtonian_state.h"
8 #include "newtonian_types.h"
9 #include "riemann_solver.h"
10 
11 // *****************************************************************************
12 // Freestream Boundary Condition
13 // *****************************************************************************
14 CEED_QFUNCTION_HELPER int Freestream(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var,
15                                      RiemannFluxType flux_type) {
16   const FreestreamContext context  = (FreestreamContext)ctx;
17   const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
18   const CeedScalar(*q_data_sur)    = in[2];
19   CeedScalar(*v)[CEED_Q_VLA]       = (CeedScalar(*)[CEED_Q_VLA])out[0];
20   CeedScalar(*jac_data_sur)        = context->newtonian_ctx.is_implicit ? out[1] : NULL;
21 
22   const NewtonianIdealGasContext newt_ctx    = &context->newtonian_ctx;
23   const bool                     is_implicit = newt_ctx->is_implicit;
24 
25   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
26     const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
27     const State      s     = StateFromQ(newt_ctx, qi, state_var);
28 
29     CeedScalar wdetJb, normal[3];
30     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
31     wdetJb *= is_implicit ? -1. : 1.;
32 
33     StateConservative flux;
34     switch (flux_type) {
35       case RIEMANN_HLL:
36         flux = RiemannFlux_HLL(newt_ctx, s, context->S_infty, normal);
37         break;
38       case RIEMANN_HLLC:
39         flux = RiemannFlux_HLLC(newt_ctx, s, context->S_infty, normal);
40         break;
41     }
42     CeedScalar Flux[5];
43     UnpackState_U(flux, Flux);
44     for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];
45 
46     if (is_implicit) {
47       CeedScalar zeros[6] = {0.};
48       StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
49       StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur);  // Every output value must be set
50     }
51   }
52   return 0;
53 }
54 
55 CEED_QFUNCTION(Freestream_Conserv_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
56   return Freestream(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLL);
57 }
58 
59 CEED_QFUNCTION(Freestream_Prim_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
60   return Freestream(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLL);
61 }
62 
63 CEED_QFUNCTION(Freestream_Entropy_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
64   return Freestream(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLL);
65 }
66 
67 CEED_QFUNCTION(Freestream_Conserv_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
68   return Freestream(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLLC);
69 }
70 
71 CEED_QFUNCTION(Freestream_Prim_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
72   return Freestream(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLLC);
73 }
74 
75 CEED_QFUNCTION(Freestream_Entropy_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
76   return Freestream(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLLC);
77 }
78 
79 CEED_QFUNCTION_HELPER int Freestream_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var,
80                                               RiemannFluxType flux_type) {
81   const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
82   const CeedScalar(*q_data_sur)     = in[2];
83   const CeedScalar(*jac_data_sur)   = in[4];
84 
85   CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
86 
87   const FreestreamContext        context     = (FreestreamContext)ctx;
88   const NewtonianIdealGasContext newt_ctx    = &context->newtonian_ctx;
89   const bool                     is_implicit = newt_ctx->is_implicit;
90   const State                    dS_infty    = {0};
91 
92   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
93     CeedScalar wdetJb, normal[3];
94     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
95     wdetJb *= is_implicit ? -1. : 1.;
96 
97     CeedScalar qi[5], dqi[5];
98     StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi);
99     for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];
100     State s  = StateFromQ(newt_ctx, qi, state_var);
101     State ds = StateFromQ_fwd(newt_ctx, s, dqi, state_var);
102 
103     StateConservative dflux;
104     switch (flux_type) {
105       case RIEMANN_HLL:
106         dflux = RiemannFlux_HLL_fwd(newt_ctx, s, ds, context->S_infty, dS_infty, normal);
107         break;
108       case RIEMANN_HLLC:
109         dflux = RiemannFlux_HLLC_fwd(newt_ctx, s, ds, context->S_infty, dS_infty, normal);
110         break;
111     }
112     CeedScalar dFlux[5];
113     UnpackState_U(dflux, dFlux);
114     for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j];
115   }
116   return 0;
117 }
118 
119 CEED_QFUNCTION(Freestream_Jacobian_Conserv_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
120   return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLL);
121 }
122 
123 CEED_QFUNCTION(Freestream_Jacobian_Prim_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
124   return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLL);
125 }
126 
127 CEED_QFUNCTION(Freestream_Jacobian_Entropy_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
128   return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLL);
129 }
130 
131 CEED_QFUNCTION(Freestream_Jacobian_Conserv_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
132   return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLLC);
133 }
134 
135 CEED_QFUNCTION(Freestream_Jacobian_Prim_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
136   return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLLC);
137 }
138 
139 CEED_QFUNCTION(Freestream_Jacobian_Entropy_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
140   return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLLC);
141 }
142 
143 // Note the identity
144 //
145 // softplus(x) - x = log(1 + exp(x)) - x
146 //                 = log(1 + exp(x)) + log(exp(-x))
147 //                 = log((1 + exp(x)) * exp(-x))
148 //                 = log(exp(-x) + 1)
149 //                 = softplus(-x)
150 CEED_QFUNCTION_HELPER CeedScalar Softplus(CeedScalar x, CeedScalar width) {
151   if (x > 40 * width) return x;
152   return width * log1p(exp(x / width));
153 }
154 
155 CEED_QFUNCTION_HELPER CeedScalar Softplus_fwd(CeedScalar x, CeedScalar dx, CeedScalar width) {
156   if (x > 40 * width) return 1;
157   const CeedScalar t = exp(x / width);
158   return t / (1 + t);
159 }
160 
161 // Viscous Outflow boundary condition, setting a constant exterior pressure and
162 // temperature as input for a Riemann solve. This condition is stable even in
163 // recirculating flow so long as the exterior temperature is sensible.
164 //
165 // The velocity in the exterior state has optional softplus regularization to
166 // keep it outflow. These parameters have been finnicky in practice and provide
167 // little or no benefit in the tests we've run thus far, thus we recommend
168 // skipping this feature and just allowing recirculation.
169 CEED_QFUNCTION_HELPER int RiemannOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
170   const OutflowContext outflow     = (OutflowContext)ctx;
171   const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
172   const CeedScalar(*Grad_q)        = in[1];
173   const CeedScalar(*q_data_sur)    = in[2];
174   CeedScalar(*v)[CEED_Q_VLA]       = (CeedScalar(*)[CEED_Q_VLA])out[0];
175   CeedScalar(*jac_data_sur)        = outflow->gas.is_implicit ? out[1] : NULL;
176 
177   const NewtonianIdealGasContext gas         = &outflow->gas;
178   const bool                     is_implicit = gas->is_implicit;
179 
180   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
181     CeedScalar wdetJb, dXdx[2][3], normal[3];
182     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal);
183     wdetJb *= is_implicit ? -1. : 1.;
184     const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
185     const State      s_int = StateFromQ(gas, qi, state_var);
186 
187     StatePrimitive y_ext      = s_int.Y;
188     y_ext.pressure            = outflow->pressure;
189     y_ext.temperature         = outflow->temperature;
190     const CeedScalar u_normal = Dot3(y_ext.velocity, normal);
191     const CeedScalar proj     = (1 - outflow->recirc) * Softplus(-u_normal, outflow->softplus_velocity);
192     for (CeedInt j = 0; j < 3; j++) {
193       y_ext.velocity[j] += normal[j] * proj;  // (I - n n^T) projects into the plane tangent to the normal
194     }
195     State s_ext = StateFromPrimitive(gas, y_ext);
196 
197     State grad_s[3];
198     StatePhysicalGradientFromReference_Boundary(Q, i, gas, s_int, 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(gas, strain_rate, kmstress);
203     KMUnpack(kmstress, stress);
204     ViscousEnergyFlux(gas, s_int.Y, grad_s, stress, Fe);
205 
206     StateConservative F_inviscid_normal = RiemannFlux_HLLC(gas, s_int, s_ext, normal);
207 
208     CeedScalar Flux[5];
209     FluxTotal_RiemannBoundary(F_inviscid_normal, stress, Fe, normal, Flux);
210 
211     for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];
212 
213     // Save values for Jacobian
214     if (is_implicit) {
215       StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
216       StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
217     }
218   }
219   return 0;
220 }
221 
222 CEED_QFUNCTION(RiemannOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
223   return RiemannOutflow(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
224 }
225 
226 CEED_QFUNCTION(RiemannOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
227   return RiemannOutflow(ctx, Q, in, out, STATEVAR_PRIMITIVE);
228 }
229 
230 CEED_QFUNCTION(RiemannOutflow_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
231   return RiemannOutflow(ctx, Q, in, out, STATEVAR_ENTROPY);
232 }
233 
234 // *****************************************************************************
235 // Jacobian for Riemann pressure/temperature outflow boundary condition
236 // *****************************************************************************
237 CEED_QFUNCTION_HELPER int RiemannOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
238                                                   StateVariable state_var) {
239   const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
240   const CeedScalar(*Grad_dq)        = in[1];
241   const CeedScalar(*q_data_sur)     = in[2];
242   const CeedScalar(*jac_data_sur)   = in[4];
243   CeedScalar(*v)[CEED_Q_VLA]        = (CeedScalar(*)[CEED_Q_VLA])out[0];
244 
245   const OutflowContext           outflow     = (OutflowContext)ctx;
246   const NewtonianIdealGasContext gas         = &outflow->gas;
247   const bool                     is_implicit = gas->is_implicit;
248 
249   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
250     CeedScalar wdetJb, dXdx[2][3], normal[3];
251     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal);
252     wdetJb *= is_implicit ? -1. : 1.;
253 
254     CeedScalar qi[5], kmstress[6], dqi[5];
255     StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi);
256     StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress);
257     for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];
258 
259     State          s_int  = StateFromQ(gas, qi, state_var);
260     const State    ds_int = StateFromQ_fwd(gas, s_int, dqi, state_var);
261     StatePrimitive y_ext = s_int.Y, dy_ext = ds_int.Y;
262     y_ext.pressure             = outflow->pressure;
263     y_ext.temperature          = outflow->temperature;
264     dy_ext.pressure            = 0;
265     dy_ext.temperature         = 0;
266     const CeedScalar u_normal  = Dot3(s_int.Y.velocity, normal);
267     const CeedScalar du_normal = Dot3(ds_int.Y.velocity, normal);
268     const CeedScalar proj      = (1 - outflow->recirc) * Softplus(-u_normal, outflow->softplus_velocity);
269     const CeedScalar dproj     = (1 - outflow->recirc) * Softplus_fwd(-u_normal, -du_normal, outflow->softplus_velocity);
270     for (CeedInt j = 0; j < 3; j++) {
271       y_ext.velocity[j] += normal[j] * proj;
272       dy_ext.velocity[j] += normal[j] * dproj;
273     }
274 
275     State s_ext  = StateFromPrimitive(gas, y_ext);
276     State ds_ext = StateFromPrimitive_fwd(gas, s_ext, dy_ext);
277 
278     State grad_ds[3];
279     StatePhysicalGradientFromReference_Boundary(Q, i, gas, s_int, state_var, Grad_dq, dXdx, grad_ds);
280 
281     CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3];
282     KMStrainRate_State(grad_ds, dstrain_rate);
283     NewtonianStress(gas, dstrain_rate, dkmstress);
284     KMUnpack(dkmstress, dstress);
285     KMUnpack(kmstress, stress);
286     ViscousEnergyFlux_fwd(gas, s_int.Y, ds_int.Y, grad_ds, stress, dstress, dFe);
287 
288     StateConservative dF_inviscid_normal = RiemannFlux_HLLC_fwd(gas, s_int, ds_int, s_ext, ds_ext, normal);
289 
290     CeedScalar dFlux[5];
291     FluxTotal_RiemannBoundary(dF_inviscid_normal, dstress, dFe, normal, dFlux);
292 
293     for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j];
294   }
295   return 0;
296 }
297 
298 CEED_QFUNCTION(RiemannOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
299   return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
300 }
301 
302 CEED_QFUNCTION(RiemannOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
303   return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
304 }
305 
306 CEED_QFUNCTION(RiemannOutflow_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
307   return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY);
308 }
309 
310 // *****************************************************************************
311 // Outflow boundary condition, weakly setting a constant pressure. This is the
312 // classic outflow condition used by PHASTA-C and retained largely for
313 // comparison. In our experiments, it is never better than RiemannOutflow, and
314 // will crash if outflow ever becomes an inflow, as occurs with strong
315 // acoustics, vortices, etc.
316 // *****************************************************************************
317 CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
318   const OutflowContext outflow     = (OutflowContext)ctx;
319   const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
320   const CeedScalar(*Grad_q)        = in[1];
321   const CeedScalar(*q_data_sur)    = in[2];
322   CeedScalar(*v)[CEED_Q_VLA]       = (CeedScalar(*)[CEED_Q_VLA])out[0];
323   CeedScalar(*jac_data_sur)        = outflow->gas.is_implicit ? out[1] : NULL;
324 
325   const NewtonianIdealGasContext gas         = &outflow->gas;
326   const bool                     is_implicit = gas->is_implicit;
327 
328   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
329     const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
330     State            s     = StateFromQ(gas, qi, state_var);
331     s.Y.pressure           = outflow->pressure;
332 
333     CeedScalar wdetJb, dXdx[2][3], normal[3];
334     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal);
335     wdetJb *= is_implicit ? -1. : 1.;
336 
337     State grad_s[3];
338     StatePhysicalGradientFromReference_Boundary(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s);
339 
340     CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
341     KMStrainRate_State(grad_s, strain_rate);
342     NewtonianStress(gas, strain_rate, kmstress);
343     KMUnpack(kmstress, stress);
344     ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe);
345 
346     StateConservative F_inviscid[3];
347     FluxInviscid(gas, s, F_inviscid);
348 
349     CeedScalar Flux[5];
350     FluxTotal_Boundary(F_inviscid, stress, Fe, normal, Flux);
351 
352     for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];
353 
354     // Save values for Jacobian
355     if (is_implicit) {
356       StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
357       StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
358     }
359   }
360   return 0;
361 }
362 
363 CEED_QFUNCTION(PressureOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
364   return PressureOutflow(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
365 }
366 
367 CEED_QFUNCTION(PressureOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
368   return PressureOutflow(ctx, Q, in, out, STATEVAR_PRIMITIVE);
369 }
370 
371 CEED_QFUNCTION(PressureOutflow_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
372   return PressureOutflow(ctx, Q, in, out, STATEVAR_ENTROPY);
373 }
374 
375 // *****************************************************************************
376 // Jacobian for weak-pressure outflow boundary condition
377 // *****************************************************************************
378 CEED_QFUNCTION_HELPER int PressureOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
379                                                    StateVariable state_var) {
380   const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
381   const CeedScalar(*Grad_dq)        = in[1];
382   const CeedScalar(*q_data_sur)     = in[2];
383   const CeedScalar(*jac_data_sur)   = in[4];
384   CeedScalar(*v)[CEED_Q_VLA]        = (CeedScalar(*)[CEED_Q_VLA])out[0];
385 
386   const OutflowContext           outflow     = (OutflowContext)ctx;
387   const NewtonianIdealGasContext gas         = &outflow->gas;
388   const bool                     is_implicit = gas->is_implicit;
389 
390   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
391     CeedScalar wdetJb, dXdx[2][3], normal[3];
392     QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal);
393     wdetJb *= is_implicit ? -1. : 1.;
394 
395     CeedScalar qi[5], kmstress[6], dqi[5];
396     StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi);
397     StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress);
398     for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];
399 
400     State s       = StateFromQ(gas, qi, state_var);
401     State ds      = StateFromQ_fwd(gas, s, dqi, state_var);
402     s.Y.pressure  = outflow->pressure;
403     ds.Y.pressure = 0.;
404 
405     State grad_ds[3];
406     StatePhysicalGradientFromReference_Boundary(Q, i, gas, 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(gas, dstrain_rate, dkmstress);
411     KMUnpack(dkmstress, dstress);
412     KMUnpack(kmstress, stress);
413     ViscousEnergyFlux_fwd(gas, s.Y, ds.Y, grad_ds, stress, dstress, dFe);
414 
415     StateConservative dF_inviscid[3];
416     FluxInviscid_fwd(gas, s, ds, dF_inviscid);
417 
418     CeedScalar dFlux[5];
419     FluxTotal_Boundary(dF_inviscid, dstress, dFe, normal, 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(PressureOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
427   return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
428 }
429 
430 CEED_QFUNCTION(PressureOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
431   return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
432 }
433 
434 CEED_QFUNCTION(PressureOutflow_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
435   return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY);
436 }
437