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