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
5 ///
6 /// "Parent" refers to the 2D plane on which statistics are collected *onto*.
7 /// "Child" refers to the 3D domain where statistics are gathered *from*.
8 /// Each quadrature point on the parent plane has several children in the child domain that it performs spanwise averaging with.
9
10 #include <ceed.h>
11 #include <petscdmplex.h>
12 #include <petscsf.h>
13 #include <spanstats.h>
14
15 #include <navierstokes.h>
16 #include <petsc_ops.h>
17
SpanStatsCtxDestroy(SpanStatsCtx * spanstats)18 PetscErrorCode SpanStatsCtxDestroy(SpanStatsCtx *spanstats) {
19 SpanStatsCtx spanstats_ = *spanstats;
20
21 PetscFunctionBeginUser;
22 if (spanstats_ == NULL) PetscFunctionReturn(PETSC_SUCCESS);
23 PetscCall(VecDestroy(&spanstats_->Child_Stats_loc));
24 PetscCall(VecDestroy(&spanstats_->Parent_Stats_loc));
25
26 PetscCall(OperatorApplyContextDestroy(spanstats_->op_stats_collect_ctx));
27 PetscCall(OperatorApplyContextDestroy(spanstats_->op_proj_rhs_ctx));
28 PetscCall(OperatorApplyContextDestroy(spanstats_->mms_error_ctx));
29
30 PetscCall(KSPDestroy(&spanstats_->ksp));
31 PetscCall(PetscSFDestroy(&spanstats_->sf));
32 PetscCall(DMDestroy(&spanstats_->dm));
33 PetscCall(PetscFree(spanstats_->prefix));
34 PetscCall(PetscFree(spanstats_));
35 *spanstats = NULL;
36 PetscFunctionReturn(PETSC_SUCCESS);
37 }
38
SpanwiseStatisticssCreateDM(Honee honee,SpanStatsCtx spanstats,PetscInt degree,PetscInt num_comps)39 static PetscErrorCode SpanwiseStatisticssCreateDM(Honee honee, SpanStatsCtx spanstats, PetscInt degree, PetscInt num_comps) {
40 PetscReal domain_min[3], domain_max[3];
41 PetscLogStage stage_stats_setup;
42 MPI_Comm comm = PetscObjectComm((PetscObject)honee->dm);
43
44 PetscFunctionBeginUser;
45 PetscCall(PetscLogStageGetId("Stats Setup", &stage_stats_setup));
46 if (stage_stats_setup == -1) PetscCall(PetscLogStageRegister("Stats Setup", &stage_stats_setup));
47 PetscCall(PetscLogStagePush(stage_stats_setup));
48
49 spanstats->num_comp_stats = num_comps;
50
51 // Get spanwise length
52 PetscCall(DMGetBoundingBox(honee->dm, domain_min, domain_max));
53 spanstats->span_width = domain_max[2] - domain_min[2];
54
55 { // Get DM from surface
56 DM parent_distributed_dm;
57 const PetscSF *isoperiodicface;
58 PetscInt num_isoperiodicface;
59 DMLabel label;
60 PetscMPIInt size;
61
62 PetscCall(DMPlexGetIsoperiodicFaceSF(honee->dm, &num_isoperiodicface, &isoperiodicface));
63
64 if (isoperiodicface) {
65 PetscSF inv_isoperiodicface;
66 PetscInt nleaves, isoperiodicface_index = -1;
67 const PetscInt *ilocal;
68 char isoperiodicface_key[1024];
69
70 PetscCall(PetscSNPrintf(isoperiodicface_key, sizeof isoperiodicface_key, "-%sisoperiodic_facesf", spanstats->prefix));
71 PetscCall(PetscOptionsGetInt(NULL, NULL, isoperiodicface_key, &isoperiodicface_index, NULL));
72 isoperiodicface_index = isoperiodicface_index == -1 ? num_isoperiodicface - 1 : isoperiodicface_index;
73 PetscCall(PetscSFCreateInverseSF(isoperiodicface[isoperiodicface_index], &inv_isoperiodicface));
74 PetscCall(PetscSFGetGraph(inv_isoperiodicface, NULL, &nleaves, &ilocal, NULL));
75 PetscCall(DMCreateLabel(honee->dm, "Periodic Face"));
76 PetscCall(DMGetLabel(honee->dm, "Periodic Face", &label));
77 for (PetscInt i = 0; i < nleaves; i++) {
78 PetscCall(DMLabelSetValue(label, ilocal[i], 1));
79 }
80 PetscCall(PetscSFDestroy(&inv_isoperiodicface));
81 } else {
82 PetscCall(DMGetLabel(honee->dm, "Face Sets", &label));
83 }
84
85 PetscCall(DMPlexLabelComplete(honee->dm, label));
86 PetscCall(DMPlexFilter(honee->dm, label, 1, PETSC_FALSE, PETSC_FALSE, PetscObjectComm((PetscObject)honee->dm), NULL, &spanstats->dm));
87 PetscCall(DMSetCoordinateDisc(spanstats->dm, NULL, PETSC_FALSE, PETSC_TRUE)); // Ensure that a coordinate FE exists
88
89 PetscCall(DMPlexDistribute(spanstats->dm, 0, NULL, &parent_distributed_dm));
90 PetscCallMPI(MPI_Comm_size(comm, &size));
91 if (parent_distributed_dm) {
92 PetscCall(DMDestroy(&spanstats->dm));
93 spanstats->dm = parent_distributed_dm;
94 } else if (size > 1) {
95 PetscCall(PetscPrintf(comm, "WARNING: Spanwise statistics: parent DM could not be distributed accross %d ranks.\n", size));
96 }
97 }
98 {
99 PetscBool is_simplex = PETSC_FALSE;
100 PetscCall(DMPlexIsSimplex(spanstats->dm, &is_simplex));
101 PetscCheck(is_simplex != PETSC_TRUE, comm, PETSC_ERR_ARG_WRONGSTATE, "Spanwise statistics is not implemented for non-tensor product grids");
102 }
103
104 PetscCall(PetscObjectSetName((PetscObject)spanstats->dm, "Spanwise_Stats"));
105 PetscCall(DMSetOptionsPrefix(spanstats->dm, spanstats->prefix));
106 PetscCall(DMSetFromOptions(spanstats->dm));
107 PetscCall(DMViewFromOptions(spanstats->dm, NULL, "-dm_view"));
108
109 // Create FE space for parent DM
110 PetscCall(DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, honee->app_ctx->degree, 1, honee->app_ctx->q_extra, 1, &spanstats->num_comp_stats,
111 spanstats->dm));
112
113 PetscCall(PetscLogStagePop());
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
117 /** @brief Create CeedElemRestriction for collocated data in component-major order.
118 a. Sets the strides of the restriction to component-major order
119 Number of quadrature points is used from the CeedBasis, and number of elements is used from the CeedElemRestriction.
120 */
CreateElemRestrColloc_CompMajor(Ceed ceed,CeedInt num_comp,CeedBasis basis,CeedElemRestriction elem_restr_base,CeedElemRestriction * elem_restr_collocated)121 static PetscErrorCode CreateElemRestrColloc_CompMajor(Ceed ceed, CeedInt num_comp, CeedBasis basis, CeedElemRestriction elem_restr_base,
122 CeedElemRestriction *elem_restr_collocated) {
123 CeedInt num_elem_qpts, loc_num_elem;
124
125 PetscFunctionBeginUser;
126 PetscCallCeed(ceed, CeedBasisGetNumQuadraturePoints(basis, &num_elem_qpts));
127 PetscCallCeed(ceed, CeedElemRestrictionGetNumElements(elem_restr_base, &loc_num_elem));
128
129 const CeedInt strides[] = {num_comp, 1, num_elem_qpts * num_comp};
130 PetscCallCeed(ceed, CeedElemRestrictionCreateStrided(ceed, loc_num_elem, num_elem_qpts, num_comp, num_comp * loc_num_elem * num_elem_qpts, strides,
131 elem_restr_collocated));
132 PetscFunctionReturn(PETSC_SUCCESS);
133 }
134
135 // Get coordinates of quadrature points
GetQuadratureCoords(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,Vec * Qx_coords,PetscInt * total_nqpnts)136 static PetscErrorCode GetQuadratureCoords(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, Vec *Qx_coords,
137 PetscInt *total_nqpnts) {
138 CeedElemRestriction elem_restr_qx, elem_restr_x;
139 CeedBasis basis_x;
140 CeedVector x_coords;
141 CeedQFunction qf_quad_coords;
142 CeedOperator op_quad_coords;
143 CeedInt num_comp_x;
144 CeedSize l_vec_size;
145 OperatorApplyContext op_quad_coords_ctx;
146
147 PetscFunctionBeginUser;
148 PetscCall(DMPlexCeedCoordinateCreateField(ceed, dm, domain_label, label_value, height, &elem_restr_x, &basis_x, &x_coords));
149 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
150 PetscCall(CreateElemRestrColloc_CompMajor(ceed, num_comp_x, basis_x, elem_restr_x, &elem_restr_qx));
151 PetscCallCeed(ceed, CeedElemRestrictionGetLVectorSize(elem_restr_qx, &l_vec_size));
152 *total_nqpnts = l_vec_size / num_comp_x;
153
154 // Create QFunction
155 PetscCallCeed(ceed, CeedQFunctionCreateIdentity(ceed, num_comp_x, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_quad_coords));
156
157 // Create Operator
158 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_quad_coords, NULL, NULL, &op_quad_coords));
159 PetscCallCeed(ceed, CeedOperatorSetField(op_quad_coords, "input", elem_restr_x, basis_x, x_coords));
160 PetscCallCeed(ceed, CeedOperatorSetField(op_quad_coords, "output", elem_restr_qx, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
161
162 PetscCall(CeedOperatorCreateLocalVecs(op_quad_coords, DMReturnVecType(dm), PETSC_COMM_SELF, NULL, Qx_coords));
163 PetscCall(OperatorApplyContextCreate(NULL, NULL, ceed, op_quad_coords, CEED_VECTOR_NONE, NULL, NULL, NULL, &op_quad_coords_ctx));
164
165 PetscCall(ApplyCeedOperatorLocalToLocal(NULL, *Qx_coords, op_quad_coords_ctx));
166
167 PetscCall(OperatorApplyContextDestroy(op_quad_coords_ctx));
168 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qx));
169 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x));
170 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
171 PetscCallCeed(ceed, CeedVectorDestroy(&x_coords));
172 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_quad_coords));
173 PetscCallCeed(ceed, CeedOperatorDestroy(&op_quad_coords));
174 PetscFunctionReturn(PETSC_SUCCESS);
175 }
176
SpanwiseStatisticsSetupDataCreate(Honee honee,SpanStatsCtx spanstats,SpanStatsSetupData * stats_setup_data)177 static PetscErrorCode SpanwiseStatisticsSetupDataCreate(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData *stats_setup_data) {
178 Ceed ceed = honee->ceed;
179 DM dm = spanstats->dm;
180 Vec X_loc;
181 CeedInt num_comp_x, num_comp_stats = spanstats->num_comp_stats;
182 PetscInt height = 0, dm_field = 0;
183 CeedElemRestriction elem_restr_q;
184 CeedBasis basis_q;
185
186 PetscFunctionBeginUser;
187 PetscCall(PetscNew(stats_setup_data));
188
189 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field,
190 &(*stats_setup_data)->elem_restr_parent_stats));
191 PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height,
192 &(*stats_setup_data)->elem_restr_parent_x));
193 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents((*stats_setup_data)->elem_restr_parent_x, &num_comp_x));
194 PetscCallCeed(ceed, CeedElemRestrictionCreateVector((*stats_setup_data)->elem_restr_parent_x, &(*stats_setup_data)->x_coord, NULL));
195 PetscCall(DMPlexCeedBasisCoordinateCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, &(*stats_setup_data)->basis_x));
196 PetscCall(DMPlexCeedBasisCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field, &(*stats_setup_data)->basis_stats));
197
198 PetscCall(CreateElemRestrColloc_CompMajor(ceed, num_comp_stats, (*stats_setup_data)->basis_stats, (*stats_setup_data)->elem_restr_parent_stats,
199 &(*stats_setup_data)->elem_restr_parent_colloc));
200 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &elem_restr_q));
201 PetscCall(DMPlexCeedBasisCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &basis_q));
202 PetscCall(CreateElemRestrColloc_CompMajor(ceed, num_comp_stats, basis_q, elem_restr_q, &(*stats_setup_data)->elem_restr_child_colloc));
203 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
204 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
205
206 { // -- Copy DM coordinates into CeedVector
207 DM cdm;
208 PetscCall(DMGetCellCoordinateDM(dm, &cdm));
209 if (cdm) {
210 PetscCall(DMGetCellCoordinatesLocal(dm, &X_loc));
211 } else {
212 PetscCall(DMGetCoordinatesLocal(dm, &X_loc));
213 }
214 }
215 // Do not need to scale the coordinates. The DMClone operation has already copied the scaled coordinates
216 PetscCall(VecCopyPetscToCeed(X_loc, (*stats_setup_data)->x_coord));
217 PetscFunctionReturn(PETSC_SUCCESS);
218 }
219
SpanwiseStatisticsSetupDataDestroy(SpanStatsSetupData data)220 static PetscErrorCode SpanwiseStatisticsSetupDataDestroy(SpanStatsSetupData data) {
221 Ceed ceed;
222
223 PetscFunctionBeginUser;
224 PetscCall(CeedElemRestrictionGetCeed(data->elem_restr_parent_x, &ceed));
225 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_parent_x));
226 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_parent_stats));
227 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_parent_colloc));
228 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_child_colloc));
229 PetscCallCeed(ceed, CeedBasisDestroy(&data->basis_x));
230 PetscCallCeed(ceed, CeedBasisDestroy(&data->basis_stats));
231 PetscCallCeed(ceed, CeedVectorDestroy(&data->x_coord));
232 PetscCall(PetscFree(data));
233 PetscCheck(CeedDestroy(&ceed) == CEED_ERROR_SUCCESS, PETSC_COMM_SELF, PETSC_ERR_LIB, "Destroying Ceed object failed");
234 PetscFunctionReturn(PETSC_SUCCESS);
235 }
236
237 // Create PetscSF for child-to-parent communication
SpanwiseStatisticsCreateSF(Honee honee,SpanStatsCtx spanstats,SpanStatsSetupData stats_setup_data,DM parentdm,DM childdm,PetscSF * statssf)238 static PetscErrorCode SpanwiseStatisticsCreateSF(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_setup_data, DM parentdm, DM childdm,
239 PetscSF *statssf) {
240 Ceed ceed = honee->ceed;
241 PetscInt child_num_qpnts, parent_num_qpnts, num_comp_x;
242 Vec Child_qx_coords, Parent_qx_coords;
243
244 PetscFunctionBeginUser;
245 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)childdm), statssf));
246
247 // Assume that child and parent have the same number of components
248 PetscCall(DMGetCoordinateNumComps(childdm, &num_comp_x));
249 const PetscInt num_comp_sf = num_comp_x - 1; // Number of coord components used in the creation of the SF
250
251 // Get quad_coords for child and parent DM
252 PetscCall(GetQuadratureCoords(ceed, childdm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, &Child_qx_coords, &child_num_qpnts));
253 PetscCall(GetQuadratureCoords(ceed, parentdm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, &Parent_qx_coords, &parent_num_qpnts));
254
255 { // Remove z component of coordinates for matching
256 const PetscReal *child_quad_coords, *parent_quad_coords;
257 PetscReal *child_coords, *parent_coords;
258
259 PetscCall(VecGetArrayRead(Child_qx_coords, &child_quad_coords));
260 PetscCall(VecGetArrayRead(Parent_qx_coords, &parent_quad_coords));
261
262 PetscCall(PetscMalloc2(child_num_qpnts * 2, &child_coords, parent_num_qpnts * 2, &parent_coords));
263 for (int i = 0; i < child_num_qpnts; i++) {
264 child_coords[0 + i * num_comp_sf] = child_quad_coords[0 + i * num_comp_x];
265 child_coords[1 + i * num_comp_sf] = child_quad_coords[1 + i * num_comp_x];
266 }
267 for (int i = 0; i < parent_num_qpnts; i++) {
268 parent_coords[0 + i * num_comp_sf] = parent_quad_coords[0 + i * num_comp_x];
269 parent_coords[1 + i * num_comp_sf] = parent_quad_coords[1 + i * num_comp_x];
270 }
271 PetscCall(VecRestoreArrayRead(Child_qx_coords, &child_quad_coords));
272 PetscCall(VecRestoreArrayRead(Parent_qx_coords, &parent_quad_coords));
273
274 PetscCall(PetscSFSetGraphFromCoordinates(*statssf, parent_num_qpnts, child_num_qpnts, num_comp_sf, 1e-12, parent_coords, child_coords));
275 PetscCall(PetscFree2(child_coords, parent_coords));
276 }
277
278 {
279 char spanstats_sf_key[1024];
280 PetscCall(PetscSNPrintf(spanstats_sf_key, sizeof spanstats_sf_key, "-%ssf_view", spanstats->prefix));
281 PetscCall(PetscObjectSetName((PetscObject)*statssf, spanstats->prefix));
282 PetscCall(PetscSFViewFromOptions(*statssf, NULL, spanstats_sf_key));
283 }
284
285 PetscCall(VecDestroy(&Child_qx_coords));
286 PetscCall(VecDestroy(&Parent_qx_coords));
287 PetscFunctionReturn(PETSC_SUCCESS);
288 }
289
290 // @brief Setup RHS and LHS for L^2 projection of statistics
SpanwiseStatisticsSetupL2Projection(Honee honee,SpanStatsCtx spanstats,SpanStatsSetupData stats_setup_data)291 static PetscErrorCode SpanwiseStatisticsSetupL2Projection(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_setup_data) {
292 Ceed ceed = honee->ceed;
293 CeedOperator op_mass, op_proj_rhs;
294 CeedQFunction qf_mass, qf_stats_proj;
295 CeedInt q_data_size, num_comp_stats = spanstats->num_comp_stats;
296 CeedElemRestriction elem_restr_qd;
297 CeedVector q_data;
298
299 PetscFunctionBeginUser;
300 // -- Create Operator for RHS of L^2 projection of statistics
301 // Simply take collocated parent data (with quadrature weight already applied) and multiply by weight function.
302 // Therefore, an Identity QF is sufficient
303 PetscCallCeed(ceed, CeedQFunctionCreateIdentity(ceed, num_comp_stats, CEED_EVAL_NONE, CEED_EVAL_INTERP, &qf_stats_proj));
304
305 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stats_proj, NULL, NULL, &op_proj_rhs));
306 PetscCallCeed(ceed, CeedOperatorSetField(op_proj_rhs, "input", stats_setup_data->elem_restr_parent_colloc, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
307 PetscCallCeed(ceed, CeedOperatorSetField(op_proj_rhs, "output", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats,
308 CEED_VECTOR_ACTIVE));
309
310 PetscCall(OperatorApplyContextCreate(NULL, spanstats->dm, ceed, op_proj_rhs, NULL, NULL, NULL, NULL, &spanstats->op_proj_rhs_ctx));
311 PetscCall(CeedOperatorCreateLocalVecs(op_proj_rhs, DMReturnVecType(spanstats->dm), PETSC_COMM_SELF, &spanstats->Parent_Stats_loc, NULL));
312 PetscCall(QDataGet(ceed, spanstats->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size));
313
314 // Create Mass CeedOperator
315 PetscCall(HoneeMassQFunctionCreate(ceed, num_comp_stats, q_data_size, &qf_mass));
316 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass));
317 PetscCallCeed(ceed,
318 CeedOperatorSetField(op_mass, "u", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, CEED_VECTOR_ACTIVE));
319 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
320 PetscCallCeed(ceed,
321 CeedOperatorSetField(op_mass, "v", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, CEED_VECTOR_ACTIVE));
322
323 { // Setup KSP for L^2 projection
324 Mat mat_mass;
325 KSP ksp;
326
327 PetscCall(MatCreateCeed(spanstats->dm, spanstats->dm, op_mass, NULL, &mat_mass));
328
329 PetscCall(KSPCreate(PetscObjectComm((PetscObject)spanstats->dm), &ksp));
330 PetscCall(KSPSetOptionsPrefix(ksp, spanstats->prefix));
331 {
332 PC pc;
333 PetscCall(KSPGetPC(ksp, &pc));
334 PetscCall(PCSetType(pc, PCJACOBI));
335 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL));
336 PetscCall(KSPSetType(ksp, KSPCG));
337 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL));
338 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
339 }
340 PetscCall(KSPSetFromOptions_WithMatCeed(ksp, mat_mass));
341 spanstats->ksp = ksp;
342 PetscCall(MatDestroy(&mat_mass));
343 }
344
345 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
346 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
347 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
348 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stats_proj));
349 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
350 PetscCallCeed(ceed, CeedOperatorDestroy(&op_proj_rhs));
351 PetscFunctionReturn(PETSC_SUCCESS);
352 }
353
SpanwiseStatisticsSetupInitialize(Honee honee,PetscInt degree,PetscInt num_comps,const char * prefix,SpanStatsSetupData * stats_setup_data,SpanStatsCtx * spanstats)354 PetscErrorCode SpanwiseStatisticsSetupInitialize(Honee honee, PetscInt degree, PetscInt num_comps, const char *prefix,
355 SpanStatsSetupData *stats_setup_data, SpanStatsCtx *spanstats) {
356 PetscLogStage stage_stats_setup;
357 SpanStatsCtx spanstats_;
358
359 PetscFunctionBeginUser;
360 PetscCall(PetscLogStageGetId("Stats Setup", &stage_stats_setup));
361 if (stage_stats_setup == -1) PetscCall(PetscLogStageRegister("Stats Setup", &stage_stats_setup));
362 PetscCall(PetscLogStagePush(stage_stats_setup));
363
364 PetscCall(PetscNew(&spanstats_));
365 PetscCall(PetscStrallocpy(prefix, &spanstats_->prefix));
366
367 spanstats_->collect_interval = 1;
368 PetscCall(PetscOptionsGetInt(NULL, prefix, "-collect_interval", &spanstats_->collect_interval, NULL));
369
370 // Create parent DM
371 PetscCall(SpanwiseStatisticssCreateDM(honee, spanstats_, degree, num_comps));
372
373 // Create necessary CeedObjects for setting up statistics
374 PetscCall(SpanwiseStatisticsSetupDataCreate(honee, spanstats_, stats_setup_data));
375 //
376 // Create SF for communicating child data back their respective parents
377 PetscCall(SpanwiseStatisticsCreateSF(honee, spanstats_, *stats_setup_data, spanstats_->dm, honee->dm, &spanstats_->sf));
378
379 *spanstats = spanstats_;
380 PetscFunctionReturn(PETSC_SUCCESS);
381 }
382
SpanwiseStatisticsSetupFinalize(TS ts,Honee honee,SpanStatsCtx spanstats,PetscViewerAndFormat * ctx,SpanStatsSetupData * stats_setup_data)383 PetscErrorCode SpanwiseStatisticsSetupFinalize(TS ts, Honee honee, SpanStatsCtx spanstats, PetscViewerAndFormat *ctx,
384 SpanStatsSetupData *stats_setup_data) {
385 PetscFunctionBeginUser;
386 // Setup KSP and Mat for L^2 projection of statistics
387 PetscCall(SpanwiseStatisticsSetupL2Projection(honee, spanstats, *stats_setup_data));
388
389 PetscCall(PetscViewerSetOptionsPrefix(ctx->viewer, spanstats->prefix));
390 PetscCall(PetscViewerSetFromOptions(ctx->viewer));
391
392 PetscCall(TSGetTime(ts, &spanstats->initial_solution_time));
393 PetscCall(TSGetStepNumber(ts, &spanstats->initial_solution_step));
394 CeedScalar initial_solution_time = spanstats->initial_solution_time; // done for type conversion
395 PetscCallCeed(honee->ceed,
396 CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->previous_time_label, &initial_solution_time));
397
398 PetscCall(SpanwiseStatisticsSetupDataDestroy(*stats_setup_data));
399 PetscCall(PetscLogStagePop());
400 PetscFunctionReturn(PETSC_SUCCESS);
401 }
402
403 // Collect statistics based on the solution Q
SpanwiseStatisticsCollect(Honee honee,SpanStatsCtx spanstats,PetscScalar solution_time,Vec Q)404 PetscErrorCode SpanwiseStatisticsCollect(Honee honee, SpanStatsCtx spanstats, PetscScalar solution_time, Vec Q) {
405 PetscFunctionBeginUser;
406 PetscLogStage stage_stats_collect;
407 PetscCall(PetscLogStageGetId("Stats Collect", &stage_stats_collect));
408 if (stage_stats_collect == -1) PetscCall(PetscLogStageRegister("Stats Collect", &stage_stats_collect));
409 PetscCall(PetscLogStagePush(stage_stats_collect));
410
411 PetscCall(UpdateBoundaryValues(honee, honee->Q_loc, solution_time));
412 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->solution_time_label, &solution_time));
413 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, honee->Q_loc));
414 PetscCall(ApplyAddCeedOperatorLocalToLocal(honee->Q_loc, spanstats->Child_Stats_loc, spanstats->op_stats_collect_ctx));
415
416 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->previous_time_label, &solution_time));
417
418 PetscCall(PetscLogStagePop());
419 PetscFunctionReturn(PETSC_SUCCESS);
420 }
421
422 // Process the child statistics into parent statistics and project them onto stats
SpanwiseStatisticsProcess(Honee honee,SpanStatsCtx spanstats,Vec stats)423 PetscErrorCode SpanwiseStatisticsProcess(Honee honee, SpanStatsCtx spanstats, Vec stats) {
424 const PetscScalar *child_stats;
425 PetscScalar *parent_stats;
426 MPI_Datatype unit;
427 Vec RHS;
428
429 PetscFunctionBeginUser;
430 PetscLogStage stage_stats_process;
431 PetscCall(PetscLogStageGetId("Stats Process", &stage_stats_process));
432 if (stage_stats_process == -1) PetscCall(PetscLogStageRegister("Stats Process", &stage_stats_process));
433 PetscCall(PetscLogStagePush(stage_stats_process));
434
435 PetscCall(VecZeroEntries(spanstats->Parent_Stats_loc));
436
437 PetscCall(VecGetArrayRead(spanstats->Child_Stats_loc, &child_stats));
438 PetscCall(VecGetArray(spanstats->Parent_Stats_loc, &parent_stats));
439
440 if (spanstats->num_comp_stats == 1) unit = MPIU_REAL;
441 else {
442 PetscCallMPI(MPI_Type_contiguous(spanstats->num_comp_stats, MPIU_REAL, &unit));
443 PetscCallMPI(MPI_Type_commit(&unit));
444 }
445
446 PetscCall(PetscSFReduceBegin(spanstats->sf, unit, child_stats, parent_stats, MPI_SUM));
447 PetscCall(PetscSFReduceEnd(spanstats->sf, unit, child_stats, parent_stats, MPI_SUM));
448
449 PetscCall(VecRestoreArrayRead(spanstats->Child_Stats_loc, &child_stats));
450 PetscCall(VecRestoreArray(spanstats->Parent_Stats_loc, &parent_stats));
451 PetscCallMPI(MPI_Type_free(&unit));
452
453 PetscReal solution_time;
454 PetscCall(DMGetOutputSequenceNumber(spanstats->dm, NULL, &solution_time));
455 PetscReal summing_duration = solution_time - honee->app_ctx->cont_time;
456 PetscCall(VecScale(spanstats->Parent_Stats_loc, 1 / (summing_duration * spanstats->span_width)));
457
458 // L^2 projection with the parent_data
459 PetscCall(DMGetGlobalVector(spanstats->dm, &RHS));
460 PetscCall(ApplyCeedOperatorLocalToGlobal(spanstats->Parent_Stats_loc, RHS, spanstats->op_proj_rhs_ctx));
461
462 PetscCall(KSPSolve(spanstats->ksp, RHS, stats));
463
464 PetscCall(DMRestoreGlobalVector(spanstats->dm, &RHS));
465 PetscCall(PetscLogStagePop());
466 PetscFunctionReturn(PETSC_SUCCESS);
467 }
468