1 static char help[] = "Tests DMComposite routines.\n\n";
2
3 #include <petscdmredundant.h>
4 #include <petscdm.h>
5 #include <petscdmda.h>
6 #include <petscdmcomposite.h>
7 #include <petscpf.h>
8
main(int argc,char ** argv)9 int main(int argc, char **argv)
10 {
11 PetscInt nredundant1 = 5, nredundant2 = 2, i;
12 ISLocalToGlobalMapping *ltog;
13 PetscMPIInt rank, size;
14 DM packer;
15 Vec global, local1, local2, redundant1, redundant2;
16 PF pf;
17 DM da1, da2, dmred1, dmred2;
18 PetscScalar *redundant1a, *redundant2a;
19 PetscViewer sviewer;
20 PetscBool gather_add = PETSC_FALSE;
21
22 PetscFunctionBeginUser;
23 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
24 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
25 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26
27 PetscCall(PetscOptionsGetBool(NULL, NULL, "-gather_add", &gather_add, NULL));
28
29 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &packer));
30
31 PetscCall(DMRedundantCreate(PETSC_COMM_WORLD, 0, nredundant1, &dmred1));
32 PetscCall(DMCreateLocalVector(dmred1, &redundant1));
33 PetscCall(DMCompositeAddDM(packer, dmred1));
34
35 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 8, 1, 1, NULL, &da1));
36 PetscCall(DMSetFromOptions(da1));
37 PetscCall(DMSetUp(da1));
38 PetscCall(DMCreateLocalVector(da1, &local1));
39 PetscCall(DMCompositeAddDM(packer, da1));
40
41 PetscCall(DMRedundantCreate(PETSC_COMM_WORLD, 1 % size, nredundant2, &dmred2));
42 PetscCall(DMCreateLocalVector(dmred2, &redundant2));
43 PetscCall(DMCompositeAddDM(packer, dmred2));
44
45 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 6, 1, 1, NULL, &da2));
46 PetscCall(DMSetFromOptions(da2));
47 PetscCall(DMSetUp(da2));
48 PetscCall(DMCreateLocalVector(da2, &local2));
49 PetscCall(DMCompositeAddDM(packer, da2));
50
51 PetscCall(DMCreateGlobalVector(packer, &global));
52 PetscCall(PFCreate(PETSC_COMM_WORLD, 1, 1, &pf));
53 PetscCall(PFSetType(pf, PFIDENTITY, NULL));
54 PetscCall(PFApplyVec(pf, NULL, global));
55 PetscCall(PFDestroy(&pf));
56 PetscCall(VecView(global, PETSC_VIEWER_STDOUT_WORLD));
57
58 PetscCall(DMCompositeScatter(packer, global, redundant1, local1, redundant2, local2));
59 PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
60 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
61 PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] My part of redundant1 vector\n", rank));
62 PetscCall(VecView(redundant1, sviewer));
63 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
64 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
65 PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] My part of da1 vector\n", rank));
66 PetscCall(VecView(local1, sviewer));
67 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
68 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
69 PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] My part of redundant2 vector\n", rank));
70 PetscCall(VecView(redundant2, sviewer));
71 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
72 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
73 PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] My part of da2 vector\n", rank));
74 PetscCall(VecView(local2, sviewer));
75 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
76 PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
77
78 PetscCall(VecGetArray(redundant1, &redundant1a));
79 PetscCall(VecGetArray(redundant2, &redundant2a));
80 for (i = 0; i < nredundant1; i++) redundant1a[i] = (rank + 2) * i;
81 for (i = 0; i < nredundant2; i++) redundant2a[i] = (rank + 10) * i;
82 PetscCall(VecRestoreArray(redundant1, &redundant1a));
83 PetscCall(VecRestoreArray(redundant2, &redundant2a));
84
85 PetscCall(DMCompositeGather(packer, gather_add ? ADD_VALUES : INSERT_VALUES, global, redundant1, local1, redundant2, local2));
86 PetscCall(VecView(global, PETSC_VIEWER_STDOUT_WORLD));
87
88 /* get the global numbering for each subvector element */
89 PetscCall(DMCompositeGetISLocalToGlobalMappings(packer, <og));
90
91 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Local to global mapping of redundant1 vector\n"));
92 PetscCall(ISLocalToGlobalMappingView(ltog[0], PETSC_VIEWER_STDOUT_WORLD));
93 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Local to global mapping of local1 vector\n"));
94 PetscCall(ISLocalToGlobalMappingView(ltog[1], PETSC_VIEWER_STDOUT_WORLD));
95 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Local to global mapping of redundant2 vector\n"));
96 PetscCall(ISLocalToGlobalMappingView(ltog[2], PETSC_VIEWER_STDOUT_WORLD));
97 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Local to global mapping of local2 vector\n"));
98 PetscCall(ISLocalToGlobalMappingView(ltog[3], PETSC_VIEWER_STDOUT_WORLD));
99
100 for (i = 0; i < 4; i++) PetscCall(ISLocalToGlobalMappingDestroy(<og[i]));
101 PetscCall(PetscFree(ltog));
102
103 PetscCall(DMDestroy(&da1));
104 PetscCall(DMDestroy(&dmred1));
105 PetscCall(DMDestroy(&dmred2));
106 PetscCall(DMDestroy(&da2));
107 PetscCall(VecDestroy(&redundant1));
108 PetscCall(VecDestroy(&redundant2));
109 PetscCall(VecDestroy(&local1));
110 PetscCall(VecDestroy(&local2));
111 PetscCall(VecDestroy(&global));
112 PetscCall(DMDestroy(&packer));
113 PetscCall(PetscFinalize());
114 return 0;
115 }
116
117 /*TEST
118
119 build:
120 requires: !complex
121
122 test:
123 nsize: 3
124
125 test:
126 suffix: 2
127 nsize: 3
128 args: -gather_add
129
130 TEST*/
131