xref: /honee/qfunctions/taylorgreen.h (revision 475f0cac5d40259768f4556cf888e8f2448554cb)
1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 #include <ceed/types.h>
4 #ifndef CEED_RUNNING_JIT_PASS
5 #include <math.h>
6 #endif
7 
8 #include "newtonian_state.h"
9 #include "newtonian_types.h"
10 #include "utils.h"
11 
12 typedef struct TaylorGreenContext_ *TaylorGreenContext;
13 struct TaylorGreenContext_ {
14   StatePrimitive                   reference;
15   struct NewtonianIdealGasContext_ newt_ctx;
16   CeedScalar                       lx;
17   CeedScalar                       ly;
18   CeedScalar                       lz;
19   CeedScalar                       u[3];
20 };
21 
22 // @brief Set initial condition for Taylor-Green Vortex problem
ICsTaylorGreen(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)23 CEED_QFUNCTION(ICsTaylorGreen)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
24   const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
25 
26   CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
27 
28   const TaylorGreenContext    context   = (TaylorGreenContext)ctx;
29   const NewtonianIGProperties gas       = context->newt_ctx.gas;
30   CeedScalar                  R         = GasConstant(gas);
31   StatePrimitive              reference = context->reference;
32   const CeedScalar            V0        = Norm3(reference.velocity);
33   const CeedScalar            density0  = reference.pressure / (reference.temperature * R);
34 
35   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
36     CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
37     CeedScalar q[5], Y[5];
38     ScaleN(x, 2 * M_PI / context->lx, 3);
39 
40     Y[0] = reference.pressure + (density0 * Square(V0) / 16) * (cos(2 * x[0]) + cos(2 * x[1])) * (cos(2 * x[2] + 2));
41     Y[1] = context->u[0] + V0 * sin(x[0]) * cos(x[1]) * cos(x[2]);
42     Y[2] = context->u[1] - V0 * cos(x[0]) * sin(x[1]) * cos(x[2]);
43     Y[3] = context->u[2];
44     Y[4] = reference.temperature;
45 
46     State s = StateFromY(gas, Y);
47     StateToQ(gas, s, q, context->newt_ctx.state_var);
48     for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
49   }
50   return 0;
51 }
52