1 // Copyright (c) 2017-2026, 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
CreateDM(MPI_Comm comm,ProblemData problem,MatType mat_type,VecType vec_type,DM * dm)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
SetUpDM(DM dm,ProblemData problem,PetscInt degree,PetscInt q_extra,SimpleBC bc,Physics phys)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
56 for (PetscInt i = 0; i < problem->num_bc_defs; i++) {
57 BCDefinition bc_def = problem->bc_defs[i];
58 PetscInt num_essential_comps, num_label_values;
59 const PetscInt *essential_comps, *label_values;
60 const char *name;
61
62 PetscCall(BCDefinitionGetEssential(bc_def, &num_essential_comps, &essential_comps));
63 if (essential_comps > 0) {
64 PetscCall(BCDefinitionGetInfo(bc_def, &name, &num_label_values, &label_values));
65 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, name, label, num_label_values, label_values, 0, num_essential_comps, essential_comps, NULL, NULL,
66 NULL, NULL));
67 }
68 }
69 {
70 PetscBool use_strongstg = PETSC_FALSE;
71 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL));
72 if (use_strongstg) PetscCall(SetupStrongStg(dm, bc, problem, phys));
73 }
74 }
75
76 PetscCall(DMSetupByOrderEnd_FEM(PETSC_TRUE, dm));
77
78 // Empty name for conserved field (because there is only one field)
79 PetscSection section;
80 PetscCall(DMGetLocalSection(dm, §ion));
81 PetscCall(PetscSectionSetFieldName(section, 0, ""));
82 switch (phys->state_var) {
83 case STATEVAR_CONSERVATIVE:
84 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Density"));
85 PetscCall(PetscSectionSetComponentName(section, 0, 1, "MomentumX"));
86 PetscCall(PetscSectionSetComponentName(section, 0, 2, "MomentumY"));
87 PetscCall(PetscSectionSetComponentName(section, 0, 3, "MomentumZ"));
88 PetscCall(PetscSectionSetComponentName(section, 0, 4, "TotalEnergy"));
89 break;
90
91 case STATEVAR_PRIMITIVE:
92 PetscCall(PetscSectionSetComponentName(section, 0, 0, "Pressure"));
93 PetscCall(PetscSectionSetComponentName(section, 0, 1, "VelocityX"));
94 PetscCall(PetscSectionSetComponentName(section, 0, 2, "VelocityY"));
95 PetscCall(PetscSectionSetComponentName(section, 0, 3, "VelocityZ"));
96 PetscCall(PetscSectionSetComponentName(section, 0, 4, "Temperature"));
97 break;
98
99 case STATEVAR_ENTROPY:
100 PetscCall(PetscSectionSetComponentName(section, 0, 0, "EntropyDensity"));
101 PetscCall(PetscSectionSetComponentName(section, 0, 1, "EntropyMomentumX"));
102 PetscCall(PetscSectionSetComponentName(section, 0, 2, "EntropyMomentumY"));
103 PetscCall(PetscSectionSetComponentName(section, 0, 3, "EntropyMomentumZ"));
104 PetscCall(PetscSectionSetComponentName(section, 0, 4, "EntropyTotalEnergy"));
105 break;
106 }
107 PetscFunctionReturn(PETSC_SUCCESS);
108 }
109
110 // Refine DM for high-order viz
VizRefineDM(DM dm,User user,ProblemData problem,SimpleBC bc,Physics phys)111 PetscErrorCode VizRefineDM(DM dm, User user, ProblemData problem, SimpleBC bc, Physics phys) {
112 DM dm_hierarchy[user->app_ctx->viz_refine + 1];
113 VecType vec_type;
114
115 PetscFunctionBeginUser;
116 PetscCall(DMPlexSetRefinementUniform(dm, PETSC_TRUE));
117
118 dm_hierarchy[0] = dm;
119 for (PetscInt i = 0, d = user->app_ctx->degree; i < user->app_ctx->viz_refine; i++) {
120 Mat interp_next;
121 PetscCall(DMRefine(dm_hierarchy[i], MPI_COMM_NULL, &dm_hierarchy[i + 1]));
122 PetscCall(DMClearDS(dm_hierarchy[i + 1]));
123 PetscCall(DMClearFields(dm_hierarchy[i + 1]));
124 PetscCall(DMSetCoarseDM(dm_hierarchy[i + 1], dm_hierarchy[i]));
125 d = (d + 1) / 2;
126 PetscInt q_order = d + user->app_ctx->q_extra;
127 if (i + 1 == user->app_ctx->viz_refine) d = 1;
128 PetscCall(DMGetVecType(dm, &vec_type));
129 PetscCall(DMSetVecType(dm_hierarchy[i + 1], vec_type));
130 PetscCall(SetUpDM(dm_hierarchy[i + 1], problem, d, q_order, bc, phys));
131 PetscCall(DMCreateInterpolation(dm_hierarchy[i], dm_hierarchy[i + 1], &interp_next, NULL));
132 if (!i) user->interp_viz = interp_next;
133 else {
134 Mat C;
135 PetscCall(MatMatMult(interp_next, user->interp_viz, MAT_INITIAL_MATRIX, PETSC_DECIDE, &C));
136 PetscCall(MatDestroy(&interp_next));
137 PetscCall(MatDestroy(&user->interp_viz));
138 user->interp_viz = C;
139 }
140 }
141 for (PetscInt i = 1; i < user->app_ctx->viz_refine; i++) {
142 PetscCall(DMDestroy(&dm_hierarchy[i]));
143 }
144 user->dm_viz = dm_hierarchy[user->app_ctx->viz_refine];
145 PetscFunctionReturn(PETSC_SUCCESS);
146 }
147