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
CreateStatisticCollectionOperator(Honee honee,SpanStatsCtx spanstats,SpanStatsSetupData stats_data,ProblemData problem)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
SpanwiseStatisticsSetup_CflPe(TS ts,PetscViewerAndFormat * ctx)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
TSMonitor_SpanwiseStatisticsCflPe(TS ts,PetscInt steps,PetscReal solution_time,Vec Q,PetscViewerAndFormat * ctx)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