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 HONEE 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 PetscBool isCGNS; 17 char filename[PETSC_MAX_PATH_LEN] = ""; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscOptionsGetString(NULL, NULL, "-dm_plex_filename", filename, sizeof(filename), NULL)); 21 PetscCall(HoneeCheckFilenameExtension(comm, filename, ".cgns", &isCGNS)); 22 23 if (isCGNS) { 24 // Must create from file to keep the sfNatural field in tact 25 PetscCall(DMPlexCreateFromFile(comm, filename, "HONEE", PETSC_TRUE, dm)); 26 PetscCall(PetscOptionsSetValue(NULL, "-dm_plex_filename", "")); 27 } else { 28 PetscCall(DMCreate(comm, dm)); 29 PetscCall(DMSetType(*dm, DMPLEX)); 30 } 31 32 { 33 PetscBool skip = PETSC_TRUE; 34 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_mat_preallocate_skip", &skip, NULL)); 35 PetscCall(DMSetMatrixPreallocateSkip(*dm, skip)); 36 } 37 PetscCall(DMSetMatType(*dm, mat_type)); 38 PetscCall(DMSetVecType(*dm, vec_type)); 39 PetscCall(HoneeOptionsSetValueDefault(NULL, "-dm_plex_simplex", "0")); 40 PetscCall(HoneeOptionsSetValueDefault(NULL, "-dm_localize", "0")); // Delay localization until DMSetupByOrderEnd_FEM 41 { // Use Mat graph generation for partitioning. See https://gitlab.com/petsc/petsc/-/merge_requests/8336 42 PetscBool distribute = PETSC_FALSE; 43 PetscMPIInt num_ranks; 44 PetscCall(DMPlexDistributeGetDefault(*dm, &distribute)); 45 PetscCallMPI(MPI_Comm_size(comm, &num_ranks)); 46 if (distribute && num_ranks > 1) PetscCall(HoneeOptionsSetValueDefault(NULL, "-dm_plex_csr_alg", "mat")); 47 } 48 PetscCall(DMSetSparseLocalize(*dm, PETSC_FALSE)); 49 PetscCall(DMSetBlockingType(*dm, DM_BLOCKING_FIELD_NODE)); 50 PetscCall(DMPlexSetPartitionBalance(*dm, PETSC_TRUE)); 51 52 // Set CL options 53 PetscCall(DMSetFromOptions(*dm)); 54 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 55 { // Force isoperiodic point SF to be created to update sfNatural. 56 // Needs to be done before removing the field corresponding to sfNatural, but after IsoperiodicFaces have been set 57 PetscInt num_fields; 58 PetscCall(DMGetNumFields(*dm, &num_fields)); 59 if (num_fields) { 60 PetscSection dummy; 61 PetscCall(DMGetGlobalSection(*dm, &dummy)); 62 } 63 } 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 // Setup DM 68 PetscErrorCode SetUpDM(DM dm, ProblemData problem, PetscInt degree, PetscInt q_extra, Physics phys) { 69 PetscInt num_comp_q = 5; 70 71 PetscFunctionBeginUser; 72 PetscCall(DMClearFields(dm)); 73 PetscCall(DMSetLocalSection(dm, NULL)); 74 PetscCall(DMSetupByOrderBegin_FEM(PETSC_TRUE, PETSC_TRUE, degree, PETSC_DECIDE, q_extra, 1, &num_comp_q, dm)); 75 76 { 77 PetscBool use_strongstg = PETSC_FALSE; 78 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); 79 if (use_strongstg) PetscCall(SetupStrongStg(dm, problem, phys)); 80 } 81 82 { // Add strong boundary conditions to DM 83 DMLabel label; 84 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 85 PetscCheck(label, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM does not have 'Face Sets' label"); 86 PetscCall(DMPlexLabelComplete(dm, label)); 87 88 for (PetscInt i = 0; i < problem->num_bc_defs; i++) { 89 BCDefinition bc_def = problem->bc_defs[i]; 90 PetscInt num_essential_comps, num_label_values; 91 const PetscInt *essential_comps, *label_values; 92 const char *name; 93 94 PetscCall(BCDefinitionSetDM(problem->bc_defs[i], dm)); 95 PetscCall(BCDefinitionGetEssential(bc_def, &num_essential_comps, &essential_comps)); 96 if (essential_comps > 0) { 97 PetscCall(BCDefinitionGetInfo(bc_def, &name, &num_label_values, &label_values)); 98 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, name, label, num_label_values, label_values, 0, num_essential_comps, essential_comps, NULL, NULL, 99 NULL, NULL)); 100 } 101 } 102 } 103 104 PetscCall(DMSetupByOrderEnd_FEM(PETSC_TRUE, dm)); 105 106 // Empty name for conserved field (because there is only one field) 107 PetscSection section; 108 PetscCall(DMGetLocalSection(dm, §ion)); 109 PetscCall(PetscSectionSetFieldName(section, 0, "")); 110 switch (phys->state_var) { 111 case STATEVAR_CONSERVATIVE: 112 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Density")); 113 PetscCall(PetscSectionSetComponentName(section, 0, 1, "MomentumX")); 114 PetscCall(PetscSectionSetComponentName(section, 0, 2, "MomentumY")); 115 PetscCall(PetscSectionSetComponentName(section, 0, 3, "MomentumZ")); 116 PetscCall(PetscSectionSetComponentName(section, 0, 4, "TotalEnergy")); 117 break; 118 119 case STATEVAR_PRIMITIVE: 120 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Pressure")); 121 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityX")); 122 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityY")); 123 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityZ")); 124 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Temperature")); 125 break; 126 127 case STATEVAR_ENTROPY: 128 PetscCall(PetscSectionSetComponentName(section, 0, 0, "EntropyDensity")); 129 PetscCall(PetscSectionSetComponentName(section, 0, 1, "EntropyMomentumX")); 130 PetscCall(PetscSectionSetComponentName(section, 0, 2, "EntropyMomentumY")); 131 PetscCall(PetscSectionSetComponentName(section, 0, 3, "EntropyMomentumZ")); 132 PetscCall(PetscSectionSetComponentName(section, 0, 4, "EntropyTotalEnergy")); 133 break; 134 } 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 // Refine DM for high-order viz 139 PetscErrorCode VizRefineDM(DM dm, Honee honee, ProblemData problem, Physics phys) { 140 DM dm_hierarchy[honee->app_ctx->viz_refine + 1]; 141 VecType vec_type; 142 143 PetscFunctionBeginUser; 144 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE)); 145 146 dm_hierarchy[0] = dm; 147 for (PetscInt i = 0, d = honee->app_ctx->degree; i < honee->app_ctx->viz_refine; i++) { 148 Mat interp_next; 149 PetscCall(DMRefine(dm_hierarchy[i], MPI_COMM_NULL, &dm_hierarchy[i + 1])); 150 PetscCall(DMClearDS(dm_hierarchy[i + 1])); 151 PetscCall(DMClearFields(dm_hierarchy[i + 1])); 152 PetscCall(DMSetCoarseDM(dm_hierarchy[i + 1], dm_hierarchy[i])); 153 d = (d + 1) / 2; 154 PetscInt q_order = d + honee->app_ctx->q_extra; 155 if (i + 1 == honee->app_ctx->viz_refine) d = 1; 156 PetscCall(DMGetVecType(dm, &vec_type)); 157 PetscCall(DMSetVecType(dm_hierarchy[i + 1], vec_type)); 158 PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, q_order, phys)); 159 PetscCall(DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i + 1], &interp_next, NULL)); 160 if (!i) honee->interp_viz = interp_next; 161 else { 162 Mat C; 163 PetscCall(MatMatMult(interp_next, honee->interp_viz, MAT_INITIAL_MATRIX, PETSC_DECIDE, &C)); 164 PetscCall(MatDestroy(&interp_next)); 165 PetscCall(MatDestroy(&honee->interp_viz)); 166 honee->interp_viz = C; 167 } 168 } 169 for (PetscInt i = 1; i < honee->app_ctx->viz_refine; i++) { 170 PetscCall(DMDestroy(&dm_hierarchy[i])); 171 } 172 honee->dm_viz = dm_hierarchy[honee->app_ctx->viz_refine]; 173 PetscFunctionReturn(PETSC_SUCCESS); 174 } 175