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