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 14 #include "../navierstokes.h" 15 #include "../problems/stg_shur14.h" 16 17 // Create mesh 18 PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, MatType mat_type, VecType vec_type, DM *dm) { 19 PetscFunctionBeginUser; 20 // Create DMPLEX 21 PetscCall(DMCreate(comm, dm)); 22 PetscCall(DMSetType(*dm, DMPLEX)); 23 { 24 PetscBool skip = PETSC_TRUE; 25 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_mat_preallocate_skip", &skip, NULL)); 26 PetscCall(DMSetMatrixPreallocateSkip(*dm, skip)); 27 } 28 PetscCall(DMSetMatType(*dm, mat_type)); 29 PetscCall(DMSetVecType(*dm, vec_type)); 30 31 // Set Tensor elements 32 PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_simplex", "0")); 33 PetscCall(PetscOptionsSetValue(NULL, "-dm_sparse_localize", "0")); 34 // Set CL options 35 PetscCall(DMSetFromOptions(*dm)); 36 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 37 PetscFunctionReturn(0); 38 } 39 40 // Setup DM 41 PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc, Physics phys) { 42 PetscFunctionBeginUser; 43 { 44 // Configure the finite element space and boundary conditions 45 PetscFE fe; 46 PetscInt num_comp_q = 5; 47 DMLabel label; 48 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); 49 PetscCall(PetscObjectSetName((PetscObject)fe, "Q")); 50 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 51 PetscCall(DMCreateDS(dm)); 52 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 53 PetscCall(DMPlexLabelComplete(dm, label)); 54 // Set wall BCs 55 if (bc->num_wall > 0) { 56 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, bc->num_wall, bc->walls, 0, bc->num_comps, bc->wall_comps, 57 (void (*)(void))problem->bc, NULL, problem->bc_ctx, NULL)); 58 } 59 // Set slip BCs in the x direction 60 if (bc->num_slip[0] > 0) { 61 PetscInt comps[1] = {1}; 62 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", label, bc->num_slip[0], bc->slips[0], 0, 1, comps, (void (*)(void))NULL, NULL, 63 problem->bc_ctx, NULL)); 64 } 65 // Set slip BCs in the y direction 66 if (bc->num_slip[1] > 0) { 67 PetscInt comps[1] = {2}; 68 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", label, bc->num_slip[1], bc->slips[1], 0, 1, comps, (void (*)(void))NULL, NULL, 69 problem->bc_ctx, NULL)); 70 } 71 // Set slip BCs in the z direction 72 if (bc->num_slip[2] > 0) { 73 PetscInt comps[1] = {3}; 74 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", label, bc->num_slip[2], bc->slips[2], 0, 1, comps, (void (*)(void))NULL, NULL, 75 problem->bc_ctx, NULL)); 76 } 77 { 78 PetscBool use_strongstg = PETSC_FALSE; 79 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); 80 81 if (use_strongstg) { 82 PetscCall(SetupStrongSTG(dm, bc, problem, phys)); 83 } 84 } 85 86 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL)); 87 PetscCall(PetscFEDestroy(&fe)); 88 } 89 90 // Empty name for conserved field (because there is only one field) 91 PetscSection section; 92 PetscCall(DMGetLocalSection(dm, §ion)); 93 PetscCall(PetscSectionSetFieldName(section, 0, "")); 94 switch (phys->state_var) { 95 case STATEVAR_CONSERVATIVE: 96 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Density")); 97 PetscCall(PetscSectionSetComponentName(section, 0, 1, "Momentum X")); 98 PetscCall(PetscSectionSetComponentName(section, 0, 2, "Momentum Y")); 99 PetscCall(PetscSectionSetComponentName(section, 0, 3, "Momentum Z")); 100 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Energy Density")); 101 break; 102 103 case STATEVAR_PRIMITIVE: 104 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Pressure")); 105 PetscCall(PetscSectionSetComponentName(section, 0, 1, "Velocity X")); 106 PetscCall(PetscSectionSetComponentName(section, 0, 2, "Velocity Y")); 107 PetscCall(PetscSectionSetComponentName(section, 0, 3, "Velocity Z")); 108 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Temperature")); 109 break; 110 } 111 PetscFunctionReturn(0); 112 } 113 114 // Refine DM for high-order viz 115 PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, SimpleBC bc, Physics phys) { 116 DM dm_hierarchy[user->app_ctx->viz_refine + 1]; 117 VecType vec_type; 118 PetscFunctionBeginUser; 119 120 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 121 122 dm_hierarchy[0] = dm; 123 for (PetscInt i = 0, d = user->app_ctx->degree; i < user->app_ctx->viz_refine; i++) { 124 Mat interp_next; 125 PetscCall(DMRefine(dm_hierarchy[i], MPI_COMM_NULL, &dm_hierarchy[i + 1])); 126 PetscCall(DMClearDS(dm_hierarchy[i + 1])); 127 PetscCall(DMClearFields(dm_hierarchy[i + 1])); 128 PetscCall(DMSetCoarseDM(dm_hierarchy[i + 1], dm_hierarchy[i])); 129 d = (d + 1) / 2; 130 if (i + 1 == user->app_ctx->viz_refine) d = 1; 131 PetscCall(DMGetVecType(dm, &vec_type)); 132 PetscCall(DMSetVecType(dm_hierarchy[i + 1], vec_type)); 133 PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, bc, phys)); 134 PetscCall(DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i + 1], &interp_next, NULL)); 135 if (!i) user->interp_viz = interp_next; 136 else { 137 Mat C; 138 PetscCall(MatMatMult(interp_next, user->interp_viz, MAT_INITIAL_MATRIX, PETSC_DECIDE, &C)); 139 PetscCall(MatDestroy(&interp_next)); 140 PetscCall(MatDestroy(&user->interp_viz)); 141 user->interp_viz = C; 142 } 143 } 144 for (PetscInt i = 1; i < user->app_ctx->viz_refine; i++) { 145 PetscCall(DMDestroy(&dm_hierarchy[i])); 146 } 147 user->dm_viz = dm_hierarchy[user->app_ctx->viz_refine]; 148 149 PetscFunctionReturn(0); 150 } 151