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