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