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