Actual source code: ex3.c
2: static char help[] = "Solves the 1-dimensional wave equation.\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
6: #include <petscdraw.h>
8: int main(int argc,char **argv)
9: {
10: PetscMPIInt rank,size;
11: PetscInt M = 60,time_steps = 100, localsize,j,i,mybase,myend,width,xbase,*localnodes = NULL;
12: DM da;
13: PetscViewer viewer,viewer_private;
14: PetscDraw draw;
15: Vec local,global;
16: PetscScalar *localptr,*globalptr;
17: PetscReal a,h,k;
18: PetscBool flg = PETSC_FALSE;
20: PetscInitialize(&argc,&argv,(char*)0,help);
21: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
22: MPI_Comm_size(PETSC_COMM_WORLD,&size);
24: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
25: PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);
26: /*
27: Test putting two nodes on each processor, exact last processor gets the rest
28: */
29: PetscOptionsGetBool(NULL,NULL,"-distribute",&flg,NULL);
30: if (flg) {
31: PetscMalloc1(size,&localnodes);
32: for (i=0; i<size-1; i++) localnodes[i] = 2;
33: localnodes[size-1] = M - 2*(size-1);
34: }
36: /* Set up the array */
37: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,1,localnodes,&da);
38: DMSetFromOptions(da);
39: DMSetUp(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,width,200,&viewer_private);
56: PetscViewerDrawGetDraw(viewer_private,0,&draw);
57: PetscDrawSetDoubleBuffer(draw);
59: /* Initialize the array */
60: VecGetLocalSize(local,&localsize);
61: VecGetArray(global,&globalptr);
63: for (i=1; i<localsize-1; i++) {
64: j = (i-1)+mybase;
65: globalptr[i-1] = PetscSinReal((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI*j*2)/((PetscReal)M))) * 2;
66: }
68: VecRestoreArray(global,&globalptr);
70: /* Assign Parameters */
71: a= 1.0;
72: h= 1.0/M;
73: k= h;
75: for (j=0; j<time_steps; j++) {
77: /* Global to Local */
78: DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);
79: DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);
81: /*Extract local array */
82: VecGetArray(local,&localptr);
83: VecGetArray(global,&globalptr);
85: /* Update Locally - Make array of new values */
86: /* Note: I don't do anything for the first and last entry */
87: for (i=1; i< localsize-1; i++) {
88: globalptr[i-1] = .5*(localptr[i+1]+localptr[i-1]) - (k / (2.0*a*h)) * (localptr[i+1] - localptr[i-1]);
89: }
90: VecRestoreArray(global,&globalptr);
91: VecRestoreArray(local,&localptr);
93: /* View my part of Wave */
94: VecView(global,viewer_private);
96: /* View global Wave */
97: VecView(global,viewer);
98: }
100: DMDestroy(&da);
101: PetscViewerDestroy(&viewer);
102: PetscViewerDestroy(&viewer_private);
103: VecDestroy(&local);
104: VecDestroy(&global);
106: PetscFinalize();
107: return 0;
108: }
110: /*TEST
112: test:
113: nsize: 3
114: args: -time 50 -nox
115: requires: x
117: TEST*/