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