Actual source code: ex5.c
2: /* This file created by Peter Mell 6/30/95 */
4: static char help[] = "Solves the one dimensional heat equation.\n\n";
6: #include petscda.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: DA 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: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,M,w,s,PETSC_NULL,&da);
29: DACreateGlobalVector(da,&global);
30: DACreateLocalVector(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: DALocalToGlobal(da,local,INSERT_VALUES,global);
61: /* Assign Parameters */
62: h= 1.0/M;
63: k= h*h/2.2;
65: for (j=0; j<time_steps; j++) {
67: /* Global to Local */
68: DAGlobalToLocalBegin(da,global,INSERT_VALUES,local);
69: DAGlobalToLocalEnd(da,global,INSERT_VALUES,local);
71: /*Extract local array */
72: VecGetArray(local,&localptr);
73: VecGetArray (copy,©ptr);
75: /* Update Locally - Make array of new values */
76: /* Note: I don't do anything for the first and last entry */
77: for (i=1; i< localsize-1; i++) {
78: copyptr[i] = localptr[i] + (k/(h*h)) *
79: (localptr[i+1]-2.0*localptr[i]+localptr[i-1]);
80: }
81:
82: VecRestoreArray(copy,©ptr);
83: VecRestoreArray(local,&localptr);
85: /* Local to Global */
86: DALocalToGlobal(da,copy,INSERT_VALUES,global);
87:
88: /* View Wave */
89: VecView(global,viewer);
91: }
93: PetscViewerDestroy(viewer);
94: VecDestroy(copy);
95: VecDestroy(local);
96: VecDestroy(global);
97: DADestroy(da);
98: PetscFinalize();
99: return 0;
100: }
101: