1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Shock tube initial condition and Euler equation operator for HONEE - modified from eulervortex.h
6
7 // Model from:
8 // On the Order of Accuracy and Numerical Performance of Two Classes of Finite Volume WENO Schemes, Zhang, Zhang, and Shu (2011).
9 #include <ceed/types.h>
10 #ifndef CEED_RUNNING_JIT_PASS
11 #include <stdbool.h>
12 #endif
13
14 #include "utils.h"
15
16 typedef struct SetupContextShock_ *SetupContextShock;
17 struct SetupContextShock_ {
18 CeedScalar theta0;
19 CeedScalar thetaC;
20 CeedScalar P0;
21 CeedScalar N;
22 CeedScalar cv;
23 CeedScalar cp;
24 CeedScalar time;
25 CeedScalar mid_point;
26 CeedScalar P_high;
27 CeedScalar rho_high;
28 CeedScalar P_low;
29 CeedScalar rho_low;
30 };
31
32 typedef struct ShockTubeContext_ *ShockTubeContext;
33 struct ShockTubeContext_ {
34 CeedScalar Cyzb;
35 CeedScalar Byzb;
36 CeedScalar c_tau;
37 bool implicit;
38 bool yzb;
39 int stabilization;
40 };
41
42 // *****************************************************************************
43 // This function sets the initial conditions
44 //
45 // Temperature:
46 // T = P / (rho * R)
47 // Density:
48 // rho = 1.0 if x <= mid_point
49 // = 0.125 if x > mid_point
50 // Pressure:
51 // P = 1.0 if x <= mid_point
52 // = 0.1 if x > mid_point
53 // Velocity:
54 // u = 0
55 // Velocity/Momentum Density:
56 // Ui = rho ui
57 // Total Energy:
58 // E = P / (gamma - 1) + rho (u u)/2
59 //
60 // Constants:
61 // cv , Specific heat, constant volume
62 // cp , Specific heat, constant pressure
63 // mid_point , Location of initial domain mid_point
64 // gamma = cp / cv, Specific heat ratio
65 //
66 // *****************************************************************************
67
68 // *****************************************************************************
69 // This helper function provides support for the exact, time-dependent solution (currently not implemented) and IC formulation for Euler traveling
70 // vortex
71 // *****************************************************************************
Exact_ShockTube(CeedInt dim,CeedScalar time,const CeedScalar X[],CeedInt Nf,CeedScalar q[],void * ctx)72 CEED_QFUNCTION_HELPER CeedInt Exact_ShockTube(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) {
73 // Context
74 const SetupContextShock context = (SetupContextShock)ctx;
75 const CeedScalar mid_point = context->mid_point; // Midpoint of the domain
76 const CeedScalar P_high = context->P_high; // Driver section pressure
77 const CeedScalar rho_high = context->rho_high; // Driver section density
78 const CeedScalar P_low = context->P_low; // Driven section pressure
79 const CeedScalar rho_low = context->rho_low; // Driven section density
80
81 // Setup
82 const CeedScalar gamma = 1.4; // ratio of specific heats
83 const CeedScalar x = X[0]; // Coordinates
84
85 CeedScalar rho, P, u[3] = {0.};
86
87 // Initial Conditions
88 if (x <= mid_point + 200 * CEED_EPSILON) {
89 rho = rho_high;
90 P = P_high;
91 } else {
92 rho = rho_low;
93 P = P_low;
94 }
95
96 // Assign exact solution
97 q[0] = rho;
98 q[1] = rho * u[0];
99 q[2] = rho * u[1];
100 q[3] = rho * u[2];
101 q[4] = P / (gamma - 1.0) + rho * (u[0] * u[0]) / 2.;
102
103 return 0;
104 }
105
106 // *****************************************************************************
107 // Helper function for computing flux Jacobian
108 // *****************************************************************************
ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5],const CeedScalar rho,const CeedScalar u[3],const CeedScalar E,const CeedScalar gamma)109 CEED_QFUNCTION_HELPER void ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5], const CeedScalar rho, const CeedScalar u[3], const CeedScalar E,
110 const CeedScalar gamma) {
111 CeedScalar u_sq = u[0] * u[0] + u[1] * u[1] + u[2] * u[2]; // Velocity square
112 for (CeedInt i = 0; i < 3; i++) { // Jacobian matrices for 3 directions
113 for (CeedInt j = 0; j < 3; j++) { // Rows of each Jacobian matrix
114 dF[i][j + 1][0] = ((i == j) ? ((gamma - 1.) * (u_sq / 2.)) : 0.) - u[i] * u[j];
115 for (CeedInt k = 0; k < 3; k++) { // Columns of each Jacobian matrix
116 dF[i][0][k + 1] = ((i == k) ? 1. : 0.);
117 dF[i][j + 1][k + 1] = ((j == k) ? u[i] : 0.) + ((i == k) ? u[j] : 0.) - ((i == j) ? u[k] : 0.) * (gamma - 1.);
118 dF[i][4][k + 1] = ((i == k) ? (E * gamma / rho - (gamma - 1.) * u_sq / 2.) : 0.) - (gamma - 1.) * u[i] * u[k];
119 }
120 dF[i][j + 1][4] = ((i == j) ? (gamma - 1.) : 0.);
121 }
122 dF[i][4][0] = u[i] * ((gamma - 1.) * u_sq - E * gamma / rho);
123 dF[i][4][4] = u[i] * gamma;
124 }
125 }
126
127 // *****************************************************************************
128 // Helper function for calculating the covariant length scale in the direction of some 3 element input vector
129 //
130 // Where
131 // vec = vector that length is measured in the direction of
132 // h = covariant element length along vec
133 // *****************************************************************************
Covariant_length_along_vector(CeedScalar vec[3],const CeedScalar dXdx[3][3])134 CEED_QFUNCTION_HELPER CeedScalar Covariant_length_along_vector(CeedScalar vec[3], const CeedScalar dXdx[3][3]) {
135 CeedScalar vec_dot_jacobian[3] = {0.0};
136
137 MatVec3(dXdx, vec, CEED_TRANSPOSE, vec_dot_jacobian);
138 return 2.0 * Norm3(vec) / Norm3(vec_dot_jacobian);
139 }
140
141 // *****************************************************************************
142 // Helper function for computing Tau elements (stabilization constant)
143 // Model from:
144 // Stabilized Methods for Compressible Flows, Hughes et al 2010
145 //
146 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix
147 // Tau[i] = c_tau h[i] Xi(Pe) / rho(A[i]) (no sum)
148 //
149 // Where
150 // c_tau = stabilization constant (0.5 is reported as "optimal")
151 // h[i] = 2 length(dxdX[i])
152 // Pe = Peclet number ( Pe = sqrt(u u) / dot(dXdx,u) diffusivity )
153 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number )
154 // rho(A[i]) = spectral radius of the convective flux Jacobian i, wave speed in direction i
155 // *****************************************************************************
Tau_spatial(CeedScalar Tau_x[3],const CeedScalar dXdx[3][3],const CeedScalar u[3],const CeedScalar sound_speed,const CeedScalar c_tau)156 CEED_QFUNCTION_HELPER void Tau_spatial(CeedScalar Tau_x[3], const CeedScalar dXdx[3][3], const CeedScalar u[3], const CeedScalar sound_speed,
157 const CeedScalar c_tau) {
158 for (CeedInt i = 0; i < 3; i++) {
159 // length of element in direction i
160 CeedScalar h = 2 / sqrt(Square(dXdx[0][i]) + Square(dXdx[1][i]) + Square(dXdx[2][i]));
161 // fastest wave in direction i
162 CeedScalar fastest_wave = fabs(u[i]) + sound_speed;
163 Tau_x[i] = c_tau * h / fastest_wave;
164 }
165 }
166
167 // *****************************************************************************
168 // This QFunction sets the initial conditions for shock tube
169 // *****************************************************************************
ICsShockTube(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)170 CEED_QFUNCTION(ICsShockTube)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
171 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
172 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
173
174 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
175 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
176 CeedScalar q[5];
177
178 Exact_ShockTube(3, 0., x, 5, q, ctx);
179
180 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
181 }
182 return 0;
183 }
184
185 // *****************************************************************************
186 // This QFunction implements the following formulation of Euler equations with explicit time stepping method
187 //
188 // This is 3D Euler for compressible gas dynamics in conservation form with state variables of density, momentum density, and total energy density.
189 //
190 // State Variables: q = ( rho, U1, U2, U3, E )
191 // rho - Mass Density
192 // Ui - Momentum Density, Ui = rho ui
193 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2
194 //
195 // Euler Equations:
196 // drho/dt + div( U ) = 0
197 // dU/dt + div( rho (u x u) + P I3 ) = 0
198 // dE/dt + div( (E + P) u ) = 0
199 //
200 // Equation of State:
201 // P = (gamma - 1) (E - rho (u u) / 2)
202 //
203 // Constants:
204 // cv , Specific heat, constant volume
205 // cp , Specific heat, constant pressure
206 // g , Gravity
207 // gamma = cp / cv, Specific heat ratio
208 // *****************************************************************************
EulerShockTube(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)209 CEED_QFUNCTION(EulerShockTube)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
210 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
211 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1];
212 const CeedScalar(*q_data) = in[2];
213 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
214 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
215
216 const CeedScalar gamma = 1.4;
217
218 ShockTubeContext context = (ShockTubeContext)ctx;
219 const CeedScalar Cyzb = context->Cyzb;
220 const CeedScalar Byzb = context->Byzb;
221 const CeedScalar c_tau = context->c_tau;
222
223 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
224 // Setup
225 // -- Interp in
226 const CeedScalar rho = q[0][i];
227 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho};
228 const CeedScalar E = q[4][i];
229 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]};
230 const CeedScalar dU[3][3] = {
231 {dq[0][1][i], dq[1][1][i], dq[2][1][i]},
232 {dq[0][2][i], dq[1][2][i], dq[2][2][i]},
233 {dq[0][3][i], dq[1][3][i], dq[2][3][i]}
234 };
235 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]};
236 CeedScalar wdetJ, dXdx[3][3];
237 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
238 // dU/dx
239 CeedScalar du[3][3] = {{0}};
240 CeedScalar drhodx[3] = {0};
241 CeedScalar dEdx[3] = {0};
242 CeedScalar dUdx[3][3] = {{0}};
243 CeedScalar dXdxdXdxT[3][3] = {{0}};
244 for (CeedInt j = 0; j < 3; j++) {
245 for (CeedInt k = 0; k < 3; k++) {
246 du[j][k] = (dU[j][k] - drho[k] * u[j]) / rho;
247 drhodx[j] += drho[k] * dXdx[k][j];
248 dEdx[j] += dE[k] * dXdx[k][j];
249 for (CeedInt l = 0; l < 3; l++) {
250 dUdx[j][k] += dU[j][l] * dXdx[l][k];
251 dXdxdXdxT[j][k] += dXdx[j][l] * dXdx[k][l]; // dXdx_j,k * dXdx_k,j
252 }
253 }
254 }
255
256 const CeedScalar E_kinetic = 0.5 * rho * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]), E_internal = E - E_kinetic,
257 P = E_internal * (gamma - 1); // P = pressure
258
259 // The Physics
260 // Zero v and dv so all future terms can safely sum into it
261 for (CeedInt j = 0; j < 5; j++) {
262 v[j][i] = 0;
263 for (CeedInt k = 0; k < 3; k++) dv[k][j][i] = 0;
264 }
265
266 // -- Density
267 // ---- u rho
268 for (CeedInt j = 0; j < 3; j++) dv[j][0][i] += wdetJ * (rho * u[0] * dXdx[j][0] + rho * u[1] * dXdx[j][1] + rho * u[2] * dXdx[j][2]);
269 // -- Momentum
270 // ---- rho (u x u) + P I3
271 for (CeedInt j = 0; j < 3; j++) {
272 for (CeedInt k = 0; k < 3; k++) {
273 dv[k][j + 1][i] += wdetJ * ((rho * u[j] * u[0] + (j == 0 ? P : 0)) * dXdx[k][0] + (rho * u[j] * u[1] + (j == 1 ? P : 0)) * dXdx[k][1] +
274 (rho * u[j] * u[2] + (j == 2 ? P : 0)) * dXdx[k][2]);
275 }
276 }
277 // -- Total Energy Density
278 // ---- (E + P) u
279 for (CeedInt j = 0; j < 3; j++) dv[j][4][i] += wdetJ * (E + P) * (u[0] * dXdx[j][0] + u[1] * dXdx[j][1] + u[2] * dXdx[j][2]);
280
281 // -- YZB stabilization
282 if (context->yzb) {
283 CeedScalar drho_norm = 0.0; // magnitude of the density gradient
284 CeedScalar j_vec[3] = {0.0}; // unit vector aligned with the density gradient
285 CeedScalar h_shock = 0.0; // element lengthscale
286 CeedScalar acoustic_vel = 0.0; // characteristic velocity, acoustic speed
287 CeedScalar tau_shock = 0.0; // timescale
288 CeedScalar nu_shock = 0.0; // artificial diffusion
289
290 // Unit vector aligned with the density gradient
291 drho_norm = Norm3(drhodx);
292 for (CeedInt j = 0; j < 3; j++) j_vec[j] = drhodx[j] / (drho_norm + 1e-20);
293
294 if (drho_norm == 0.0) {
295 nu_shock = 0.0;
296 } else {
297 h_shock = Covariant_length_along_vector(j_vec, dXdx);
298 h_shock /= Cyzb;
299 acoustic_vel = sqrt(gamma * P / rho);
300 tau_shock = h_shock / (2 * acoustic_vel) * pow(drho_norm * h_shock / rho, Byzb);
301 nu_shock = fabs(tau_shock * acoustic_vel * acoustic_vel);
302 }
303
304 for (CeedInt j = 0; j < 3; j++) dv[j][0][i] -= wdetJ * nu_shock * drhodx[j];
305
306 for (CeedInt k = 0; k < 3; k++) {
307 for (CeedInt j = 0; j < 3; j++) dv[j][k][i] -= wdetJ * nu_shock * du[k][j];
308 }
309
310 for (CeedInt j = 0; j < 3; j++) dv[j][4][i] -= wdetJ * nu_shock * dEdx[j];
311 }
312
313 // Stabilization
314 // Need the Jacobian for the advective fluxes for stabilization
315 // indexed as: jacob_F_conv[direction][flux component][solution component]
316 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}};
317 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma);
318
319 // dqdx collects drhodx, dUdx and dEdx in one vector
320 CeedScalar dqdx[5][3];
321 for (CeedInt j = 0; j < 3; j++) {
322 dqdx[0][j] = drhodx[j];
323 dqdx[4][j] = dEdx[j];
324 for (CeedInt k = 0; k < 3; k++) dqdx[k + 1][j] = dUdx[k][j];
325 }
326
327 // strong_conv = dF/dq * dq/dx (Strong convection)
328 CeedScalar strong_conv[5] = {0};
329 for (CeedInt j = 0; j < 3; j++) {
330 for (CeedInt k = 0; k < 5; k++) {
331 for (CeedInt l = 0; l < 5; l++) strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j];
332 }
333 }
334
335 // Stabilization
336 // -- Tau elements
337 const CeedScalar sound_speed = sqrt(gamma * P / rho);
338 CeedScalar Tau_x[3] = {0.};
339 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau);
340
341 CeedScalar stab[5][3] = {0};
342 switch (context->stabilization) {
343 case 0: // Galerkin
344 break;
345 case 1: // SU
346 for (CeedInt j = 0; j < 3; j++) {
347 for (CeedInt k = 0; k < 5; k++) {
348 for (CeedInt l = 0; l < 5; l++) {
349 stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l];
350 }
351 }
352 }
353 for (CeedInt j = 0; j < 5; j++) {
354 for (CeedInt k = 0; k < 3; k++) dv[k][j][i] -= wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]);
355 }
356 break;
357 }
358 }
359 return 0;
360 }
361