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