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