Actual source code: ex15.c
petsc-3.8.4 2018-03-24
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 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 };
16: PetscErrorCode MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
17: {
18: MPI_Comm comm;
19: PetscViewer viewer;
20: PetscBool ismpiio,isskip;
24: PetscObjectGetComm((PetscObject)x,&comm);
26: PetscViewerCreate(comm,&viewer);
27: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
28: if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
29: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
30: if (usempiio) { PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE); }
31: PetscViewerFileSetName(viewer,fname);
33: VecView(x,viewer);
35: PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
36: if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n"); }
37: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
38: if (isskip) { PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n"); }
40: PetscViewerDestroy(&viewer);
41: return(0);
42: }
44: PetscErrorCode MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
45: {
46: MPI_Comm comm;
47: PetscViewer viewer;
48: PetscBool ismpiio,isskip;
52: PetscObjectGetComm((PetscObject)x,&comm);
54: PetscViewerCreate(comm,&viewer);
55: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
56: if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
57: PetscViewerFileSetMode(viewer,FILE_MODE_READ);
58: if (usempiio) { PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE); }
59: PetscViewerFileSetName(viewer,fname);
61: VecLoad(x,viewer);
63: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
64: if (isskip) { PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n"); }
65: PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
66: if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n"); }
68: PetscViewerDestroy(&viewer);
69: return(0);
70: }
72: PetscErrorCode DMDAVecGenerateEntries(DM dm,Vec a)
73: {
74: PetscScalar ****LA_v;
75: PetscInt i,j,k,si,sj,sk,ni,nj,nk,M,N;
79: DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
80: DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk);
81: DMDAVecGetArrayDOF(dm,a,&LA_v);
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: LA_v[k][j][i][0] = 3.0 * test_value_s;
89: LA_v[k][j][i][1] = 3.0 * test_value_s + 1.0;
90: LA_v[k][j][i][2] = 3.0 * test_value_s + 2.0;
91: }
92: }
93: }
94: DMDAVecRestoreArrayDOF(dm,a,&LA_v);
95: return(0);
96: }
98: PetscErrorCode HeaderlessBinaryReadCheck(DM dm,const char name[])
99: {
101: int fdes;
102: PetscScalar buffer[DMDA_I*DMDA_J*DMDA_K*3];
103: PetscInt len,d,i,j,k,M,N;
104: PetscMPIInt rank;
105: PetscBool dataverified = PETSC_TRUE;
108: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
109: DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
110: len = DMDA_I*DMDA_J*DMDA_K*3;
111: if (!rank) {
112: PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
113: PetscBinaryRead(fdes,buffer,len,PETSC_SCALAR);
114: PetscBinaryClose(fdes);
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<3; d++) {
120: PetscScalar v,test_value_s,test_value;
121: PetscInt index;
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 = 3.0 * test_value_s + (PetscScalar)d;
126: index = 3*(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: 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);
131: dataverified = PETSC_FALSE;
132: }
133: #else
134: if (PetscRealPart(v) > 1.0e-10) {
135: 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);
136: dataverified = PETSC_FALSE;
137: }
138: #endif
139: }
140: }
141: }
142: }
143: if (dataverified) {
144: PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified for: %s\n",name);
145: }
146: }
147: return(0);
148: }
150: PetscErrorCode VecCompare(Vec a,Vec b)
151: {
152: PetscInt locmin[2],locmax[2];
153: PetscReal min[2],max[2];
154: Vec ref;
158: VecMin(a,&locmin[0],&min[0]);
159: VecMax(a,&locmax[0],&max[0]);
161: VecMin(b,&locmin[1],&min[1]);
162: VecMax(b,&locmax[1],&max[1]);
164: PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n");
165: PetscPrintf(PETSC_COMM_WORLD," min(a) = %+1.2e [loc %D]\n",(double)min[0],locmin[0]);
166: PetscPrintf(PETSC_COMM_WORLD," max(a) = %+1.2e [loc %D]\n",(double)max[0],locmax[0]);
168: PetscPrintf(PETSC_COMM_WORLD," min(b) = %+1.2e [loc %D]\n",(double)min[1],locmin[1]);
169: PetscPrintf(PETSC_COMM_WORLD," max(b) = %+1.2e [loc %D]\n",(double)max[1],locmax[1]);
171: VecDuplicate(a,&ref);
172: VecCopy(a,ref);
173: VecAXPY(ref,-1.0,b);
174: VecMin(ref,&locmin[0],&min[0]);
175: if (PetscAbsReal(min[0]) > 1.0e-10) {
176: PetscPrintf(PETSC_COMM_WORLD," ERROR: min(a-b) > 1.0e-10\n");
177: PetscPrintf(PETSC_COMM_WORLD," min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
178: } else {
179: PetscPrintf(PETSC_COMM_WORLD," min(a-b) < 1.0e-10\n");
180: }
181: VecDestroy(&ref);
182: return(0);
183: }
185: PetscErrorCode TestDMDAVec(PetscBool usempiio)
186: {
187: DM dm;
188: Vec x_ref,x_test;
189: PetscBool skipheader = PETSC_TRUE;
193: if (!usempiio) { PetscPrintf(PETSC_COMM_WORLD,"%s\n",PETSC_FUNCTION_NAME); }
194: else { PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",PETSC_FUNCTION_NAME); }
195: 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);
197: DMSetFromOptions(dm);
198: DMSetUp(dm);
200: DMCreateGlobalVector(dm,&x_ref);
201: DMDAVecGenerateEntries(dm,x_ref);
203: if (!usempiio) {
204: MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref);
205: } else {
206: MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref);
207: }
209: DMCreateGlobalVector(dm,&x_test);
211: if (!usempiio) {
212: MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test);
213: } else {
214: MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test);
215: }
217: VecCompare(x_ref,x_test);
219: if (!usempiio) {
220: HeaderlessBinaryReadCheck(dm,"dmda.pbvec");
221: } else {
222: HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec");
223: }
224: VecDestroy(&x_ref);
225: VecDestroy(&x_test);
226: DMDestroy(&dm);
227: return(0);
228: }
230: int main(int argc,char **args)
231: {
233: PetscBool usempiio = PETSC_FALSE;
235: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
236: PetscOptionsGetBool(NULL,NULL,"-usempiio",&usempiio,NULL);
237: if (!usempiio) {
238: TestDMDAVec(PETSC_FALSE);
239: } else {
240: #if defined(PETSC_HAVE_MPIIO)
241: TestDMDAVec(PETSC_TRUE);
242: #else
243: PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");
244: #endif
245: }
246: PetscFinalize();
247: return ierr;
248: }
250: /*TEST
252: test:
253: suffix: 3
254: nsize: 12
255: requires: define(PETSC_HAVE_MPIIO) !define(PETSC_HAVE_LIBMSMPI)
256: args: -usempiio
258: TEST*/