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