Actual source code: ex5.c
petsc-3.3-p6 2013-02-11
2: /* This file created by Peter Mell 6/30/95 */
4: static char help[] = "Solves the one dimensional heat equation.\n\n";
6: #include <petscdmda.h>
10: int main(int argc,char **argv)
11: {
12: PetscMPIInt rank,size;
14: PetscInt M = 14,time_steps = 1000,w=1,s=1,localsize,j,i,mybase,myend;
15: DM da;
16: PetscViewer viewer;
17: PetscDraw draw;
18: Vec local,global,copy;
19: PetscScalar *localptr,*copyptr;
20: PetscReal h,k;
21:
22: PetscInitialize(&argc,&argv,(char*)0,help);
24: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
25: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
26:
27: /* Set up the array */
28: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,M,w,s,PETSC_NULL,&da);
29: DMCreateGlobalVector(da,&global);
30: DMCreateLocalVector(da,&local);
31: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
32: MPI_Comm_size(PETSC_COMM_WORLD,&size);
34: /* Make copy of local array for doing updates */
35: VecDuplicate(local,©);
37: /* Set Up Display to Show Heat Graph */
38: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,480,500,160,&viewer);
39: PetscViewerDrawGetDraw(viewer,0,&draw);
40: PetscDrawSetDoubleBuffer(draw);
42: /* determine starting point of each processor */
43: VecGetOwnershipRange(global,&mybase,&myend);
45: /* Initialize the Array */
46: VecGetLocalSize (local,&localsize);
47: VecGetArray (local,&localptr);
48: VecGetArray (copy,©ptr);
49: localptr[0] = copyptr[0] = 0.0;
50: localptr[localsize-1] = copyptr[localsize-1] = 1.0;
51: for (i=1; i<localsize-1; i++) {
52: j=(i-1)+mybase;
53: localptr[i] = sin((PETSC_PI*j*6)/((PetscReal)M)
54: + 1.2 * sin((PETSC_PI*j*2)/((PetscReal)M))) * 4+4;
55: }
57: VecRestoreArray(local,&localptr);
58: VecRestoreArray(copy,©ptr);
59: DMLocalToGlobalBegin(da,local,INSERT_VALUES,global);
60: DMLocalToGlobalEnd(da,local,INSERT_VALUES,global);
62: /* Assign Parameters */
63: h= 1.0/M;
64: k= h*h/2.2;
66: for (j=0; j<time_steps; j++) {
68: /* Global to Local */
69: DMGlobalToLocalBegin(da,global,INSERT_VALUES,local);
70: DMGlobalToLocalEnd(da,global,INSERT_VALUES,local);
72: /*Extract local array */
73: VecGetArray(local,&localptr);
74: VecGetArray (copy,©ptr);
76: /* Update Locally - Make array of new values */
77: /* Note: I don't do anything for the first and last entry */
78: for (i=1; i< localsize-1; i++) {
79: copyptr[i] = localptr[i] + (k/(h*h)) *
80: (localptr[i+1]-2.0*localptr[i]+localptr[i-1]);
81: }
82:
83: VecRestoreArray(copy,©ptr);
84: VecRestoreArray(local,&localptr);
86: /* Local to Global */
87: DMLocalToGlobalBegin(da,copy,INSERT_VALUES,global);
88: DMLocalToGlobalEnd(da,copy,INSERT_VALUES,global);
89:
90: /* View Wave */
91: VecView(global,viewer);
93: }
95: PetscViewerDestroy(&viewer);
96: VecDestroy(©);
97: VecDestroy(&local);
98: VecDestroy(&global);
99: DMDestroy(&da);
100: PetscFinalize();
101: return 0;
102: }
103: