1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 /// @file 8 /// Functions for setting up and performing statistics collection 9 10 #include "../qfunctions/turb_spanstats.h" 11 12 #include <petscsf.h> 13 14 #include "../include/matops.h" 15 #include "../navierstokes.h" 16 #include "ceed/ceed.h" 17 #include "petscerror.h" 18 #include "petscmat.h" 19 #include "petscsys.h" 20 #include "petscvec.h" 21 22 PetscErrorCode CreateStatsDM(User user, ProblemData *problem, PetscInt degree, SimpleBC bc) { 23 user->spanstats.num_comp_stats = 22; 24 PetscReal domain_min[3], domain_max[3]; 25 PetscFunctionBeginUser; 26 27 // Get spanwise length 28 PetscCall(DMGetBoundingBox(user->dm, domain_min, domain_max)); 29 user->spanstats.span_width = domain_max[2] - domain_min[1]; 30 31 // Get DM from surface 32 { 33 DMLabel label; 34 PetscCall(DMGetLabel(user->dm, "Face Sets", &label)); 35 PetscCall(DMPlexLabelComplete(user->dm, label)); 36 PetscCall(DMPlexFilter(user->dm, label, 1, &user->spanstats.dm)); 37 PetscCall(DMProjectCoordinates(user->spanstats.dm, NULL)); // Ensure that a coordinate FE exists 38 } 39 40 PetscCall(PetscObjectSetName((PetscObject)user->spanstats.dm, "Spanwise_Stats")); 41 PetscCall(DMSetOptionsPrefix(user->spanstats.dm, "spanstats_")); 42 PetscCall(DMSetFromOptions(user->spanstats.dm)); 43 PetscCall(DMViewFromOptions(user->spanstats.dm, NULL, "-dm_view")); // -spanstats_dm_view 44 { 45 PetscFE fe; 46 DMLabel label; 47 48 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim - 1, user->spanstats.num_comp_stats, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); 49 PetscCall(PetscObjectSetName((PetscObject)fe, "stats")); 50 PetscCall(DMAddField(user->spanstats.dm, NULL, (PetscObject)fe)); 51 PetscCall(DMCreateDS(user->spanstats.dm)); 52 PetscCall(DMGetLabel(user->spanstats.dm, "Face Sets", &label)); 53 54 PetscCall(DMPlexSetClosurePermutationTensor(user->spanstats.dm, PETSC_DETERMINE, NULL)); 55 PetscCall(PetscFEDestroy(&fe)); 56 } 57 58 PetscSection section; 59 PetscCall(DMGetLocalSection(user->spanstats.dm, §ion)); 60 PetscCall(PetscSectionSetFieldName(section, 0, "")); 61 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Mean Density")); 62 PetscCall(PetscSectionSetComponentName(section, 0, 1, "Mean Pressure")); 63 PetscCall(PetscSectionSetComponentName(section, 0, 2, "Mean Pressure Squared")); 64 PetscCall(PetscSectionSetComponentName(section, 0, 3, "Mean Pressure Velocity X")); 65 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Mean Pressure Velocity Y")); 66 PetscCall(PetscSectionSetComponentName(section, 0, 5, "Mean Pressure Velocity Z")); 67 PetscCall(PetscSectionSetComponentName(section, 0, 6, "Mean Density Temperature")); 68 PetscCall(PetscSectionSetComponentName(section, 0, 7, "Mean Density Temperature Flux X")); 69 PetscCall(PetscSectionSetComponentName(section, 0, 8, "Mean Density Temperature Flux Y")); 70 PetscCall(PetscSectionSetComponentName(section, 0, 9, "Mean Density Temperature Flux Z")); 71 PetscCall(PetscSectionSetComponentName(section, 0, 10, "Mean Momentum X")); 72 PetscCall(PetscSectionSetComponentName(section, 0, 11, "Mean Momentum Y")); 73 PetscCall(PetscSectionSetComponentName(section, 0, 12, "Mean Momentum Z")); 74 PetscCall(PetscSectionSetComponentName(section, 0, 13, "Mean Momentum Flux XX")); 75 PetscCall(PetscSectionSetComponentName(section, 0, 14, "Mean Momentum Flux YY")); 76 PetscCall(PetscSectionSetComponentName(section, 0, 15, "Mean Momentum Flux ZZ")); 77 PetscCall(PetscSectionSetComponentName(section, 0, 16, "Mean Momentum Flux YZ")); 78 PetscCall(PetscSectionSetComponentName(section, 0, 17, "Mean Momentum Flux XZ")); 79 PetscCall(PetscSectionSetComponentName(section, 0, 18, "Mean Momentum Flux XY")); 80 PetscCall(PetscSectionSetComponentName(section, 0, 19, "Mean Velocity X")); 81 PetscCall(PetscSectionSetComponentName(section, 0, 20, "Mean Velocity Y")); 82 PetscCall(PetscSectionSetComponentName(section, 0, 21, "Mean Velocity Z")); 83 84 PetscFunctionReturn(0); 85 } 86 87 // Create CeedElemRestriction for collocated data based on associated CeedBasis and CeedElemRestriction 88 // Number of quadrature points is used from the CeedBasis, and number of elements is used from the CeedElemRestriction 89 PetscErrorCode CreateElemRestrColloc(Ceed ceed, CeedInt num_comp, CeedBasis basis, CeedElemRestriction elem_restr_base, 90 CeedElemRestriction *elem_restr_collocated, CeedVector *l_vec, CeedVector *e_vec) { 91 CeedInt num_elem_qpts, loc_num_elem; 92 PetscFunctionBeginUser; 93 94 CeedBasisGetNumQuadraturePoints(basis, &num_elem_qpts); 95 CeedElemRestrictionGetNumElements(elem_restr_base, &loc_num_elem); 96 97 const CeedInt strides[] = {num_comp, 1, num_elem_qpts * num_comp}; 98 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, num_elem_qpts, num_comp, num_comp * loc_num_elem * num_elem_qpts, strides, 99 elem_restr_collocated); 100 CeedElemRestrictionCreateVector(*elem_restr_collocated, l_vec, e_vec); 101 PetscFunctionReturn(0); 102 } 103 104 // Get coordinates of quadrature points 105 PetscErrorCode GetQuadratureCoords(Ceed ceed, DM dm, CeedElemRestriction elem_restr_x, CeedBasis basis_x, CeedVector x_coords, CeedVector *qx_coords, 106 PetscInt *total_nqpnts) { 107 CeedQFunction qf_quad_coords; 108 CeedOperator op_quad_coords; 109 PetscInt num_comp_x, loc_num_elem, num_elem_qpts; 110 CeedElemRestriction elem_restr_qx; 111 PetscFunctionBeginUser; 112 113 // Create Element Restriction and CeedVector for quadrature coordinates 114 CeedBasisGetNumQuadraturePoints(basis_x, &num_elem_qpts); 115 CeedElemRestrictionGetNumElements(elem_restr_x, &loc_num_elem); 116 CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x); 117 *total_nqpnts = num_elem_qpts * loc_num_elem; 118 PetscCall(CreateElemRestrColloc(ceed, num_comp_x, basis_x, elem_restr_x, &elem_restr_qx, qx_coords, NULL)); 119 120 // Create QFunction 121 CeedQFunctionCreateIdentity(ceed, num_comp_x, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_quad_coords); 122 123 // Create Operator 124 CeedOperatorCreate(ceed, qf_quad_coords, NULL, NULL, &op_quad_coords); 125 CeedOperatorSetField(op_quad_coords, "input", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 126 CeedOperatorSetField(op_quad_coords, "output", elem_restr_qx, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 127 128 CeedOperatorApply(op_quad_coords, x_coords, *qx_coords, CEED_REQUEST_IMMEDIATE); 129 130 CeedQFunctionDestroy(&qf_quad_coords); 131 CeedOperatorDestroy(&op_quad_coords); 132 PetscFunctionReturn(0); 133 } 134 135 // Create PetscSF for child-to-parent communication 136 PetscErrorCode CreateStatsSF(Ceed ceed, CeedData ceed_data, DM parentdm, DM childdm, PetscSF statssf) { 137 PetscInt child_num_qpnts, parent_num_qpnts, num_comp_x; 138 CeedVector child_qx_coords, parent_qx_coords; 139 PetscReal *child_coords, *parent_coords; 140 PetscFunctionBeginUser; 141 142 // Assume that child and parent have the same number of components 143 CeedBasisGetNumComponents(ceed_data->basis_x, &num_comp_x); 144 const PetscInt num_comp_sf = num_comp_x - 1; // Number of coord components used in the creation of the SF 145 146 // Get quad_coords for child DM 147 PetscCall(GetQuadratureCoords(ceed, childdm, ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord, &child_qx_coords, &child_num_qpnts)); 148 149 // Get quad_coords for parent DM 150 PetscCall(GetQuadratureCoords(ceed, parentdm, ceed_data->spanstats.elem_restr_parent_x, ceed_data->spanstats.basis_x, ceed_data->spanstats.x_coord, 151 &parent_qx_coords, &parent_num_qpnts)); 152 153 // Remove z component of coordinates for matching 154 { 155 const PetscReal *child_quad_coords, *parent_quad_coords; 156 157 CeedVectorGetArrayRead(child_qx_coords, CEED_MEM_HOST, &child_quad_coords); 158 CeedVectorGetArrayRead(parent_qx_coords, CEED_MEM_HOST, &parent_quad_coords); 159 160 PetscCall(PetscMalloc2(child_num_qpnts * 2, &child_coords, parent_num_qpnts * 2, &parent_coords)); 161 for (int i = 0; i < child_num_qpnts; i++) { 162 child_coords[0 + i * num_comp_sf] = child_quad_coords[0 + i * num_comp_x]; 163 child_coords[1 + i * num_comp_sf] = child_quad_coords[1 + i * num_comp_x]; 164 } 165 for (int i = 0; i < parent_num_qpnts; i++) { 166 parent_coords[0 + i * num_comp_sf] = parent_quad_coords[0 + i * num_comp_x]; 167 parent_coords[1 + i * num_comp_sf] = parent_quad_coords[1 + i * num_comp_x]; 168 } 169 CeedVectorRestoreArrayRead(child_qx_coords, &child_quad_coords); 170 CeedVectorRestoreArrayRead(parent_qx_coords, &parent_quad_coords); 171 } 172 173 PetscCall(PetscSFSetGraphFromCoordinates(statssf, parent_num_qpnts, child_num_qpnts, num_comp_sf, 1e-12, parent_coords, child_coords)); 174 175 PetscCall(PetscSFViewFromOptions(statssf, NULL, "-spanstats_sf_view")); 176 177 PetscCall(PetscFree2(child_coords, parent_coords)); 178 CeedVectorDestroy(&child_qx_coords); 179 CeedVectorDestroy(&parent_qx_coords); 180 PetscFunctionReturn(0); 181 } 182 183 // Compute mass matrix for statistics projection 184 PetscErrorCode SetupL2ProjectionStats(Ceed ceed, User user, CeedData ceed_data) { 185 CeedQFunction qf_mass; 186 CeedOperator op_mass; 187 CeedInt num_comp_q, q_data_size; 188 PetscFunctionBeginUser; 189 190 // CEED Restriction 191 CeedElemRestrictionGetNumComponents(ceed_data->spanstats.elem_restr_parent_stats, &num_comp_q); 192 CeedElemRestrictionGetNumComponents(ceed_data->spanstats.elem_restr_parent_qd, &q_data_size); 193 194 // Create Mass CeedOperator 195 PetscCall(CreateMassQFunction(ceed, num_comp_q, q_data_size, &qf_mass)); 196 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 197 CeedOperatorSetField(op_mass, "q", ceed_data->spanstats.elem_restr_parent_stats, ceed_data->spanstats.basis_stats, CEED_VECTOR_ACTIVE); 198 CeedOperatorSetField(op_mass, "qdata", ceed_data->spanstats.elem_restr_parent_qd, CEED_BASIS_COLLOCATED, ceed_data->spanstats.q_data); 199 CeedOperatorSetField(op_mass, "v", ceed_data->spanstats.elem_restr_parent_stats, ceed_data->spanstats.basis_stats, CEED_VECTOR_ACTIVE); 200 201 // Setup KSP for L^2 projection 202 { 203 MatopApplyContext M_ctx; 204 PetscInt l_size, g_size; 205 Mat mat_mass; 206 VecType vec_type; 207 KSP ksp; 208 Vec ones, M_inv; 209 CeedVector x_ceed, y_ceed; 210 211 PetscCall(DMCreateGlobalVector(user->spanstats.dm, &M_inv)); 212 PetscCall(VecGetLocalSize(M_inv, &l_size)); 213 PetscCall(VecGetSize(M_inv, &g_size)); 214 PetscCall(VecGetType(M_inv, &vec_type)); 215 216 PetscCall(PetscMalloc1(1, &M_ctx)); 217 PetscCall(MatCreateShell(PETSC_COMM_WORLD, l_size, l_size, g_size, g_size, M_ctx, &mat_mass)); 218 PetscCall(MatShellSetOperation(mat_mass, MATOP_MULT, (void (*)(void))MatMult_Ceed)); 219 PetscCall(MatShellSetOperation(mat_mass, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiag_Ceed)); 220 PetscCall(MatShellSetVecType(mat_mass, vec_type)); 221 222 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_stats, &x_ceed, NULL); 223 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_stats, &y_ceed, NULL); 224 225 PetscCall(SetupMatopApplyCtx(PETSC_COMM_WORLD, user->spanstats.dm, user->ceed, op_mass, x_ceed, y_ceed, NULL, M_ctx)); 226 user->spanstats.M_ctx = M_ctx; 227 228 // Create lumped mass matrix inverse 229 PetscCall(DMGetGlobalVector(user->spanstats.dm, &ones)); 230 PetscCall(VecZeroEntries(M_inv)); 231 PetscCall(VecSet(ones, 1)); 232 PetscCall(MatMult(mat_mass, ones, M_inv)); 233 PetscCall(VecReciprocal(M_inv)); 234 user->spanstats.M_inv = M_inv; 235 PetscCall(DMRestoreGlobalVector(user->spanstats.dm, &ones)); 236 237 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 238 PetscCall(KSPSetOptionsPrefix(ksp, "spanstats_")); 239 { 240 PC pc; 241 PetscCall(KSPGetPC(ksp, &pc)); 242 PetscCall(PCSetType(pc, PCJACOBI)); 243 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 244 PetscCall(KSPSetType(ksp, KSPCG)); 245 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 246 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 247 } 248 PetscCall(KSPSetOperators(ksp, mat_mass, mat_mass)); 249 PetscCall(KSPSetFromOptions(ksp)); 250 user->spanstats.ksp = ksp; 251 } 252 253 // Cleanup 254 CeedQFunctionDestroy(&qf_mass); 255 PetscFunctionReturn(0); 256 } 257 258 // Create CeedOperators and KSP for the statistics collection and processing 259 PetscErrorCode CreateStatisticsOperators(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 260 CeedInt num_comp_stats = user->spanstats.num_comp_stats, num_comp_x = problem->dim, num_comp_q; 261 CeedOperator op_setup_sur; 262 PetscFunctionBeginUser; 263 CeedBasisGetNumComponents(ceed_data->basis_q, &num_comp_q); 264 265 // Create Operator for statistics collection 266 switch (user->phys->state_var) { 267 case STATEVAR_PRIMITIVE: 268 CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_Prim, ChildStatsCollection_Prim_loc, &ceed_data->spanstats.qf_stats_collect); 269 break; 270 case STATEVAR_CONSERVATIVE: 271 CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollection_Conserv, ChildStatsCollection_Conserv_loc, &ceed_data->spanstats.qf_stats_collect); 272 break; 273 } 274 275 if (user->app_ctx->stats_test) { 276 CeedQFunctionDestroy(&ceed_data->spanstats.qf_stats_collect); 277 CeedQFunctionCreateInterior(ceed, 1, ChildStatsCollectionTest, ChildStatsCollectionTest_loc, &ceed_data->spanstats.qf_stats_collect); 278 } 279 280 CeedQFunctionSetContext(ceed_data->spanstats.qf_stats_collect, problem->apply_vol_ifunction.qfunction_context); 281 CeedQFunctionAddInput(ceed_data->spanstats.qf_stats_collect, "q", num_comp_q, CEED_EVAL_INTERP); 282 CeedQFunctionAddInput(ceed_data->spanstats.qf_stats_collect, "q_data", problem->q_data_size_vol, CEED_EVAL_NONE); 283 CeedQFunctionAddInput(ceed_data->spanstats.qf_stats_collect, "x", num_comp_x, CEED_EVAL_INTERP); 284 CeedQFunctionAddOutput(ceed_data->spanstats.qf_stats_collect, "v", num_comp_stats, CEED_EVAL_NONE); 285 286 CeedOperatorCreate(ceed, ceed_data->spanstats.qf_stats_collect, NULL, NULL, &user->spanstats.op_stats_collect); 287 CeedOperatorSetField(user->spanstats.op_stats_collect, "q", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 288 CeedOperatorSetField(user->spanstats.op_stats_collect, "q_data", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 289 CeedOperatorSetField(user->spanstats.op_stats_collect, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord); 290 CeedOperatorSetField(user->spanstats.op_stats_collect, "v", ceed_data->spanstats.elem_restr_child_colloc, CEED_BASIS_COLLOCATED, 291 CEED_VECTOR_ACTIVE); 292 293 // Create Operator for statistics projection / processing 294 // Simply take collocated parent data (with quadrature weight already applied) and multiply by weight function. 295 // Therefore, an Identity QF is sufficient 296 CeedQFunctionCreateIdentity(ceed, num_comp_stats, CEED_EVAL_NONE, CEED_EVAL_INTERP, &ceed_data->spanstats.qf_stats_proj); 297 298 CeedOperatorCreate(ceed, ceed_data->spanstats.qf_stats_proj, NULL, NULL, &user->spanstats.op_stats_proj); 299 CeedOperatorSetField(user->spanstats.op_stats_proj, "input", ceed_data->spanstats.elem_restr_parent_colloc, CEED_BASIS_COLLOCATED, 300 CEED_VECTOR_ACTIVE); 301 CeedOperatorSetField(user->spanstats.op_stats_proj, "output", ceed_data->spanstats.elem_restr_parent_stats, ceed_data->spanstats.basis_stats, 302 CEED_VECTOR_ACTIVE); 303 304 // Get q_data for lumped mass matrix formation 305 CeedOperatorCreate(ceed, ceed_data->qf_setup_sur, NULL, NULL, &op_setup_sur); 306 CeedOperatorSetField(op_setup_sur, "dx", ceed_data->spanstats.elem_restr_parent_x, ceed_data->spanstats.basis_x, CEED_VECTOR_ACTIVE); 307 CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, ceed_data->spanstats.basis_x, CEED_VECTOR_NONE); 308 CeedOperatorSetField(op_setup_sur, "surface qdata", ceed_data->spanstats.elem_restr_parent_qd, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 309 CeedOperatorApply(op_setup_sur, ceed_data->spanstats.x_coord, ceed_data->spanstats.q_data, CEED_REQUEST_IMMEDIATE); 310 311 CeedOperatorDestroy(&op_setup_sur); 312 PetscFunctionReturn(0); 313 } 314 315 // Setup for statistics collection 316 PetscErrorCode SetupStatsCollection(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 317 DM dm = user->spanstats.dm; 318 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 319 CeedInt dim, P, Q, num_comp_x; 320 Vec X_loc; 321 PetscMemType X_loc_memtype; 322 const PetscScalar *X_loc_array; 323 PetscFunctionBeginUser; 324 325 PetscCall(DMGetDimension(dm, &dim)); 326 CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &Q); 327 CeedBasisGetNumNodes1D(ceed_data->basis_q, &P); 328 329 PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, Q, problem->q_data_size_sur, &ceed_data->spanstats.elem_restr_parent_stats, 330 &ceed_data->spanstats.elem_restr_parent_x, &ceed_data->spanstats.elem_restr_parent_qd)); 331 CeedElemRestrictionGetNumComponents(ceed_data->spanstats.elem_restr_parent_x, &num_comp_x); 332 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_x, &ceed_data->spanstats.x_coord, NULL); 333 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_stats, &user->spanstats.rhs_ceed, NULL); 334 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_qd, &ceed_data->spanstats.q_data, NULL); 335 336 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS, &ceed_data->spanstats.basis_x); 337 CeedBasisCreateTensorH1Lagrange(ceed, dim, user->spanstats.num_comp_stats, P, Q, CEED_GAUSS, &ceed_data->spanstats.basis_stats); 338 339 PetscCall(CreateElemRestrColloc(ceed, user->spanstats.num_comp_stats, ceed_data->spanstats.basis_stats, 340 ceed_data->spanstats.elem_restr_parent_stats, &ceed_data->spanstats.elem_restr_parent_colloc, 341 &user->spanstats.parent_stats, NULL)); 342 PetscCall(CreateElemRestrColloc(ceed, user->spanstats.num_comp_stats, ceed_data->basis_q, ceed_data->elem_restr_q, 343 &ceed_data->spanstats.elem_restr_child_colloc, &user->spanstats.child_stats, NULL)); 344 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_child_colloc, &user->spanstats.child_inst_stats, NULL); 345 CeedVectorSetValue(user->spanstats.child_stats, 0); 346 347 // -- Copy DM coordinates into CeedVector 348 { 349 DM cdm; 350 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 351 if (cdm) { 352 PetscCall(DMGetCellCoordinatesLocal(dm, &X_loc)); 353 } else { 354 PetscCall(DMGetCoordinatesLocal(dm, &X_loc)); 355 } 356 } 357 PetscCall(VecScale(X_loc, problem->dm_scale)); 358 PetscCall(VecGetArrayReadAndMemType(X_loc, &X_loc_array, &X_loc_memtype)); 359 CeedVectorSetArray(ceed_data->spanstats.x_coord, MemTypeP2C(X_loc_memtype), CEED_COPY_VALUES, (PetscScalar *)X_loc_array); 360 PetscCall(VecRestoreArrayRead(X_loc, &X_loc_array)); 361 362 // Create SF for communicating child data back their respective parents 363 PetscCall(PetscSFCreate(comm, &user->spanstats.sf)); 364 PetscCall(CreateStatsSF(ceed, ceed_data, user->dm, user->spanstats.dm, user->spanstats.sf)); 365 366 // Create CeedOperators for statistics collection 367 PetscCall(CreateStatisticsOperators(ceed, user, ceed_data, problem)); 368 369 // Setup KSP and Mat for L^2 projection of statistics 370 PetscCall(SetupL2ProjectionStats(ceed, user, ceed_data)); 371 372 PetscFunctionReturn(0); 373 } 374 375 // Collect statistics based on the solution Q 376 PetscErrorCode CollectStatistics(User user, PetscScalar solution_time, Vec Q) { 377 PetscMemType q_mem_type; 378 const PetscScalar *q_arr; 379 Vec Q_loc; 380 PetscFunctionBeginUser; 381 382 PetscCall(DMGetLocalVector(user->dm, &Q_loc)); 383 PetscCall(VecZeroEntries(Q_loc)); 384 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 385 386 PetscCall(VecGetArrayReadAndMemType(Q_loc, &q_arr, &q_mem_type)); 387 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, (PetscScalar *)q_arr); 388 389 CeedOperatorApply(user->spanstats.op_stats_collect, user->q_ceed, user->spanstats.child_inst_stats, CEED_REQUEST_IMMEDIATE); 390 391 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 392 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, &q_arr)); 393 PetscCall(DMRestoreLocalVector(user->dm, &Q_loc)); 394 395 // Record averaging using left rectangle rule 396 PetscScalar delta_t = solution_time - user->spanstats.prev_time; 397 PetscScalar prev_timeinterval = user->spanstats.prev_time - user->app_ctx->cont_time; 398 CeedVectorScale(user->spanstats.child_stats, prev_timeinterval / (prev_timeinterval + delta_t)); 399 CeedVectorAXPY(user->spanstats.child_stats, delta_t / (prev_timeinterval + delta_t), user->spanstats.child_inst_stats); 400 user->spanstats.prev_time = solution_time; 401 402 PetscFunctionReturn(0); 403 } 404 405 // Process the child statistics into parent statistics and project them onto stats 406 PetscErrorCode ProcessStatistics(User user, Vec *stats) { 407 Span_Stats user_stats = user->spanstats; 408 const PetscScalar *child_stats; 409 PetscScalar *parent_stats; 410 MPI_Datatype unit; 411 Vec rhs_loc, rhs; 412 PetscMemType rhs_mem_type; 413 CeedScalar *rhs_arr; 414 CeedMemType ceed_mem_type; 415 PetscFunctionBeginUser; 416 417 CeedGetPreferredMemType(user->ceed, &ceed_mem_type); 418 CeedVectorSetValue(user_stats.parent_stats, 0); 419 420 CeedVectorGetArrayRead(user_stats.child_stats, ceed_mem_type, &child_stats); 421 CeedVectorGetArray(user_stats.parent_stats, ceed_mem_type, &parent_stats); 422 423 if (user_stats.num_comp_stats == 1) unit = MPIU_REAL; 424 else { 425 PetscCallMPI(MPI_Type_contiguous(user_stats.num_comp_stats, MPIU_REAL, &unit)); 426 PetscCallMPI(MPI_Type_commit(&unit)); 427 } 428 429 PetscCall(PetscSFReduceBegin(user_stats.sf, unit, child_stats, parent_stats, MPI_SUM)); 430 PetscCall(PetscSFReduceEnd(user_stats.sf, unit, child_stats, parent_stats, MPI_SUM)); 431 432 CeedVectorRestoreArrayRead(user_stats.child_stats, &child_stats); 433 CeedVectorRestoreArray(user_stats.parent_stats, &parent_stats); 434 PetscCallMPI(MPI_Type_free(&unit)); 435 436 CeedVectorScale(user_stats.parent_stats, 1 / user_stats.span_width); 437 438 // L^2 projection with the parent_data 439 PetscCall(DMGetGlobalVector(user_stats.dm, &rhs)); 440 PetscCall(DMGetLocalVector(user_stats.dm, &rhs_loc)); 441 PetscCall(VecZeroEntries(rhs)); 442 PetscCall(VecZeroEntries(rhs_loc)); 443 PetscCall(VecGetArrayWriteAndMemType(rhs_loc, &rhs_arr, &rhs_mem_type)); 444 CeedVectorSetArray(user_stats.rhs_ceed, MemTypeP2C(rhs_mem_type), CEED_USE_POINTER, (PetscScalar *)rhs_arr); 445 446 CeedOperatorApply(user_stats.op_stats_proj, user_stats.parent_stats, user_stats.rhs_ceed, CEED_REQUEST_IMMEDIATE); 447 448 CeedVectorTakeArray(user_stats.rhs_ceed, MemTypeP2C(rhs_mem_type), &rhs_arr); 449 PetscCall(VecRestoreArrayAndMemType(rhs_loc, &rhs_arr)); 450 PetscCall(DMLocalToGlobal(user_stats.dm, rhs_loc, ADD_VALUES, rhs)); 451 452 PetscCall(VecDuplicate(rhs, stats)); 453 PetscCall(VecPointwiseMult(*stats, rhs, user_stats.M_inv)); 454 455 PetscCall(KSPSolve(user_stats.ksp, rhs, *stats)); 456 457 PetscFunctionReturn(0); 458 } 459 460 // TSMonitor for the statistics collection and processing 461 PetscErrorCode TSMonitor_Statistics(TS ts, PetscInt steps, PetscReal solution_time, Vec Q, void *ctx) { 462 User user = (User)ctx; 463 Vec stats; 464 PetscFunctionBeginUser; 465 466 // Do not collect or process on the first step of the run (ie. on the initial condition) 467 if (steps == user->app_ctx->cont_steps) PetscFunctionReturn(0); 468 469 if (steps % user->app_ctx->stats_collect_interval == 0) { 470 PetscCall(CollectStatistics(user, solution_time, Q)); 471 } 472 473 if (steps % user->app_ctx->stats_write_interval == 0 && user->app_ctx->stats_write_interval != -1) { 474 PetscCall(DMGetGlobalVector(user->spanstats.dm, &stats)); 475 PetscCall(ProcessStatistics(user, &stats)); 476 PetscCall(VecViewFromOptions(stats, NULL, "-stats_write_view")); 477 PetscCall(DMRestoreGlobalVector(user->spanstats.dm, &stats)); 478 } 479 PetscFunctionReturn(0); 480 } 481 482 // Function to be called at the end of a simulation 483 PetscErrorCode StatsCollectFinalCall(User user, PetscReal solution_time, Vec Q) { 484 Vec stats; 485 PetscFunctionBeginUser; 486 487 PetscCall(CollectStatistics(user, solution_time, Q)); 488 489 PetscCall(DMGetGlobalVector(user->spanstats.dm, &stats)); 490 PetscCall(ProcessStatistics(user, &stats)); 491 PetscCall(VecViewFromOptions(stats, NULL, "-stats_write_view")); 492 493 PetscScalar *stats_arr; 494 PetscCall(VecGetArray(stats, &stats_arr)); 495 PetscCall(VecRestoreArray(stats, &stats_arr)); 496 PetscCall(DMRestoreGlobalVector(user->spanstats.dm, &stats)); 497 498 PetscFunctionReturn(0); 499 } 500