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 CFL and Peclet number 5 6 #include "../qfunctions/spanstats/cflpe.h" 7 #include "../qfunctions/advection_types.h" 8 9 #include <ceed.h> 10 #include <petscdmplex.h> 11 #include <petscsf.h> 12 #include <spanstats.h> 13 14 #include <navierstokes.h> 15 #include <petsc_ops.h> 16 17 // @brief Create CeedOperator for statistics collection 18 static PetscErrorCode CreateStatisticCollectionOperator(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_data, ProblemData problem) { 19 Ceed ceed = honee->ceed; 20 CeedInt num_comp_stats = spanstats->num_comp_stats, q_data_size; 21 PetscInt dim, num_comp_q, num_comp_x; 22 CflPe_SpanStatsContext collect_ctx; 23 CeedQFunctionContext collect_qfctx; 24 CeedQFunction qf_stats_collect = NULL; 25 CeedOperator op_stats_collect; 26 CeedVector q_data, x_coord; 27 CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x; 28 CeedBasis basis_q, basis_x; 29 30 PetscFunctionBeginUser; 31 PetscCall(DMGetDimension(honee->dm, &dim)); 32 PetscCall(DMGetCoordinateNumComps(honee->dm, &num_comp_x)); 33 PetscCall(DMGetFieldNumComps(honee->dm, 0, &num_comp_q)); 34 35 // Create Operator for statistics collection 36 switch (dim) { 37 case 2: 38 switch (honee->phys->state_var) { 39 case STATEVAR_PRIMITIVE: 40 PetscCallCeed(ceed, 41 CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_2D_Prim, ChildStatsCollection_2D_Prim_loc, &qf_stats_collect)); 42 break; 43 case STATEVAR_CONSERVATIVE: 44 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_2D_Conserv, ChildStatsCollection_2D_Conserv_loc, 45 &qf_stats_collect)); 46 break; 47 case STATEVAR_ENTROPY: 48 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_2D_Entropy, ChildStatsCollection_2D_Entropy_loc, 49 &qf_stats_collect)); 50 break; 51 } 52 break; 53 case 3: 54 switch (honee->phys->state_var) { 55 case STATEVAR_PRIMITIVE: 56 PetscCallCeed(ceed, 57 CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_3D_Prim, ChildStatsCollection_3D_Prim_loc, &qf_stats_collect)); 58 break; 59 case STATEVAR_CONSERVATIVE: 60 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_3D_Conserv, ChildStatsCollection_3D_Conserv_loc, 61 &qf_stats_collect)); 62 break; 63 case STATEVAR_ENTROPY: 64 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_3D_Entropy, ChildStatsCollection_3D_Entropy_loc, 65 &qf_stats_collect)); 66 break; 67 } 68 break; 69 } 70 PetscCheck(qf_stats_collect, PetscObjectComm((PetscObject)honee->dm), PETSC_ERR_SUP, 71 "Could not create CFL/Pe spanstats collection QFunction for dim %" PetscInt_FMT " and state variable %s", dim, 72 StateVariables[honee->phys->state_var]); 73 74 PetscCall(DMPlexCeedCoordinateCreateField(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, &elem_restr_x, &basis_x, &x_coord)); 75 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &elem_restr_q)); 76 PetscCall(DMPlexCeedBasisCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &basis_q)); 77 PetscCall(QDataGet(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size)); 78 { // Setup Collection Context 79 PetscBool is_advection; 80 81 PetscCall(PetscNew(&collect_ctx)); 82 PetscCall(PetscStrcmp("advection", honee->app_ctx->problem_name, &is_advection)); 83 if (is_advection) { 84 AdvectionContext advection_ctx; 85 86 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &advection_ctx)); 87 collect_ctx->newt_ctx = (struct NewtonianIdealGasContext_){0}; 88 collect_ctx->diffusion_coeff = advection_ctx->diffusion_coeff; 89 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &advection_ctx)); 90 } else { 91 NewtonianIdealGasContext newtonian_ig_ctx; 92 93 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 94 collect_ctx->newt_ctx = *newtonian_ig_ctx; 95 collect_ctx->diffusion_coeff = newtonian_ig_ctx->gas.mu; 96 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 97 } 98 99 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &collect_qfctx)); 100 PetscCallCeed(ceed, CeedQFunctionContextSetData(collect_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*collect_ctx), collect_ctx)); 101 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(collect_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 102 103 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(collect_qfctx, "solution time", offsetof(struct CflPe_SpanStatsContext_, solution_time), 1, 104 "Current solution time")); 105 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(collect_qfctx, "previous time", offsetof(struct CflPe_SpanStatsContext_, previous_time), 1, 106 "Previous time statistics collection was done")); 107 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(collect_qfctx, "timestep", offsetof(struct CflPe_SpanStatsContext_, timestep), 1, 108 "Current timestep taken by TS")); 109 } 110 111 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_stats_collect, collect_qfctx)); 112 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&collect_qfctx)); 113 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_stats_collect, "q", num_comp_q, CEED_EVAL_INTERP)); 114 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_stats_collect, "q_data", q_data_size, CEED_EVAL_NONE)); 115 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_stats_collect, "x", num_comp_x, CEED_EVAL_INTERP)); 116 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_stats_collect, "v", num_comp_stats, CEED_EVAL_NONE)); 117 118 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stats_collect, NULL, NULL, &op_stats_collect)); 119 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 120 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "q_data", elem_restr_qd, CEED_BASIS_NONE, q_data)); 121 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "x", elem_restr_x, basis_x, x_coord)); 122 PetscCallCeed(ceed, CeedOperatorSetField(op_stats_collect, "v", stats_data->elem_restr_child_colloc, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 123 124 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_stats_collect, "solution time", &spanstats->solution_time_label)); 125 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_stats_collect, "previous time", &spanstats->previous_time_label)); 126 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_stats_collect, "timestep", &spanstats->timestep_label)); 127 128 PetscCall(OperatorApplyContextCreate(honee->dm, spanstats->dm, honee->ceed, op_stats_collect, honee->q_ceed, NULL, NULL, NULL, 129 &spanstats->op_stats_collect_ctx)); 130 131 PetscCall(CeedOperatorCreateLocalVecs(op_stats_collect, DMReturnVecType(spanstats->dm), PETSC_COMM_SELF, NULL, &spanstats->Child_Stats_loc)); 132 PetscCall(VecZeroEntries(spanstats->Child_Stats_loc)); 133 134 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 135 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 136 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stats_collect)); 137 PetscCallCeed(ceed, CeedOperatorDestroy(&op_stats_collect)); 138 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 139 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 140 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x)); 141 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x)); 142 PetscCallCeed(ceed, CeedVectorDestroy(&x_coord)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 // @brief Setup for statistics collection of CFL and Peclet number 147 PetscErrorCode SpanwiseStatisticsSetup_CflPe(TS ts, PetscViewerAndFormat *ctx) { 148 Honee honee; 149 SpanStatsSetupData stats_setup_data; 150 SpanStatsCtx spanstats; 151 const char *prefix = "ts_monitor_spanstats_cflpe_"; 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, 6, 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, 0, "MeanCFL")); 168 PetscCall(PetscSectionSetComponentName(section, 0, 1, "MeanCFLSquared")); 169 PetscCall(PetscSectionSetComponentName(section, 0, 2, "MeanCFLCubed")); 170 PetscCall(PetscSectionSetComponentName(section, 0, 3, "MeanPe")); 171 PetscCall(PetscSectionSetComponentName(section, 0, 4, "MeanPeSquared")); 172 PetscCall(PetscSectionSetComponentName(section, 0, 5, "MeanPeCubed")); 173 } 174 175 // Create CeedOperators for statistics collection 176 PetscCall(CreateStatisticCollectionOperator(honee, spanstats, stats_setup_data, honee->problem_data)); 177 178 ctx->data = spanstats; 179 ctx->data_destroy = (PetscCtxDestroyFn *)SpanStatsCtxDestroy; 180 181 PetscCall(SpanwiseStatisticsSetupFinalize(ts, honee, spanstats, ctx, &stats_setup_data)); 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184 185 // @brief TSMonitor for collection and processing of CFL and Peclet number statistics 186 PetscErrorCode TSMonitor_SpanwiseStatisticsCflPe(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, PetscViewerAndFormat *ctx) { 187 Honee honee; 188 Vec stats; 189 TSConvergedReason reason; 190 SpanStatsCtx spanstats = ctx->data; 191 PetscInt collect_interval = spanstats->collect_interval, viewer_interval = ctx->view_interval; 192 PetscReal timestep; 193 194 PetscFunctionBeginUser; 195 PetscCall(TSGetApplicationContext(ts, &honee)); 196 PetscCall(TSGetConvergedReason(ts, &reason)); 197 // Do not collect or process on the first step of the run (ie. on the initial condition) 198 if (steps == honee->app_ctx->cont_steps) PetscFunctionReturn(PETSC_SUCCESS); 199 200 PetscBool run_processing_and_viewer = (steps % viewer_interval == 0 && viewer_interval != -1) || reason != TS_CONVERGED_ITERATING; 201 202 if (steps % collect_interval == 0 || run_processing_and_viewer) { 203 PetscCall(TSGetTimeStep(ts, ×tep)); 204 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->timestep_label, ×tep)); 205 PetscCall(SpanwiseStatisticsCollect(honee, spanstats, solution_time, Q)); 206 207 if (run_processing_and_viewer) { 208 PetscCall(DMSetOutputSequenceNumber(spanstats->dm, steps, solution_time)); 209 PetscCall(DMGetGlobalVector(spanstats->dm, &stats)); 210 PetscCall(SpanwiseStatisticsProcess(honee, spanstats, stats)); 211 212 if (honee->app_ctx->test_type == TESTTYPE_NONE) { 213 PetscCall(PetscViewerPushFormat(ctx->viewer, ctx->format)); 214 PetscCall(VecView(stats, ctx->viewer)); 215 PetscCall(PetscViewerPopFormat(ctx->viewer)); 216 { 217 PetscInt tab_level; 218 PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 219 CeedScalar second = honee->units->second; 220 const char *filename; 221 222 PetscCall(PetscViewerFileGetName(ctx->viewer, &filename)); 223 PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tab_level)); 224 PetscCall(PetscViewerASCIIAddTab(viewer, tab_level + 1)); 225 if (filename) { 226 PetscCall(PetscViewerASCIIPrintf( 227 viewer, "Spanwise cflpe statistics written to %s for time span (%0.12e,%0.12e] and step span [%" PetscInt_FMT ",%" PetscInt_FMT "]\n", 228 filename, spanstats->initial_solution_time / second, solution_time / second, spanstats->initial_solution_step, steps)); 229 } else { 230 PetscCall(PetscViewerASCIIPrintf( 231 viewer, "Spanwise statistics (%s) file written for time span (%0.12e,%0.12e] and step span [%" PetscInt_FMT ",%" PetscInt_FMT "]\n", 232 spanstats->prefix, spanstats->initial_solution_time / second, solution_time / second, spanstats->initial_solution_step, steps)); 233 } 234 PetscCall(PetscViewerASCIISubtractTab(viewer, tab_level + 1)); 235 } 236 } else if (honee->app_ctx->test_type == TESTTYPE_SPANSTATS && reason != TS_CONVERGED_ITERATING) { 237 PetscCall(RegressionTest(honee->app_ctx, stats)); 238 } 239 240 PetscCall(DMRestoreGlobalVector(spanstats->dm, &stats)); 241 } 242 } 243 PetscFunctionReturn(PETSC_SUCCESS); 244 } 245