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