1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Newtonian fluids operator for HONEE
6 #include <ceed/types.h>
7
8 #include "newtonian_state.h"
9 #include "newtonian_types.h"
10 #include "stabilization.h"
11 #include "utils.h"
12
13 // *****************************************************************************
14 // This QFunction sets a "still" initial condition for generic Newtonian IG problems
15 // *****************************************************************************
ICsNewtonianIG(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)16 CEED_QFUNCTION_HELPER int ICsNewtonianIG(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
17 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
18
19 const SetupContext context = (SetupContext)ctx;
20 NewtonianIGProperties gas = context->newt_ctx.gas;
21
22 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
23 CeedScalar q[5];
24 State s = StateFromPrimitive(gas, context->reference);
25 StateToQ(gas, s, q, state_var);
26 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
27 }
28 return 0;
29 }
30
ICsNewtonianIG_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)31 CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
32 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
33 }
34
ICsNewtonianIG_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)35 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
36 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_PRIMITIVE);
37 }
38
ICsNewtonianIG_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)39 CEED_QFUNCTION(ICsNewtonianIG_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
40 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_ENTROPY);
41 }
42
MassFunction_Newtonian(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)43 CEED_QFUNCTION_HELPER int MassFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
44 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
45 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
46 const CeedScalar(*q_data) = in[2];
47 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
48 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
49
50 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
51 NewtonianIGProperties gas = context->gas;
52
53 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
54 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
55 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]};
56 const State s = StateFromQ(gas, qi, state_var);
57 const State s_dot = StateFromQ(gas, qi_dot, state_var);
58 CeedScalar wdetJ, dXdx[3][3];
59 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
60
61 // Standard mass matrix term
62 for (CeedInt f = 0; f < 5; f++) {
63 v[f][i] = wdetJ * qi_dot[f];
64 }
65
66 // Stabilization method: none (Galerkin), SU, or SUPG
67 State grad_s[3] = {{{0.}}};
68 CeedScalar Tau_d[3], stab[5][3], body_force[5] = {0.}, divFdiff[5] = {0.}, U_dot[5];
69 UnpackState_U(s_dot.U, U_dot);
70 Tau_diagPrim(context->tau_coeffs, gas, s, dXdx, context->dt, Tau_d);
71 Stabilization(context->stabilization, gas, s, Tau_d, grad_s, U_dot, body_force, divFdiff, stab);
72
73 // Stabilized mass term
74 for (CeedInt j = 0; j < 5; j++) {
75 for (CeedInt k = 0; k < 3; k++) {
76 Grad_v[k][j][i] = wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]);
77 }
78 }
79 }
80 return 0;
81 }
82
MassFunction_Newtonian_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)83 CEED_QFUNCTION(MassFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
84 return MassFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
85 }
86
87 // @brief Computes the residual created by IDL
InternalDampingLayer_Residual(const NewtonianIGProperties gas,const State s,const CeedScalar sigma,CeedScalar damp_Y[5],CeedScalar damp_residual[5])88 CEED_QFUNCTION_HELPER void InternalDampingLayer_Residual(const NewtonianIGProperties gas, const State s, const CeedScalar sigma, CeedScalar damp_Y[5],
89 CeedScalar damp_residual[5]) {
90 ScaleN(damp_Y, sigma, 5);
91 State damp_s = StateFromY_fwd(gas, s, damp_Y);
92
93 CeedScalar U[5];
94 UnpackState_U(damp_s.U, U);
95 for (int i = 0; i < 5; i++) damp_residual[i] += U[i];
96 }
97
98 /**
99 @brief IFunction integrand for Internal Damping Layer
100
101 `location` refers to whatever scalar distance is desired for IDL to ramp from.
102 See `LinearRampCoefficient()` for details on the `amplitude`, `length`, `start`, and `location` arguments.
103
104 @param[in] s Solution `State`
105 @param[in] gas Newtonian ideal gas properties
106 @param[in] amplitude Amplitude of the IDL ramp
107 @param[in] length Length of the IDL ramp
108 @param[in] start Start of the IDL ramp
109 @param[in] location Quadrature point location (relative to IDL ramp specification)
110 @param[in] pressure Pressure used to damp to
111 @param[inout] v_i Output to be multiplied by weight function, summed into
112 @param[out] sigma IDL ramp coefficient
113 **/
InternalDampingLayer_IFunction_Integrand(const State s,const NewtonianIGProperties gas,CeedScalar amplitude,CeedScalar length,CeedScalar start,CeedScalar location,CeedScalar pressure,CeedScalar v_i[5],CeedScalar * sigma)114 CEED_QFUNCTION_HELPER void InternalDampingLayer_IFunction_Integrand(const State s, const NewtonianIGProperties gas, CeedScalar amplitude,
115 CeedScalar length, CeedScalar start, CeedScalar location, CeedScalar pressure,
116 CeedScalar v_i[5], CeedScalar *sigma) {
117 const CeedScalar sigma_ = LinearRampCoefficient(amplitude, length, start, location);
118 CeedScalar damp_state[5] = {s.Y.pressure - pressure, 0, 0, 0, 0}, idl_residual[5] = {0.};
119 InternalDampingLayer_Residual(gas, s, sigma_, damp_state, idl_residual);
120 AXPY(1, idl_residual, v_i, 5);
121 *sigma = sigma_;
122 }
123
124 /**
125 @brief IJacobian integrand for Internal Damping Layer
126
127 @note This uses a Picard-type linearization of the damping and could be replaced by an `InternalDampingLayer_fwd` that uses s and ds.
128
129 @param[in] s Solution `State`
130 @param[in] ds Change in `State` of solution
131 @param[in] gas Newtonian ideal gas properties
132 @param[in] sigma IDL ramp coefficient
133 @param[inout] v_i Output to be multiplied by weight function, summed into
134 **/
InternalDampingLayer_IJacobian_Integrand(const State s,const State ds,const NewtonianIGProperties gas,CeedScalar sigma,CeedScalar v_i[5])135 CEED_QFUNCTION_HELPER void InternalDampingLayer_IJacobian_Integrand(const State s, const State ds, const NewtonianIGProperties gas, CeedScalar sigma,
136 CeedScalar v_i[5]) {
137 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.};
138 InternalDampingLayer_Residual(gas, s, sigma, damp_state, idl_residual);
139 AXPY(1, idl_residual, v_i, 5);
140 }
141
142 // *****************************************************************************
143 // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method
144 //
145 // This is 3D compressible Navier-Stokes in conservation form with state variables of density, momentum density, and total energy density.
146 //
147 // State Variables: q = ( rho, U1, U2, U3, E )
148 // rho - Mass Density
149 // Ui - Momentum Density, Ui = rho ui
150 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z)
151 //
152 // Navier-Stokes Equations:
153 // drho/dt + div( U ) = 0
154 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu )
155 // dE/dt + div( (E + P) u ) = div( Fe )
156 //
157 // Viscous Stress:
158 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3)
159 //
160 // Thermal Stress:
161 // Fe = u Fu + k grad( T )
162 // Equation of State
163 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z)
164 //
165 // Stabilization:
166 // Tau = diag(TauC, TauM, TauM, TauM, TauE)
167 // f1 = rho sqrt(ui uj gij)
168 // gij = dXi/dX * dXi/dX
169 // TauC = Cc f1 / (8 gii)
170 // TauM = min( 1 , 1 / f1 )
171 // TauE = TauM / (Ce cv)
172 //
173 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) )
174 //
175 // Constants:
176 // lambda = - 2 / 3, From Stokes hypothesis
177 // mu , Dynamic viscosity
178 // k , Thermal conductivity
179 // cv , Specific heat, constant volume
180 // cp , Specific heat, constant pressure
181 // g , Gravity
182 // gamma = cp / cv, Specific heat ratio
183 //
184 // 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
185 // gradu )
186 // *****************************************************************************
RHSFunction_Newtonian(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)187 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
188 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
189 const bool use_divFdiff = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE;
190
191 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
192 const CeedScalar(*Grad_q) = in[1];
193 const CeedScalar(*q_data) = in[2];
194 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3];
195 const CeedScalar(*divFdiff)[CEED_Q_VLA] = use_divFdiff ? (const CeedScalar(*)[CEED_Q_VLA])in[4] : NULL;
196 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
197 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
198
199 const CeedScalar *g = context->g;
200 const CeedScalar dt = context->dt;
201 const NewtonianIGProperties gas = context->gas;
202
203 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
204 CeedScalar U[5], wdetJ, dXdx[3][3];
205 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
206 for (int j = 0; j < 5; j++) U[j] = q[j][i];
207 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
208 State s = StateFromU(gas, U);
209
210 State grad_s[3];
211 StatePhysicalGradientFromReference(Q, i, gas, s, STATEVAR_CONSERVATIVE, Grad_q, dXdx, grad_s);
212
213 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
214 KMStrainRate_State(grad_s, strain_rate);
215 NewtonianStress(gas, strain_rate, kmstress);
216 KMUnpack(kmstress, stress);
217 ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe);
218
219 StateConservative F_inviscid[3];
220 FluxInviscid(gas, s, F_inviscid);
221
222 // Total flux
223 CeedScalar Flux[5][3];
224 FluxTotal(F_inviscid, stress, Fe, Flux);
225
226 for (CeedInt j = 0; j < 5; j++) {
227 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]);
228 }
229
230 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)};
231 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j];
232
233 if (context->idl_enable) {
234 const CeedScalar idl_pressure = context->idl_pressure;
235 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]);
236 CeedScalar damp_state[5] = {s.Y.pressure - idl_pressure, 0, 0, 0, 0}, idl_residual[5] = {0.};
237 InternalDampingLayer_Residual(gas, s, sigma, damp_state, idl_residual);
238 for (int j = 0; j < 5; j++) v[j][i] -= wdetJ * idl_residual[j];
239 }
240
241 CeedScalar divFdiff_i[5] = {0.};
242 if (use_divFdiff)
243 for (int j = 1; j < 5; j++) divFdiff_i[j] = divFdiff[j - 1][i];
244
245 // -- Stabilization method: none (Galerkin), SU, or SUPG
246 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0};
247 Tau_diagPrim(context->tau_coeffs, gas, s, dXdx, dt, Tau_d);
248 Stabilization(context->stabilization, gas, s, Tau_d, grad_s, U_dot, body_force, divFdiff_i, stab);
249
250 for (CeedInt j = 0; j < 5; j++) {
251 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]);
252 }
253 }
254 return 0;
255 }
256
257 /**
258 @brief IFunction integrand of Navier-Stokes for Newtonian ideal gas
259
260 This is used in the quadrature point loop within a larger QFunction.
261 `v_i` and `dv_i` are summed into (meaning they must be some initialized value).
262 `kmstress` and `Tau_d` are given to be included as Jacobian data.
263
264 @param[in] s `State` of solution
265 @param[in] grad_s Physical gradient of solution
266 @param[in] s_dot Time derivative of solution
267 @param[in] divFdiff_i Divergence of diffusive flux
268 @param[in] x_i Coordinate location of quadrature point
269 @param[in] gas Ideal gas properties
270 @param[in] context Newtonian context
271 @param[in] dXdx Inverse of element mapping Jacobian (d\xi / dx)
272 @param[inout] v_i Output to be multiplied by weight function, summed into
273 @param[inout] grad_v_i Output to be multiplied by gradient of weight function, summed into
274 @param[out] kmstress Viscous stress, in Kelvin-Mandel ordering
275 @param[out] Tau_d Diagonal Tau coefficients
276 **/
IFunction_Newtonian_Integrand(const State s,const State grad_s[3],const State s_dot,const CeedScalar divFdiff_i[5],const CeedScalar x_i[3],const NewtonianIGProperties gas,const NewtonianIdealGasContext context,const CeedScalar dXdx[3][3],CeedScalar v_i[5],CeedScalar grad_v_i[5][3],CeedScalar kmstress[6],CeedScalar Tau_d[3])277 CEED_QFUNCTION_HELPER void IFunction_Newtonian_Integrand(const State s, const State grad_s[3], const State s_dot, const CeedScalar divFdiff_i[5],
278 const CeedScalar x_i[3], const NewtonianIGProperties gas,
279 const NewtonianIdealGasContext context, const CeedScalar dXdx[3][3], CeedScalar v_i[5],
280 CeedScalar grad_v_i[5][3], CeedScalar kmstress[6], CeedScalar Tau_d[3]) {
281 CeedScalar strain_rate[6], stress[3][3], F_visc_energy[3], F_total[5][3];
282 StateConservative F_inviscid[3];
283 const CeedScalar *g = context->g, dt = context->dt;
284
285 // Advective and viscous fluxes
286 KMStrainRate_State(grad_s, strain_rate);
287 NewtonianStress(gas, strain_rate, kmstress);
288 KMUnpack(kmstress, stress);
289 ViscousEnergyFlux(gas, s.Y, grad_s, stress, F_visc_energy);
290 FluxInviscid(gas, s, F_inviscid);
291 FluxTotal(F_inviscid, stress, F_visc_energy, F_total);
292 AXPY(-1, (CeedScalar *)F_total, (CeedScalar *)grad_v_i, 15);
293
294 // Body force and time derivative
295 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)};
296 CeedScalar U_dot[5];
297 UnpackState_U(s_dot.U, U_dot);
298 for (CeedInt j = 0; j < 5; j++) v_i[j] += U_dot[j] - body_force[j];
299
300 // Stabilization
301 CeedScalar stab[5][3];
302 Tau_diagPrim(context->tau_coeffs, gas, s, dXdx, dt, Tau_d);
303 Stabilization(context->stabilization, gas, s, Tau_d, grad_s, U_dot, body_force, divFdiff_i, stab);
304 AXPY(1, (CeedScalar *)stab, (CeedScalar *)grad_v_i, 15);
305 }
306
307 // @brief State-independent IFunction of Navier-Stokes for Newtonian ideal gas
IFunction_Newtonian(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)308 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
309 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
310 const bool use_divFdiff = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE;
311
312 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
313 const CeedScalar(*grad_q) = in[1];
314 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
315 const CeedScalar(*q_data) = in[3];
316 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4];
317 const CeedScalar(*divFdiff)[CEED_Q_VLA] = use_divFdiff ? (const CeedScalar(*)[CEED_Q_VLA])in[5] : NULL;
318 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
319 CeedScalar(*grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
320 CeedScalar(*jac_data) = out[2];
321
322 const NewtonianIGProperties gas = context->gas;
323
324 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
325 const CeedScalar q_i[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
326 const CeedScalar q_i_dot[5] = {q_dot[0][i], q_dot[1][i], q_dot[2][i], q_dot[3][i], q_dot[4][i]};
327 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
328 const State s = StateFromQ(gas, q_i, state_var);
329 const State s_dot = StateFromQ_fwd(gas, s, q_i_dot, state_var);
330
331 CeedScalar wdetJ, dXdx[3][3];
332 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
333 State grad_s[3];
334 StatePhysicalGradientFromReference(Q, i, gas, s, state_var, grad_q, dXdx, grad_s);
335 CeedScalar divFdiff_i[5] = {0.};
336 if (use_divFdiff)
337 for (int j = 1; j < 5; j++) divFdiff_i[j] = divFdiff[j - 1][i];
338
339 CeedScalar v_i[5] = {0.}, grad_v_i[5][3] = {{0.}}, kmstress[6], Tau_d[3], sigma = 0;
340 IFunction_Newtonian_Integrand(s, grad_s, s_dot, divFdiff_i, x_i, gas, context, dXdx, v_i, grad_v_i, kmstress, Tau_d);
341 if (context->idl_enable)
342 InternalDampingLayer_IFunction_Integrand(s, gas, context->idl_amplitude, context->idl_length, context->idl_start, x_i[0], context->idl_pressure,
343 v_i, &sigma);
344
345 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * v_i[j];
346 for (CeedInt j = 0; j < 5; j++) {
347 for (CeedInt k = 0; k < 3; k++)
348 grad_v[k][j][i] = wdetJ * (grad_v_i[j][0] * dXdx[k][0] + grad_v_i[j][1] * dXdx[k][1] + grad_v_i[j][2] * dXdx[k][2]);
349 }
350
351 StoredValuesPack(Q, i, 0, 5, q_i, jac_data);
352 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data);
353 StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data);
354 if (context->idl_enable) StoredValuesPack(Q, i, 14, 1, &sigma, jac_data);
355 }
356 return 0;
357 }
358
IFunction_Newtonian_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)359 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
360 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
361 }
362
IFunction_Newtonian_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)363 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
364 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
365 }
366
IFunction_Newtonian_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)367 CEED_QFUNCTION(IFunction_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
368 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY);
369 }
370
371 /**
372 @brief IJacobian integrand of Navier-Stokes for Newtonian ideal gas
373
374 This is used in the quadrature point loop within a larger QFunction.
375 `v_i` and `dv_i` are summed into (meaning they must be some initialized value).
376 `kmstress` and `Tau_d` are (generally) calculated and stored by the IFunction.
377
378 @param[in] s `State` of solution
379 @param[in] ds Change in `State` of solution
380 @param[in] grad_ds Physical gradient of change in `State` of solution
381 @param[in] gas Ideal gas properties
382 @param[in] context Newtonian context
383 @param[in] kmstress Viscous stress, in Kelvin-Mandel ordering
384 @param[in] Tau_d Diagonal Tau coefficients
385 @param[inout] v_i Output to be multiplied by weight function, summed into
386 @param[inout] grad_v_i Output to be multiplied by gradient of weight function, summed into
387 **/
IJacobian_Newtonian_Integrand(const State s,const State ds,const State grad_ds[3],const NewtonianIGProperties gas,const NewtonianIdealGasContext context,const CeedScalar kmstress[6],const CeedScalar Tau_d[3],CeedScalar v_i[5],CeedScalar grad_v_i[5][3])388 CEED_QFUNCTION_HELPER void IJacobian_Newtonian_Integrand(const State s, const State ds, const State grad_ds[3], const NewtonianIGProperties gas,
389 const NewtonianIdealGasContext context, const CeedScalar kmstress[6],
390 const CeedScalar Tau_d[3], CeedScalar v_i[5], CeedScalar grad_v_i[5][3]) {
391 const CeedScalar *g = context->g;
392 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dF_visc_energy[3], dF_total[5][3];
393 StateConservative dF_inviscid[3];
394
395 // Advective and viscous fluxes
396 KMStrainRate_State(grad_ds, dstrain_rate);
397 NewtonianStress(gas, dstrain_rate, dkmstress);
398 KMUnpack(dkmstress, dstress);
399 KMUnpack(kmstress, stress);
400 ViscousEnergyFlux_fwd(gas, s.Y, ds.Y, grad_ds, stress, dstress, dF_visc_energy);
401 FluxInviscid_fwd(gas, s, ds, dF_inviscid);
402 FluxTotal(dF_inviscid, dstress, dF_visc_energy, dF_total);
403 AXPY(-1, (CeedScalar *)dF_total, (CeedScalar *)grad_v_i, 15);
404
405 // Body force and time derivative
406 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)};
407 CeedScalar dU[5], dU_dot[5];
408 UnpackState_U(ds.U, dU);
409 for (CeedInt j = 0; j < 5; j++) {
410 dU_dot[j] = context->ijacobian_time_shift * dU[j];
411 v_i[j] = dU_dot[j] - dbody_force[j];
412 }
413
414 // Stabilization
415 CeedScalar dstab[5][3];
416 const CeedScalar zeroFlux[5] = {0.};
417 Stabilization(context->stabilization, gas, s, Tau_d, grad_ds, dU_dot, dbody_force, zeroFlux, dstab);
418 AXPY(1, (CeedScalar *)dstab, (CeedScalar *)grad_v_i, 15);
419 }
420
421 // @brief State-independent IJacobian of Navier-Stokes for Newtonian ideal gas
IJacobian_Newtonian(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)422 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
423 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
424 const CeedScalar(*grad_dq) = in[1];
425 const CeedScalar(*q_data) = in[2];
426 const CeedScalar(*jac_data) = in[3];
427 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
428 CeedScalar(*grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
429
430 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
431 const NewtonianIGProperties gas = context->gas;
432
433 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
434 const CeedScalar dq_i[5] = {dq[0][i], dq[1][i], dq[2][i], dq[3][i], dq[4][i]};
435 CeedScalar qi[5], kmstress[6], Tau_d[3];
436 StoredValuesUnpack(Q, i, 0, 5, jac_data, qi);
437 StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress);
438 StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d);
439 const State s = StateFromQ(gas, qi, state_var);
440 const State ds = StateFromQ_fwd(gas, s, dq_i, state_var);
441
442 CeedScalar wdetJ, dXdx[3][3];
443 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
444 State grad_ds[3];
445 StatePhysicalGradientFromReference(Q, i, gas, s, state_var, grad_dq, dXdx, grad_ds);
446
447 CeedScalar v_i[5] = {0.}, grad_v_i[5][3] = {{0.}};
448 IJacobian_Newtonian_Integrand(s, ds, grad_ds, gas, context, kmstress, Tau_d, v_i, grad_v_i);
449 if (context->idl_enable) {
450 CeedScalar sigma;
451 StoredValuesUnpack(Q, i, 14, 1, jac_data, &sigma);
452 InternalDampingLayer_IJacobian_Integrand(s, ds, gas, sigma, v_i);
453 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * v_i[j];
454 }
455
456 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * v_i[j];
457 for (int j = 0; j < 5; j++) {
458 for (int k = 0; k < 3; k++) grad_v[k][j][i] = wdetJ * (grad_v_i[j][0] * dXdx[k][0] + grad_v_i[j][1] * dXdx[k][1] + grad_v_i[j][2] * dXdx[k][2]);
459 }
460 }
461 return 0;
462 }
463
IJacobian_Newtonian_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)464 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
465 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
466 }
467
IJacobian_Newtonian_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)468 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
469 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
470 }
471
IJacobian_Newtonian_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)472 CEED_QFUNCTION(IJacobian_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
473 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY);
474 }
475
476 // *****************************************************************************
477 // Compute boundary integral (ie. for strongly set inflows)
478 // *****************************************************************************
BoundaryIntegral(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)479 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
480 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
481 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
482 const CeedScalar(*Grad_q) = in[1];
483 const CeedScalar(*q_data_sur) = in[2];
484 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
485 CeedScalar(*jac_data_sur) = context->is_implicit ? out[1] : NULL;
486
487 const bool is_implicit = context->is_implicit;
488 const NewtonianIGProperties gas = context->gas;
489
490 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
491 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
492 State s = StateFromQ(gas, qi, state_var);
493
494 CeedScalar wdetJb, dXdx[2][3], normal[3];
495 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal);
496 wdetJb *= is_implicit ? -1. : 1.;
497
498 State grad_s[3];
499 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s);
500
501 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
502 KMStrainRate_State(grad_s, strain_rate);
503 NewtonianStress(gas, strain_rate, kmstress);
504 KMUnpack(kmstress, stress);
505 ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe);
506
507 StateConservative F_inviscid[3];
508 FluxInviscid(gas, s, F_inviscid);
509
510 CeedScalar Flux[5];
511 FluxTotal_Boundary(F_inviscid, stress, Fe, normal, Flux);
512
513 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];
514
515 if (is_implicit) {
516 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
517 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
518 }
519 }
520 return 0;
521 }
522
BoundaryIntegral_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)523 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
524 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
525 }
526
BoundaryIntegral_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)527 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
528 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE);
529 }
530
BoundaryIntegral_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)531 CEED_QFUNCTION(BoundaryIntegral_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
532 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_ENTROPY);
533 }
534
535 // *****************************************************************************
536 // Jacobian for "set nothing" boundary integral
537 // *****************************************************************************
BoundaryIntegral_Jacobian(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)538 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
539 StateVariable state_var) {
540 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
541 const CeedScalar(*Grad_dq) = in[1];
542 const CeedScalar(*q_data_sur) = in[2];
543 const CeedScalar(*jac_data_sur) = in[4];
544 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
545
546 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
547 const NewtonianIGProperties gas = context->gas;
548 const bool is_implicit = context->is_implicit;
549
550 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
551 CeedScalar wdetJb, dXdx[2][3], normal[3];
552 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal);
553 wdetJb *= is_implicit ? -1. : 1.;
554
555 CeedScalar qi[5], kmstress[6], dqi[5];
556 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi);
557 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress);
558 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];
559
560 State s = StateFromQ(gas, qi, state_var);
561 State ds = StateFromQ_fwd(gas, s, dqi, state_var);
562
563 State grad_ds[3];
564 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s, state_var, Grad_dq, dXdx, grad_ds);
565
566 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3];
567 KMStrainRate_State(grad_ds, dstrain_rate);
568 NewtonianStress(gas, dstrain_rate, dkmstress);
569 KMUnpack(dkmstress, dstress);
570 KMUnpack(kmstress, stress);
571 ViscousEnergyFlux_fwd(gas, s.Y, ds.Y, grad_ds, stress, dstress, dFe);
572
573 StateConservative dF_inviscid[3];
574 FluxInviscid_fwd(gas, s, ds, dF_inviscid);
575
576 CeedScalar dFlux[5];
577 FluxTotal_Boundary(dF_inviscid, dstress, dFe, normal, dFlux);
578
579 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j];
580 }
581 return 0;
582 }
583
BoundaryIntegral_Jacobian_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)584 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
585 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
586 }
587
BoundaryIntegral_Jacobian_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)588 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
589 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
590 }
591
BoundaryIntegral_Jacobian_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)592 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
593 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY);
594 }
595
596 // @brief Volume integral for RHS of divergence of diffusive flux direct projection
DivDiffusiveFluxVolumeRHS_NS(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)597 CEED_QFUNCTION_HELPER int DivDiffusiveFluxVolumeRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
598 StateVariable state_var) {
599 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
600 const CeedScalar(*Grad_q) = in[1];
601 const CeedScalar(*q_data) = in[2];
602 CeedScalar(*Grad_v)[4][CEED_Q_VLA] = (CeedScalar(*)[4][CEED_Q_VLA])out[0];
603
604 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
605 const NewtonianIGProperties gas = context->gas;
606 const StateConservative ZeroInviscidFluxes[3] = {{0}};
607
608 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
609 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
610 const State s = StateFromQ(gas, qi, state_var);
611 CeedScalar wdetJ, dXdx[3][3];
612 CeedScalar stress[3][3], Fe[3], Fdiff[5][3];
613
614 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
615 { // Get stress and Fe
616 State grad_s[3];
617 CeedScalar strain_rate[6], kmstress[6];
618
619 StatePhysicalGradientFromReference(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s);
620 KMStrainRate_State(grad_s, strain_rate);
621 NewtonianStress(gas, strain_rate, kmstress);
622 KMUnpack(kmstress, stress);
623 ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe);
624 }
625
626 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff);
627
628 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip
629 for (CeedInt k = 0; k < 3; k++) {
630 Grad_v[k][j - 1][i] = -wdetJ * Dot3(dXdx[k], Fdiff[j]);
631 }
632 }
633 }
634 return 0;
635 }
636
DivDiffusiveFluxVolumeRHS_NS_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)637 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
638 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
639 }
640
DivDiffusiveFluxVolumeRHS_NS_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)641 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
642 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE);
643 }
644
DivDiffusiveFluxVolumeRHS_NS_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)645 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
646 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY);
647 }
648
649 // @brief Boundary integral for RHS of divergence of diffusive flux direct projection
DivDiffusiveFluxBoundaryRHS_NS(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)650 CEED_QFUNCTION_HELPER int DivDiffusiveFluxBoundaryRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
651 StateVariable state_var) {
652 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
653 const CeedScalar(*Grad_q) = in[1];
654 const CeedScalar(*q_data) = in[2];
655 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
656
657 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
658 const NewtonianIGProperties gas = context->gas;
659 const StateConservative ZeroInviscidFluxes[3] = {{0}};
660
661 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
662 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
663 const State s = StateFromQ(gas, qi, state_var);
664 CeedScalar wdetJ, dXdx[3][3], normal[3];
665 CeedScalar stress[3][3], Fe[3], Fdiff[5];
666
667 QdataBoundaryGradientUnpack_3D(Q, i, q_data, &wdetJ, dXdx, normal);
668 { // Get stress and Fe
669 State grad_s[3];
670 CeedScalar strain_rate[6], kmstress[6];
671
672 StatePhysicalGradientFromReference(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s);
673 KMStrainRate_State(grad_s, strain_rate);
674 NewtonianStress(gas, strain_rate, kmstress);
675 KMUnpack(kmstress, stress);
676 ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe);
677 }
678
679 FluxTotal_Boundary(ZeroInviscidFluxes, stress, Fe, normal, Fdiff);
680
681 // Continuity has no diffusive flux, therefore skip
682 for (CeedInt j = 1; j < 5; j++) v[j - 1][i] = wdetJ * Fdiff[j];
683 }
684 return 0;
685 }
686
DivDiffusiveFluxBoundaryRHS_NS_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)687 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
688 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
689 }
690
DivDiffusiveFluxBoundaryRHS_NS_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)691 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
692 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE);
693 }
694
DivDiffusiveFluxBoundaryRHS_NS_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)695 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
696 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY);
697 }
698
699 // @brief Integral for RHS of diffusive flux indirect projection
DiffusiveFluxRHS_NS(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)700 CEED_QFUNCTION_HELPER int DiffusiveFluxRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
701 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
702 const CeedScalar(*Grad_q) = in[1];
703 const CeedScalar(*q_data) = in[2];
704 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
705
706 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
707 const NewtonianIGProperties gas = context->gas;
708 const StateConservative ZeroInviscidFluxes[3] = {{0}};
709
710 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
711 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
712 const State s = StateFromQ(gas, qi, state_var);
713 CeedScalar wdetJ, dXdx[3][3];
714 CeedScalar stress[3][3], Fe[3], Fdiff[5][3];
715
716 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
717 { // Get stress and Fe
718 State grad_s[3];
719 CeedScalar strain_rate[6], kmstress[6];
720
721 StatePhysicalGradientFromReference(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s);
722 KMStrainRate_State(grad_s, strain_rate);
723 NewtonianStress(gas, strain_rate, kmstress);
724 KMUnpack(kmstress, stress);
725 ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe);
726 }
727
728 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff);
729
730 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip
731 for (CeedInt k = 0; k < 3; k++) {
732 v[(j - 1) * 3 + k][i] = wdetJ * Fdiff[j][k];
733 }
734 }
735 }
736 return 0;
737 }
738
DiffusiveFluxRHS_NS_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)739 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
740 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
741 }
742
DiffusiveFluxRHS_NS_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)743 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
744 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE);
745 }
746
DiffusiveFluxRHS_NS_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)747 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
748 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY);
749 }
750