Actual source code: ex3.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Solves the 1-dimensional wave equation.\n\n";
4: #include <petscdmda.h>
5: #include <petscdraw.h>
9: int main(int argc,char **argv)
10: {
11: PetscMPIInt rank,size;
13: PetscInt M = 60,time_steps = 100, localsize,j,i,mybase,myend,width,xbase,*localnodes = NULL;
14: DM da;
15: PetscViewer viewer,viewer_private;
16: PetscDraw draw;
17: Vec local,global,copy;
18: PetscScalar *localptr,*copyptr;
19: PetscReal a,h,k;
20: PetscBool flg = PETSC_FALSE;
22: PetscInitialize(&argc,&argv,(char*)0,help);
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: PetscOptionsGetInt(NULL,"-M",&M,NULL);
27: PetscOptionsGetInt(NULL,"-time",&time_steps,NULL);
28: /*
29: Test putting two nodes on each processor, exact last processor gets the rest
30: */
31: PetscOptionsGetBool(NULL,"-distribute",&flg,NULL);
32: if (flg) {
33: PetscMalloc(size*sizeof(PetscInt),&localnodes);
34: for (i=0; i<size-1; i++) localnodes[i] = 2;
35: localnodes[size-1] = M - 2*(size-1);
36: }
38: /* Set up the array */
39: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,M,1,1,localnodes,&da);
40: PetscFree(localnodes);
41: DMCreateGlobalVector(da,&global);
42: DMCreateLocalVector(da,&local);
44: /* Set up display to show combined wave graph */
45: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"Entire Solution",20,480,800,200,&viewer);
46: PetscViewerDrawGetDraw(viewer,0,&draw);
47: PetscDrawSetDoubleBuffer(draw);
49: /* determine starting point of each processor */
50: VecGetOwnershipRange(global,&mybase,&myend);
52: /* set up display to show my portion of the wave */
53: xbase = (int)((mybase)*((800.0 - 4.0*size)/M) + 4.0*rank);
54: width = (int)((myend-mybase)*800./M);
55: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Local Portion of Solution",xbase,200,
56: width,200,&viewer_private);
57: PetscViewerDrawGetDraw(viewer_private,0,&draw);
58: PetscDrawSetDoubleBuffer(draw);
62: /* Initialize the array */
63: VecGetLocalSize(local,&localsize);
64: VecGetArray(local,&localptr);
66: localptr[0] = 0.0;
67: localptr[localsize-1] = 0.0;
68: for (i=1; i<localsize-1; i++) {
69: j = (i-1)+mybase;
70: localptr[i] = sin((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * sin((PETSC_PI*j*2)/((PetscReal)M))) * 2;
71: }
73: VecRestoreArray(local,&localptr);
74: DMLocalToGlobalBegin(da,local,INSERT_VALUES,global);
75: DMLocalToGlobalEnd(da,local,INSERT_VALUES,global);
77: /* Make copy of local array for doing updates */
78: VecDuplicate(local,©);
80: /* Assign Parameters */
81: a= 1.0;
82: h= 1.0/M;
83: k= h;
85: for (j=0; j<time_steps; j++) {
87: /* Global to Local */
88: DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);
89: DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);
91: /*Extract local array */
92: VecGetArray(local,&localptr);
93: VecGetArray(copy,©ptr);
95: /* Update Locally - Make array of new values */
96: /* Note: I don't do anything for the first and last entry */
97: for (i=1; i< localsize-1; i++) {
98: copyptr[i] = .5*(localptr[i+1]+localptr[i-1]) - (k / (2.0*a*h)) * (localptr[i+1] - localptr[i-1]);
99: }
100: VecRestoreArray(copy,©ptr);
101: VecRestoreArray(local,&localptr);
103: /* Local to Global */
104: DMLocalToGlobalBegin(da,copy,INSERT_VALUES,global);
105: DMLocalToGlobalEnd(da,copy,INSERT_VALUES,global);
107: /* View my part of Wave */
108: VecView(copy,viewer_private);
110: /* View global Wave */
111: VecView(global,viewer);
112: }
114: DMDestroy(&da);
115: PetscViewerDestroy(&viewer);
116: PetscViewerDestroy(&viewer_private);
117: VecDestroy(©);
118: VecDestroy(&local);
119: VecDestroy(&global);
121: PetscFinalize();
122: return 0;
123: }