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