Actual source code: ex15.c

petsc-master 2019-09-15
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>


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

 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 };

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

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

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

 34:   VecView(x,viewer);

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

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

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

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

 55:   PetscViewerCreate(comm,&viewer);
 56:   PetscViewerSetType(viewer,PETSCVIEWERBINARY);
 57:   if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
 58:   PetscViewerFileSetMode(viewer,FILE_MODE_READ);
 59:   if (usempiio) { PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE); }
 60:   PetscViewerFileSetName(viewer,fname);

 62:   VecLoad(x,viewer);

 64:   PetscViewerBinaryGetSkipHeader(viewer,&isskip);
 65:   if (isskip) { PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n"); }
 66:   PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
 67:   if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n"); }

 69:   PetscViewerDestroy(&viewer);
 70:   return(0);
 71: }

 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;

 80:   DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
 81:   DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk);
 82:   DMDAVecGetArrayDOF(dm,a,&LA_v);
 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;

 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:   DMDAVecRestoreArrayDOF(dm,a,&LA_v);
 96:   return(0);
 97: }

 99: PetscErrorCode HeaderlessBinaryReadCheck(DM dm,const char name[])
100: {
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;

109:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
110:   DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
111:   len = DMDA_I*DMDA_J*DMDA_K*dof;
112:   if (!rank) {
113:     PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
114:     PetscBinaryRead(fdes,buffer,len,NULL,PETSC_SCALAR);
115:     PetscBinaryClose(fdes);

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;

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;

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:               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);
132:               dataverified = PETSC_FALSE;
133:             }
134: #else
135:             if (PetscRealPart(v) > 1.0e-10) {
136:               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);
137:               dataverified = PETSC_FALSE;
138:             }
139: #endif
140:           }
141:         }
142:       }
143:     }
144:     if (dataverified) {
145:       PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified for: %s\n",name);
146:     }
147:   }
148:   return(0);
149: }

151: PetscErrorCode VecCompare(Vec a,Vec b)
152: {
153:   PetscInt       locmin[2],locmax[2];
154:   PetscReal      min[2],max[2];
155:   Vec            ref;

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

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

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

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

172:   VecDuplicate(a,&ref);
173:   VecCopy(a,ref);
174:   VecAXPY(ref,-1.0,b);
175:   VecMin(ref,&locmin[0],&min[0]);
176:   if (PetscAbsReal(min[0]) > 1.0e-10) {
177:     PetscPrintf(PETSC_COMM_WORLD,"  ERROR: min(a-b) > 1.0e-10\n");
178:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
179:   } else {
180:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) < 1.0e-10\n");
181:   }
182:   VecDestroy(&ref);
183:   return(0);
184: }

186: PetscErrorCode TestDMDAVec(PetscBool usempiio)
187: {
188:   DM             dm;
189:   Vec            x_ref,x_test;
190:   PetscBool      skipheader = PETSC_TRUE;

194:   if (!usempiio) { PetscPrintf(PETSC_COMM_WORLD,"%s\n",PETSC_FUNCTION_NAME); }
195:   else { PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",PETSC_FUNCTION_NAME); }
196:   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);
198:   DMSetFromOptions(dm);
199:   DMSetUp(dm);

201:   DMCreateGlobalVector(dm,&x_ref);
202:   DMDAVecGenerateEntries(dm,x_ref);

204:   if (!usempiio) {
205:     MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref);
206:   } else {
207:     MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref);
208:   }

210:   DMCreateGlobalVector(dm,&x_test);

212:   if (!usempiio) {
213:     MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test);
214:   } else {
215:     MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test);
216:   }

218:   VecCompare(x_ref,x_test);

220:   if (!usempiio) {
221:     HeaderlessBinaryReadCheck(dm,"dmda.pbvec");
222:   } else {
223:     HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec");
224:   }
225:   VecDestroy(&x_ref);
226:   VecDestroy(&x_test);
227:   DMDestroy(&dm);
228:   return(0);
229: }

231: int main(int argc,char **args)
232: {
234:   PetscBool      usempiio = PETSC_FALSE;

236:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
237:   PetscOptionsGetBool(NULL,NULL,"-usempiio",&usempiio,NULL);
238:   if (!usempiio) {
239:     TestDMDAVec(PETSC_FALSE);
240:   } else {
241: #if defined(PETSC_HAVE_MPIIO)
242:     TestDMDAVec(PETSC_TRUE);
243: #else
244:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");
245: #endif
246:   }
247:   PetscFinalize();
248:   return ierr;
249: }


252: /*TEST

254:    test:

256:    test:
257:       suffix: 2
258:       nsize: 12

260:    test:
261:       suffix: 3
262:       nsize: 12
263:       requires: define(PETSC_HAVE_MPIIO)
264:       args: -usempiio

266: TEST*/