1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Euler traveling vortex initial condition and operator for HONEE
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 EulerContext_ *EulerContext;
17 struct EulerContext_ {
18 CeedScalar center[3];
19 CeedScalar curr_time;
20 CeedScalar vortex_strength;
21 CeedScalar c_tau;
22 CeedScalar mean_velocity[3];
23 bool implicit;
24 int euler_test;
25 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG
26 };
27
28 // *****************************************************************************
29 // This function sets the initial conditions
30 //
31 // Temperature:
32 // T = 1 - (gamma - 1) vortex_strength**2 exp(1 - r**2) / (8 gamma pi**2)
33 // Density:
34 // rho = (T/S_vortex)^(1 / (gamma - 1))
35 // Pressure:
36 // P = rho * T
37 // Velocity:
38 // ui = 1 + vortex_strength exp((1 - r**2)/2.) [yc - y, x - xc] / (2 pi)
39 // r = sqrt( (x - xc)**2 + (y - yc)**2 )
40 // Velocity/Momentum Density:
41 // Ui = rho ui
42 // Total Energy:
43 // E = P / (gamma - 1) + rho (u u)/2
44 //
45 // Constants:
46 // cv , Specific heat, constant volume
47 // cp , Specific heat, constant pressure
48 // vortex_strength , Strength of vortex
49 // center , Location of bubble center
50 // gamma = cp / cv, Specific heat ratio
51 //
52 // *****************************************************************************
53
54 // *****************************************************************************
55 // This helper function provides support for the exact, time-dependent solution (currently not implemented) and IC formulation for Euler traveling
56 // vortex
57 // *****************************************************************************
Exact_Euler(CeedInt dim,CeedScalar time,const CeedScalar X[],CeedInt Nf,CeedScalar q[],void * ctx)58 CEED_QFUNCTION_HELPER int Exact_Euler(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) {
59 // Context
60 const EulerContext context = (EulerContext)ctx;
61 const CeedScalar vortex_strength = context->vortex_strength;
62 const CeedScalar *center = context->center; // Center of the domain
63 const CeedScalar *mean_velocity = context->mean_velocity;
64
65 // Setup
66 const CeedScalar gamma = 1.4;
67 const CeedScalar cv = 2.5;
68 const CeedScalar R = 1.;
69 const CeedScalar x = X[0], y = X[1]; // Coordinates
70 // Vortex center
71 const CeedScalar xc = center[0] + mean_velocity[0] * time;
72 const CeedScalar yc = center[1] + mean_velocity[1] * time;
73
74 const CeedScalar x0 = x - xc;
75 const CeedScalar y0 = y - yc;
76 const CeedScalar r = sqrt(Square(x0) + Square(y0));
77 const CeedScalar C = vortex_strength * exp((1. - r * r) / 2.) / (2. * M_PI);
78 const CeedScalar delta_T = -(gamma - 1.) * vortex_strength * vortex_strength * exp(1 - r * r) / (8. * gamma * M_PI * M_PI);
79 const CeedScalar S_vortex = 1; // no perturbation in the entropy P / rho^gamma
80 const CeedScalar S_bubble = (gamma - 1.) * vortex_strength * vortex_strength / (8. * gamma * M_PI * M_PI);
81 CeedScalar rho, P, T, E, u[3] = {0.};
82
83 // Initial Conditions
84 switch (context->euler_test) {
85 case 0: // Traveling vortex
86 T = 1 + delta_T;
87 // P = rho * T
88 // P = S * rho^gamma
89 // Solve for rho, then substitute for P
90 rho = pow(T / S_vortex, 1 / (gamma - 1.));
91 P = rho * T;
92 u[0] = mean_velocity[0] - C * y0;
93 u[1] = mean_velocity[1] + C * x0;
94
95 // Assign exact solution
96 q[0] = rho;
97 q[1] = rho * u[0];
98 q[2] = rho * u[1];
99 q[3] = rho * u[2];
100 q[4] = P / (gamma - 1.) + rho * (u[0] * u[0] + u[1] * u[1]) / 2.;
101 break;
102 case 1: // Constant zero velocity, density constant, total energy constant
103 rho = 1.;
104 E = 2.;
105
106 // Assign exact solution
107 q[0] = rho;
108 q[1] = rho * u[0];
109 q[2] = rho * u[1];
110 q[3] = rho * u[2];
111 q[4] = E;
112 break;
113 case 2: // Constant nonzero velocity, density constant, total energy constant
114 rho = 1.;
115 E = 2.;
116 u[0] = mean_velocity[0];
117 u[1] = mean_velocity[1];
118
119 // Assign exact solution
120 q[0] = rho;
121 q[1] = rho * u[0];
122 q[2] = rho * u[1];
123 q[3] = rho * u[2];
124 q[4] = E;
125 break;
126 case 3: // Velocity zero, pressure constant (so density and internal energy will be non-constant), but the velocity should stay zero and the
127 // bubble won't diffuse
128 // (for Euler, where there is no thermal conductivity)
129 P = 1.;
130 T = 1. - S_bubble * exp(1. - r * r);
131 rho = P / (R * T);
132
133 // Assign exact solution
134 q[0] = rho;
135 q[1] = rho * u[0];
136 q[2] = rho * u[1];
137 q[3] = rho * u[2];
138 q[4] = rho * (cv * T + (u[0] * u[0] + u[1] * u[1]) / 2.);
139 break;
140 case 4: // Constant nonzero velocity, pressure constant (so density and internal energy will be non-constant),
141 // It should be transported across the domain, but velocity stays constant
142 P = 1.;
143 T = 1. - S_bubble * exp(1. - r * r);
144 rho = P / (R * T);
145 u[0] = mean_velocity[0];
146 u[1] = mean_velocity[1];
147
148 // Assign exact solution
149 q[0] = rho;
150 q[1] = rho * u[0];
151 q[2] = rho * u[1];
152 q[3] = rho * u[2];
153 q[4] = rho * (cv * T + (u[0] * u[0] + u[1] * u[1]) / 2.);
154 break;
155 case 5: // non-smooth thermal bubble - cylinder
156 P = 1.;
157 T = 1. - (r < 1. ? S_bubble : 0.);
158 rho = P / (R * T);
159 u[0] = mean_velocity[0];
160 u[1] = mean_velocity[1];
161
162 // Assign exact solution
163 q[0] = rho;
164 q[1] = rho * u[0];
165 q[2] = rho * u[1];
166 q[3] = rho * u[2];
167 q[4] = rho * (cv * T + (u[0] * u[0] + u[1] * u[1]) / 2.);
168 break;
169 }
170 return 0;
171 }
172
173 // *****************************************************************************
174 // Helper function for computing flux Jacobian
175 // *****************************************************************************
ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5],const CeedScalar rho,const CeedScalar u[3],const CeedScalar E,const CeedScalar gamma)176 CEED_QFUNCTION_HELPER void ConvectiveFluxJacobian_Euler(CeedScalar dF[3][5][5], const CeedScalar rho, const CeedScalar u[3], const CeedScalar E,
177 const CeedScalar gamma) {
178 CeedScalar u_sq = u[0] * u[0] + u[1] * u[1] + u[2] * u[2]; // Velocity square
179 for (CeedInt i = 0; i < 3; i++) { // Jacobian matrices for 3 directions
180 for (CeedInt j = 0; j < 3; j++) { // Rows of each Jacobian matrix
181 dF[i][j + 1][0] = ((i == j) ? ((gamma - 1.) * (u_sq / 2.)) : 0.) - u[i] * u[j];
182 for (CeedInt k = 0; k < 3; k++) { // Columns of each Jacobian matrix
183 dF[i][0][k + 1] = ((i == k) ? 1. : 0.);
184 dF[i][j + 1][k + 1] = ((j == k) ? u[i] : 0.) + ((i == k) ? u[j] : 0.) - ((i == j) ? u[k] : 0.) * (gamma - 1.);
185 dF[i][4][k + 1] = ((i == k) ? (E * gamma / rho - (gamma - 1.) * u_sq / 2.) : 0.) - (gamma - 1.) * u[i] * u[k];
186 }
187 dF[i][j + 1][4] = ((i == j) ? (gamma - 1.) : 0.);
188 }
189 dF[i][4][0] = u[i] * ((gamma - 1.) * u_sq - E * gamma / rho);
190 dF[i][4][4] = u[i] * gamma;
191 }
192 }
193
194 // *****************************************************************************
195 // Helper function for computing Tau elements (stabilization constant)
196 // Model from:
197 // Stabilized Methods for Compressible Flows, Hughes et al 2010
198 //
199 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix
200 // Tau[i] = c_tau h[i] Xi(Pe) / rho(A[i]) (no sum)
201 //
202 // Where
203 // c_tau = stabilization constant (0.5 is reported as "optimal")
204 // h[i] = 2 length(dxdX[i])
205 // Pe = Peclet number ( Pe = sqrt(u u) / dot(dXdx,u) diffusivity )
206 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number )
207 // rho(A[i]) = spectral radius of the convective flux Jacobian i, wave speed in direction i
208 // *****************************************************************************
Tau_spatial(CeedScalar Tau_x[3],const CeedScalar dXdx[3][3],const CeedScalar u[3],const CeedScalar sound_speed,const CeedScalar c_tau)209 CEED_QFUNCTION_HELPER void Tau_spatial(CeedScalar Tau_x[3], const CeedScalar dXdx[3][3], const CeedScalar u[3], const CeedScalar sound_speed,
210 const CeedScalar c_tau) {
211 for (CeedInt i = 0; i < 3; i++) {
212 // length of element in direction i
213 CeedScalar h = 2 / sqrt(Square(dXdx[0][i]) + Square(dXdx[1][i]) + Square(dXdx[2][i]));
214 // fastest wave in direction i
215 CeedScalar fastest_wave = fabs(u[i]) + sound_speed;
216 Tau_x[i] = c_tau * h / fastest_wave;
217 }
218 }
219
220 // *****************************************************************************
221 // This QFunction sets the initial conditions for Euler traveling vortex
222 // *****************************************************************************
ICsEuler(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)223 CEED_QFUNCTION(ICsEuler)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
224 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
225 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
226
227 const EulerContext context = (EulerContext)ctx;
228
229 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
230 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
231 CeedScalar q[5] = {0.};
232
233 Exact_Euler(3, context->curr_time, x, 5, q, ctx);
234 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
235 }
236 return 0;
237 }
238
239 // *****************************************************************************
240 // This QFunction implements the following formulation of Euler equations with explicit time stepping method
241 //
242 // This is 3D Euler for compressible gas dynamics in conservation form with state variables of density, momentum density, and total energy density.
243 //
244 // State Variables: q = ( rho, U1, U2, U3, E )
245 // rho - Mass Density
246 // Ui - Momentum Density, Ui = rho ui
247 // E - Total Energy Density, E = P / (gamma - 1) + rho (u u)/2
248 //
249 // Euler Equations:
250 // drho/dt + div( U ) = 0
251 // dU/dt + div( rho (u x u) + P I3 ) = 0
252 // dE/dt + div( (E + P) u ) = 0
253 //
254 // Equation of State:
255 // P = (gamma - 1) (E - rho (u u) / 2)
256 //
257 // Constants:
258 // cv , Specific heat, constant volume
259 // cp , Specific heat, constant pressure
260 // g , Gravity
261 // gamma = cp / cv, Specific heat ratio
262 // *****************************************************************************
Euler(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)263 CEED_QFUNCTION(Euler)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
264 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
265 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1];
266 const CeedScalar(*q_data) = in[2];
267 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
268 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
269
270 EulerContext context = (EulerContext)ctx;
271 const CeedScalar c_tau = context->c_tau;
272 const CeedScalar gamma = 1.4;
273
274 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
275 // Setup
276 // -- Interp in
277 const CeedScalar rho = q[0][i];
278 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho};
279 const CeedScalar E = q[4][i];
280 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]};
281 const CeedScalar dU[3][3] = {
282 {dq[0][1][i], dq[1][1][i], dq[2][1][i]},
283 {dq[0][2][i], dq[1][2][i], dq[2][2][i]},
284 {dq[0][3][i], dq[1][3][i], dq[2][3][i]}
285 };
286 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]};
287 CeedScalar wdetJ, dXdx[3][3];
288 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
289 // dU/dx
290 CeedScalar drhodx[3] = {0.};
291 CeedScalar dEdx[3] = {0.};
292 CeedScalar dUdx[3][3] = {{0.}};
293 CeedScalar dXdxdXdxT[3][3] = {{0.}};
294 for (CeedInt j = 0; j < 3; j++) {
295 for (CeedInt k = 0; k < 3; k++) {
296 drhodx[j] += drho[k] * dXdx[k][j];
297 dEdx[j] += dE[k] * dXdx[k][j];
298 for (CeedInt l = 0; l < 3; l++) {
299 dUdx[j][k] += dU[j][l] * dXdx[l][k];
300 dXdxdXdxT[j][k] += dXdx[j][l] * dXdx[k][l]; // dXdx_j,k * dXdx_k,j
301 }
302 }
303 }
304 // Pressure
305 const CeedScalar E_kinetic = 0.5 * rho * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]), E_internal = E - E_kinetic,
306 P = E_internal * (gamma - 1.); // P = pressure
307
308 // The Physics
309 // Zero v and dv so all future terms can safely sum into it
310 for (CeedInt j = 0; j < 5; j++) {
311 v[j][i] = 0.;
312 for (CeedInt k = 0; k < 3; k++) dv[k][j][i] = 0.;
313 }
314
315 // -- Density
316 // ---- u rho
317 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]);
318 // -- Momentum
319 // ---- rho (u x u) + P I3
320 for (CeedInt j = 0; j < 3; j++) {
321 for (CeedInt k = 0; k < 3; k++) {
322 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] +
323 (rho * u[j] * u[2] + (j == 2 ? P : 0.)) * dXdx[k][2]);
324 }
325 }
326 // -- Total Energy Density
327 // ---- (E + P) u
328 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]);
329
330 // --Stabilization terms
331 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction
332 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}};
333 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma);
334
335 // ---- dqdx collects drhodx, dUdx and dEdx in one vector
336 CeedScalar dqdx[5][3];
337 for (CeedInt j = 0; j < 3; j++) {
338 dqdx[0][j] = drhodx[j];
339 dqdx[4][j] = dEdx[j];
340 for (CeedInt k = 0; k < 3; k++) dqdx[k + 1][j] = dUdx[k][j];
341 }
342
343 // ---- strong_conv = dF/dq * dq/dx (Strong convection)
344 CeedScalar strong_conv[5] = {0.};
345 for (CeedInt j = 0; j < 3; j++) {
346 for (CeedInt k = 0; k < 5; k++) {
347 for (CeedInt l = 0; l < 5; l++) strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j];
348 }
349 }
350
351 // Stabilization
352 // -- Tau elements
353 const CeedScalar sound_speed = sqrt(gamma * P / rho);
354 CeedScalar Tau_x[3] = {0.};
355 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau);
356
357 // -- Stabilization method: none or SU
358 CeedScalar stab[5][3] = {{0.}};
359 switch (context->stabilization) {
360 case 0: // Galerkin
361 break;
362 case 1: // SU
363 for (CeedInt j = 0; j < 3; j++) {
364 for (CeedInt k = 0; k < 5; k++) {
365 for (CeedInt l = 0; l < 5; l++) stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l];
366 }
367 }
368
369 for (CeedInt j = 0; j < 5; j++) {
370 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]);
371 }
372 break;
373 case 2: // SUPG is not implemented for explicit scheme
374 break;
375 }
376 }
377 return 0;
378 }
379
380 // *****************************************************************************
381 // This QFunction implements the Euler equations with (mentioned above) with implicit time stepping method
382 // *****************************************************************************
IFunction_Euler(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)383 CEED_QFUNCTION(IFunction_Euler)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
384 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
385 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1];
386 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
387 const CeedScalar(*q_data) = in[3];
388 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
389 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1];
390
391 EulerContext context = (EulerContext)ctx;
392 const CeedScalar c_tau = context->c_tau;
393 const CeedScalar gamma = 1.4;
394
395 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
396 // Setup
397 // -- Interp in
398 const CeedScalar rho = q[0][i];
399 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho};
400 const CeedScalar E = q[4][i];
401 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]};
402 const CeedScalar dU[3][3] = {
403 {dq[0][1][i], dq[1][1][i], dq[2][1][i]},
404 {dq[0][2][i], dq[1][2][i], dq[2][2][i]},
405 {dq[0][3][i], dq[1][3][i], dq[2][3][i]}
406 };
407 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]};
408 CeedScalar wdetJ, dXdx[3][3];
409 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
410 // dU/dx
411 CeedScalar drhodx[3] = {0.};
412 CeedScalar dEdx[3] = {0.};
413 CeedScalar dUdx[3][3] = {{0.}};
414 CeedScalar dXdxdXdxT[3][3] = {{0.}};
415 for (CeedInt j = 0; j < 3; j++) {
416 for (CeedInt k = 0; k < 3; k++) {
417 drhodx[j] += drho[k] * dXdx[k][j];
418 dEdx[j] += dE[k] * dXdx[k][j];
419 for (CeedInt l = 0; l < 3; l++) {
420 dUdx[j][k] += dU[j][l] * dXdx[l][k];
421 dXdxdXdxT[j][k] += dXdx[j][l] * dXdx[k][l]; // dXdx_j,k * dXdx_k,j
422 }
423 }
424 }
425 const CeedScalar E_kinetic = 0.5 * rho * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]), E_internal = E - E_kinetic,
426 P = E_internal * (gamma - 1.); // P = pressure
427
428 // The Physics
429 // Zero v and dv so all future terms can safely sum into it
430 for (CeedInt j = 0; j < 5; j++) {
431 v[j][i] = 0.;
432 for (CeedInt k = 0; k < 3; k++) dv[k][j][i] = 0.;
433 }
434 //-----mass matrix
435 for (CeedInt j = 0; j < 5; j++) v[j][i] += wdetJ * q_dot[j][i];
436
437 // -- Density
438 // ---- u rho
439 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]);
440 // -- Momentum
441 // ---- rho (u x u) + P I3
442 for (CeedInt j = 0; j < 3; j++) {
443 for (CeedInt k = 0; k < 3; k++) {
444 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] +
445 (rho * u[j] * u[2] + (j == 2 ? P : 0.)) * dXdx[k][2]);
446 }
447 }
448 // -- Total Energy Density
449 // ---- (E + P) u
450 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]);
451
452 // -- Stabilization terms
453 // ---- jacob_F_conv[3][5][5] = dF(convective)/dq at each direction
454 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}};
455 ConvectiveFluxJacobian_Euler(jacob_F_conv, rho, u, E, gamma);
456
457 // ---- dqdx collects drhodx, dUdx and dEdx in one vector
458 CeedScalar dqdx[5][3];
459 for (CeedInt j = 0; j < 3; j++) {
460 dqdx[0][j] = drhodx[j];
461 dqdx[4][j] = dEdx[j];
462 for (CeedInt k = 0; k < 3; k++) dqdx[k + 1][j] = dUdx[k][j];
463 }
464
465 // ---- strong_conv = dF/dq * dq/dx (Strong convection)
466 CeedScalar strong_conv[5] = {0.};
467 for (CeedInt j = 0; j < 3; j++) {
468 for (CeedInt k = 0; k < 5; k++) {
469 for (CeedInt l = 0; l < 5; l++) strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j];
470 }
471 }
472
473 // ---- Strong residual
474 CeedScalar strong_res[5];
475 for (CeedInt j = 0; j < 5; j++) strong_res[j] = q_dot[j][i] + strong_conv[j];
476
477 // Stabilization
478 // -- Tau elements
479 const CeedScalar sound_speed = sqrt(gamma * P / rho);
480 CeedScalar Tau_x[3] = {0.};
481 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau);
482
483 // -- Stabilization method: none, SU, or SUPG
484 CeedScalar stab[5][3] = {{0.}};
485 switch (context->stabilization) {
486 case 0: // Galerkin
487 break;
488 case 1: // SU
489 for (CeedInt j = 0; j < 3; j++) {
490 for (CeedInt k = 0; k < 5; k++) {
491 for (CeedInt l = 0; l < 5; l++) stab[k][j] += jacob_F_conv[j][k][l] * Tau_x[j] * strong_conv[l];
492 }
493 }
494
495 for (CeedInt j = 0; j < 5; j++) {
496 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]);
497 }
498 break;
499 case 2: // SUPG
500 for (CeedInt j = 0; j < 3; j++) {
501 for (CeedInt k = 0; k < 5; k++) {
502 for (CeedInt l = 0; l < 5; l++) stab[k][j] = jacob_F_conv[j][k][l] * Tau_x[j] * strong_res[l];
503 }
504 }
505
506 for (CeedInt j = 0; j < 5; j++) {
507 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]);
508 }
509 break;
510 }
511 }
512 return 0;
513 }
514 // *****************************************************************************
515 // This QFunction sets the inflow boundary conditions for the traveling vortex problem.
516 //
517 // Prescribed T_inlet and P_inlet are converted to conservative variables and applied weakly.
518 // *****************************************************************************
TravelingVortex_Inflow(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)519 CEED_QFUNCTION(TravelingVortex_Inflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
520 const CeedScalar(*q_data_sur) = in[2];
521 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
522
523 EulerContext context = (EulerContext)ctx;
524 const int euler_test = context->euler_test;
525 const bool is_implicit = context->implicit;
526 CeedScalar *mean_velocity = context->mean_velocity;
527 const CeedScalar cv = 2.5;
528 const CeedScalar R = 1.;
529 CeedScalar T_inlet;
530 CeedScalar P_inlet;
531
532 // For test cases 1 and 3 the background velocity is zero
533 if (euler_test == 1 || euler_test == 3) {
534 for (CeedInt i = 0; i < 3; i++) mean_velocity[i] = 0.;
535 }
536
537 // For test cases 1 and 2, T_inlet = T_inlet = 0.4
538 if (euler_test == 1 || euler_test == 2) T_inlet = P_inlet = .4;
539 else T_inlet = P_inlet = 1.;
540
541 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
542 CeedScalar wdetJb, normal[3];
543 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
544 wdetJb *= is_implicit ? -1. : 1.;
545
546 // face_normal = Normal vector of the face
547 const CeedScalar face_normal = Dot3(normal, mean_velocity);
548 // The Physics
549 // Zero v so all future terms can safely sum into it
550 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.;
551
552 // Implementing in/outflow BCs
553 if (face_normal > 0) {
554 } else { // inflow
555 const CeedScalar rho_inlet = P_inlet / (R * T_inlet);
556 const CeedScalar E_kinetic_inlet = (mean_velocity[0] * mean_velocity[0] + mean_velocity[1] * mean_velocity[1]) / 2.;
557 // incoming total energy
558 const CeedScalar E_inlet = rho_inlet * (cv * T_inlet + E_kinetic_inlet);
559
560 // The Physics
561 // -- Density
562 v[0][i] -= wdetJb * rho_inlet * face_normal;
563
564 // -- Momentum
565 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho_inlet * face_normal * mean_velocity[j] + normal[j] * P_inlet);
566
567 // -- Total Energy Density
568 v[4][i] -= wdetJb * face_normal * (E_inlet + P_inlet);
569 }
570 }
571 return 0;
572 }
573
574 // *****************************************************************************
575 // This QFunction sets the outflow boundary conditions for the Euler solver.
576 //
577 // Outflow BCs:
578 // The validity of the weak form of the governing equations is extended to the outflow.
579 // *****************************************************************************
Euler_Outflow(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)580 CEED_QFUNCTION(Euler_Outflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
581 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
582 const CeedScalar(*q_data_sur) = in[2];
583 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
584
585 EulerContext context = (EulerContext)ctx;
586 const bool is_implicit = context->implicit;
587 CeedScalar *mean_velocity = context->mean_velocity;
588
589 const CeedScalar gamma = 1.4;
590
591 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
592 // Setup
593 // -- Interp in
594 const CeedScalar rho = q[0][i];
595 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho};
596 const CeedScalar E = q[4][i];
597
598 CeedScalar wdetJb, normal[3];
599 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
600 wdetJb *= is_implicit ? -1. : 1.;
601
602 // face_normal = Normal vector of the face
603 const CeedScalar face_normal = Dot3(normal, mean_velocity);
604 // The Physics
605 // Zero v so all future terms can safely sum into it
606 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0;
607
608 // Implementing in/outflow BCs
609 if (face_normal > 0) { // outflow
610 const CeedScalar E_kinetic = (u[0] * u[0] + u[1] * u[1]) / 2.;
611 const CeedScalar P = (E - E_kinetic * rho) * (gamma - 1.); // pressure
612 const CeedScalar u_normal = Dot3(normal, u); // Normal velocity
613 // The Physics
614 // -- Density
615 v[0][i] -= wdetJb * rho * u_normal;
616
617 // -- Momentum
618 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho * u_normal * u[j] + normal[j] * P);
619
620 // -- Total Energy Density
621 v[4][i] -= wdetJb * u_normal * (E + P);
622 }
623 }
624 return 0;
625 }
626