1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Setup DM for Navier-Stokes example using PETSc 19 20 #include "../navierstokes.h" 21 22 // Create mesh 23 PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, DM *dm) { 24 PetscErrorCode ierr; 25 PetscFunctionBeginUser; 26 // Create DMPLEX 27 ierr = DMCreate(comm, dm); CHKERRQ(ierr); 28 ierr = DMSetType(*dm, DMPLEX); CHKERRQ(ierr); 29 // Set Tensor elements 30 ierr = PetscOptionsSetValue(NULL, "-dm_plex_simplex", "0"); CHKERRQ(ierr); 31 // Set CL options 32 ierr = DMSetFromOptions(*dm); CHKERRQ(ierr); 33 ierr = DMViewFromOptions(*dm, NULL, "-dm_view"); CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 // Setup DM 38 PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, 39 SimpleBC bc, Physics phys, void *setup_ctx) { 40 PetscErrorCode ierr; 41 PetscFunctionBeginUser; 42 { 43 // Configure the finite element space and boundary conditions 44 PetscFE fe; 45 PetscInt num_comp_q = 5; 46 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, 47 PETSC_FALSE, degree, PETSC_DECIDE, 48 &fe); CHKERRQ(ierr); 49 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 50 ierr = DMAddField(dm, NULL,(PetscObject)fe); CHKERRQ(ierr); 51 ierr = DMCreateDS(dm); CHKERRQ(ierr); 52 { 53 /* create FE field for coordinates */ 54 PetscFE fe_coords; 55 PetscInt num_comp_coord; 56 ierr = DMGetCoordinateDim(dm, &num_comp_coord); CHKERRQ(ierr); 57 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_coord, 58 PETSC_FALSE, 1, 1, &fe_coords); CHKERRQ(ierr); 59 ierr = DMProjectCoordinates(dm, fe_coords); CHKERRQ(ierr); 60 ierr = PetscFEDestroy(&fe_coords); CHKERRQ(ierr); 61 } 62 // Set wall BCs 63 if (bc->num_wall > 0) { 64 DMLabel label; 65 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 66 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 67 bc->num_wall, bc->walls, 0, bc->num_comps, 68 bc->wall_comps, (void(*)(void))problem->bc, 69 NULL, setup_ctx, NULL); CHKERRQ(ierr); 70 } 71 // Set slip BCs in the x direction 72 if (bc->num_slip[0] > 0) { 73 DMLabel label; 74 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 75 PetscInt comps[1] = {1}; 76 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", label, 77 bc->num_slip[0], bc->slips[0], 0, 1, comps, 78 (void(*)(void))NULL, NULL, setup_ctx, NULL); CHKERRQ(ierr); 79 } 80 // Set slip BCs in the y direction 81 if (bc->num_slip[1] > 0) { 82 DMLabel label; 83 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 84 PetscInt comps[1] = {2}; 85 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", label, 86 bc->num_slip[1], bc->slips[1], 0, 1, comps, 87 (void(*)(void))NULL, NULL, setup_ctx, NULL); CHKERRQ(ierr); 88 } 89 // Set slip BCs in the z direction 90 if (bc->num_slip[2] > 0) { 91 DMLabel label; 92 ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); 93 PetscInt comps[1] = {3}; 94 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", label, 95 bc->num_slip[2], bc->slips[2], 0, 1, comps, 96 (void(*)(void))NULL, NULL, setup_ctx, NULL); CHKERRQ(ierr); 97 } 98 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 99 CHKERRQ(ierr); 100 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 101 } 102 { 103 // Empty name for conserved field (because there is only one field) 104 PetscSection section; 105 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 106 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 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 } 118 PetscFunctionReturn(0); 119 } 120 121 // Refine DM for high-order viz 122 PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, 123 SimpleBC bc, Physics phys, void *setup_ctx) { 124 PetscErrorCode ierr; 125 DM dm_hierarchy[user->app_ctx->viz_refine + 1]; 126 VecType vec_type; 127 PetscFunctionBeginUser; 128 129 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 130 131 dm_hierarchy[0] = dm; 132 for (PetscInt i = 0, d = user->app_ctx->degree; 133 i < user->app_ctx->viz_refine; i++) { 134 Mat interp_next; 135 ierr = DMRefine(dm_hierarchy[i], MPI_COMM_NULL, &dm_hierarchy[i+1]); 136 CHKERRQ(ierr); 137 ierr = DMClearDS(dm_hierarchy[i+1]); CHKERRQ(ierr); 138 ierr = DMClearFields(dm_hierarchy[i+1]); CHKERRQ(ierr); 139 ierr = DMSetCoarseDM(dm_hierarchy[i+1], dm_hierarchy[i]); CHKERRQ(ierr); 140 d = (d + 1) / 2; 141 if (i + 1 == user->app_ctx->viz_refine) d = 1; 142 ierr = DMGetVecType(dm, &vec_type); CHKERRQ(ierr); 143 ierr = DMSetVecType(dm_hierarchy[i+1], vec_type); CHKERRQ(ierr); 144 ierr = SetUpDM(dm_hierarchy[i+1], problem, d, bc, phys, setup_ctx); 145 CHKERRQ(ierr); 146 ierr = DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i+1], &interp_next, 147 NULL); CHKERRQ(ierr); 148 if (!i) user->interp_viz = interp_next; 149 else { 150 Mat C; 151 ierr = MatMatMult(interp_next, user->interp_viz, MAT_INITIAL_MATRIX, 152 PETSC_DECIDE, &C); CHKERRQ(ierr); 153 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 154 ierr = MatDestroy(&user->interp_viz); CHKERRQ(ierr); 155 user->interp_viz = C; 156 } 157 } 158 for (PetscInt i=1; i<user->app_ctx->viz_refine; i++) { 159 ierr = DMDestroy(&dm_hierarchy[i]); CHKERRQ(ierr); 160 } 161 user->dm_viz = dm_hierarchy[user->app_ctx->viz_refine]; 162 163 PetscFunctionReturn(0); 164 } 165