1 // Copyright (c) 2017-2026, 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 #include <ceed/types.h>
8
9 #include "newtonian_state.h"
10 #include "turb_stats_types.h"
11 #include "utils.h"
12
ChildStatsCollection(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,StateVariable state_var)13 CEED_QFUNCTION_HELPER int ChildStatsCollection(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
14 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
15 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
16 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
17
18 Turbulence_SpanStatsContext context = (Turbulence_SpanStatsContext)ctx;
19 NewtonianIdealGasContext gas = &context->gas;
20 CeedScalar delta_t = context->solution_time - context->previous_time;
21
22 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
23 const CeedScalar wdetJ = q_data[0][i] * delta_t;
24
25 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
26 const State s = StateFromQ(gas, qi, state_var);
27
28 v[TURB_MEAN_DENSITY][i] = wdetJ * s.U.density;
29 v[TURB_MEAN_PRESSURE][i] = wdetJ * s.Y.pressure;
30 v[TURB_MEAN_PRESSURE_SQUARED][i] = wdetJ * Square(s.Y.pressure);
31 v[TURB_MEAN_PRESSURE_VELOCITY_X][i] = wdetJ * s.Y.pressure * s.Y.velocity[0];
32 v[TURB_MEAN_PRESSURE_VELOCITY_Y][i] = wdetJ * s.Y.pressure * s.Y.velocity[1];
33 v[TURB_MEAN_PRESSURE_VELOCITY_Z][i] = wdetJ * s.Y.pressure * s.Y.velocity[2];
34 v[TURB_MEAN_DENSITY_TEMPERATURE][i] = wdetJ * s.U.density * s.Y.temperature;
35 v[TURB_MEAN_DENSITY_TEMPERATURE_FLUX_X][i] = wdetJ * s.U.density * s.Y.temperature * s.Y.velocity[0];
36 v[TURB_MEAN_DENSITY_TEMPERATURE_FLUX_Y][i] = wdetJ * s.U.density * s.Y.temperature * s.Y.velocity[1];
37 v[TURB_MEAN_DENSITY_TEMPERATURE_FLUX_Z][i] = wdetJ * s.U.density * s.Y.temperature * s.Y.velocity[2];
38 v[TURB_MEAN_MOMENTUM_X][i] = wdetJ * s.U.momentum[0];
39 v[TURB_MEAN_MOMENTUM_Y][i] = wdetJ * s.U.momentum[1];
40 v[TURB_MEAN_MOMENTUM_Z][i] = wdetJ * s.U.momentum[2];
41 v[TURB_MEAN_MOMENTUMFLUX_XX][i] = wdetJ * s.U.momentum[0] * s.Y.velocity[0];
42 v[TURB_MEAN_MOMENTUMFLUX_YY][i] = wdetJ * s.U.momentum[1] * s.Y.velocity[1];
43 v[TURB_MEAN_MOMENTUMFLUX_ZZ][i] = wdetJ * s.U.momentum[2] * s.Y.velocity[2];
44 v[TURB_MEAN_MOMENTUMFLUX_YZ][i] = wdetJ * s.U.momentum[1] * s.Y.velocity[2];
45 v[TURB_MEAN_MOMENTUMFLUX_XZ][i] = wdetJ * s.U.momentum[0] * s.Y.velocity[2];
46 v[TURB_MEAN_MOMENTUMFLUX_XY][i] = wdetJ * s.U.momentum[0] * s.Y.velocity[1];
47 v[TURB_MEAN_VELOCITY_X][i] = wdetJ * s.Y.velocity[0];
48 v[TURB_MEAN_VELOCITY_Y][i] = wdetJ * s.Y.velocity[1];
49 v[TURB_MEAN_VELOCITY_Z][i] = wdetJ * s.Y.velocity[2];
50 }
51 return 0;
52 }
53
ChildStatsCollection_Conserv(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)54 CEED_QFUNCTION(ChildStatsCollection_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
55 return ChildStatsCollection(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
56 }
57
ChildStatsCollection_Prim(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)58 CEED_QFUNCTION(ChildStatsCollection_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
59 return ChildStatsCollection(ctx, Q, in, out, STATEVAR_PRIMITIVE);
60 }
61
ChildStatsCollection_Entropy(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)62 CEED_QFUNCTION(ChildStatsCollection_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
63 return ChildStatsCollection(ctx, Q, in, out, STATEVAR_ENTROPY);
64 }
65
66 // QFunctions for testing
ChildStatsCollectionTest_Exact(const CeedScalar x_i[3])67 CEED_QFUNCTION_HELPER CeedScalar ChildStatsCollectionTest_Exact(const CeedScalar x_i[3]) { return x_i[0] + Square(x_i[1]); }
68
ChildStatsCollectionMMSTest(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)69 CEED_QFUNCTION(ChildStatsCollectionMMSTest)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
70 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
71 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
72 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
73
74 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
75 const CeedScalar t = context->time;
76
77 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
78 const CeedScalar wdetJ = q_data[0][i];
79 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
80
81 // set spanwise domain to [0,1] and integrate from t \in [0,1] to recover exact solution
82 v[0][i] = wdetJ * (ChildStatsCollectionTest_Exact(x_i) + t - 0.5) * 4 * Cube(x_i[2]);
83 for (int j = 1; j < 22; j++) v[j][i] = 0;
84 }
85 return 0;
86 }
87
ChildStatsCollectionMMSTest_Error(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)88 CEED_QFUNCTION(ChildStatsCollectionMMSTest_Error)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
89 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
90 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
91 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
92 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
93
94 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
95 const CeedScalar wdetJ = q_data[0][i];
96 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
97
98 v[0][i] = wdetJ * Square(ChildStatsCollectionTest_Exact(x_i) - q[0][i]);
99 }
100 return 0;
101 }
102