Actual source code: ex15.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  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";

  5: #include <petscdm.h>
  6: #include <petscdmda.h>

  8: #define DMDA_I 5
  9: #define DMDA_J 4
 10: #define DMDA_K 6

 12: const PetscScalar dmda_i_val[] = { 1.10, 2.3006, 2.32444, 3.44006, 66.9009 };
 13: const PetscScalar dmda_j_val[] = { 0.0, 0.25, 0.5, 0.75 };
 14: const PetscScalar dmda_k_val[] = { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 };

 18: PetscErrorCode MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
 19: {
 20:   MPI_Comm       comm;
 21:   PetscViewer    viewer;
 22:   PetscBool      ismpiio,isskip;
 24: 
 26:   PetscObjectGetComm((PetscObject)x,&comm);

 28:   PetscViewerCreate(comm,&viewer);
 29:   PetscViewerSetType(viewer,PETSCVIEWERBINARY);
 30:   if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
 31:   PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
 32:   if (usempiio) { PetscViewerBinarySetMPIIO(viewer); }
 33:   PetscViewerFileSetName(viewer,fname);

 35:   VecView(x,viewer);

 37:   PetscViewerBinaryGetMPIIO(viewer,&ismpiio);
 38:   if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n"); }
 39:   PetscViewerBinaryGetSkipHeader(viewer,&isskip);
 40:   if (isskip) { PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n"); }

 42:   PetscViewerDestroy(&viewer);
 43:   return(0);
 44: }

 48: PetscErrorCode MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
 49: {
 50:   MPI_Comm       comm;
 51:   PetscViewer    viewer;
 52:   PetscBool      ismpiio,isskip;

 56:   PetscObjectGetComm((PetscObject)x,&comm);

 58:   PetscViewerCreate(comm,&viewer);
 59:   PetscViewerSetType(viewer,PETSCVIEWERBINARY);
 60:   if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
 61:   PetscViewerFileSetMode(viewer,FILE_MODE_READ);
 62:   if (usempiio) { PetscViewerBinarySetMPIIO(viewer); }
 63:   PetscViewerFileSetName(viewer,fname);

 65:   VecLoad(x,viewer);

 67:   PetscViewerBinaryGetSkipHeader(viewer,&isskip);
 68:   if (isskip) { PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n"); }
 69:   PetscViewerBinaryGetMPIIO(viewer,&ismpiio);
 70:   if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n"); }

 72:   PetscViewerDestroy(&viewer);
 73:   return(0);
 74: }

 78: PetscErrorCode DMDAVecGenerateEntries(DM dm,Vec a)
 79: {
 80:   PetscScalar    ****LA_v;
 81:   PetscInt       i,j,k,si,sj,sk,ni,nj,nk,M,N;

 85:   DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
 86:   DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk);
 87:   DMDAVecGetArrayDOF(dm,a,&LA_v);
 88:   for (k=sk; k<sk+nk; k++) {
 89:     for (j=sj; j<sj+nj; j++) {
 90:       for (i=si; i<si+ni; i++) {
 91:         PetscScalar test_value_s;
 92: 
 93:         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));
 94:         LA_v[k][j][i][0] = 3.0 * test_value_s;
 95:         LA_v[k][j][i][1] = 3.0 * test_value_s + 1.0;
 96:         LA_v[k][j][i][2] = 3.0 * test_value_s + 2.0;
 97:       }
 98:     }
 99:   }
100:   DMDAVecRestoreArrayDOF(dm,a,&LA_v);
101:   return(0);
102: }

106: PetscErrorCode HeaderlessBinaryReadCheck(DM dm,const char name[])
107: {
109:   int            fdes;
110:   PetscScalar    buffer[DMDA_I*DMDA_J*DMDA_K*3];
111:   PetscInt       len,d,i,j,k,M,N;
112:   PetscMPIInt    rank;
113:   PetscBool      dataverified = PETSC_TRUE;

116:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
117:   DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
118:   len = DMDA_I*DMDA_J*DMDA_K*3;
119:   if (!rank) {
120:     PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
121:     PetscBinaryRead(fdes,buffer,len,PETSC_SCALAR);
122:     PetscBinaryClose(fdes);

124:     for (k=0; k<DMDA_K; k++) {
125:       for (j=0; j<DMDA_J; j++) {
126:         for (i=0; i<DMDA_I; i++) {
127:           for (d=0; d<3; d++) {
128:             PetscScalar v,test_value_s,test_value;
129:             PetscInt    index;

131:             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));
132:             test_value = 3.0 * test_value_s + (PetscScalar)d;

134:             index = 3*(i + j*M + k*M*N) + d;
135:             v = PetscAbsScalar(test_value-buffer[index]);
136: #if defined(PETSC_USE_COMPLEX)
137:             if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
138:               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);
139:               dataverified = PETSC_FALSE;
140:             }
141: #else
142:             if (PetscRealPart(v) > 1.0e-10) {
143:               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);
144:               dataverified = PETSC_FALSE;
145:             }
146: #endif
147:           }
148:         }
149:       }
150:     }
151:     if (dataverified) {
152:       PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified for: %s\n",name);
153:     }
154:   }
155:   return(0);
156: }

160: PetscErrorCode VecCompare(Vec a,Vec b)
161: {
162:   PetscInt       locmin[2],locmax[2];
163:   PetscReal      min[2],max[2];
164:   Vec            ref;

168:   VecMin(a,&locmin[0],&min[0]);
169:   VecMax(a,&locmax[0],&max[0]);

171:   VecMin(b,&locmin[1],&min[1]);
172:   VecMax(b,&locmax[1],&max[1]);

174:   PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n");
175:   PetscPrintf(PETSC_COMM_WORLD,"  min(a)   = %+1.2e [loc %D]\n",(double)min[0],locmin[0]);
176:   PetscPrintf(PETSC_COMM_WORLD,"  max(a)   = %+1.2e [loc %D]\n",(double)max[0],locmax[0]);

178:   PetscPrintf(PETSC_COMM_WORLD,"  min(b)   = %+1.2e [loc %D]\n",(double)min[1],locmin[1]);
179:   PetscPrintf(PETSC_COMM_WORLD,"  max(b)   = %+1.2e [loc %D]\n",(double)max[1],locmax[1]);

181:   VecDuplicate(a,&ref);
182:   VecCopy(a,ref);
183:   VecAXPY(ref,-1.0,b);
184:   VecMin(ref,&locmin[0],&min[0]);
185:   if (PetscAbsReal(min[0]) > 1.0e-10) {
186:     PetscPrintf(PETSC_COMM_WORLD,"  ERROR: min(a-b) > 1.0e-10\n");
187:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
188:   } else {
189:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) < 1.0e-10\n");
190:   }
191:   VecDestroy(&ref);
192:   return(0);
193: }

197: PetscErrorCode TestDMDAVec(PetscBool usempiio)
198: {
199:   DM             dm;
200:   Vec            x_ref,x_test;
201:   PetscBool      skipheader = PETSC_TRUE;

205:   if (!usempiio) { PetscPrintf(PETSC_COMM_WORLD,"%s\n",__FUNCT__); }
206:   else { PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",__FUNCT__); }
207:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
208:                       DMDA_STENCIL_BOX,DMDA_I,DMDA_J,DMDA_K,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
209:                         3,2,NULL,NULL,NULL,&dm);

211:   DMCreateGlobalVector(dm,&x_ref);
212:   DMDAVecGenerateEntries(dm,x_ref);

214:   if (!usempiio) {
215:     MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref);
216:   } else {
217:     MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref);
218:   }

220:   DMCreateGlobalVector(dm,&x_test);

222:   if (!usempiio) {
223:     MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test);
224:   } else {
225:     MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test);
226:   }

228:   VecCompare(x_ref,x_test);

230:   if (!usempiio) {
231:     HeaderlessBinaryReadCheck(dm,"dmda.pbvec");
232:   } else {
233:     HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec");
234:   }
235:   VecDestroy(&x_ref);
236:   VecDestroy(&x_test);
237:   DMDestroy(&dm);
238:   return(0);
239: }

243: int main(int argc,char **args)
244: {
246:   PetscBool      usempiio = PETSC_FALSE;

248:   PetscInitialize(&argc,&args,(char*)0,help);

250:   PetscOptionsGetBool(NULL,"-usempiio",&usempiio,NULL);
251:   if (!usempiio) {
252:     TestDMDAVec(PETSC_FALSE);
253:   } else {
254: #if defined(PETSC_HAVE_MPIIO)
255:     TestDMDAVec(PETSC_TRUE);
256: #else
257:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");
258: #endif
259:   }
260:   PetscFinalize();
261:   return 0;
262: }