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 8 /// @file 9 /// Utility functions for setting up statistics collection 10 11 #include <petscsf.h> 12 13 #include "../navierstokes.h" 14 #include "petscsys.h" 15 16 PetscErrorCode CreateStatsDM(User user, ProblemData *problem, PetscInt degree, SimpleBC bc) { 17 user->spanstats.num_comp_stats = 1; 18 PetscFunctionBeginUser; 19 20 // Get DM from surface 21 { 22 DMLabel label; 23 PetscCall(DMGetLabel(user->dm, "Face Sets", &label)); 24 PetscCall(DMPlexLabelComplete(user->dm, label)); 25 PetscCall(DMPlexFilter(user->dm, label, 1, &user->spanstats.dm)); 26 PetscCall(DMProjectCoordinates(user->spanstats.dm, NULL)); // Ensure that a coordinate FE exists 27 } 28 29 PetscCall(PetscObjectSetName((PetscObject)user->spanstats.dm, "Spanwise_Stats")); 30 PetscCall(DMSetOptionsPrefix(user->spanstats.dm, "spanstats_")); 31 PetscCall(PetscOptionsSetValue(NULL, "-spanstats_dm_sparse_localize", "0")); // [Jed] Not relevant because not periodic in this direction 32 33 PetscCall(DMSetFromOptions(user->spanstats.dm)); 34 PetscCall(DMViewFromOptions(user->spanstats.dm, NULL, "-dm_view")); // -spanstats_dm_view (option includes prefix) 35 { 36 PetscFE fe; 37 DMLabel label; 38 39 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim - 1, user->spanstats.num_comp_stats, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); 40 PetscCall(PetscObjectSetName((PetscObject)fe, "stats")); 41 PetscCall(DMAddField(user->spanstats.dm, NULL, (PetscObject)fe)); 42 PetscCall(DMCreateDS(user->spanstats.dm)); 43 PetscCall(DMGetLabel(user->spanstats.dm, "Face Sets", &label)); 44 45 // // Set wall BCs 46 // if (bc->num_wall > 0) { 47 // PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, bc->num_wall, bc->walls, 0, bc->num_comps, bc->wall_comps, 48 // (void (*)(void))problem->bc, NULL, problem->bc_ctx, NULL)); 49 // } 50 // // Set slip BCs in the x direction 51 // if (bc->num_slip[0] > 0) { 52 // PetscInt comps[1] = {1}; 53 // PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", label, bc->num_slip[0], bc->slips[0], 0, 1, comps, (void (*)(void))NULL, NULL, 54 // problem->bc_ctx, NULL)); 55 // } 56 // // Set slip BCs in the y direction 57 // if (bc->num_slip[1] > 0) { 58 // PetscInt comps[1] = {2}; 59 // PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", label, bc->num_slip[1], bc->slips[1], 0, 1, comps, (void (*)(void))NULL, NULL, 60 // problem->bc_ctx, NULL)); 61 // } 62 // // Set slip BCs in the z direction 63 // if (bc->num_slip[2] > 0) { 64 // PetscInt comps[1] = {3}; 65 // PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", label, bc->num_slip[2], bc->slips[2], 0, 1, comps, (void (*)(void))NULL, NULL, 66 // problem->bc_ctx, NULL)); 67 // } 68 69 PetscCall(DMPlexSetClosurePermutationTensor(user->spanstats.dm, PETSC_DETERMINE, NULL)); 70 PetscCall(PetscFEDestroy(&fe)); 71 } 72 73 PetscSection section; 74 PetscCall(DMGetLocalSection(user->spanstats.dm, §ion)); 75 PetscCall(PetscSectionSetFieldName(section, 0, "")); 76 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Test")); 77 // PetscCall(PetscSectionSetComponentName(section, 0, 0, "Mean Velocity Products XX")); 78 // PetscCall(PetscSectionSetComponentName(section, 0, 1, "Mean Velocity Products YY")); 79 // PetscCall(PetscSectionSetComponentName(section, 0, 2, "Mean Velocity Products ZZ")); 80 // PetscCall(PetscSectionSetComponentName(section, 0, 3, "Mean Velocity Products YZ")); 81 // PetscCall(PetscSectionSetComponentName(section, 0, 4, "Mean Velocity Products XZ")); 82 // PetscCall(PetscSectionSetComponentName(section, 0, 5, "Mean Velocity Products XY")); 83 84 Vec test; 85 PetscCall(DMCreateLocalVector(user->spanstats.dm, &test)); 86 PetscCall(VecZeroEntries(test)); 87 PetscCall(VecViewFromOptions(test, NULL, "-test_view")); 88 89 PetscFunctionReturn(0); 90 } 91 92 PetscErrorCode GetQuadratureCoords(Ceed ceed, DM dm, CeedElemRestriction elem_restr_x, CeedBasis basis_x, CeedVector x_coords, CeedVector *qx_coords, 93 PetscInt *total_nqpnts) { 94 CeedQFunction qf_quad_coords; 95 CeedOperator op_quad_coords; 96 PetscInt num_comp_x, loc_num_elem, num_elem_qpts; 97 CeedElemRestriction elem_restr_qx; 98 PetscFunctionBeginUser; 99 100 // Create Element Restriction and CeedVector for quadrature coordinates 101 CeedBasisGetNumQuadraturePoints(basis_x, &num_elem_qpts); 102 CeedElemRestrictionGetNumElements(elem_restr_x, &loc_num_elem); 103 CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x); 104 *total_nqpnts = num_elem_qpts * loc_num_elem; 105 const CeedInt strides[] = {num_comp_x, 1, num_elem_qpts * num_comp_x}; 106 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, num_elem_qpts, num_comp_x, num_comp_x * loc_num_elem * num_elem_qpts, strides, &elem_restr_qx); 107 CeedElemRestrictionCreateVector(elem_restr_qx, qx_coords, NULL); 108 109 // Create QFunction 110 CeedQFunctionCreateIdentity(ceed, num_comp_x, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_quad_coords); 111 112 // Create Operator 113 CeedOperatorCreate(ceed, qf_quad_coords, NULL, NULL, &op_quad_coords); 114 CeedOperatorSetField(op_quad_coords, "input", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 115 CeedOperatorSetField(op_quad_coords, "output", elem_restr_qx, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 116 117 CeedOperatorApply(op_quad_coords, x_coords, *qx_coords, CEED_REQUEST_IMMEDIATE); 118 119 CeedQFunctionDestroy(&qf_quad_coords); 120 CeedOperatorDestroy(&op_quad_coords); 121 PetscFunctionReturn(0); 122 } 123 124 PetscErrorCode CreateStatsSF(Ceed ceed, CeedData ceed_data, DM parentdm, DM childdm, PetscSF statssf) { 125 PetscInt child_num_qpnts, parent_num_qpnts, num_comp_x; 126 CeedVector child_qx_coords, parent_qx_coords; 127 PetscReal *child_coords, *parent_coords; 128 PetscFunctionBeginUser; 129 130 // Assume that child and parent have the same number of components 131 CeedBasisGetNumComponents(ceed_data->basis_x, &num_comp_x); 132 const PetscInt num_comp_sf = num_comp_x - 1; // Number of coord components used in the creation of the SF 133 134 // Get quad_coords for child DM 135 PetscCall(GetQuadratureCoords(ceed, childdm, ceed_data->elem_restr_x, ceed_data->basis_xc, ceed_data->x_coord, &child_qx_coords, &child_num_qpnts)); 136 137 // Get quad_coords for parent DM 138 PetscCall(GetQuadratureCoords(ceed, parentdm, ceed_data->spanstats.elem_restr_parent_x, ceed_data->spanstats.basis_x, ceed_data->spanstats.x_coord, 139 &parent_qx_coords, &parent_num_qpnts)); 140 141 // Remove z component of coordinates for matching 142 { 143 const PetscReal *child_quad_coords, *parent_quad_coords; 144 145 CeedVectorGetArrayRead(child_qx_coords, CEED_MEM_HOST, &child_quad_coords); 146 CeedVectorGetArrayRead(parent_qx_coords, CEED_MEM_HOST, &parent_quad_coords); 147 148 PetscCall(PetscMalloc2(child_num_qpnts * 2, &child_coords, parent_num_qpnts * 2, &parent_coords)); 149 for (int i = 0; i < child_num_qpnts; i++) { 150 child_coords[0 + i * num_comp_sf] = child_quad_coords[0 + i * num_comp_x]; 151 child_coords[1 + i * num_comp_sf] = child_quad_coords[1 + i * num_comp_x]; 152 } 153 for (int i = 0; i < parent_num_qpnts; i++) { 154 parent_coords[0 + i * num_comp_sf] = parent_quad_coords[0 + i * num_comp_x]; 155 parent_coords[1 + i * num_comp_sf] = parent_quad_coords[1 + i * num_comp_x]; 156 } 157 CeedVectorRestoreArrayRead(child_qx_coords, &child_quad_coords); 158 CeedVectorRestoreArrayRead(parent_qx_coords, &parent_quad_coords); 159 } 160 161 // Only check the first two components of the coordinates 162 PetscCall(PetscSFSetGraphFromCoordinates(statssf, parent_num_qpnts, child_num_qpnts, num_comp_sf, 1e-12, parent_coords, child_coords)); 163 164 PetscCall(PetscFree2(child_coords, parent_coords)); 165 CeedVectorDestroy(&ceed_data->spanstats.x_coord); 166 CeedVectorDestroy(&child_qx_coords); 167 CeedVectorDestroy(&parent_qx_coords); 168 PetscFunctionReturn(0); 169 } 170 171 PetscErrorCode SetupStatsCollection(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 172 DM dm = user->spanstats.dm; 173 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 174 CeedInt dim, P, Q, num_comp_x; 175 Vec X_loc; 176 PetscMemType X_loc_memtype; 177 const PetscScalar *X_loc_array; 178 PetscFunctionBeginUser; 179 180 PetscCall(DMGetDimension(dm, &dim)); 181 CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_x, &Q); 182 CeedBasisGetNumNodes1D(ceed_data->basis_x, &P); 183 184 // TODO: Possibly need to create a elem_restr_qcollocated for the global domain as well 185 PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, Q, user->spanstats.num_comp_stats, &ceed_data->spanstats.elem_restr_parent_stats, 186 &ceed_data->spanstats.elem_restr_parent_x, &ceed_data->spanstats.elem_restr_parent_qd)); 187 CeedElemRestrictionGetNumComponents(ceed_data->spanstats.elem_restr_parent_x, &num_comp_x); 188 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_x, &ceed_data->spanstats.x_coord, NULL); 189 190 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS_LOBATTO, &ceed_data->spanstats.basis_x); 191 CeedBasisCreateTensorH1Lagrange(ceed, dim, user->spanstats.num_comp_stats, P, Q, CEED_GAUSS, &ceed_data->spanstats.basis_stats); 192 193 // -- Copy DM coordinates into CeedVector 194 { 195 DM cdm; 196 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 197 if (cdm) { 198 PetscCall(DMGetCellCoordinatesLocal(dm, &X_loc)); 199 } else { 200 PetscCall(DMGetCoordinatesLocal(dm, &X_loc)); 201 } 202 } 203 PetscCall(VecScale(X_loc, problem->dm_scale)); 204 PetscCall(VecGetArrayReadAndMemType(X_loc, &X_loc_array, &X_loc_memtype)); 205 CeedVectorSetArray(ceed_data->spanstats.x_coord, MemTypeP2C(X_loc_memtype), CEED_COPY_VALUES, (PetscScalar *)X_loc_array); 206 PetscCall(VecRestoreArrayRead(X_loc, &X_loc_array)); 207 208 // Create SF for communicating child data back their respective parents 209 PetscCall(PetscSFCreate(comm, &user->spanstats.sf)); 210 PetscCall(CreateStatsSF(ceed, ceed_data, user->dm, user->spanstats.dm, user->spanstats.sf)); 211 212 PetscFunctionReturn(0); 213 } 214