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(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 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, 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( 111 DMSetupByOrder_FEM(PETSC_TRUE, PETSC_TRUE, honee->app_ctx->degree, 1, honee->app_ctx->q_extra, 1, &spanstats->num_comp_stats, 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 */ 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 136 static PetscErrorCode GetQuadratureCoords(Ceed ceed, DM dm, CeedElemRestriction elem_restr_x, CeedBasis basis_x, CeedVector x_coords, Vec *Qx_coords, 137 PetscInt *total_nqpnts) { 138 CeedElemRestriction elem_restr_qx; 139 CeedQFunction qf_quad_coords; 140 CeedOperator op_quad_coords; 141 CeedInt num_comp_x; 142 CeedSize l_vec_size; 143 OperatorApplyContext op_quad_coords_ctx; 144 145 PetscFunctionBeginUser; 146 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x)); 147 PetscCall(CreateElemRestrColloc_CompMajor(ceed, num_comp_x, basis_x, elem_restr_x, &elem_restr_qx)); 148 PetscCallCeed(ceed, CeedElemRestrictionGetLVectorSize(elem_restr_qx, &l_vec_size)); 149 *total_nqpnts = l_vec_size / num_comp_x; 150 151 // Create QFunction 152 PetscCallCeed(ceed, CeedQFunctionCreateIdentity(ceed, num_comp_x, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_quad_coords)); 153 154 // Create Operator 155 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_quad_coords, NULL, NULL, &op_quad_coords)); 156 PetscCallCeed(ceed, CeedOperatorSetField(op_quad_coords, "input", elem_restr_x, basis_x, x_coords)); 157 PetscCallCeed(ceed, CeedOperatorSetField(op_quad_coords, "output", elem_restr_qx, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 158 159 PetscCall(CeedOperatorCreateLocalVecs(op_quad_coords, DMReturnVecType(dm), PETSC_COMM_SELF, NULL, Qx_coords)); 160 PetscCall(OperatorApplyContextCreate(NULL, NULL, ceed, op_quad_coords, CEED_VECTOR_NONE, NULL, NULL, NULL, &op_quad_coords_ctx)); 161 162 PetscCall(ApplyCeedOperatorLocalToLocal(NULL, *Qx_coords, op_quad_coords_ctx)); 163 164 PetscCall(OperatorApplyContextDestroy(op_quad_coords_ctx)); 165 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qx)); 166 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_quad_coords)); 167 PetscCallCeed(ceed, CeedOperatorDestroy(&op_quad_coords)); 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170 171 static PetscErrorCode SpanwiseStatisticsSetupDataCreate(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData *stats_setup_data) { 172 Ceed ceed = honee->ceed; 173 DM dm = spanstats->dm; 174 Vec X_loc; 175 DMLabel domain_label = NULL; 176 CeedInt num_comp_x, num_comp_stats = spanstats->num_comp_stats; 177 PetscInt label_value = 0, height = 0, dm_field = 0; 178 179 PetscFunctionBeginUser; 180 PetscCall(PetscNew(stats_setup_data)); 181 182 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, label_value, height, dm_field, &(*stats_setup_data)->elem_restr_parent_stats)); 183 PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height, &(*stats_setup_data)->elem_restr_parent_x)); 184 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents((*stats_setup_data)->elem_restr_parent_x, &num_comp_x)); 185 PetscCallCeed(ceed, CeedElemRestrictionCreateVector((*stats_setup_data)->elem_restr_parent_x, &(*stats_setup_data)->x_coord, NULL)); 186 187 { 188 DM dm_coord; 189 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 190 PetscCall(CreateBasisFromPlex(ceed, dm_coord, domain_label, label_value, height, dm_field, &(*stats_setup_data)->basis_x)); 191 PetscCall(CreateBasisFromPlex(ceed, dm, domain_label, label_value, height, dm_field, &(*stats_setup_data)->basis_stats)); 192 } 193 194 PetscCall(CreateElemRestrColloc_CompMajor(ceed, num_comp_stats, (*stats_setup_data)->basis_stats, (*stats_setup_data)->elem_restr_parent_stats, 195 &(*stats_setup_data)->elem_restr_parent_colloc)); 196 PetscCall( 197 CreateElemRestrColloc_CompMajor(ceed, num_comp_stats, honee->basis_q, honee->elem_restr_q, &(*stats_setup_data)->elem_restr_child_colloc)); 198 199 { // -- Copy DM coordinates into CeedVector 200 DM cdm; 201 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 202 if (cdm) { 203 PetscCall(DMGetCellCoordinatesLocal(dm, &X_loc)); 204 } else { 205 PetscCall(DMGetCoordinatesLocal(dm, &X_loc)); 206 } 207 } 208 // Do not need to scale the coordinates. The DMClone operation has already copied the scaled coordinates 209 PetscCall(VecCopyPetscToCeed(X_loc, (*stats_setup_data)->x_coord)); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212 213 static PetscErrorCode SpanwiseStatisticsSetupDataDestroy(SpanStatsSetupData data) { 214 Ceed ceed; 215 216 PetscFunctionBeginUser; 217 PetscCall(CeedElemRestrictionGetCeed(data->elem_restr_parent_x, &ceed)); 218 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_parent_x)); 219 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_parent_stats)); 220 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_parent_colloc)); 221 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&data->elem_restr_child_colloc)); 222 PetscCallCeed(ceed, CeedBasisDestroy(&data->basis_x)); 223 PetscCallCeed(ceed, CeedBasisDestroy(&data->basis_stats)); 224 PetscCallCeed(ceed, CeedVectorDestroy(&data->x_coord)); 225 PetscCall(PetscFree(data)); 226 PetscCheck(CeedDestroy(&ceed) == CEED_ERROR_SUCCESS, PETSC_COMM_SELF, PETSC_ERR_LIB, "Destroying Ceed object failed"); 227 PetscFunctionReturn(PETSC_SUCCESS); 228 } 229 230 // Create PetscSF for child-to-parent communication 231 static PetscErrorCode SpanwiseStatisticsCreateSF(Honee honee, SpanStatsSetupData stats_setup_data, DM parentdm, DM childdm, PetscSF *statssf) { 232 Ceed ceed = honee->ceed; 233 PetscInt child_num_qpnts, parent_num_qpnts; 234 CeedInt num_comp_x; 235 Vec Child_qx_coords, Parent_qx_coords; 236 237 PetscFunctionBeginUser; 238 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)childdm), statssf)); 239 240 // Assume that child and parent have the same number of components 241 PetscCallCeed(ceed, CeedBasisGetNumComponents(honee->basis_x, &num_comp_x)); 242 const PetscInt num_comp_sf = num_comp_x - 1; // Number of coord components used in the creation of the SF 243 244 // Get quad_coords for child and parent DM 245 PetscCall(GetQuadratureCoords(ceed, childdm, honee->elem_restr_x, honee->basis_x, honee->x_coord, &Child_qx_coords, &child_num_qpnts)); 246 PetscCall(GetQuadratureCoords(ceed, parentdm, stats_setup_data->elem_restr_parent_x, stats_setup_data->basis_x, stats_setup_data->x_coord, 247 &Parent_qx_coords, &parent_num_qpnts)); 248 249 { // Remove z component of coordinates for matching 250 const PetscReal *child_quad_coords, *parent_quad_coords; 251 PetscReal *child_coords, *parent_coords; 252 253 PetscCall(VecGetArrayRead(Child_qx_coords, &child_quad_coords)); 254 PetscCall(VecGetArrayRead(Parent_qx_coords, &parent_quad_coords)); 255 256 PetscCall(PetscMalloc2(child_num_qpnts * 2, &child_coords, parent_num_qpnts * 2, &parent_coords)); 257 for (int i = 0; i < child_num_qpnts; i++) { 258 child_coords[0 + i * num_comp_sf] = child_quad_coords[0 + i * num_comp_x]; 259 child_coords[1 + i * num_comp_sf] = child_quad_coords[1 + i * num_comp_x]; 260 } 261 for (int i = 0; i < parent_num_qpnts; i++) { 262 parent_coords[0 + i * num_comp_sf] = parent_quad_coords[0 + i * num_comp_x]; 263 parent_coords[1 + i * num_comp_sf] = parent_quad_coords[1 + i * num_comp_x]; 264 } 265 PetscCall(VecRestoreArrayRead(Child_qx_coords, &child_quad_coords)); 266 PetscCall(VecRestoreArrayRead(Parent_qx_coords, &parent_quad_coords)); 267 268 PetscCall(PetscSFSetGraphFromCoordinates(*statssf, parent_num_qpnts, child_num_qpnts, num_comp_sf, 1e-12, parent_coords, child_coords)); 269 PetscCall(PetscFree2(child_coords, parent_coords)); 270 } 271 272 PetscCall(PetscSFViewFromOptions(*statssf, NULL, "-spanstats_sf_view")); 273 274 PetscCall(VecDestroy(&Child_qx_coords)); 275 PetscCall(VecDestroy(&Parent_qx_coords)); 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 // @brief Setup RHS and LHS for L^2 projection of statistics 280 static PetscErrorCode SpanwiseStatisticsSetupL2Projection(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_setup_data) { 281 Ceed ceed = honee->ceed; 282 CeedOperator op_mass, op_proj_rhs; 283 CeedQFunction qf_mass, qf_stats_proj; 284 CeedInt q_data_size, num_comp_stats = spanstats->num_comp_stats; 285 CeedElemRestriction elem_restr_qd; 286 CeedVector q_data; 287 DMLabel domain_label = NULL; 288 PetscInt label_value = 0; 289 290 PetscFunctionBeginUser; 291 // -- Create Operator for RHS of L^2 projection of statistics 292 // Simply take collocated parent data (with quadrature weight already applied) and multiply by weight function. 293 // Therefore, an Identity QF is sufficient 294 PetscCallCeed(ceed, CeedQFunctionCreateIdentity(ceed, num_comp_stats, CEED_EVAL_NONE, CEED_EVAL_INTERP, &qf_stats_proj)); 295 296 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stats_proj, NULL, NULL, &op_proj_rhs)); 297 PetscCallCeed(ceed, CeedOperatorSetField(op_proj_rhs, "input", stats_setup_data->elem_restr_parent_colloc, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 298 PetscCallCeed(ceed, CeedOperatorSetField(op_proj_rhs, "output", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, 299 CEED_VECTOR_ACTIVE)); 300 301 PetscCall(OperatorApplyContextCreate(NULL, spanstats->dm, ceed, op_proj_rhs, NULL, NULL, NULL, NULL, &spanstats->op_proj_rhs_ctx)); 302 PetscCall(CeedOperatorCreateLocalVecs(op_proj_rhs, DMReturnVecType(spanstats->dm), PETSC_COMM_SELF, &spanstats->Parent_Stats_loc, NULL)); 303 PetscCall(QDataGet(ceed, spanstats->dm, domain_label, label_value, stats_setup_data->elem_restr_parent_x, stats_setup_data->basis_x, 304 stats_setup_data->x_coord, &elem_restr_qd, &q_data, &q_data_size)); 305 306 // Create Mass CeedOperator 307 PetscCall(HoneeMassQFunctionCreate(ceed, num_comp_stats, q_data_size, &qf_mass)); 308 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 309 PetscCallCeed(ceed, 310 CeedOperatorSetField(op_mass, "u", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, CEED_VECTOR_ACTIVE)); 311 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 312 PetscCallCeed(ceed, 313 CeedOperatorSetField(op_mass, "v", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, CEED_VECTOR_ACTIVE)); 314 315 { // Setup KSP for L^2 projection 316 Mat mat_mass; 317 KSP ksp; 318 319 PetscCall(MatCreateCeed(spanstats->dm, spanstats->dm, op_mass, NULL, &mat_mass)); 320 321 PetscCall(KSPCreate(PetscObjectComm((PetscObject)spanstats->dm), &ksp)); 322 PetscCall(KSPSetOptionsPrefix(ksp, spanstats->prefix)); 323 { 324 PC pc; 325 PetscCall(KSPGetPC(ksp, &pc)); 326 PetscCall(PCSetType(pc, PCJACOBI)); 327 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 328 PetscCall(KSPSetType(ksp, KSPCG)); 329 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 330 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 331 } 332 PetscCall(KSPSetFromOptions_WithMatCeed(ksp, mat_mass)); 333 spanstats->ksp = ksp; 334 PetscCall(MatDestroy(&mat_mass)); 335 } 336 337 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 338 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 339 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 340 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stats_proj)); 341 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 342 PetscCallCeed(ceed, CeedOperatorDestroy(&op_proj_rhs)); 343 PetscFunctionReturn(PETSC_SUCCESS); 344 } 345 346 PetscErrorCode SpanwiseStatisticsSetupInitialize(Honee honee, PetscInt degree, PetscInt num_comps, const char *prefix, 347 SpanStatsSetupData *stats_setup_data, SpanStatsCtx *spanstats) { 348 PetscLogStage stage_stats_setup; 349 SpanStatsCtx spanstats_; 350 351 PetscFunctionBeginUser; 352 PetscCall(PetscLogStageGetId("Stats Setup", &stage_stats_setup)); 353 if (stage_stats_setup == -1) PetscCall(PetscLogStageRegister("Stats Setup", &stage_stats_setup)); 354 PetscCall(PetscLogStagePush(stage_stats_setup)); 355 356 PetscCall(PetscNew(&spanstats_)); 357 PetscCall(PetscStrallocpy(prefix, &spanstats_->prefix)); 358 359 spanstats_->collect_interval = 1; 360 PetscCall(PetscOptionsGetInt(NULL, prefix, "-collect_interval", &spanstats_->collect_interval, NULL)); 361 362 // Create parent DM 363 PetscCall(SpanwiseStatisticssCreateDM(honee, spanstats_, degree, num_comps)); 364 365 // Create necessary CeedObjects for setting up statistics 366 PetscCall(SpanwiseStatisticsSetupDataCreate(honee, spanstats_, stats_setup_data)); 367 // 368 // Create SF for communicating child data back their respective parents 369 PetscCall(SpanwiseStatisticsCreateSF(honee, *stats_setup_data, honee->dm, spanstats_->dm, &spanstats_->sf)); 370 371 *spanstats = spanstats_; 372 PetscFunctionReturn(PETSC_SUCCESS); 373 } 374 375 PetscErrorCode SpanwiseStatisticsSetupFinalize(TS ts, Honee honee, SpanStatsCtx spanstats, PetscViewerAndFormat *ctx, 376 SpanStatsSetupData *stats_setup_data) { 377 PetscFunctionBeginUser; 378 // Setup KSP and Mat for L^2 projection of statistics 379 PetscCall(SpanwiseStatisticsSetupL2Projection(honee, spanstats, *stats_setup_data)); 380 381 PetscCall(PetscViewerSetOptionsPrefix(ctx->viewer, spanstats->prefix)); 382 PetscCall(PetscViewerSetFromOptions(ctx->viewer)); 383 384 PetscCall(TSGetTime(ts, &spanstats->initial_solution_time)); 385 PetscCall(TSGetStepNumber(ts, &spanstats->initial_solution_step)); 386 CeedScalar initial_solution_time = spanstats->initial_solution_time; // done for type conversion 387 PetscCallCeed(honee->ceed, 388 CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->previous_time_label, &initial_solution_time)); 389 390 PetscCall(SpanwiseStatisticsSetupDataDestroy(*stats_setup_data)); 391 PetscCall(PetscLogStagePop()); 392 PetscFunctionReturn(PETSC_SUCCESS); 393 } 394 395 // Collect statistics based on the solution Q 396 PetscErrorCode SpanwiseStatisticsCollect(Honee honee, SpanStatsCtx spanstats, PetscScalar solution_time, Vec Q) { 397 PetscFunctionBeginUser; 398 PetscLogStage stage_stats_collect; 399 PetscCall(PetscLogStageGetId("Stats Collect", &stage_stats_collect)); 400 if (stage_stats_collect == -1) PetscCall(PetscLogStageRegister("Stats Collect", &stage_stats_collect)); 401 PetscCall(PetscLogStagePush(stage_stats_collect)); 402 403 PetscCall(UpdateBoundaryValues(honee, honee->Q_loc, solution_time)); 404 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->solution_time_label, &solution_time)); 405 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, honee->Q_loc)); 406 PetscCall(ApplyAddCeedOperatorLocalToLocal(honee->Q_loc, spanstats->Child_Stats_loc, spanstats->op_stats_collect_ctx)); 407 408 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->previous_time_label, &solution_time)); 409 410 PetscCall(PetscLogStagePop()); 411 PetscFunctionReturn(PETSC_SUCCESS); 412 } 413 414 // Process the child statistics into parent statistics and project them onto stats 415 PetscErrorCode SpanwiseStatisticsProcess(Honee honee, SpanStatsCtx spanstats, Vec stats) { 416 const PetscScalar *child_stats; 417 PetscScalar *parent_stats; 418 MPI_Datatype unit; 419 Vec RHS; 420 421 PetscFunctionBeginUser; 422 PetscLogStage stage_stats_process; 423 PetscCall(PetscLogStageGetId("Stats Process", &stage_stats_process)); 424 if (stage_stats_process == -1) PetscCall(PetscLogStageRegister("Stats Process", &stage_stats_process)); 425 PetscCall(PetscLogStagePush(stage_stats_process)); 426 427 PetscCall(VecZeroEntries(spanstats->Parent_Stats_loc)); 428 429 PetscCall(VecGetArrayRead(spanstats->Child_Stats_loc, &child_stats)); 430 PetscCall(VecGetArray(spanstats->Parent_Stats_loc, &parent_stats)); 431 432 if (spanstats->num_comp_stats == 1) unit = MPIU_REAL; 433 else { 434 PetscCallMPI(MPI_Type_contiguous(spanstats->num_comp_stats, MPIU_REAL, &unit)); 435 PetscCallMPI(MPI_Type_commit(&unit)); 436 } 437 438 PetscCall(PetscSFReduceBegin(spanstats->sf, unit, child_stats, parent_stats, MPI_SUM)); 439 PetscCall(PetscSFReduceEnd(spanstats->sf, unit, child_stats, parent_stats, MPI_SUM)); 440 441 PetscCall(VecRestoreArrayRead(spanstats->Child_Stats_loc, &child_stats)); 442 PetscCall(VecRestoreArray(spanstats->Parent_Stats_loc, &parent_stats)); 443 PetscCallMPI(MPI_Type_free(&unit)); 444 445 PetscReal solution_time; 446 PetscCall(DMGetOutputSequenceNumber(spanstats->dm, NULL, &solution_time)); 447 PetscReal summing_duration = solution_time - honee->app_ctx->cont_time; 448 PetscCall(VecScale(spanstats->Parent_Stats_loc, 1 / (summing_duration * spanstats->span_width))); 449 450 // L^2 projection with the parent_data 451 PetscCall(DMGetGlobalVector(spanstats->dm, &RHS)); 452 PetscCall(ApplyCeedOperatorLocalToGlobal(spanstats->Parent_Stats_loc, RHS, spanstats->op_proj_rhs_ctx)); 453 454 PetscCall(KSPSolve(spanstats->ksp, RHS, stats)); 455 456 PetscCall(DMRestoreGlobalVector(spanstats->dm, &RHS)); 457 PetscCall(PetscLogStagePop()); 458 PetscFunctionReturn(PETSC_SUCCESS); 459 } 460