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, &degree, 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*/