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