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