Actual source code: ex15.c
petsc-3.13.6 2020-09-29
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*/