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