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