1*9ba83ac0SJeremy L Thompson // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
214712a6bSJames Wright // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
314712a6bSJames Wright //
414712a6bSJames Wright // SPDX-License-Identifier: BSD-2-Clause
514712a6bSJames Wright //
614712a6bSJames Wright // This file is part of CEED: http://github.com/ceed
7c0b5abf0SJeremy L Thompson #include <ceed/types.h>
8c0b5abf0SJeremy L Thompson #ifndef CEED_RUNNING_JIT_PASS
914712a6bSJames Wright #include <math.h>
10c0b5abf0SJeremy L Thompson #endif
1114712a6bSJames Wright
1214712a6bSJames Wright #include "newtonian_state.h"
1314712a6bSJames Wright #include "newtonian_types.h"
1414712a6bSJames Wright #include "utils.h"
1514712a6bSJames Wright
1614712a6bSJames Wright // @brief Set initial condition for Taylor-Green Vortex problem
ICsTaylorGreen(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)1714712a6bSJames Wright CEED_QFUNCTION(ICsTaylorGreen)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
1814712a6bSJames Wright const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
1914712a6bSJames Wright
2014712a6bSJames Wright CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
2114712a6bSJames Wright
2214712a6bSJames Wright const SetupContext context = (SetupContext)ctx;
23a2d72b6fSJames Wright const NewtonianIdealGasContext gas = &context->gas;
2414712a6bSJames Wright CeedScalar R = GasConstant(gas);
2514712a6bSJames Wright StatePrimitive reference = context->reference;
2614712a6bSJames Wright const CeedScalar V0 = sqrt(Dot3(reference.velocity, reference.velocity));
2714712a6bSJames Wright const CeedScalar density0 = reference.pressure / (reference.temperature * R);
2814712a6bSJames Wright
2914712a6bSJames Wright CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
3014712a6bSJames Wright CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
3114712a6bSJames Wright CeedScalar q[5] = {0}, Y[5];
3214712a6bSJames Wright ScaleN(x, 2 * M_PI / context->lx, 3);
3314712a6bSJames Wright
3414712a6bSJames Wright Y[0] = reference.pressure + (density0 * Square(V0) / 16) * (cos(2 * x[0]) + cos(2 * x[1])) * (cos(2 * x[2] + 2));
3514712a6bSJames Wright Y[1] = V0 * sin(x[0]) * cos(x[1]) * cos(x[2]);
3614712a6bSJames Wright Y[2] = -V0 * cos(x[0]) * sin(x[1]) * cos(x[2]);
3714712a6bSJames Wright Y[3] = 0;
3814712a6bSJames Wright Y[4] = reference.temperature;
3914712a6bSJames Wright
403bd61617SKenneth E. Jansen State s = StateFromY(gas, Y);
41a2d72b6fSJames Wright StateToQ(gas, s, q, gas->state_var);
4214712a6bSJames Wright for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
4314712a6bSJames Wright }
4414712a6bSJames Wright return 0;
4514712a6bSJames Wright }
46