xref: /honee/src/spanstats/turbulence.c (revision 00359db47665a79ecb0241f6ccbf886b649022df) 
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
CreateStatisticCollectionOperator(Honee honee,SpanStatsCtx spanstats,SpanStatsSetupData stats_data,ProblemData problem)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, x_coord;
27   CeedElemRestriction         elem_restr_qd, elem_restr_q, elem_restr_x;
28   CeedBasis                   basis_q, basis_x;
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(DMPlexCeedCoordinateCreateField(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, &elem_restr_x, &basis_x, &x_coord));
53   PetscCall(DMPlexCeedElemRestrictionCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &elem_restr_q));
54   PetscCall(DMPlexCeedBasisCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &basis_q));
55   PetscCall(QDataGet(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &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", elem_restr_x, basis_x, 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, CeedVectorDestroy(&x_coord));
99   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
100   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
101   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x));
102   PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
103   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
104   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stats_collect));
105   PetscCallCeed(ceed, CeedOperatorDestroy(&op_stats_collect));
106   PetscFunctionReturn(PETSC_SUCCESS);
107 }
108 
109 // @brief Creates operator for calculating error of method of manufactured solution (MMS) test
SetupMMSErrorChecking(Honee honee,SpanStatsCtx spanstats,SpanStatsSetupData stats_data)110 static PetscErrorCode SetupMMSErrorChecking(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_data) {
111   Ceed                ceed           = honee->ceed;
112   CeedInt             num_comp_stats = spanstats->num_comp_stats, num_comp_x, q_data_size;
113   CeedQFunction       qf_error;
114   CeedOperator        op_error;
115   CeedVector          x_ceed, y_ceed;
116   CeedVector          q_data;
117   CeedElemRestriction elem_restr_parent_qd;
118 
119   PetscFunctionBeginUser;
120   PetscCall(QDataGet(ceed, spanstats->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_parent_qd, &q_data, &q_data_size));
121   PetscCallCeed(ceed, CeedBasisGetNumComponents(stats_data->basis_x, &num_comp_x));
122 
123   PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollectionMMSTest_Error, ChildStatsCollectionMMSTest_Error_loc, &qf_error));
124   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_error, "q", num_comp_stats, CEED_EVAL_INTERP));
125   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_error, "qdata", q_data_size, CEED_EVAL_NONE));
126   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_error, "x", num_comp_x, CEED_EVAL_INTERP));
127   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_error, "v", num_comp_stats, CEED_EVAL_INTERP));
128 
129   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_error, NULL, NULL, &op_error));
130   PetscCallCeed(ceed, CeedOperatorSetField(op_error, "q", stats_data->elem_restr_parent_stats, stats_data->basis_stats, CEED_VECTOR_ACTIVE));
131   PetscCallCeed(ceed, CeedOperatorSetField(op_error, "qdata", elem_restr_parent_qd, CEED_BASIS_NONE, q_data));
132   PetscCallCeed(ceed, CeedOperatorSetField(op_error, "x", stats_data->elem_restr_parent_x, stats_data->basis_x, stats_data->x_coord));
133   PetscCallCeed(ceed, CeedOperatorSetField(op_error, "v", stats_data->elem_restr_parent_stats, stats_data->basis_stats, CEED_VECTOR_ACTIVE));
134 
135   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(stats_data->elem_restr_parent_stats, &x_ceed, NULL));
136   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(stats_data->elem_restr_parent_stats, &y_ceed, NULL));
137   PetscCall(OperatorApplyContextCreate(spanstats->dm, spanstats->dm, honee->ceed, op_error, x_ceed, y_ceed, NULL, NULL, &spanstats->mms_error_ctx));
138 
139   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
140   PetscCallCeed(ceed, CeedVectorDestroy(&x_ceed));
141   PetscCallCeed(ceed, CeedVectorDestroy(&y_ceed));
142   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_parent_qd));
143   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_error));
144   PetscCallCeed(ceed, CeedOperatorDestroy(&op_error));
145   PetscFunctionReturn(PETSC_SUCCESS);
146 }
147 
148 // @brief Setup for statistics collection for turbulence quantities
SpanwiseStatisticsSetup_Turbulence(TS ts,PetscViewerAndFormat * ctx)149 PetscErrorCode SpanwiseStatisticsSetup_Turbulence(TS ts, PetscViewerAndFormat *ctx) {
150   Honee              honee;
151   SpanStatsSetupData stats_setup_data;
152   SpanStatsCtx       spanstats;
153   const char        *prefix = "ts_monitor_spanstats_turbulence_";
154   PetscBool          is_ascii;
155 
156   PetscFunctionBeginUser;
157   PetscCall(TSGetApplicationContext(ts, &honee));
158   PetscCall(PetscObjectTypeCompare((PetscObject)ctx->viewer, PETSCVIEWERASCII, &is_ascii));
159   PetscCheck(!is_ascii || honee->app_ctx->test_type != TESTTYPE_NONE, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP,
160              "Turbulence spanwise statistics does not support ASCII viewers");
161 
162   PetscCall(SpanwiseStatisticsSetupInitialize(honee, honee->app_ctx->degree, TURB_NUM_COMPONENTS, prefix, &stats_setup_data, &spanstats));
163 
164   {  // Create Section for data
165     PetscSection section;
166 
167     PetscCall(DMGetLocalSection(spanstats->dm, &section));
168     PetscCall(PetscSectionSetFieldName(section, 0, ""));
169     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY, "MeanDensity"));
170     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE, "MeanPressure"));
171     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_SQUARED, "MeanPressureSquared"));
172     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_VELOCITY_X, "MeanPressureVelocityX"));
173     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_VELOCITY_Y, "MeanPressureVelocityY"));
174     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_PRESSURE_VELOCITY_Z, "MeanPressureVelocityZ"));
175     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE, "MeanDensityTemperature"));
176     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE_FLUX_X, "MeanDensityTemperatureFluxX"));
177     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE_FLUX_Y, "MeanDensityTemperatureFluxY"));
178     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_DENSITY_TEMPERATURE_FLUX_Z, "MeanDensityTemperatureFluxZ"));
179     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUM_X, "MeanMomentumX"));
180     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUM_Y, "MeanMomentumY"));
181     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUM_Z, "MeanMomentumZ"));
182     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_XX, "MeanMomentumFluxXX"));
183     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_YY, "MeanMomentumFluxYY"));
184     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_ZZ, "MeanMomentumFluxZZ"));
185     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_YZ, "MeanMomentumFluxYZ"));
186     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_XZ, "MeanMomentumFluxXZ"));
187     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_MOMENTUMFLUX_XY, "MeanMomentumFluxXY"));
188     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_VELOCITY_X, "MeanVelocityX"));
189     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_VELOCITY_Y, "MeanVelocityY"));
190     PetscCall(PetscSectionSetComponentName(section, 0, TURB_MEAN_VELOCITY_Z, "MeanVelocityZ"));
191   }
192 
193   // Create CeedOperators for statistics collection
194   PetscCall(CreateStatisticCollectionOperator(honee, spanstats, stats_setup_data, honee->problem_data));
195 
196   PetscCall(PetscOptionsGetBool(NULL, NULL, "-ts_monitor_spanstats_turbulence_mms", &spanstats->do_mms_test, NULL));
197   if (spanstats->do_mms_test) {
198     PetscCall(SetupMMSErrorChecking(honee, spanstats, stats_setup_data));
199   }
200 
201   ctx->data         = spanstats;
202   ctx->data_destroy = (PetscCtxDestroyFn *)SpanStatsCtxDestroy;
203 
204   PetscCall(SpanwiseStatisticsSetupFinalize(ts, honee, spanstats, ctx, &stats_setup_data));
205   PetscFunctionReturn(PETSC_SUCCESS);
206 }
207 
208 // @brief TSMonitor for collection and processing of turbulence statistics
TSMonitor_SpanwiseStatisticsTurbulence(TS ts,PetscInt steps,PetscReal solution_time,Vec Q,PetscViewerAndFormat * ctx)209 PetscErrorCode TSMonitor_SpanwiseStatisticsTurbulence(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, PetscViewerAndFormat *ctx) {
210   Honee             honee;
211   Vec               stats;
212   TSConvergedReason reason;
213   SpanStatsCtx      spanstats        = ctx->data;
214   PetscInt          collect_interval = spanstats->collect_interval, viewer_interval = ctx->view_interval;
215 
216   PetscFunctionBeginUser;
217   PetscCall(TSGetApplicationContext(ts, &honee));
218   PetscCall(TSGetConvergedReason(ts, &reason));
219   // Do not collect or process on the first step of the run (ie. on the initial condition)
220   if (steps == honee->app_ctx->cont_steps) PetscFunctionReturn(PETSC_SUCCESS);
221 
222   PetscBool run_processing_and_viewer = (steps % viewer_interval == 0 && viewer_interval != -1) || reason != TS_CONVERGED_ITERATING;
223 
224   if (steps % collect_interval == 0 || run_processing_and_viewer) {
225     PetscCall(SpanwiseStatisticsCollect(honee, spanstats, solution_time, Q));
226 
227     if (run_processing_and_viewer) {
228       PetscCall(DMSetOutputSequenceNumber(spanstats->dm, steps, solution_time));
229       PetscCall(DMGetGlobalVector(spanstats->dm, &stats));
230       PetscCall(SpanwiseStatisticsProcess(honee, spanstats, stats));
231       if (honee->app_ctx->test_type == TESTTYPE_NONE) {
232         PetscCall(PetscViewerPushFormat(ctx->viewer, ctx->format));
233         PetscCall(VecView(stats, ctx->viewer));
234         PetscCall(PetscViewerPopFormat(ctx->viewer));
235         {
236           PetscInt    tab_level;
237           PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts));
238           CeedScalar  second = honee->units->second;
239           const char *filename;
240 
241           PetscCall(PetscViewerFileGetName(ctx->viewer, &filename));
242           PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tab_level));
243           PetscCall(PetscViewerASCIIAddTab(viewer, tab_level + 1));
244           if (filename) {
245             PetscCall(PetscViewerASCIIPrintf(viewer,
246                                              "Spanwise turbulent statistics written to %s for time span (%0.12e,%0.12e] and step span [%" PetscInt_FMT
247                                              ",%" PetscInt_FMT "]\n",
248                                              filename, spanstats->initial_solution_time / second, solution_time / second,
249                                              spanstats->initial_solution_step, steps));
250           } else {
251             PetscCall(PetscViewerASCIIPrintf(
252                 viewer, "Spanwise statistics (%s) file written for time span (%0.12e,%0.12e] and step span [%" PetscInt_FMT ",%" PetscInt_FMT "]\n",
253                 spanstats->prefix, spanstats->initial_solution_time / second, solution_time / second, spanstats->initial_solution_step, steps));
254           }
255           PetscCall(PetscViewerASCIISubtractTab(viewer, tab_level + 1));
256         }
257       }
258       if (honee->app_ctx->test_type == TESTTYPE_TURB_SPANSTATS && reason != TS_CONVERGED_ITERATING) {
259         PetscCall(RegressionTest(honee->app_ctx, stats));
260       }
261       if (spanstats->do_mms_test && reason != TS_CONVERGED_ITERATING) {
262         Vec error;
263         PetscCall(VecDuplicate(stats, &error));
264         PetscCall(ApplyCeedOperatorGlobalToGlobal(stats, error, spanstats->mms_error_ctx));
265         PetscScalar error_sq = 0;
266         PetscCall(VecSum(error, &error_sq));
267         PetscScalar l2_error = sqrt(error_sq);
268         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error: %.5e\n", l2_error));
269       }
270       PetscCall(DMRestoreGlobalVector(spanstats->dm, &stats));
271     }
272   }
273   PetscFunctionReturn(PETSC_SUCCESS);
274 }
275