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 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 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 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