Actual source code: ex101.c
1: static char help[] = "Verify isoperiodic cone corrections";
3: #include <petscdmplex.h>
4: #include <petscsf.h>
5: #define EX "ex101.c"
7: // Creates periodic solution on a [0,1] x D domain for D dimension
8: static PetscErrorCode project_function(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
9: {
10: PetscReal x_tot = 0;
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: }
22: PetscErrorCode CreateFEField(DM dm, PetscBool use_natural_fe, PetscInt num_comps)
23: {
24: PetscInt degree;
26: PetscFunctionBeginUser;
27: { // Get degree of the coords section
28: PetscFE fe;
29: PetscSpace basis_space;
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: }
42: PetscCall(DMClearFields(dm));
43: PetscCall(DMSetLocalSection(dm, NULL)); // See https://gitlab.com/petsc/petsc/-/issues/1669
45: { // Setup fe to load in the initial condition data
46: PetscFE fe;
47: PetscInt dim;
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: }
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;
67: PetscErrorCode (*funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) = {project_function};
69: PetscFunctionBeginUser;
70: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
71: comm = PETSC_COMM_WORLD;
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();
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));
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));
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"));
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));
100: if (test_cgns_load) {
101: #ifndef 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";
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));
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"));
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 dummy_section;
122: PetscCall(DMGetGlobalSection(dm_read, &dummy_section));
123: }
124: PetscCall(CreateFEField(dm_read, PETSC_TRUE, num_comps));
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"));
133: PetscCall(DMProjectFunctionLocal(dm_read_output, 0, &funcs, NULL, INSERT_VALUES, V_read_project2local));
134: PetscCall(VecViewFromOptions(V_read_project2local, NULL, "-project2local_view"));
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"));
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));
144: PetscCall(PetscViewerDestroy(&viewer));
145: PetscCall(DMRestoreGlobalVector(dm_read_output, &V_read));
146: }
147: PetscCall(VecViewFromOptions(V_read_output2local, NULL, "-output2local_view"));
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));
153: PetscCall(DMRestoreLocalVector(dm_read, &V_read_project2local));
154: PetscCall(DMRestoreLocalVector(dm_read, &V_read_output2local));
155: PetscCall(DMDestroy(&dm_read));
156: #endif
157: }
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: }
167: /*TEST
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
174: test:
175: requires: cgns
176: suffix: cgns
177: nsize: 3
178: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/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
180: test:
181: requires: cgns parmetis
182: suffix: cgns_parmetis
183: nsize: 3
184: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/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
186: TEST*/