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