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(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 */ 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(CreateElemRestrColloc_CompMajor(ceed, num_comp_stats, honee->basis_q, honee->elem_restr_q, 197 &(*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, SpanStatsCtx spanstats, SpanStatsSetupData stats_setup_data, DM parentdm, DM childdm, 232 PetscSF *statssf) { 233 Ceed ceed = honee->ceed; 234 PetscInt child_num_qpnts, parent_num_qpnts; 235 CeedInt num_comp_x; 236 Vec Child_qx_coords, Parent_qx_coords; 237 238 PetscFunctionBeginUser; 239 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)childdm), statssf)); 240 241 // Assume that child and parent have the same number of components 242 PetscCallCeed(ceed, CeedBasisGetNumComponents(honee->basis_x, &num_comp_x)); 243 const PetscInt num_comp_sf = num_comp_x - 1; // Number of coord components used in the creation of the SF 244 245 // Get quad_coords for child and parent DM 246 PetscCall(GetQuadratureCoords(ceed, childdm, honee->elem_restr_x, honee->basis_x, honee->x_coord, &Child_qx_coords, &child_num_qpnts)); 247 PetscCall(GetQuadratureCoords(ceed, parentdm, stats_setup_data->elem_restr_parent_x, stats_setup_data->basis_x, stats_setup_data->x_coord, 248 &Parent_qx_coords, &parent_num_qpnts)); 249 250 { // Remove z component of coordinates for matching 251 const PetscReal *child_quad_coords, *parent_quad_coords; 252 PetscReal *child_coords, *parent_coords; 253 254 PetscCall(VecGetArrayRead(Child_qx_coords, &child_quad_coords)); 255 PetscCall(VecGetArrayRead(Parent_qx_coords, &parent_quad_coords)); 256 257 PetscCall(PetscMalloc2(child_num_qpnts * 2, &child_coords, parent_num_qpnts * 2, &parent_coords)); 258 for (int i = 0; i < child_num_qpnts; i++) { 259 child_coords[0 + i * num_comp_sf] = child_quad_coords[0 + i * num_comp_x]; 260 child_coords[1 + i * num_comp_sf] = child_quad_coords[1 + i * num_comp_x]; 261 } 262 for (int i = 0; i < parent_num_qpnts; i++) { 263 parent_coords[0 + i * num_comp_sf] = parent_quad_coords[0 + i * num_comp_x]; 264 parent_coords[1 + i * num_comp_sf] = parent_quad_coords[1 + i * num_comp_x]; 265 } 266 PetscCall(VecRestoreArrayRead(Child_qx_coords, &child_quad_coords)); 267 PetscCall(VecRestoreArrayRead(Parent_qx_coords, &parent_quad_coords)); 268 269 PetscCall(PetscSFSetGraphFromCoordinates(*statssf, parent_num_qpnts, child_num_qpnts, num_comp_sf, 1e-12, parent_coords, child_coords)); 270 PetscCall(PetscFree2(child_coords, parent_coords)); 271 } 272 273 { 274 char spanstats_sf_key[1024]; 275 PetscCall(PetscSNPrintf(spanstats_sf_key, sizeof spanstats_sf_key, "-%ssf_view", spanstats->prefix)); 276 PetscCall(PetscObjectSetName((PetscObject)*statssf, spanstats->prefix)); 277 PetscCall(PetscSFViewFromOptions(*statssf, NULL, spanstats_sf_key)); 278 } 279 280 PetscCall(VecDestroy(&Child_qx_coords)); 281 PetscCall(VecDestroy(&Parent_qx_coords)); 282 PetscFunctionReturn(PETSC_SUCCESS); 283 } 284 285 // @brief Setup RHS and LHS for L^2 projection of statistics 286 static PetscErrorCode SpanwiseStatisticsSetupL2Projection(Honee honee, SpanStatsCtx spanstats, SpanStatsSetupData stats_setup_data) { 287 Ceed ceed = honee->ceed; 288 CeedOperator op_mass, op_proj_rhs; 289 CeedQFunction qf_mass, qf_stats_proj; 290 CeedInt q_data_size, num_comp_stats = spanstats->num_comp_stats; 291 CeedElemRestriction elem_restr_qd; 292 CeedVector q_data; 293 DMLabel domain_label = NULL; 294 PetscInt label_value = 0; 295 296 PetscFunctionBeginUser; 297 // -- Create Operator for RHS of L^2 projection of statistics 298 // Simply take collocated parent data (with quadrature weight already applied) and multiply by weight function. 299 // Therefore, an Identity QF is sufficient 300 PetscCallCeed(ceed, CeedQFunctionCreateIdentity(ceed, num_comp_stats, CEED_EVAL_NONE, CEED_EVAL_INTERP, &qf_stats_proj)); 301 302 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stats_proj, NULL, NULL, &op_proj_rhs)); 303 PetscCallCeed(ceed, CeedOperatorSetField(op_proj_rhs, "input", stats_setup_data->elem_restr_parent_colloc, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 304 PetscCallCeed(ceed, CeedOperatorSetField(op_proj_rhs, "output", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, 305 CEED_VECTOR_ACTIVE)); 306 307 PetscCall(OperatorApplyContextCreate(NULL, spanstats->dm, ceed, op_proj_rhs, NULL, NULL, NULL, NULL, &spanstats->op_proj_rhs_ctx)); 308 PetscCall(CeedOperatorCreateLocalVecs(op_proj_rhs, DMReturnVecType(spanstats->dm), PETSC_COMM_SELF, &spanstats->Parent_Stats_loc, NULL)); 309 PetscCall(QDataGet(ceed, spanstats->dm, domain_label, label_value, stats_setup_data->elem_restr_parent_x, stats_setup_data->basis_x, 310 stats_setup_data->x_coord, &elem_restr_qd, &q_data, &q_data_size)); 311 312 // Create Mass CeedOperator 313 PetscCall(HoneeMassQFunctionCreate(ceed, num_comp_stats, q_data_size, &qf_mass)); 314 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass)); 315 PetscCallCeed(ceed, 316 CeedOperatorSetField(op_mass, "u", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, CEED_VECTOR_ACTIVE)); 317 PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data)); 318 PetscCallCeed(ceed, 319 CeedOperatorSetField(op_mass, "v", stats_setup_data->elem_restr_parent_stats, stats_setup_data->basis_stats, CEED_VECTOR_ACTIVE)); 320 321 { // Setup KSP for L^2 projection 322 Mat mat_mass; 323 KSP ksp; 324 325 PetscCall(MatCreateCeed(spanstats->dm, spanstats->dm, op_mass, NULL, &mat_mass)); 326 327 PetscCall(KSPCreate(PetscObjectComm((PetscObject)spanstats->dm), &ksp)); 328 PetscCall(KSPSetOptionsPrefix(ksp, spanstats->prefix)); 329 { 330 PC pc; 331 PetscCall(KSPGetPC(ksp, &pc)); 332 PetscCall(PCSetType(pc, PCJACOBI)); 333 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 334 PetscCall(KSPSetType(ksp, KSPCG)); 335 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 336 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 337 } 338 PetscCall(KSPSetFromOptions_WithMatCeed(ksp, mat_mass)); 339 spanstats->ksp = ksp; 340 PetscCall(MatDestroy(&mat_mass)); 341 } 342 343 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 344 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 345 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass)); 346 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stats_proj)); 347 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass)); 348 PetscCallCeed(ceed, CeedOperatorDestroy(&op_proj_rhs)); 349 PetscFunctionReturn(PETSC_SUCCESS); 350 } 351 352 PetscErrorCode SpanwiseStatisticsSetupInitialize(Honee honee, PetscInt degree, PetscInt num_comps, const char *prefix, 353 SpanStatsSetupData *stats_setup_data, SpanStatsCtx *spanstats) { 354 PetscLogStage stage_stats_setup; 355 SpanStatsCtx spanstats_; 356 357 PetscFunctionBeginUser; 358 PetscCall(PetscLogStageGetId("Stats Setup", &stage_stats_setup)); 359 if (stage_stats_setup == -1) PetscCall(PetscLogStageRegister("Stats Setup", &stage_stats_setup)); 360 PetscCall(PetscLogStagePush(stage_stats_setup)); 361 362 PetscCall(PetscNew(&spanstats_)); 363 PetscCall(PetscStrallocpy(prefix, &spanstats_->prefix)); 364 365 spanstats_->collect_interval = 1; 366 PetscCall(PetscOptionsGetInt(NULL, prefix, "-collect_interval", &spanstats_->collect_interval, NULL)); 367 368 // Create parent DM 369 PetscCall(SpanwiseStatisticssCreateDM(honee, spanstats_, degree, num_comps)); 370 371 // Create necessary CeedObjects for setting up statistics 372 PetscCall(SpanwiseStatisticsSetupDataCreate(honee, spanstats_, stats_setup_data)); 373 // 374 // Create SF for communicating child data back their respective parents 375 PetscCall(SpanwiseStatisticsCreateSF(honee, spanstats_, *stats_setup_data, honee->dm, spanstats_->dm, &spanstats_->sf)); 376 377 *spanstats = spanstats_; 378 PetscFunctionReturn(PETSC_SUCCESS); 379 } 380 381 PetscErrorCode SpanwiseStatisticsSetupFinalize(TS ts, Honee honee, SpanStatsCtx spanstats, PetscViewerAndFormat *ctx, 382 SpanStatsSetupData *stats_setup_data) { 383 PetscFunctionBeginUser; 384 // Setup KSP and Mat for L^2 projection of statistics 385 PetscCall(SpanwiseStatisticsSetupL2Projection(honee, spanstats, *stats_setup_data)); 386 387 PetscCall(PetscViewerSetOptionsPrefix(ctx->viewer, spanstats->prefix)); 388 PetscCall(PetscViewerSetFromOptions(ctx->viewer)); 389 390 PetscCall(TSGetTime(ts, &spanstats->initial_solution_time)); 391 PetscCall(TSGetStepNumber(ts, &spanstats->initial_solution_step)); 392 CeedScalar initial_solution_time = spanstats->initial_solution_time; // done for type conversion 393 PetscCallCeed(honee->ceed, 394 CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->previous_time_label, &initial_solution_time)); 395 396 PetscCall(SpanwiseStatisticsSetupDataDestroy(*stats_setup_data)); 397 PetscCall(PetscLogStagePop()); 398 PetscFunctionReturn(PETSC_SUCCESS); 399 } 400 401 // Collect statistics based on the solution Q 402 PetscErrorCode SpanwiseStatisticsCollect(Honee honee, SpanStatsCtx spanstats, PetscScalar solution_time, Vec Q) { 403 PetscFunctionBeginUser; 404 PetscLogStage stage_stats_collect; 405 PetscCall(PetscLogStageGetId("Stats Collect", &stage_stats_collect)); 406 if (stage_stats_collect == -1) PetscCall(PetscLogStageRegister("Stats Collect", &stage_stats_collect)); 407 PetscCall(PetscLogStagePush(stage_stats_collect)); 408 409 PetscCall(UpdateBoundaryValues(honee, honee->Q_loc, solution_time)); 410 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->solution_time_label, &solution_time)); 411 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, honee->Q_loc)); 412 PetscCall(ApplyAddCeedOperatorLocalToLocal(honee->Q_loc, spanstats->Child_Stats_loc, spanstats->op_stats_collect_ctx)); 413 414 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(spanstats->op_stats_collect_ctx->op, spanstats->previous_time_label, &solution_time)); 415 416 PetscCall(PetscLogStagePop()); 417 PetscFunctionReturn(PETSC_SUCCESS); 418 } 419 420 // Process the child statistics into parent statistics and project them onto stats 421 PetscErrorCode SpanwiseStatisticsProcess(Honee honee, SpanStatsCtx spanstats, Vec stats) { 422 const PetscScalar *child_stats; 423 PetscScalar *parent_stats; 424 MPI_Datatype unit; 425 Vec RHS; 426 427 PetscFunctionBeginUser; 428 PetscLogStage stage_stats_process; 429 PetscCall(PetscLogStageGetId("Stats Process", &stage_stats_process)); 430 if (stage_stats_process == -1) PetscCall(PetscLogStageRegister("Stats Process", &stage_stats_process)); 431 PetscCall(PetscLogStagePush(stage_stats_process)); 432 433 PetscCall(VecZeroEntries(spanstats->Parent_Stats_loc)); 434 435 PetscCall(VecGetArrayRead(spanstats->Child_Stats_loc, &child_stats)); 436 PetscCall(VecGetArray(spanstats->Parent_Stats_loc, &parent_stats)); 437 438 if (spanstats->num_comp_stats == 1) unit = MPIU_REAL; 439 else { 440 PetscCallMPI(MPI_Type_contiguous(spanstats->num_comp_stats, MPIU_REAL, &unit)); 441 PetscCallMPI(MPI_Type_commit(&unit)); 442 } 443 444 PetscCall(PetscSFReduceBegin(spanstats->sf, unit, child_stats, parent_stats, MPI_SUM)); 445 PetscCall(PetscSFReduceEnd(spanstats->sf, unit, child_stats, parent_stats, MPI_SUM)); 446 447 PetscCall(VecRestoreArrayRead(spanstats->Child_Stats_loc, &child_stats)); 448 PetscCall(VecRestoreArray(spanstats->Parent_Stats_loc, &parent_stats)); 449 PetscCallMPI(MPI_Type_free(&unit)); 450 451 PetscReal solution_time; 452 PetscCall(DMGetOutputSequenceNumber(spanstats->dm, NULL, &solution_time)); 453 PetscReal summing_duration = solution_time - honee->app_ctx->cont_time; 454 PetscCall(VecScale(spanstats->Parent_Stats_loc, 1 / (summing_duration * spanstats->span_width))); 455 456 // L^2 projection with the parent_data 457 PetscCall(DMGetGlobalVector(spanstats->dm, &RHS)); 458 PetscCall(ApplyCeedOperatorLocalToGlobal(spanstats->Parent_Stats_loc, RHS, spanstats->op_proj_rhs_ctx)); 459 460 PetscCall(KSPSolve(spanstats->ksp, RHS, stats)); 461 462 PetscCall(DMRestoreGlobalVector(spanstats->dm, &RHS)); 463 PetscCall(PetscLogStagePop()); 464 PetscFunctionReturn(PETSC_SUCCESS); 465 } 466