xref: /petsc/src/dm/tutorials/ex15.c (revision 6a98f8dc3f2c9149905a87dc2e9d0fedaf64e09a) ! 
1 
2 static char help[] = "Tests VecView() functionality with DMDA objects when using:"\
3 "(i) a PetscViewer binary with MPI-IO support; and (ii) when the binary header is skipped.\n\n";
4 
5 #include <petscdm.h>
6 #include <petscdmda.h>
7 
8 
9 #define DMDA_I 5
10 #define DMDA_J 4
11 #define DMDA_K 6
12 
13 const PetscReal dmda_i_val[] = { 1.10, 2.3006, 2.32444, 3.44006, 66.9009 };
14 const PetscReal dmda_j_val[] = { 0.0, 0.25, 0.5, 0.75 };
15 const PetscReal dmda_k_val[] = { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 };
16 
17 PetscErrorCode MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
18 {
19   MPI_Comm       comm;
20   PetscViewer    viewer;
21   PetscBool      ismpiio,isskip;
22   PetscErrorCode ierr;
23 
24   PetscFunctionBeginUser;
25   ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr);
26 
27   ierr = PetscViewerCreate(comm,&viewer);CHKERRQ(ierr);
28   ierr = PetscViewerSetType(viewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
29   if (skippheader) { ierr = PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE);CHKERRQ(ierr); }
30   ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
31   if (usempiio) { ierr = PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE);CHKERRQ(ierr); }
32   ierr = PetscViewerFileSetName(viewer,fname);CHKERRQ(ierr);
33 
34   ierr = VecView(x,viewer);CHKERRQ(ierr);
35 
36   ierr = PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);CHKERRQ(ierr);
37   if (ismpiio) { ierr = PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n");CHKERRQ(ierr); }
38   ierr = PetscViewerBinaryGetSkipHeader(viewer,&isskip);CHKERRQ(ierr);
39   if (isskip) { ierr = PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n");CHKERRQ(ierr); }
40 
41   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
42   PetscFunctionReturn(0);
43 }
44 
45 PetscErrorCode MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
46 {
47   MPI_Comm       comm;
48   PetscViewer    viewer;
49   PetscBool      ismpiio,isskip;
50   PetscErrorCode ierr;
51 
52   PetscFunctionBeginUser;
53   ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr);
54 
55   ierr = PetscViewerCreate(comm,&viewer);CHKERRQ(ierr);
56   ierr = PetscViewerSetType(viewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
57   if (skippheader) { ierr = PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE);CHKERRQ(ierr); }
58   ierr = PetscViewerFileSetMode(viewer,FILE_MODE_READ);CHKERRQ(ierr);
59   if (usempiio) { ierr = PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE);CHKERRQ(ierr); }
60   ierr = PetscViewerFileSetName(viewer,fname);CHKERRQ(ierr);
61 
62   ierr = VecLoad(x,viewer);CHKERRQ(ierr);
63 
64   ierr = PetscViewerBinaryGetSkipHeader(viewer,&isskip);CHKERRQ(ierr);
65   if (isskip) { ierr = PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n");CHKERRQ(ierr); }
66   ierr = PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);CHKERRQ(ierr);
67   if (ismpiio) { ierr = PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n");CHKERRQ(ierr); }
68 
69   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
70   PetscFunctionReturn(0);
71 }
72 
73 PetscErrorCode DMDAVecGenerateEntries(DM dm,Vec a)
74 {
75   PetscScalar    ****LA_v;
76   PetscInt       i,j,k,l,si,sj,sk,ni,nj,nk,M,N,dof;
77   PetscErrorCode ierr;
78 
79   PetscFunctionBeginUser;
80   ierr = DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
81   ierr = DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk);CHKERRQ(ierr);
82   ierr = DMDAVecGetArrayDOF(dm,a,&LA_v);CHKERRQ(ierr);
83   for (k=sk; k<sk+nk; k++) {
84     for (j=sj; j<sj+nj; j++) {
85       for (i=si; i<si+ni; i++) {
86         PetscScalar test_value_s;
87 
88         test_value_s = dmda_i_val[i]*((PetscScalar)i) + dmda_j_val[j]*((PetscScalar)(i+j*M)) + dmda_k_val[k]*((PetscScalar)(i + j*M + k*M*N));
89         for (l=0; l<dof; l++) {
90           LA_v[k][j][i][l] = (PetscScalar)dof * test_value_s + (PetscScalar)l;
91         }
92       }
93     }
94   }
95   ierr = DMDAVecRestoreArrayDOF(dm,a,&LA_v);CHKERRQ(ierr);
96   PetscFunctionReturn(0);
97 }
98 
99 PetscErrorCode HeaderlessBinaryReadCheck(DM dm,const char name[])
100 {
101   PetscErrorCode ierr;
102   int            fdes;
103   PetscScalar    buffer[DMDA_I*DMDA_J*DMDA_K*10];
104   PetscInt       len,d,i,j,k,M,N,dof;
105   PetscMPIInt    rank;
106   PetscBool      dataverified = PETSC_TRUE;
107 
108   PetscFunctionBeginUser;
109   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
110   ierr = DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
111   len = DMDA_I*DMDA_J*DMDA_K*dof;
112   if (!rank) {
113     ierr = PetscBinaryOpen(name,FILE_MODE_READ,&fdes);CHKERRQ(ierr);
114     ierr = PetscBinaryRead(fdes,buffer,len,NULL,PETSC_SCALAR);CHKERRQ(ierr);
115     ierr = PetscBinaryClose(fdes);CHKERRQ(ierr);
116 
117     for (k=0; k<DMDA_K; k++) {
118       for (j=0; j<DMDA_J; j++) {
119         for (i=0; i<DMDA_I; i++) {
120           for (d=0; d<dof; d++) {
121             PetscScalar v,test_value_s,test_value;
122             PetscInt    index;
123 
124             test_value_s = dmda_i_val[i]*((PetscScalar)i) + dmda_j_val[j]*((PetscScalar)(i+j*M)) + dmda_k_val[k]*((PetscScalar)(i + j*M + k*M*N));
125             test_value = (PetscScalar)dof * test_value_s + (PetscScalar)d;
126 
127             index = dof*(i + j*M + k*M*N) + d;
128             v = PetscAbsScalar(test_value-buffer[index]);
129 #if defined(PETSC_USE_COMPLEX)
130             if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
131               ierr = PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %D,%D,%D(%D)])\n",(double)PetscRealPart(test_value),(double)PetscImaginaryPart(test_value),i,j,k,d);CHKERRQ(ierr);
132               dataverified = PETSC_FALSE;
133             }
134 #else
135             if (PetscRealPart(v) > 1.0e-10) {
136               ierr = PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %D,%D,%D(%D)])\n",(double)PetscRealPart(test_value),i,j,k,d);CHKERRQ(ierr);
137               dataverified = PETSC_FALSE;
138             }
139 #endif
140           }
141         }
142       }
143     }
144     if (dataverified) {
145       ierr = PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified for: %s\n",name);CHKERRQ(ierr);
146     }
147   }
148   PetscFunctionReturn(0);
149 }
150 
151 PetscErrorCode VecCompare(Vec a,Vec b)
152 {
153   PetscInt       locmin[2],locmax[2];
154   PetscReal      min[2],max[2];
155   Vec            ref;
156   PetscErrorCode ierr;
157 
158   PetscFunctionBeginUser;
159   ierr = VecMin(a,&locmin[0],&min[0]);CHKERRQ(ierr);
160   ierr = VecMax(a,&locmax[0],&max[0]);CHKERRQ(ierr);
161 
162   ierr = VecMin(b,&locmin[1],&min[1]);CHKERRQ(ierr);
163   ierr = VecMax(b,&locmax[1],&max[1]);CHKERRQ(ierr);
164 
165   ierr = PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n");CHKERRQ(ierr);
166   ierr = PetscPrintf(PETSC_COMM_WORLD,"  min(a)   = %+1.2e [loc %D]\n",(double)min[0],locmin[0]);CHKERRQ(ierr);
167   ierr = PetscPrintf(PETSC_COMM_WORLD,"  max(a)   = %+1.2e [loc %D]\n",(double)max[0],locmax[0]);CHKERRQ(ierr);
168 
169   ierr = PetscPrintf(PETSC_COMM_WORLD,"  min(b)   = %+1.2e [loc %D]\n",(double)min[1],locmin[1]);CHKERRQ(ierr);
170   ierr = PetscPrintf(PETSC_COMM_WORLD,"  max(b)   = %+1.2e [loc %D]\n",(double)max[1],locmax[1]);CHKERRQ(ierr);
171 
172   ierr = VecDuplicate(a,&ref);CHKERRQ(ierr);
173   ierr = VecCopy(a,ref);CHKERRQ(ierr);
174   ierr = VecAXPY(ref,-1.0,b);CHKERRQ(ierr);
175   ierr = VecMin(ref,&locmin[0],&min[0]);CHKERRQ(ierr);
176   if (PetscAbsReal(min[0]) > 1.0e-10) {
177     ierr = PetscPrintf(PETSC_COMM_WORLD,"  ERROR: min(a-b) > 1.0e-10\n");CHKERRQ(ierr);
178     ierr = PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));CHKERRQ(ierr);
179   } else {
180     ierr = PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) < 1.0e-10\n");CHKERRQ(ierr);
181   }
182   ierr = VecDestroy(&ref);CHKERRQ(ierr);
183   PetscFunctionReturn(0);
184 }
185 
186 PetscErrorCode TestDMDAVec(PetscBool usempiio)
187 {
188   DM             dm;
189   Vec            x_ref,x_test;
190   PetscBool      skipheader = PETSC_TRUE;
191   PetscErrorCode ierr;
192 
193   PetscFunctionBeginUser;
194   if (!usempiio) { ierr = PetscPrintf(PETSC_COMM_WORLD,"%s\n",PETSC_FUNCTION_NAME);CHKERRQ(ierr); }
195   else { ierr = PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",PETSC_FUNCTION_NAME);CHKERRQ(ierr); }
196   ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,DMDA_I,DMDA_J,DMDA_K,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
197                         3,2,NULL,NULL,NULL,&dm);CHKERRQ(ierr);
198   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
199   ierr = DMSetUp(dm);CHKERRQ(ierr);
200 
201   ierr = DMCreateGlobalVector(dm,&x_ref);CHKERRQ(ierr);
202   ierr = DMDAVecGenerateEntries(dm,x_ref);CHKERRQ(ierr);
203 
204   if (!usempiio) {
205     ierr = MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref);CHKERRQ(ierr);
206   } else {
207     ierr = MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref);CHKERRQ(ierr);
208   }
209 
210   ierr = DMCreateGlobalVector(dm,&x_test);CHKERRQ(ierr);
211 
212   if (!usempiio) {
213     ierr = MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test);CHKERRQ(ierr);
214   } else {
215     ierr = MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test);CHKERRQ(ierr);
216   }
217 
218   ierr = VecCompare(x_ref,x_test);CHKERRQ(ierr);
219 
220   if (!usempiio) {
221     ierr = HeaderlessBinaryReadCheck(dm,"dmda.pbvec");CHKERRQ(ierr);
222   } else {
223     ierr = HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec");CHKERRQ(ierr);
224   }
225   ierr = VecDestroy(&x_ref);CHKERRQ(ierr);
226   ierr = VecDestroy(&x_test);CHKERRQ(ierr);
227   ierr = DMDestroy(&dm);CHKERRQ(ierr);
228   PetscFunctionReturn(0);
229 }
230 
231 int main(int argc,char **args)
232 {
233   PetscErrorCode ierr;
234   PetscBool      usempiio = PETSC_FALSE;
235 
236   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
237   ierr = PetscOptionsGetBool(NULL,NULL,"-usempiio",&usempiio,NULL);CHKERRQ(ierr);
238   if (!usempiio) {
239     ierr = TestDMDAVec(PETSC_FALSE);CHKERRQ(ierr);
240   } else {
241 #if defined(PETSC_HAVE_MPIIO)
242     ierr = TestDMDAVec(PETSC_TRUE);CHKERRQ(ierr);
243 #else
244     ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");CHKERRQ(ierr);
245 #endif
246   }
247   ierr = PetscFinalize();
248   return ierr;
249 }
250 
251 
252 /*TEST
253 
254    test:
255 
256    test:
257       suffix: 2
258       nsize: 12
259 
260    test:
261       suffix: 3
262       nsize: 12
263       requires: define(PETSC_HAVE_MPIIO)
264       args: -usempiio
265 
266 TEST*/
267 
268