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;
12: PetscInt M = 60,time_steps = 100, localsize,j,i,mybase,myend,width,xbase,*localnodes = NULL;
13: DM da;
14: PetscViewer viewer,viewer_private;
15: PetscDraw draw;
16: Vec local,global;
17: PetscScalar *localptr,*globalptr;
18: PetscReal a,h,k;
19: PetscBool flg = PETSC_FALSE;
21: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
22: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
26: PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);
27: /*
28: Test putting two nodes on each processor, exact last processor gets the rest
29: */
30: PetscOptionsGetBool(NULL,NULL,"-distribute",&flg,NULL);
31: if (flg) {
32: PetscMalloc1(size,&localnodes);
33: for (i=0; i<size-1; i++) localnodes[i] = 2;
34: localnodes[size-1] = M - 2*(size-1);
35: }
37: /* Set up the array */
38: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,1,localnodes,&da);
39: DMSetFromOptions(da);
40: DMSetUp(da);
41: PetscFree(localnodes);
42: DMCreateGlobalVector(da,&global);
43: DMCreateLocalVector(da,&local);
45: /* Set up display to show combined wave graph */
46: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"Entire Solution",20,480,800,200,&viewer);
47: PetscViewerDrawGetDraw(viewer,0,&draw);
48: PetscDrawSetDoubleBuffer(draw);
50: /* determine starting point of each processor */
51: VecGetOwnershipRange(global,&mybase,&myend);
53: /* set up display to show my portion of the wave */
54: xbase = (int)((mybase)*((800.0 - 4.0*size)/M) + 4.0*rank);
55: width = (int)((myend-mybase)*800./M);
56: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Local Portion of Solution",xbase,200,width,200,&viewer_private);
57: PetscViewerDrawGetDraw(viewer_private,0,&draw);
58: PetscDrawSetDoubleBuffer(draw);
60: /* Initialize the array */
61: VecGetLocalSize(local,&localsize);
62: VecGetArray(global,&globalptr);
64: for (i=1; i<localsize-1; i++) {
65: j = (i-1)+mybase;
66: globalptr[i-1] = PetscSinReal((PETSC_PI*j*6)/((PetscReal)M) + 1.2 * PetscSinReal((PETSC_PI*j*2)/((PetscReal)M))) * 2;
67: }
69: VecRestoreArray(global,&globalptr);
71: /* Assign Parameters */
72: a= 1.0;
73: h= 1.0/M;
74: k= h;
76: for (j=0; j<time_steps; j++) {
78: /* Global to Local */
79: DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);
80: DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);
82: /*Extract local array */
83: VecGetArray(local,&localptr);
84: VecGetArray(global,&globalptr);
86: /* Update Locally - Make array of new values */
87: /* Note: I don't do anything for the first and last entry */
88: for (i=1; i< localsize-1; i++) {
89: globalptr[i-1] = .5*(localptr[i+1]+localptr[i-1]) - (k / (2.0*a*h)) * (localptr[i+1] - localptr[i-1]);
90: }
91: VecRestoreArray(global,&globalptr);
92: VecRestoreArray(local,&localptr);
94: /* View my part of Wave */
95: VecView(global,viewer_private);
97: /* View global Wave */
98: VecView(global,viewer);
99: }
101: DMDestroy(&da);
102: PetscViewerDestroy(&viewer);
103: PetscViewerDestroy(&viewer_private);
104: VecDestroy(&local);
105: VecDestroy(&global);
107: PetscFinalize();
108: return ierr;
109: }
111: /*TEST
113: test:
114: nsize: 3
115: args: -time 50 -nox
116: requires: x
118: TEST*/