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