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