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 // Create CeedElemRestriction for collocated data based on associated CeedBasis and CeedElemRestriction 93 // Number of quadrature points is used from the CeedBasis, and number of elements is used from the CeedElemRestriction 94 PetscErrorCode CreateElemRestrColloc(Ceed ceed, CeedInt num_comp, CeedBasis basis, CeedElemRestriction elem_restr_base, 95 CeedElemRestriction *elem_restr_collocated, CeedVector *l_vec, CeedVector *e_vec) { 96 CeedInt num_elem_qpts, loc_num_elem; 97 PetscFunctionBeginUser; 98 99 CeedBasisGetNumQuadraturePoints(basis, &num_elem_qpts); 100 CeedElemRestrictionGetNumElements(elem_restr_base, &loc_num_elem); 101 102 const CeedInt strides[] = {num_comp, 1, num_elem_qpts * num_comp}; 103 CeedElemRestrictionCreateStrided(ceed, loc_num_elem, num_elem_qpts, num_comp, num_comp * loc_num_elem * num_elem_qpts, strides, 104 elem_restr_collocated); 105 CeedElemRestrictionCreateVector(*elem_restr_collocated, l_vec, e_vec); 106 PetscFunctionReturn(0); 107 } 108 109 PetscErrorCode GetQuadratureCoords(Ceed ceed, DM dm, CeedElemRestriction elem_restr_x, CeedBasis basis_x, CeedVector x_coords, CeedVector *qx_coords, 110 PetscInt *total_nqpnts) { 111 CeedQFunction qf_quad_coords; 112 CeedOperator op_quad_coords; 113 PetscInt num_comp_x, loc_num_elem, num_elem_qpts; 114 CeedElemRestriction elem_restr_qx; 115 PetscFunctionBeginUser; 116 117 // Create Element Restriction and CeedVector for quadrature coordinates 118 CeedBasisGetNumQuadraturePoints(basis_x, &num_elem_qpts); 119 CeedElemRestrictionGetNumElements(elem_restr_x, &loc_num_elem); 120 CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x); 121 *total_nqpnts = num_elem_qpts * loc_num_elem; 122 PetscCall(CreateElemRestrColloc(ceed, num_comp_x, basis_x, elem_restr_x, &elem_restr_qx, qx_coords, NULL)); 123 124 // Create QFunction 125 CeedQFunctionCreateIdentity(ceed, num_comp_x, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_quad_coords); 126 127 // Create Operator 128 CeedOperatorCreate(ceed, qf_quad_coords, NULL, NULL, &op_quad_coords); 129 CeedOperatorSetField(op_quad_coords, "input", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 130 CeedOperatorSetField(op_quad_coords, "output", elem_restr_qx, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 131 132 CeedOperatorApply(op_quad_coords, x_coords, *qx_coords, CEED_REQUEST_IMMEDIATE); 133 134 CeedQFunctionDestroy(&qf_quad_coords); 135 CeedOperatorDestroy(&op_quad_coords); 136 PetscFunctionReturn(0); 137 } 138 139 PetscErrorCode CreateStatsSF(Ceed ceed, CeedData ceed_data, DM parentdm, DM childdm, PetscSF statssf) { 140 PetscInt child_num_qpnts, parent_num_qpnts, num_comp_x; 141 CeedVector child_qx_coords, parent_qx_coords; 142 PetscReal *child_coords, *parent_coords; 143 PetscFunctionBeginUser; 144 145 // Assume that child and parent have the same number of components 146 CeedBasisGetNumComponents(ceed_data->basis_x, &num_comp_x); 147 const PetscInt num_comp_sf = num_comp_x - 1; // Number of coord components used in the creation of the SF 148 149 // Get quad_coords for child DM 150 PetscCall(GetQuadratureCoords(ceed, childdm, ceed_data->elem_restr_x, ceed_data->basis_xc, ceed_data->x_coord, &child_qx_coords, &child_num_qpnts)); 151 152 // Get quad_coords for parent DM 153 PetscCall(GetQuadratureCoords(ceed, parentdm, ceed_data->spanstats.elem_restr_parent_x, ceed_data->spanstats.basis_x, ceed_data->spanstats.x_coord, 154 &parent_qx_coords, &parent_num_qpnts)); 155 156 // Remove z component of coordinates for matching 157 { 158 const PetscReal *child_quad_coords, *parent_quad_coords; 159 160 CeedVectorGetArrayRead(child_qx_coords, CEED_MEM_HOST, &child_quad_coords); 161 CeedVectorGetArrayRead(parent_qx_coords, CEED_MEM_HOST, &parent_quad_coords); 162 163 PetscCall(PetscMalloc2(child_num_qpnts * 2, &child_coords, parent_num_qpnts * 2, &parent_coords)); 164 for (int i = 0; i < child_num_qpnts; i++) { 165 child_coords[0 + i * num_comp_sf] = child_quad_coords[0 + i * num_comp_x]; 166 child_coords[1 + i * num_comp_sf] = child_quad_coords[1 + i * num_comp_x]; 167 } 168 for (int i = 0; i < parent_num_qpnts; i++) { 169 parent_coords[0 + i * num_comp_sf] = parent_quad_coords[0 + i * num_comp_x]; 170 parent_coords[1 + i * num_comp_sf] = parent_quad_coords[1 + i * num_comp_x]; 171 } 172 CeedVectorRestoreArrayRead(child_qx_coords, &child_quad_coords); 173 CeedVectorRestoreArrayRead(parent_qx_coords, &parent_quad_coords); 174 } 175 176 // Only check the first two components of the coordinates 177 PetscCall(PetscSFSetGraphFromCoordinates(statssf, parent_num_qpnts, child_num_qpnts, num_comp_sf, 1e-12, parent_coords, child_coords)); 178 179 PetscCall(PetscFree2(child_coords, parent_coords)); 180 CeedVectorDestroy(&ceed_data->spanstats.x_coord); 181 CeedVectorDestroy(&child_qx_coords); 182 CeedVectorDestroy(&parent_qx_coords); 183 PetscFunctionReturn(0); 184 } 185 186 PetscErrorCode SetupStatsCollection(Ceed ceed, User user, CeedData ceed_data, ProblemData *problem) { 187 DM dm = user->spanstats.dm; 188 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 189 CeedInt dim, P, Q, num_comp_x; 190 Vec X_loc; 191 PetscMemType X_loc_memtype; 192 const PetscScalar *X_loc_array; 193 PetscFunctionBeginUser; 194 195 PetscCall(DMGetDimension(dm, &dim)); 196 CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_x, &Q); 197 CeedBasisGetNumNodes1D(ceed_data->basis_x, &P); 198 199 PetscCall(GetRestrictionForDomain(ceed, dm, 0, 0, 0, Q, problem->q_data_size_sur, &ceed_data->spanstats.elem_restr_parent_stats, 200 &ceed_data->spanstats.elem_restr_parent_x, &ceed_data->spanstats.elem_restr_parent_qd)); 201 CeedElemRestrictionGetNumComponents(ceed_data->spanstats.elem_restr_parent_x, &num_comp_x); 202 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_x, &ceed_data->spanstats.x_coord, NULL); 203 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_stats, &user->spanstats.stats_ceed, NULL); 204 CeedElemRestrictionCreateVector(ceed_data->spanstats.elem_restr_parent_qd, &ceed_data->spanstats.q_data, NULL); 205 206 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS_LOBATTO, &ceed_data->spanstats.basis_x); 207 CeedBasisCreateTensorH1Lagrange(ceed, dim, user->spanstats.num_comp_stats, P, Q, CEED_GAUSS, &ceed_data->spanstats.basis_stats); 208 209 PetscCall(CreateElemRestrColloc(ceed, user->spanstats.num_comp_stats, ceed_data->spanstats.basis_stats, 210 ceed_data->spanstats.elem_restr_parent_stats, &ceed_data->spanstats.elem_restr_parent_colloc, 211 &user->spanstats.parent_stats, NULL)); 212 PetscCall(CreateElemRestrColloc(ceed, user->spanstats.num_comp_stats, ceed_data->basis_q, ceed_data->elem_restr_q, 213 &ceed_data->spanstats.elem_restr_child_colloc, &user->spanstats.child_stats, NULL)); 214 215 // -- Copy DM coordinates into CeedVector 216 { 217 DM cdm; 218 PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 219 if (cdm) { 220 PetscCall(DMGetCellCoordinatesLocal(dm, &X_loc)); 221 } else { 222 PetscCall(DMGetCoordinatesLocal(dm, &X_loc)); 223 } 224 } 225 PetscCall(VecScale(X_loc, problem->dm_scale)); 226 PetscCall(VecGetArrayReadAndMemType(X_loc, &X_loc_array, &X_loc_memtype)); 227 CeedVectorSetArray(ceed_data->spanstats.x_coord, MemTypeP2C(X_loc_memtype), CEED_COPY_VALUES, (PetscScalar *)X_loc_array); 228 PetscCall(VecRestoreArrayRead(X_loc, &X_loc_array)); 229 230 // Create SF for communicating child data back their respective parents 231 PetscCall(PetscSFCreate(comm, &user->spanstats.sf)); 232 PetscCall(CreateStatsSF(ceed, ceed_data, user->dm, user->spanstats.dm, user->spanstats.sf)); 233 234 PetscFunctionReturn(0); 235 } 236