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_ gas; 39 }; 40 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 NewtonianIdealGasContext gas = &context->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 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 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 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 95 CEED_QFUNCTION_HELPER CeedScalar ChildStatsCollectionTest_Exact(const CeedScalar x_i[3]) { return x_i[0] + Square(x_i[1]); } 96 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 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