xref: /petsc/src/dm/tutorials/ex15.c (revision f97672e55eacc8688507b9471cd7ec2664d7f203) 
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   PetscCall(PetscObjectGetComm((PetscObject)x,&comm));
24 
25   PetscCall(PetscViewerCreate(comm,&viewer));
26   PetscCall(PetscViewerSetType(viewer,PETSCVIEWERBINARY));
27   if (skippheader) PetscCall(PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE));
28   PetscCall(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE));
29   if (usempiio) PetscCall(PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE));
30   PetscCall(PetscViewerFileSetName(viewer,fname));
31 
32   PetscCall(VecView(x,viewer));
33 
34   PetscCall(PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio));
35   if (ismpiio) PetscCall(PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n"));
36   PetscCall(PetscViewerBinaryGetSkipHeader(viewer,&isskip));
37   if (isskip) PetscCall(PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n"));
38 
39   PetscCall(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   PetscCall(PetscObjectGetComm((PetscObject)x,&comm));
51 
52   PetscCall(PetscViewerCreate(comm,&viewer));
53   PetscCall(PetscViewerSetType(viewer,PETSCVIEWERBINARY));
54   if (skippheader) PetscCall(PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE));
55   PetscCall(PetscViewerFileSetMode(viewer,FILE_MODE_READ));
56   if (usempiio) PetscCall(PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE));
57   PetscCall(PetscViewerFileSetName(viewer,fname));
58 
59   PetscCall(VecLoad(x,viewer));
60 
61   PetscCall(PetscViewerBinaryGetSkipHeader(viewer,&isskip));
62   if (isskip) PetscCall(PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n"));
63   PetscCall(PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio));
64   if (ismpiio) PetscCall(PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n"));
65 
66   PetscCall(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   PetscCall(DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL));
77   PetscCall(DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk));
78   PetscCall(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   PetscCall(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   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
105   PetscCall(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     PetscCall(PetscBinaryOpen(name,FILE_MODE_READ,&fdes));
109     PetscCall(PetscBinaryRead(fdes,buffer,len,NULL,PETSC_SCALAR));
110     PetscCall(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               PetscCall(PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "(%" PetscInt_FMT ")])\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               PetscCall(PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "(%" PetscInt_FMT ")])\n",(double)PetscRealPart(test_value),i,j,k,d));
132               dataverified = PETSC_FALSE;
133             }
134 #endif
135           }
136         }
137       }
138     }
139     if (dataverified) {
140       PetscCall(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   PetscCall(VecMin(a,&locmin[0],&min[0]));
154   PetscCall(VecMax(a,&locmax[0],&max[0]));
155 
156   PetscCall(VecMin(b,&locmin[1],&min[1]));
157   PetscCall(VecMax(b,&locmax[1],&max[1]));
158 
159   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n"));
160   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  min(a)   = %+1.2e [loc %" PetscInt_FMT "]\n",(double)min[0],locmin[0]));
161   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  max(a)   = %+1.2e [loc %" PetscInt_FMT "]\n",(double)max[0],locmax[0]));
162 
163   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  min(b)   = %+1.2e [loc %" PetscInt_FMT "]\n",(double)min[1],locmin[1]));
164   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  max(b)   = %+1.2e [loc %" PetscInt_FMT "]\n",(double)max[1],locmax[1]));
165 
166   PetscCall(VecDuplicate(a,&ref));
167   PetscCall(VecCopy(a,ref));
168   PetscCall(VecAXPY(ref,-1.0,b));
169   PetscCall(VecMin(ref,&locmin[0],&min[0]));
170   if (PetscAbsReal(min[0]) > 1.0e-10) {
171     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  ERROR: min(a-b) > 1.0e-10\n"));
172     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0])));
173   } else {
174     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) < 1.0e-10\n"));
175   }
176   PetscCall(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 
186   PetscFunctionBeginUser;
187   if (!usempiio) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s\n",PETSC_FUNCTION_NAME));
188   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",PETSC_FUNCTION_NAME));
189   PetscCall(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,
190                          3,2,NULL,NULL,NULL,&dm));
191   PetscCall(DMSetFromOptions(dm));
192   PetscCall(DMSetUp(dm));
193 
194   PetscCall(DMCreateGlobalVector(dm,&x_ref));
195   PetscCall(DMDAVecGenerateEntries(dm,x_ref));
196 
197   if (!usempiio) {
198     PetscCall(MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref));
199   } else {
200     PetscCall(MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref));
201   }
202 
203   PetscCall(DMCreateGlobalVector(dm,&x_test));
204 
205   if (!usempiio) {
206     PetscCall(MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test));
207   } else {
208     PetscCall(MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test));
209   }
210 
211   PetscCall(VecCompare(x_ref,x_test));
212 
213   if (!usempiio) {
214     PetscCall(HeaderlessBinaryReadCheck(dm,"dmda.pbvec"));
215   } else {
216     PetscCall(HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec"));
217   }
218   PetscCall(VecDestroy(&x_ref));
219   PetscCall(VecDestroy(&x_test));
220   PetscCall(DMDestroy(&dm));
221   PetscFunctionReturn(0);
222 }
223 
224 int main(int argc,char **args)
225 {
226   PetscBool      usempiio = PETSC_FALSE;
227 
228   PetscCall(PetscInitialize(&argc,&args,(char*)0,help));
229   PetscCall(PetscOptionsGetBool(NULL,NULL,"-usempiio",&usempiio,NULL));
230   if (!usempiio) {
231     PetscCall(TestDMDAVec(PETSC_FALSE));
232   } else {
233 #if defined(PETSC_HAVE_MPIIO)
234     PetscCall(TestDMDAVec(PETSC_TRUE));
235 #else
236     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n"));
237 #endif
238   }
239   PetscCall(PetscFinalize());
240   return 0;
241 }
242 
243 /*TEST
244 
245    test:
246 
247    test:
248       suffix: 2
249       nsize: 12
250 
251    test:
252       suffix: 3
253       nsize: 12
254       requires: defined(PETSC_HAVE_MPIIO)
255       args: -usempiio
256 
257 TEST*/
258