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