1 // Copyright (c) 2017-2024, 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 <ceed.h> 12 #include <petscdmplex.h> 13 #include <petscds.h> 14 15 #include "../navierstokes.h" 16 #include "../problems/stg_shur14.h" 17 18 // Create mesh 19 PetscErrorCode CreateDM(MPI_Comm comm, ProblemData problem, MatType mat_type, VecType vec_type, DM *dm) { 20 PetscFunctionBeginUser; 21 // Create DMPLEX 22 PetscCall(DMCreate(comm, dm)); 23 PetscCall(DMSetType(*dm, DMPLEX)); 24 { 25 PetscBool skip = PETSC_TRUE; 26 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_mat_preallocate_skip", &skip, NULL)); 27 PetscCall(DMSetMatrixPreallocateSkip(*dm, skip)); 28 } 29 PetscCall(DMSetMatType(*dm, mat_type)); 30 PetscCall(DMSetVecType(*dm, vec_type)); 31 32 // Set Tensor elements 33 PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_simplex", "0")); 34 PetscCall(PetscOptionsSetValue(NULL, "-dm_sparse_localize", "0")); 35 PetscCall(PetscOptionsSetValue(NULL, "-dm_localize", "0")); // Localization done in DMSetupByOrderEnd_FEM 36 PetscCall(PetscOptionsSetValue(NULL, "-dm_blocking_type", "field_node")); 37 38 // Set CL options 39 PetscCall(DMSetFromOptions(*dm)); 40 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 41 PetscFunctionReturn(PETSC_SUCCESS); 42 } 43 44 // Setup DM 45 PetscErrorCode SetUpDM(DM dm, ProblemData problem, PetscInt degree, PetscInt q_extra, SimpleBC bc, Physics phys) { 46 PetscInt num_comp_q = 5; 47 PetscFunctionBeginUser; 48 49 PetscCall(DMSetupByOrderBegin_FEM(PETSC_TRUE, PETSC_TRUE, degree, PETSC_DECIDE, q_extra, 1, &num_comp_q, dm)); 50 51 { // Add strong boundary conditions to DM 52 DMLabel label; 53 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 54 PetscCall(DMPlexLabelComplete(dm, label)); 55 // Set wall BCs 56 if (bc->num_wall > 0) { 57 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, bc->num_wall, bc->walls, 0, bc->num_comps, bc->wall_comps, NULL, NULL, NULL, NULL)); 58 } 59 // Set symmetry BCs in the x direction 60 if (bc->num_symmetry[0] > 0) { 61 PetscInt comps[1] = {1}; 62 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "symmetry_x", label, bc->num_symmetry[0], bc->symmetries[0], 0, 1, comps, NULL, NULL, NULL, NULL)); 63 } 64 // Set symmetry BCs in the y direction 65 if (bc->num_symmetry[1] > 0) { 66 PetscInt comps[1] = {2}; 67 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "symmetry_y", label, bc->num_symmetry[1], bc->symmetries[1], 0, 1, comps, NULL, NULL, NULL, NULL)); 68 } 69 // Set symmetry BCs in the z direction 70 if (bc->num_symmetry[2] > 0) { 71 PetscInt comps[1] = {3}; 72 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "symmetry_z", label, bc->num_symmetry[2], bc->symmetries[2], 0, 1, comps, NULL, NULL, NULL, NULL)); 73 } 74 { 75 PetscBool use_strongstg = PETSC_FALSE; 76 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); 77 if (use_strongstg) PetscCall(SetupStrongStg(dm, bc, problem, phys)); 78 } 79 } 80 81 PetscCall(DMSetupByOrderEnd_FEM(PETSC_TRUE, dm)); 82 83 // Empty name for conserved field (because there is only one field) 84 PetscSection section; 85 PetscCall(DMGetLocalSection(dm, §ion)); 86 PetscCall(PetscSectionSetFieldName(section, 0, "")); 87 switch (phys->state_var) { 88 case STATEVAR_CONSERVATIVE: 89 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Density")); 90 PetscCall(PetscSectionSetComponentName(section, 0, 1, "MomentumX")); 91 PetscCall(PetscSectionSetComponentName(section, 0, 2, "MomentumY")); 92 PetscCall(PetscSectionSetComponentName(section, 0, 3, "MomentumZ")); 93 PetscCall(PetscSectionSetComponentName(section, 0, 4, "TotalEnergy")); 94 break; 95 96 case STATEVAR_PRIMITIVE: 97 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Pressure")); 98 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityX")); 99 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityY")); 100 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityZ")); 101 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Temperature")); 102 break; 103 104 case STATEVAR_ENTROPY: 105 PetscCall(PetscSectionSetComponentName(section, 0, 0, "EntropyDensity")); 106 PetscCall(PetscSectionSetComponentName(section, 0, 1, "EntropyMomentumX")); 107 PetscCall(PetscSectionSetComponentName(section, 0, 2, "EntropyMomentumY")); 108 PetscCall(PetscSectionSetComponentName(section, 0, 3, "EntropyMomentumZ")); 109 PetscCall(PetscSectionSetComponentName(section, 0, 4, "EntropyTotalEnergy")); 110 break; 111 } 112 PetscFunctionReturn(PETSC_SUCCESS); 113 } 114 115 // Refine DM for high-order viz 116 PetscErrorCode VizRefineDM(DM dm, User user, ProblemData problem, SimpleBC bc, Physics phys) { 117 DM dm_hierarchy[user->app_ctx->viz_refine + 1]; 118 VecType vec_type; 119 120 PetscFunctionBeginUser; 121 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 122 123 dm_hierarchy[0] = dm; 124 for (PetscInt i = 0, d = user->app_ctx->degree; i < user->app_ctx->viz_refine; i++) { 125 Mat interp_next; 126 PetscCall(DMRefine(dm_hierarchy[i], MPI_COMM_NULL, &dm_hierarchy[i + 1])); 127 PetscCall(DMClearDS(dm_hierarchy[i + 1])); 128 PetscCall(DMClearFields(dm_hierarchy[i + 1])); 129 PetscCall(DMSetCoarseDM(dm_hierarchy[i + 1], dm_hierarchy[i])); 130 d = (d + 1) / 2; 131 PetscInt q_order = d + user->app_ctx->q_extra; 132 if (i + 1 == user->app_ctx->viz_refine) d = 1; 133 PetscCall(DMGetVecType(dm, &vec_type)); 134 PetscCall(DMSetVecType(dm_hierarchy[i + 1], vec_type)); 135 PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, q_order, bc, phys)); 136 PetscCall(DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i + 1], &interp_next, NULL)); 137 if (!i) user->interp_viz = interp_next; 138 else { 139 Mat C; 140 PetscCall(MatMatMult(interp_next, user->interp_viz, MAT_INITIAL_MATRIX, PETSC_DECIDE, &C)); 141 PetscCall(MatDestroy(&interp_next)); 142 PetscCall(MatDestroy(&user->interp_viz)); 143 user->interp_viz = C; 144 } 145 } 146 for (PetscInt i = 1; i < user->app_ctx->viz_refine; i++) { 147 PetscCall(DMDestroy(&dm_hierarchy[i])); 148 } 149 user->dm_viz = dm_hierarchy[user->app_ctx->viz_refine]; 150 PetscFunctionReturn(PETSC_SUCCESS); 151 } 152