xref: /libCEED/examples/fluids/qfunctions/taylorgreen.h (revision e64bb3f3ed2986a0c10dec3b47522d734c6e367d)
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 ///
10 
11 #ifndef taylorgreen_h
12 #define taylorgreen_h
13 
14 #include <ceed.h>
15 #include <math.h>
16 
17 #include "newtonian_state.h"
18 #include "newtonian_types.h"
19 #include "utils.h"
20 
21 // @brief Set initial condition for Taylor-Green Vortex problem
22 CEED_QFUNCTION(ICsTaylorGreen)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
23   const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
24 
25   CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
26 
27   const SetupContext                context   = (SetupContext)ctx;
28   struct NewtonianIdealGasContext_ *gas       = &context->gas;
29   CeedScalar                        R         = GasConstant(gas);
30   StatePrimitive                    reference = context->reference;
31   const CeedScalar                  V0        = sqrt(Dot3(reference.velocity, reference.velocity));
32   const CeedScalar                  density0  = reference.pressure / (reference.temperature * R);
33   const CeedScalar                  x0[3]     = {0.};
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] = {0}, 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] = V0 * sin(x[0]) * cos(x[1]) * cos(x[2]);
42     Y[2] = -V0 * cos(x[0]) * sin(x[1]) * cos(x[2]);
43     Y[3] = 0;
44     Y[4] = reference.temperature;
45 
46     State s = StateFromY(gas, Y, x0);
47     switch (gas->state_var) {
48       case STATEVAR_CONSERVATIVE:
49         UnpackState_U(s.U, q);
50         break;
51       case STATEVAR_PRIMITIVE:
52         UnpackState_Y(s.Y, q);
53         break;
54     }
55 
56     for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
57   }
58   return 0;
59 }
60 
61 #endif  // taylorgreen_h
62