1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 /// @file 4 /// Functions for setting up and performing spanwise-statistics collection of turbulence quantities 5 6 #include "../../qfunctions/spanstats/turbulence.h" 7 8 #include <ceed.h> 9 #include <petscdmplex.h> 10 #include <petscsf.h> 11 #include <spanstats.h> 12 13 #include <navierstokes.h> 14 #include <petsc_ops.h> 15 16 // @brief Create CeedOperator for statistics collection 17 static PetscErrorCode CreateStatisticCollectionOperator(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_data, ProblemData problem) { 18 Ceed ceed = honee->ceed; 19 CeedInt num_comp_stats = spanstats->num_comp_stats, q_data_size; 20 PetscInt num_comp_q, num_comp_x; 21 Turbulence_SpanStatsContext collect_ctx; 22 NewtonianIdealGasContext newtonian_ig_ctx; 23 CeedQFunctionContext collect_qfctx; 24 CeedQFunction qf_stats_collect; 25 CeedOperator op_stats_collect; 26 CeedVector q_data; 27 CeedElemRestriction elem_restr_qd, elem_restr_q; 28 CeedBasis basis_q; 29 30 PetscFunctionBeginUser; 31 PetscCall(DMGetCoordinateNumComps(honee->dm, &num_comp_x)); 32 PetscCall(DMGetFieldNumComps(honee->dm, 0, &num_comp_q)); 33 34 // Create Operator for statistics collection 35 switch (honee->phys->state_var) { 36 case STATEVAR_PRIMITIVE: 37 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_Prim, ChildStatsCollection_Prim_loc, &qf_stats_collect)); 38 break; 39 case STATEVAR_CONSERVATIVE: 40 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_Conserv, ChildStatsCollection_Conserv_loc, &qf_stats_collect)); 41 break; 42 case STATEVAR_ENTROPY: 43 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_Entropy, ChildStatsCollection_Entropy_loc, &qf_stats_collect)); 44 break; 45 } 46 47 if (spanstats->do_mms_test) { 48 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stats_collect)); 49 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollectionMMSTest, ChildStatsCollectionMMSTest_loc, &qf_stats_collect)); 50 } 51 52 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &elem_restr_q)); 53 PetscCall(DMPlexCeedBasisCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &basis_q)); 54 PetscCall(QDataGet(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, honee->elem_restr_x, honee->basis_x, honee->x_coord, &elem_restr_qd, 55 &q_data, &q_data_size)); 56 57 { // Setup Collection Context 58 PetscCall(PetscNew(&collect_ctx)); 59 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 60 collect_ctx->newt_ctx = *newtonian_ig_ctx; 61 62 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &collect_qfctx)); 63 PetscCallCeed(ceed, CeedQFunctionContextSetData(collect_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*collect_ctx), collect_ctx)); 64 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(collect_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 65 66 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(collect_qfctx, "solution time", 67 offsetof(struct Turbulence_SpanStatsContext_, solution_time), 1, "Current solution time")); 68 PetscCallCeed(ceed, 69 CeedQFunctionContextRegisterDouble(collect_qfctx, "previous time", offsetof(struct Turbulence_SpanStatsContext_, previous_time), 1, 70 "Previous time statistics collection was done")); 71 72 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 73 } 74 75 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_stats_collect, collect_qfctx)); 76 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&collect_qfctx)); 77 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_stats_collect, "q", num_comp_q, CEED_EVAL_INTERP)); 78 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_stats_collect, "q_data", q_data_size, CEED_EVAL_NONE)); 79 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_stats_collect, "x", num_comp_x, CEED_EVAL_INTERP)); 80 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_stats_collect, "v", num_comp_stats, CEED_EVAL_NONE)); 81 82 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stats_collect, NULL, NULL, &op_stats_collect)); 83 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 84 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "q_data", elem_restr_qd, CEED_BASIS_NONE, q_data)); 85 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "x", honee->elem_restr_x, honee->basis_x, honee->x_coord)); 86 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "v", stats_data->elem_restr_child_colloc, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 87 88 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_stats_collect, "solution time", &spanstats->solution_time_label)); 89 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_stats_collect, "previous time", &spanstats->previous_time_label)); 90 91 PetscCall(OperatorApplyContextCreate(honee->dm, spanstats->dm, honee->ceed, op_stats_collect, honee->q_ceed, NULL, NULL, NULL, 92 &spanstats->op_stats_collect_ctx)); 93 94 PetscCall(CeedOperatorCreateLocalVecs(op_stats_collect, DMReturnVecType(spanstats->dm), PETSC_COMM_SELF, NULL, &spanstats->Child_Stats_loc)); 95 PetscCall(VecZeroEntries(spanstats->Child_Stats_loc)); 96 97 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 98 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 99 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 100 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 101 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stats_collect)); 102 PetscCallCeed(ceed, CeedOperatorDestroy(&op_stats_collect)); 103 PetscFunctionReturn(PETSC_SUCCESS); 104 } 105 106 // @brief Creates operator for calculating error of method of manufactured solution (MMS) test 107 static PetscErrorCode SetupMMSErrorChecking(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_data) { 108 Ceed ceed = honee->ceed; 109 CeedInt num_comp_stats = spanstats->num_comp_stats, num_comp_x, q_data_size; 110 CeedQFunction qf_error; 111 CeedOperator op_error; 112 CeedVector x_ceed, y_ceed; 113 CeedVector q_data; 114 CeedElemRestriction elem_restr_parent_qd; 115 116 PetscFunctionBeginUser; 117 PetscCall(QDataGet(ceed, spanstats->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, stats_data->elem_restr_parent_x, stats_data->basis_x, 118 stats_data->x_coord, &elem_restr_parent_qd, &q_data, &q_data_size)); 119 PetscCallCeed(ceed, CeedBasisGetNumComponents(stats_data->basis_x, &num_comp_x)); 120 121 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollectionMMSTest_Error, ChildStatsCollectionMMSTest_Error_loc, &qf_error)); 122 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_error, "q", num_comp_stats, CEED_EVAL_INTERP)); 123 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_error, "qdata", q_data_size, CEED_EVAL_NONE)); 124 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_error, "x", num_comp_x, CEED_EVAL_INTERP)); 125 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_error, "v", num_comp_stats, CEED_EVAL_INTERP)); 126 127 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_error, NULL, NULL, &op_error)); 128 PetscCallCeed(ceed, CeedOperatorSetField(op_error, "q", stats_data->elem_restr_parent_stats, stats_data->basis_stats, CEED_VECTOR_ACTIVE)); 129 PetscCallCeed(ceed, CeedOperatorSetField(op_error, "qdata", elem_restr_parent_qd, CEED_BASIS_NONE, q_data)); 130 PetscCallCeed(ceed, CeedOperatorSetField(op_error, "x", stats_data->elem_restr_parent_x, stats_data->basis_x, stats_data->x_coord)); 131 PetscCallCeed(ceed, CeedOperatorSetField(op_error, "v", stats_data->elem_restr_parent_stats, stats_data->basis_stats, CEED_VECTOR_ACTIVE)); 132 133 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(stats_data->elem_restr_parent_stats, &x_ceed, NULL)); 134 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(stats_data->elem_restr_parent_stats, &y_ceed, NULL)); 135 PetscCall(OperatorApplyContextCreate(spanstats->dm, spanstats->dm, honee->ceed, op_error, x_ceed, y_ceed, NULL, NULL, &spanstats->mms_error_ctx)); 136 137 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 138 PetscCallCeed(ceed, CeedVectorDestroy(&x_ceed)); 139 PetscCallCeed(ceed, CeedVectorDestroy(&y_ceed)); 140 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_parent_qd)); 141 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_error)); 142 PetscCallCeed(ceed, CeedOperatorDestroy(&op_error)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 // @brief Setup for statistics collection for turbulence quantities 147 PetscErrorCode SpanwiseStatisticsSetup_Turbulence(TS ts, PetscViewerAndFormat *ctx) { 148 Honee honee; 149 SpanStatsSetupData stats_setup_data; 150 SpanStatsCtx spanstats; 151 const char *prefix = "ts_monitor_spanstats_turbulence_"; 152 PetscBool is_ascii; 153 154 PetscFunctionBeginUser; 155 PetscCall(TSGetApplicationContext(ts, &honee)); 156 PetscCall(PetscObjectTypeCompare((PetscObject)ctx->viewer, PETSCVIEWERASCII, &is_ascii)); 157 PetscCheck(!is_ascii || honee->app_ctx->test_type != TESTTYPE_NONE, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, 158 "Turbulence spanwise statistics does not support ASCII viewers"); 159 160 PetscCall(SpanwiseStatisticsSetupInitialize(honee, honee->app_ctx->degree, TURB_NUM_COMPONENTS, prefix, &stats_setup_data, &spanstats)); 161 162 { // Create Section for data 163 PetscSection section; 164 165 PetscCall(DMGetLocalSection(spanstats->dm, §ion)); 166 PetscCall(PetscSectionSetFieldName(section, 0, "")); 167 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY, "MeanDensity")); 168 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE, "MeanPressure")); 169 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_SQUARED, "MeanPressureSquared")); 170 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_VELOCITY_X, "MeanPressureVelocityX")); 171 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_VELOCITY_Y, "MeanPressureVelocityY")); 172 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_VELOCITY_Z, "MeanPressureVelocityZ")); 173 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE, "MeanDensityTemperature")); 174 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE_FLUX_X, "MeanDensityTemperatureFluxX")); 175 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE_FLUX_Y, "MeanDensityTemperatureFluxY")); 176 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE_FLUX_Z, "MeanDensityTemperatureFluxZ")); 177 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUM_X, "MeanMomentumX")); 178 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUM_Y, "MeanMomentumY")); 179 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUM_Z, "MeanMomentumZ")); 180 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_XX, "MeanMomentumFluxXX")); 181 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_YY, "MeanMomentumFluxYY")); 182 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_ZZ, "MeanMomentumFluxZZ")); 183 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_YZ, "MeanMomentumFluxYZ")); 184 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_XZ, "MeanMomentumFluxXZ")); 185 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_XY, "MeanMomentumFluxXY")); 186 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_VELOCITY_X, "MeanVelocityX")); 187 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_VELOCITY_Y, "MeanVelocityY")); 188 PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_VELOCITY_Z, "MeanVelocityZ")); 189 } 190 191 // Create CeedOperators for statistics collection 192 PetscCall(CreateStatisticCollectionOperator(honee, spanstats, stats_setup_data, honee->problem_data)); 193 194 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ts_monitor_spanstats_turbulence_mms", &spanstats->do_mms_test, NULL)); 195 if (spanstats->do_mms_test) { 196 PetscCall(SetupMMSErrorChecking(honee, spanstats, stats_setup_data)); 197 } 198 199 ctx->data = spanstats; 200 ctx->data_destroy = (PetscCtxDestroyFn *)SpanStatsCtxDestroy; 201 202 PetscCall(SpanwiseStatisticsSetupFinalize(ts, honee, spanstats, ctx, &stats_setup_data)); 203 PetscFunctionReturn(PETSC_SUCCESS); 204 } 205 206 // @brief TSMonitor for collection and processing of turbulence statistics 207 PetscErrorCode TSMonitor_SpanwiseStatisticsTurbulence(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, PetscViewerAndFormat *ctx) { 208 Honee honee; 209 Vec stats; 210 TSConvergedReason reason; 211 SpanStatsCtx spanstats = ctx->data; 212 PetscInt collect_interval = spanstats->collect_interval, viewer_interval = ctx->view_interval; 213 214 PetscFunctionBeginUser; 215 PetscCall(TSGetApplicationContext(ts, &honee)); 216 PetscCall(TSGetConvergedReason(ts, &reason)); 217 // Do not collect or process on the first step of the run (ie. on the initial condition) 218 if (steps == honee->app_ctx->cont_steps) PetscFunctionReturn(PETSC_SUCCESS); 219 220 PetscBool run_processing_and_viewer = (steps % viewer_interval == 0 && viewer_interval != -1) || reason != TS_CONVERGED_ITERATING; 221 222 if (steps % collect_interval == 0 || run_processing_and_viewer) { 223 PetscCall(SpanwiseStatisticsCollect(honee, spanstats, solution_time, Q)); 224 225 if (run_processing_and_viewer) { 226 PetscCall(DMSetOutputSequenceNumber(spanstats->dm, steps, solution_time)); 227 PetscCall(DMGetGlobalVector(spanstats->dm, &stats)); 228 PetscCall(SpanwiseStatisticsProcess(honee, spanstats, stats)); 229 if (honee->app_ctx->test_type == TESTTYPE_NONE) { 230 PetscCall(PetscViewerPushFormat(ctx->viewer, ctx->format)); 231 PetscCall(VecView(stats, ctx->viewer)); 232 PetscCall(PetscViewerPopFormat(ctx->viewer)); 233 { 234 PetscInt tab_level; 235 PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 236 CeedScalar second = honee->units->second; 237 const char *filename; 238 239 PetscCall(PetscViewerFileGetName(ctx->viewer, &filename)); 240 PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tab_level)); 241 PetscCall(PetscViewerASCIIAddTab(viewer, tab_level + 1)); 242 if (filename) { 243 PetscCall(PetscViewerASCIIPrintf(viewer, 244 "Spanwise turbulent statistics written to %s for time span (%0.12e,%0.12e] and step span [%" PetscInt_FMT 245 ",%" PetscInt_FMT "]\n", 246 filename, spanstats->initial_solution_time / second, solution_time / second, 247 spanstats->initial_solution_step, steps)); 248 } else { 249 PetscCall(PetscViewerASCIIPrintf( 250 viewer, "Spanwise statistics (%s) file written for time span (%0.12e,%0.12e] and step span [%" PetscInt_FMT ",%" PetscInt_FMT "]\n", 251 spanstats->prefix, spanstats->initial_solution_time / second, solution_time / second, spanstats->initial_solution_step, steps)); 252 } 253 PetscCall(PetscViewerASCIISubtractTab(viewer, tab_level + 1)); 254 } 255 } 256 if (honee->app_ctx->test_type == TESTTYPE_TURB_SPANSTATS && reason != TS_CONVERGED_ITERATING) { 257 PetscCall(RegressionTest(honee->app_ctx, stats)); 258 } 259 if (spanstats->do_mms_test && reason != TS_CONVERGED_ITERATING) { 260 Vec error; 261 PetscCall(VecDuplicate(stats, &error)); 262 PetscCall(ApplyCeedOperatorGlobalToGlobal(stats, error, spanstats->mms_error_ctx)); 263 PetscScalar error_sq = 0; 264 PetscCall(VecSum(error, &error_sq)); 265 PetscScalar l2_error = sqrt(error_sq); 266 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error: %.5e\n", l2_error)); 267 } 268 PetscCall(DMRestoreGlobalVector(spanstats->dm, &stats)); 269 } 270 } 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273