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 /// Thermodynamic wave propogation for testing freestream/non-reflecting boundary conditions. Proposed in Mengaldo et. al. 2014
10 #include <ceed/types.h>
11 #ifndef CEED_RUNNING_JIT_PASS
12 #include <math.h>
13 #endif
14
15 #include "newtonian_state.h"
16 #include "utils.h"
17
18 typedef struct GaussianWaveContext_ *GaussianWaveContext;
19 struct GaussianWaveContext_ {
20 CeedScalar epicenter[3]; // Location of the perturbation
21 CeedScalar width; // Controls width of the perturbation
22 CeedScalar amplitude; // Amplitude of the perturbation
23 State S_infty; // Flow state at infinity
24 struct NewtonianIdealGasContext_ newt_ctx;
25 };
26
IC_GaussianWave(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)27 CEED_QFUNCTION_HELPER int IC_GaussianWave(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
28 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
29
30 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
31
32 const GaussianWaveContext context = (GaussianWaveContext)ctx;
33 const NewtonianIdealGasContext newt_ctx = &context->newt_ctx;
34
35 const CeedScalar amplitude = context->amplitude;
36 const CeedScalar width = context->width;
37 const State S_infty = context->S_infty;
38 const CeedScalar xc = context->epicenter[0];
39 const CeedScalar yc = context->epicenter[1];
40
41 const CeedScalar gamma = HeatCapacityRatio(newt_ctx);
42
43 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
44 CeedScalar U[5] = {0.}, qi[5] = {0.};
45 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]};
46 const CeedScalar x0 = x[0] - xc;
47 const CeedScalar y0 = x[1] - yc;
48 const CeedScalar e_kinetic = 0.5 * S_infty.U.density * Dot3(S_infty.Y.velocity, S_infty.Y.velocity);
49
50 const CeedScalar perturbation = 1 + amplitude * exp(-(Square(x0) + Square(y0)) / (2 * Square(width)));
51
52 U[0] = S_infty.U.density * perturbation;
53 U[1] = S_infty.Y.velocity[0] * U[0];
54 U[2] = S_infty.Y.velocity[1] * U[0];
55 U[3] = S_infty.Y.velocity[2] * U[0];
56 U[4] = S_infty.Y.pressure / (gamma - 1) * perturbation + e_kinetic;
57
58 State initCond = StateFromU(newt_ctx, U);
59 StateToQ(newt_ctx, initCond, qi, state_var);
60
61 for (CeedInt j = 0; j < 5; j++) q0[j][i] = qi[j];
62 }
63
64 return 0;
65 }
66
IC_GaussianWave_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)67 CEED_QFUNCTION(IC_GaussianWave_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
68 return IC_GaussianWave(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
69 }
70
IC_GaussianWave_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)71 CEED_QFUNCTION(IC_GaussianWave_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
72 return IC_GaussianWave(ctx, Q, in, out, STATEVAR_PRIMITIVE);
73 }
74
IC_GaussianWave_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)75 CEED_QFUNCTION(IC_GaussianWave_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
76 return IC_GaussianWave(ctx, Q, in, out, STATEVAR_ENTROPY);
77 }
78