1 static char help[] = "Verify isoperiodic cone corrections";
2
3 #include <petscdmplex.h>
4 #include <petscsf.h>
5 #define EX "ex101.c"
6
7 // Creates periodic solution on a [0,1] x D domain for D dimension
project_function(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)8 static PetscErrorCode project_function(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
9 {
10 PetscReal x_tot = 0;
11
12 PetscFunctionBeginUser;
13 for (PetscInt d = 0; d < dim; d++) x_tot += x[d];
14 for (PetscInt c = 0; c < Nc; c++) {
15 PetscScalar value = c % 2 ? PetscSinReal(2 * M_PI * x_tot) : PetscCosReal(2 * M_PI * x_tot);
16 if (PetscAbsScalar(value) < 1e-7) value = 0.;
17 u[c] = value;
18 }
19 PetscFunctionReturn(PETSC_SUCCESS);
20 }
21
CreateFEField(DM dm,PetscBool use_natural_fe,PetscInt num_comps)22 PetscErrorCode CreateFEField(DM dm, PetscBool use_natural_fe, PetscInt num_comps)
23 {
24 PetscInt degree;
25
26 PetscFunctionBeginUser;
27 { // Get degree of the coords section
28 PetscFE fe;
29 PetscSpace basis_space;
30
31 if (use_natural_fe) {
32 PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fe));
33 } else {
34 DM cdm;
35 PetscCall(DMGetCoordinateDM(dm, &cdm));
36 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
37 }
38 PetscCall(PetscFEGetBasisSpace(fe, &basis_space));
39 PetscCall(PetscSpaceGetDegree(basis_space, °ree, NULL));
40 }
41
42 PetscCall(DMClearFields(dm));
43 PetscCall(DMSetLocalSection(dm, NULL)); // See https://gitlab.com/petsc/petsc/-/issues/1669
44
45 { // Setup fe to load in the initial condition data
46 PetscFE fe;
47 PetscInt dim;
48
49 PetscCall(DMGetDimension(dm, &dim));
50 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, num_comps, PETSC_FALSE, degree, PETSC_DETERMINE, &fe));
51 PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
52 PetscCall(DMCreateDS(dm));
53 PetscCall(PetscFEDestroy(&fe));
54 }
55 PetscFunctionReturn(PETSC_SUCCESS);
56 }
57
main(int argc,char ** argv)58 int main(int argc, char **argv)
59 {
60 MPI_Comm comm;
61 DM dm = NULL;
62 Vec V, V_G2L, V_local;
63 PetscReal norm;
64 PetscBool test_cgns_load = PETSC_FALSE;
65 PetscInt num_comps = 1;
66
67 PetscErrorCode (*funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) = {project_function};
68
69 PetscFunctionBeginUser;
70 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
71 comm = PETSC_COMM_WORLD;
72
73 PetscOptionsBegin(comm, "", "ex101.c Options", "DMPLEX");
74 PetscCall(PetscOptionsBool("-test_cgns_load", "Test VecLoad using CGNS file", EX, test_cgns_load, &test_cgns_load, NULL));
75 PetscCall(PetscOptionsInt("-num_comps", "Number of components in FE field", EX, num_comps, &num_comps, NULL));
76 PetscOptionsEnd();
77
78 PetscCall(DMCreate(comm, &dm));
79 PetscCall(DMSetType(dm, DMPLEX));
80 PetscCall(PetscObjectSetName((PetscObject)dm, "ex101_dm"));
81 PetscCall(DMSetFromOptions(dm));
82 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
83 PetscCall(CreateFEField(dm, PETSC_FALSE, num_comps));
84
85 // Verify that projected function on global vector (then projected onto local vector) is equal to projected function onto a local vector
86 PetscCall(DMGetLocalVector(dm, &V_G2L));
87 PetscCall(DMGetGlobalVector(dm, &V));
88 PetscCall(DMProjectFunction(dm, 0, &funcs, NULL, INSERT_VALUES, V));
89 PetscCall(DMGlobalToLocal(dm, V, INSERT_VALUES, V_G2L));
90
91 PetscCall(DMGetLocalVector(dm, &V_local));
92 PetscCall(DMProjectFunctionLocal(dm, 0, &funcs, NULL, INSERT_VALUES, V_local));
93 PetscCall(VecViewFromOptions(V_local, NULL, "-local_view"));
94
95 PetscCall(VecAXPY(V_G2L, -1, V_local));
96 PetscCall(VecNorm(V_G2L, NORM_MAX, &norm));
97 PetscReal tol = PetscDefined(USE_REAL___FLOAT128) ? 1e-12 : 1e4 * PETSC_MACHINE_EPSILON;
98 if (norm > tol) PetscCall(PetscPrintf(comm, "Error! GlobalToLocal result does not match Local projection by norm %g\n", (double)norm));
99
100 if (test_cgns_load) {
101 #if !defined(PETSC_HAVE_CGNS)
102 SETERRQ(comm, PETSC_ERR_SUP, "PETSc not compiled with CGNS support");
103 #else
104 PetscViewer viewer;
105 DM dm_read, dm_read_output;
106 Vec V_read, V_read_project2local, V_read_output2local;
107 const char *filename = "test_file.cgns";
108
109 // Write previous solution to filename
110 PetscCall(PetscViewerCGNSOpen(comm, filename, FILE_MODE_WRITE, &viewer));
111 PetscCall(VecView(V_local, viewer));
112 PetscCall(PetscViewerDestroy(&viewer));
113
114 // Write DM from written CGNS file
115 PetscCall(DMPlexCreateFromFile(comm, filename, "ex101_dm_read", PETSC_TRUE, &dm_read));
116 PetscCall(DMSetFromOptions(dm_read));
117 PetscCall(DMViewFromOptions(dm_read, NULL, "-dm_view"));
118
119 { // Force isoperiodic point SF to be created to update sfNatural.
120 // Needs to be done before removing the field corresponding to sfNatural
121 PetscSection unused_section;
122 PetscCall(DMGetGlobalSection(dm_read, &unused_section));
123 }
124 PetscCall(CreateFEField(dm_read, PETSC_TRUE, num_comps));
125
126 // Use OutputDM as this doesn't use the isoperiodic pointSF and is compatible with sfNatural
127 PetscCall(DMGetOutputDM(dm_read, &dm_read_output));
128 PetscCall(DMGetLocalVector(dm_read, &V_read_project2local));
129 PetscCall(DMGetLocalVector(dm_read, &V_read_output2local));
130 PetscCall(PetscObjectSetName((PetscObject)V_read_output2local, "V_read_output2local"));
131 PetscCall(PetscObjectSetName((PetscObject)V_read_project2local, "V_read_project2local"));
132
133 PetscCall(DMProjectFunctionLocal(dm_read_output, 0, &funcs, NULL, INSERT_VALUES, V_read_project2local));
134 PetscCall(VecViewFromOptions(V_read_project2local, NULL, "-project2local_view"));
135
136 { // Read solution from file and communicate to local Vec
137 PetscCall(DMGetGlobalVector(dm_read_output, &V_read));
138 PetscCall(PetscObjectSetName((PetscObject)V_read, "V_read"));
139
140 PetscCall(PetscViewerCGNSOpen(comm, filename, FILE_MODE_READ, &viewer));
141 PetscCall(VecLoad(V_read, viewer));
142 PetscCall(DMGlobalToLocal(dm_read_output, V_read, INSERT_VALUES, V_read_output2local));
143
144 PetscCall(PetscViewerDestroy(&viewer));
145 PetscCall(DMRestoreGlobalVector(dm_read_output, &V_read));
146 }
147 PetscCall(VecViewFromOptions(V_read_output2local, NULL, "-output2local_view"));
148
149 PetscCall(VecAXPY(V_read_output2local, -1, V_read_project2local));
150 PetscCall(VecNorm(V_read_output2local, NORM_MAX, &norm));
151 if (norm > tol) PetscCall(PetscPrintf(comm, "Error! CGNS VecLoad result does not match Local projection by norm %g\n", (double)norm));
152
153 PetscCall(DMRestoreLocalVector(dm_read, &V_read_project2local));
154 PetscCall(DMRestoreLocalVector(dm_read, &V_read_output2local));
155 PetscCall(DMDestroy(&dm_read));
156 #endif
157 }
158
159 PetscCall(DMRestoreGlobalVector(dm, &V));
160 PetscCall(DMRestoreLocalVector(dm, &V_G2L));
161 PetscCall(DMRestoreLocalVector(dm, &V_local));
162 PetscCall(DMDestroy(&dm));
163 PetscCall(PetscFinalize());
164 return 0;
165 }
166
167 /*TEST
168
169 test:
170 nsize: 2
171 args: -dm_plex_shape zbox -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 1,2,3 -dm_coord_space -dm_coord_petscspace_degree 3
172 args: -dm_plex_box_bd periodic,periodic,periodic -dm_view ::ascii_info_detail -petscpartitioner_type simple
173
174 test:
175 requires: cgns
176 suffix: cgns
177 nsize: 3
178 args: -dm_plex_filename ${DATAFILESPATH}/meshes/2x2x2_Q3_wave.cgns -dm_plex_cgns_parallel -dm_view ::ascii_info_detail -dm_plex_box_label true -dm_plex_box_label_bd periodic,periodic,periodic -petscpartitioner_type simple -test_cgns_load -num_comps 2
179
180 test:
181 requires: cgns parmetis
182 suffix: cgns_parmetis
183 nsize: 3
184 args: -dm_plex_filename ${DATAFILESPATH}/meshes/2x2x2_Q3_wave.cgns -dm_plex_cgns_parallel -dm_plex_box_label true -dm_plex_box_label_bd periodic,periodic,periodic -petscpartitioner_type parmetis -test_cgns_load -num_comps 2
185 output_file: output/empty.out
186
187 TEST*/
188