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