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 /// Setup DM for Navier-Stokes example using PETSc 10 11 #include <ceed.h> 12 #include <petscdmplex.h> 13 #include <petscds.h> 14 15 #include "../navierstokes.h" 16 #include "../problems/stg_shur14.h" 17 18 // Create mesh 19 PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, MatType mat_type, VecType vec_type, DM *dm) { 20 PetscFunctionBeginUser; 21 // Create DMPLEX 22 PetscCall(DMCreate(comm, dm)); 23 PetscCall(DMSetType(*dm, DMPLEX)); 24 { 25 PetscBool skip = PETSC_TRUE; 26 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_mat_preallocate_skip", &skip, NULL)); 27 PetscCall(DMSetMatrixPreallocateSkip(*dm, skip)); 28 } 29 PetscCall(DMSetMatType(*dm, mat_type)); 30 PetscCall(DMSetVecType(*dm, vec_type)); 31 32 // Set Tensor elements 33 PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_simplex", "0")); 34 PetscCall(PetscOptionsSetValue(NULL, "-dm_sparse_localize", "0")); 35 // Set CL options 36 PetscCall(DMSetFromOptions(*dm)); 37 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 38 PetscFunctionReturn(PETSC_SUCCESS); 39 } 40 41 // Setup DM 42 PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_extra, SimpleBC bc, Physics phys) { 43 PetscInt q_order = degree + q_extra; 44 PetscFunctionBeginUser; 45 { 46 PetscBool is_simplex = PETSC_TRUE; 47 48 // Check if simplex or tensor-product mesh 49 PetscCall(DMPlexIsSimplex(dm, &is_simplex)); 50 // Configure the finite element space and boundary conditions 51 PetscFE fe; 52 PetscInt num_comp_q = 5; 53 DMLabel label; 54 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, is_simplex, degree, q_order, &fe)); 55 PetscCall(PetscObjectSetName((PetscObject)fe, "Q")); 56 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 57 PetscCall(DMCreateDS(dm)); 58 { // Project coordinates to enrich quadrature space 59 DM dm_coord; 60 PetscDS ds_coord; 61 PetscFE fe_coord_current, fe_coord_new, fe_coord_face_new; 62 PetscDualSpace fe_coord_dual_space; 63 PetscInt fe_coord_order, num_comp_coord; 64 65 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 66 PetscCall(DMGetCoordinateDim(dm, &num_comp_coord)); 67 PetscCall(DMGetRegionDS(dm_coord, NULL, NULL, &ds_coord, NULL)); 68 PetscCall(PetscDSGetDiscretization(ds_coord, 0, (PetscObject *)&fe_coord_current)); 69 PetscCall(PetscFEGetDualSpace(fe_coord_current, &fe_coord_dual_space)); 70 PetscCall(PetscDualSpaceGetOrder(fe_coord_dual_space, &fe_coord_order)); 71 72 // Create FE for coordinates 73 PetscCheck(fe_coord_order == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER_INPUT, 74 "Only linear mesh geometry supported. Recieved %d order", fe_coord_order); 75 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_coord, is_simplex, fe_coord_order, q_order, &fe_coord_new)); 76 PetscCall(PetscFEGetHeightSubspace(fe_coord_new, 1, &fe_coord_face_new)); 77 PetscCall(DMProjectCoordinates(dm, fe_coord_new)); 78 PetscCall(PetscFEDestroy(&fe_coord_new)); 79 } 80 81 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 82 PetscCall(DMPlexLabelComplete(dm, label)); 83 // Set wall BCs 84 if (bc->num_wall > 0) { 85 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, bc->num_wall, bc->walls, 0, bc->num_comps, bc->wall_comps, NULL, NULL, NULL, NULL)); 86 } 87 // Set slip BCs in the x direction 88 if (bc->num_slip[0] > 0) { 89 PetscInt comps[1] = {1}; 90 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", label, bc->num_slip[0], bc->slips[0], 0, 1, comps, NULL, NULL, NULL, NULL)); 91 } 92 // Set slip BCs in the y direction 93 if (bc->num_slip[1] > 0) { 94 PetscInt comps[1] = {2}; 95 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", label, bc->num_slip[1], bc->slips[1], 0, 1, comps, NULL, NULL, NULL, NULL)); 96 } 97 // Set slip BCs in the z direction 98 if (bc->num_slip[2] > 0) { 99 PetscInt comps[1] = {3}; 100 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", label, bc->num_slip[2], bc->slips[2], 0, 1, comps, NULL, NULL, NULL, NULL)); 101 } 102 { 103 PetscBool use_strongstg = PETSC_FALSE; 104 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); 105 if (use_strongstg) PetscCall(SetupStrongSTG(dm, bc, problem, phys)); 106 } 107 108 if (!is_simplex) { 109 DM dm_coord; 110 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 111 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 112 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); 113 } 114 115 PetscCall(PetscFEDestroy(&fe)); 116 } 117 118 // Empty name for conserved field (because there is only one field) 119 PetscSection section; 120 PetscCall(DMGetLocalSection(dm, §ion)); 121 PetscCall(PetscSectionSetFieldName(section, 0, "")); 122 switch (phys->state_var) { 123 case STATEVAR_CONSERVATIVE: 124 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Density")); 125 PetscCall(PetscSectionSetComponentName(section, 0, 1, "MomentumX")); 126 PetscCall(PetscSectionSetComponentName(section, 0, 2, "MomentumY")); 127 PetscCall(PetscSectionSetComponentName(section, 0, 3, "MomentumZ")); 128 PetscCall(PetscSectionSetComponentName(section, 0, 4, "TotalEnergy")); 129 break; 130 131 case STATEVAR_PRIMITIVE: 132 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Pressure")); 133 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityX")); 134 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityY")); 135 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityZ")); 136 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Temperature")); 137 break; 138 } 139 PetscFunctionReturn(PETSC_SUCCESS); 140 } 141 142 // Refine DM for high-order viz 143 PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, SimpleBC bc, Physics phys) { 144 DM dm_hierarchy[user->app_ctx->viz_refine + 1]; 145 VecType vec_type; 146 PetscFunctionBeginUser; 147 148 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 149 150 dm_hierarchy[0] = dm; 151 for (PetscInt i = 0, d = user->app_ctx->degree; i < user->app_ctx->viz_refine; i++) { 152 Mat interp_next; 153 PetscCall(DMRefine(dm_hierarchy[i], MPI_COMM_NULL, &dm_hierarchy[i + 1])); 154 PetscCall(DMClearDS(dm_hierarchy[i + 1])); 155 PetscCall(DMClearFields(dm_hierarchy[i + 1])); 156 PetscCall(DMSetCoarseDM(dm_hierarchy[i + 1], dm_hierarchy[i])); 157 d = (d + 1) / 2; 158 PetscInt q_order = d + user->app_ctx->q_extra; 159 if (i + 1 == user->app_ctx->viz_refine) d = 1; 160 PetscCall(DMGetVecType(dm, &vec_type)); 161 PetscCall(DMSetVecType(dm_hierarchy[i + 1], vec_type)); 162 PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, q_order, bc, phys)); 163 PetscCall(DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i + 1], &interp_next, NULL)); 164 if (!i) user->interp_viz = interp_next; 165 else { 166 Mat C; 167 PetscCall(MatMatMult(interp_next, user->interp_viz, MAT_INITIAL_MATRIX, PETSC_DECIDE, &C)); 168 PetscCall(MatDestroy(&interp_next)); 169 PetscCall(MatDestroy(&user->interp_viz)); 170 user->interp_viz = C; 171 } 172 } 173 for (PetscInt i = 1; i < user->app_ctx->viz_refine; i++) { 174 PetscCall(DMDestroy(&dm_hierarchy[i])); 175 } 176 user->dm_viz = dm_hierarchy[user->app_ctx->viz_refine]; 177 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180