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