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