Actual source code: ex7.c
1: /*$Id: ex7.c,v 1.7 2001/08/21 21:04:28 bsmith Exp $*/
3: /* Program usage: mpirun -np <procs> ex5 [-help] [all PETSc options] */
5: static char help[] = "Nonlinear, time-dependent PDE in 2d.n";
8: /*
9: Include "petscda.h" so that we can use distributed arrays (DAs).
10: Include "petscts.h" so that we can use SNES solvers. Note that this
11: file automatically includes:
12: petsc.h - base PETSc routines petscvec.h - vectors
13: petscsys.h - system routines petscmat.h - matrices
14: petscis.h - index sets petscksp.h - Krylov subspace methods
15: petscviewer.h - viewers petscpc.h - preconditioners
16: petscsles.h - linear solvers
17: */
18: #include petscda.h
19: #include petscts.h
22: /*
23: User-defined routines
24: */
25: extern int FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DA,Vec);
26: extern int Monitor(TS,int,PetscReal,Vec,void*);
28: #undef __FUNCT__
30: int main(int argc,char **argv)
31: {
32: TS ts; /* nonlinear solver */
33: Vec x,r; /* solution, residual vectors */
34: Mat J; /* Jacobian matrix */
35: int steps; /* iterations for convergence */
36: int ierr;
37: DA da;
38: MatFDColoring matfdcoloring;
39: ISColoring iscoloring;
40: PetscReal ftime;
42: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: Initialize program
44: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46: PetscInitialize(&argc,&argv,(char *)0,help);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create distributed array (DA) to manage parallel grid and vectors
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,
53: 1,1,PETSC_NULL,PETSC_NULL,&da);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Extract global vectors from DA; then duplicate for remaining
57: vectors that are the same types
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: DACreateGlobalVector(da,&x);
60: VecDuplicate(x,&r);
62: TSCreate(PETSC_COMM_WORLD,&ts);
63: TSSetProblemType(ts,TS_NONLINEAR);
64: TSSetRHSFunction(ts,FormFunction,da);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create matrix data structure; set Jacobian evaluation routine
70: Set Jacobian matrix data structure and default Jacobian evaluation
71: routine. User can override with:
72: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
73: (unless user explicitly sets preconditioner)
74: -snes_mf_operator : form preconditioning matrix as set by the user,
75: but use matrix-free approx for Jacobian-vector
76: products within Newton-Krylov method
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: DAGetColoring(da,IS_COLORING_LOCAL,&iscoloring);
80: DAGetMatrix(da,MATMPIAIJ,&J);
81: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
82: ISColoringDestroy(iscoloring);
83: MatFDColoringSetFunction(matfdcoloring,(int (*)(void))FormFunction,da);
84: MatFDColoringSetFromOptions(matfdcoloring);
85: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Customize nonlinear solver; set runtime options
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: TSSetInitialTimeStep(ts,0.0,.0001);
90: TSSetType(ts,TS_BEULER);
91: TSSetDuration(ts,100,1.0);
92: TSSetFromOptions(ts);
93: TSSetMonitor(ts,Monitor,0,0);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Set initial conditions
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: FormInitialSolution(da,x);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Solve nonlinear system
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: TSSetSolution(ts,x);
104: TSStep(ts,&steps,&ftime);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Free work space. All PETSc objects should be destroyed when they
109: are no longer needed.
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: MatDestroy(J);
113: MatFDColoringDestroy(matfdcoloring);
114: VecDestroy(x);
115: VecDestroy(r);
116: TSDestroy(ts);
117: DADestroy(da);
118: PetscFinalize();
120: return(0);
121: }
122: /* ------------------------------------------------------------------- */
123: #undef __FUNCT__
125: /*
126: FormFunction - Evaluates nonlinear function, F(x).
128: Input Parameters:
129: . ts - the TS context
130: . X - input vector
131: . ptr - optional user-defined context, as set by SNESSetFunction()
133: Output Parameter:
134: . F - function vector
135: */
136: int FormFunction(TS ts,PetscReal time,Vec X,Vec F,void *ptr)
137: {
138: DA da = (DA)ptr;
139: int ierr,i,j,Mx,My,xs,ys,xm,ym;
140: PetscReal two = 2.0,hx,hy,hxdhy,hydhx,sx,sy;
141: PetscScalar u,uxx,uyy,**x,**f;
142: Vec localX;
145: DAGetLocalVector(da,&localX);
146: DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
147: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
149: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
150: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
151: hxdhy = hx/hy;
152: hydhx = hy/hx;
154: /*
155: Scatter ghost points to local vector,using the 2-step process
156: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
157: By placing code between these two statements, computations can be
158: done while messages are in transition.
159: */
160: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
161: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
163: /*
164: Get pointers to vector data
165: */
166: DAVecGetArray(da,localX,(void**)&x);
167: DAVecGetArray(da,F,(void**)&f);
169: /*
170: Get local grid boundaries
171: */
172: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
174: /*
175: Compute function over the locally owned part of the grid
176: */
177: for (j=ys; j<ys+ym; j++) {
178: for (i=xs; i<xs+xm; i++) {
179: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
180: f[j][i] = x[j][i];
181: continue;
182: }
183: u = x[j][i];
184: uxx = (two*u - x[j][i-1] - x[j][i+1])*sx;
185: uyy = (two*u - x[j-1][i] - x[j+1][i])*sy;
186: /* f[j][i] = -(uxx + uyy); */
187: f[j][i] = -u*(uxx + uyy) - (4.0 - 1.0)*((x[j][i+1] - x[j][i-1])*(x[j][i+1] - x[j][i-1])*.25*sx +
188: (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
189: }
190: }
192: /*
193: Restore vectors
194: */
195: DAVecRestoreArray(da,localX,(void**)&x);
196: DAVecRestoreArray(da,F,(void**)&f);
197: DARestoreLocalVector(da,&localX);
198: PetscLogFlops(11*ym*xm);
199: return(0);
200: }
202: /* ------------------------------------------------------------------- */
203: #undef __FUNCT__
205: int FormInitialSolution(DA da,Vec U)
206: {
207: int ierr,i,j,xs,ys,xm,ym,Mx,My;
208: PetscScalar **u;
209: PetscReal hx,hy,x,y,r;
212: DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
213: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
215: hx = 1.0/(PetscReal)(Mx-1);
216: hy = 1.0/(PetscReal)(My-1);
218: /*
219: Get pointers to vector data
220: */
221: DAVecGetArray(da,U,(void**)&u);
223: /*
224: Get local grid boundaries
225: */
226: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
228: /*
229: Compute function over the locally owned part of the grid
230: */
231: for (j=ys; j<ys+ym; j++) {
232: y = j*hy;
233: for (i=xs; i<xs+xm; i++) {
234: x = i*hx;
235: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
236: if (r < .125) {
237: u[j][i] = PetscExpScalar(-30.0*r*r*r);
238: } else {
239: u[j][i] = 0.0;
240: }
241: }
242: }
244: /*
245: Restore vectors
246: */
247: DAVecRestoreArray(da,U,(void**)&u);
248: return(0);
249: }
251: #undef __FUNCT__
253: int Monitor(TS ts,int step,PetscReal ptime,Vec v,void *ctx)
254: {
255: int ierr;
256: PetscReal norm;
257: MPI_Comm comm;
260: VecNorm(v,NORM_2,&norm);
261: PetscObjectGetComm((PetscObject)ts,&comm);
262: PetscPrintf(comm,"timestep %d time %g norm %gn",step,ptime,norm);
263: return(0);
264: }