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