Actual source code: ex20.c
2: static char help[] = "Tests SDALocalToLocalxxx().\n\n";
4: #include petscda.h
5: #include petscsys.h
7: /*
8: For testing purposes this example also creates a
9: DA context. Actually codes using SDA routines will probably
10: not also work with DA contexts.
11: */
15: int main(int argc,char **argv)
16: {
17: PetscMPIInt size,rank;
18: PetscInt M=8,dof=1,stencil_width=1,i,start,end,P=5,N = 6,m=PETSC_DECIDE,n=PETSC_DECIDE,p=PETSC_DECIDE,pt = 0,st = 0;
20: PetscTruth flg2,flg3,flg;
21: DAPeriodicType periodic = DA_NONPERIODIC;
22: DAStencilType stencil_type = DA_STENCIL_STAR;
23: DA da;
24: SDA sda;
25: Vec local,global,local_copy;
26: PetscScalar value,*in,*out;
27: PetscReal norm,work;
28: PetscViewer viewer;
29: char filename[PETSC_MAX_PATH_LEN];
30: FILE *file;
33: PetscInitialize(&argc,&argv,(char*)0,help);
35: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
36: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
37: PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);
38: PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
39: PetscOptionsGetInt(PETSC_NULL,"-stencil_width",&stencil_width,PETSC_NULL);
40: PetscOptionsGetInt(PETSC_NULL,"-periodic",&pt,PETSC_NULL);
41: periodic = (DAPeriodicType) pt;
42: PetscOptionsGetInt(PETSC_NULL,"-stencil_type",&st,PETSC_NULL);
43: stencil_type = (DAStencilType) st;
45: PetscOptionsHasName(PETSC_NULL,"-1d",&flg2);
46: PetscOptionsHasName(PETSC_NULL,"-2d",&flg2);
47: PetscOptionsHasName(PETSC_NULL,"-3d",&flg3);
48: if (flg2) {
49: DACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&da);
50: SDACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&sda);
51: } else if (flg3) {
52: DACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,0,0,0,&da);
53: SDACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,0,0,0,&sda);
54: }
55: else {
56: DACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&da);
57: SDACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&sda);
58: }
60: DACreateGlobalVector(da,&global);
61: DACreateLocalVector(da,&local);
62: VecDuplicate(local,&local_copy);
63:
64: /* zero out vectors so that ghostpoints are zero */
65: value = 0;
66: VecSet(local,value);
67: VecSet(local_copy,value);
69: VecGetOwnershipRange(global,&start,&end);
70: for (i=start; i<end; i++) {
71: value = i + 1;
72: VecSetValues(global,1,&i,&value,INSERT_VALUES);
73: }
74: VecAssemblyBegin(global);
75: VecAssemblyEnd(global);
77: DAGlobalToLocalBegin(da,global,INSERT_VALUES,local);
78: DAGlobalToLocalEnd(da,global,INSERT_VALUES,local);
81: PetscOptionsHasName(PETSC_NULL,"-same_array",&flg);
82: if (flg) {
83: /* test the case where the input and output array is the same */
84: VecCopy(local,local_copy);
85: VecGetArray(local_copy,&in);
86: VecRestoreArray(local_copy,PETSC_NULL);
87: SDALocalToLocalBegin(sda,in,INSERT_VALUES,in);
88: SDALocalToLocalEnd(sda,in,INSERT_VALUES,in);
89: } else {
90: VecGetArray(local,&out);
91: VecRestoreArray(local,PETSC_NULL);
92: VecGetArray(local_copy,&in);
93: VecRestoreArray(local_copy,PETSC_NULL);
94: SDALocalToLocalBegin(sda,out,INSERT_VALUES,in);
95: SDALocalToLocalEnd(sda,out,INSERT_VALUES,in);
96: }
98: PetscOptionsHasName(PETSC_NULL,"-save",&flg);
99: if (flg) {
100: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
101: sprintf(filename,"local.%d",rank);
102: PetscViewerASCIIOpen(PETSC_COMM_SELF,filename,&viewer);
103: PetscViewerASCIIGetPointer(viewer,&file);
104: VecView(local,viewer);
105: fprintf(file,"Vector with correct ghost points\n");
106: VecView(local_copy,viewer);
107: PetscViewerDestroy(viewer);
108: }
110: VecAXPY(local_copy,-1.0,local);
111: VecNorm(local_copy,NORM_MAX,&work);
112: MPI_Allreduce(&work,&norm,1,MPIU_REAL,MPI_MAX,PETSC_COMM_WORLD);
113: if (norm != 0.0) {
114: MPI_Comm_size(PETSC_COMM_WORLD,&size);
115: PetscPrintf(PETSC_COMM_WORLD,"Norm of difference %G should be zero\n",norm);
116: PetscPrintf(PETSC_COMM_WORLD," Number of processors %d\n",size);
117: PetscPrintf(PETSC_COMM_WORLD," M,N,P,dof %D %D %D %D\n",M,N,P,dof);
118: PetscPrintf(PETSC_COMM_WORLD," stencil_width %D stencil_type %d periodic %d\n",stencil_width,(int)stencil_type,(int)periodic);
119: PetscPrintf(PETSC_COMM_WORLD," dimension %d\n",1 + (int) flg2 + (int) flg3);
120: }
121:
122: DADestroy(da);
123: SDADestroy(sda);
124: VecDestroy(local_copy);
125: VecDestroy(local);
126: VecDestroy(global);
127: PetscFinalize();
128: return 0;
129: }