Actual source code: ex3.c
2: static char help[] = "Solves the 1-dimensional wave equation.\n\n";
4: #include petscda.h
8: int main(int argc,char **argv)
9: {
10: PetscMPIInt rank,size;
12: PetscInt M = 60,time_steps = 100, localsize,j,i,mybase,myend,width,xbase,*localnodes = PETSC_NULL;
13: DA da;
14: PetscViewer viewer,viewer_private;
15: PetscDraw draw;
16: Vec local,global,copy;
17: PetscScalar *localptr,*copyptr;
18: PetscReal a,h,k;
19: PetscTruth flg = PETSC_FALSE;
20:
21: PetscInitialize(&argc,&argv,(char*)0,help);
22: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
26: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
27: /*
28: Test putting two nodes on each processor, exact last processor gets the rest
29: */
30: PetscOptionsGetTruth(PETSC_NULL,"-distribute",&flg,PETSC_NULL);
31: if (flg) {
32: PetscMalloc(size*sizeof(PetscInt),&localnodes);
33: for (i=0; i<size-1; i++) { localnodes[i] = 2;}
34: localnodes[size-1] = M - 2*(size-1);
35: }
36:
37: /* Set up the array */
38: DACreate1d(PETSC_COMM_WORLD,DA_XPERIODIC,M,1,1,localnodes,&da);
39: PetscFree(localnodes);
40: DACreateGlobalVector(da,&global);
41: DACreateLocalVector(da,&local);
43: /* Set up display to show combined wave graph */
44: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"Entire Solution",20,480,800,200,&viewer);
45: PetscViewerDrawGetDraw(viewer,0,&draw);
46: PetscDrawSetDoubleBuffer(draw);
48: /* determine starting point of each processor */
49: VecGetOwnershipRange(global,&mybase,&myend);
51: /* set up display to show my portion of the wave */
52: xbase = (int)((mybase)*((800.0 - 4.0*size)/M) + 4.0*rank);
53: width = (int)((myend-mybase)*800./M);
54: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Local Portion of Solution",xbase,200,
55: width,200,&viewer_private);
56: PetscViewerDrawGetDraw(viewer_private,0,&draw);
57: PetscDrawSetDoubleBuffer(draw);
61: /* Initialize the array */
62: VecGetLocalSize(local,&localsize);
63: VecGetArray(local,&localptr);
64: localptr[0] = 0.0;
65: localptr[localsize-1] = 0.0;
66: for (i=1; i<localsize-1; i++) {
67: j=(i-1)+mybase;
68: localptr[i] = sin((PETSC_PI*j*6)/((PetscReal)M)
69: + 1.2 * sin((PETSC_PI*j*2)/((PetscReal)M))) * 2;
70: }
72: VecRestoreArray(local,&localptr);
73: DALocalToGlobal(da,local,INSERT_VALUES,global);
75: /* Make copy of local array for doing updates */
76: VecDuplicate(local,©);
78: /* Assign Parameters */
79: a= 1.0;
80: h= 1.0/M;
81: k= h;
83: for (j=0; j<time_steps; j++) {
85: /* Global to Local */
86: DAGlobalToLocalBegin(da,global,INSERT_VALUES,local);
87: DAGlobalToLocalEnd(da,global,INSERT_VALUES,local);
89: /*Extract local array */
90: VecGetArray(local,&localptr);
91: VecGetArray(copy,©ptr);
93: /* Update Locally - Make array of new values */
94: /* Note: I don't do anything for the first and last entry */
95: for (i=1; i< localsize-1; i++) {
96: copyptr[i] = .5*(localptr[i+1]+localptr[i-1]) -
97: (k / (2.0*a*h)) * (localptr[i+1] - localptr[i-1]);
98: }
99: VecRestoreArray(copy,©ptr);
100: VecRestoreArray(local,&localptr);
102: /* Local to Global */
103: DALocalToGlobal(da,copy,INSERT_VALUES,global);
104:
105: /* View my part of Wave */
106: VecView(copy,viewer_private);
108: /* View global Wave */
109: VecView(global,viewer);
110: }
112: DADestroy(da);
113: PetscViewerDestroy(viewer);
114: PetscViewerDestroy(viewer_private);
115: VecDestroy(copy);
116: VecDestroy(local);
117: VecDestroy(global);
119: PetscFinalize();
120: return 0;
121: }
122: