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, 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 */ 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, 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 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 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 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 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 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 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 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 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