Actual source code: ex7.c
2: /* Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options] */
4: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
7: /*
8: Include "petscda.h" so that we can use distributed arrays (DAs).
9: Include "petscts.h" so that we can use SNES solvers. Note that this
10: file automatically includes:
11: petscsys.h - base PETSc routines petscvec.h - vectors
12: petscmat.h - matrices
13: petscis.h - index sets petscksp.h - Krylov subspace methods
14: petscviewer.h - viewers petscpc.h - preconditioners
15: petscksp.h - linear solvers
16: */
17: #include petscda.h
18: #include petscts.h
21: /*
22: User-defined routines
23: */
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: PetscInt steps,maxsteps = 100; /* iterations for convergence */
36: PetscErrorCode ierr;
37: DA da;
38: MatFDColoring matfdcoloring;
39: ISColoring iscoloring;
40: PetscReal ftime;
41: SNES ts_snes;
43: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Initialize program
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46: PetscInitialize(&argc,&argv,(char *)0,help);
47: PetscOptionsGetInt(PETSC_NULL,"-max_steps",&maxsteps,PETSC_NULL);
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: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create timestepping solver context
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: TSCreate(PETSC_COMM_WORLD,&ts);
66: TSSetProblemType(ts,TS_NONLINEAR);
67: TSSetRHSFunction(ts,FormFunction,da);
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Create matrix data structure; set Jacobian evaluation routine
72: Set Jacobian matrix data structure and default Jacobian evaluation
73: routine. User can override with:
74: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
75: (unless user explicitly sets preconditioner)
76: -snes_mf_operator : form preconditioning matrix as set by the user,
77: but use matrix-free approx for Jacobian-vector
78: products within Newton-Krylov method
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: DAGetColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
82: DAGetMatrix(da,MATAIJ,&J);
83: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
84: ISColoringDestroy(iscoloring);
85: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,da);
86: MatFDColoringSetFromOptions(matfdcoloring);
87: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
89: TSSetDuration(ts,maxsteps,1.0);
90: TSMonitorSet(ts,MyTSMonitor,0,0);
91:
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Customize nonlinear solver
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: TSSetType(ts,TSBEULER);
96: TSGetSNES(ts,&ts_snes);
97: SNESMonitorSet(ts_snes,MySNESMonitor,PETSC_NULL,PETSC_NULL);
98:
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Set initial conditions
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: FormInitialSolution(da,x);
103: TSSetInitialTimeStep(ts,0.0,.0001);
104: TSSetSolution(ts,x);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Set runtime options
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: TSSetFromOptions(ts);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Solve nonlinear system
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: TSStep(ts,&steps,&ftime);
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Free work space. All PETSc objects should be destroyed when they
118: are no longer needed.
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: MatDestroy(J);
121: MatFDColoringDestroy(matfdcoloring);
122: VecDestroy(x);
123: VecDestroy(r);
124: TSDestroy(ts);
125: DADestroy(da);
127: PetscFinalize();
128: return(0);
129: }
130: /* ------------------------------------------------------------------- */
133: /*
134: FormFunction - Evaluates nonlinear function, F(x).
136: Input Parameters:
137: . ts - the TS context
138: . X - input vector
139: . ptr - optional user-defined context, as set by SNESSetFunction()
141: Output Parameter:
142: . F - function vector
143: */
144: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
145: {
146: DA da = (DA)ptr;
148: PetscInt i,j,Mx,My,xs,ys,xm,ym;
149: PetscReal two = 2.0,hx,hy,hxdhy,hydhx,sx,sy;
150: PetscScalar u,uxx,uyy,**x,**f;
151: Vec localX;
154: DAGetLocalVector(da,&localX);
155: DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
156: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
158: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
159: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
160: hxdhy = hx/hy;
161: hydhx = hy/hx;
163: /*
164: Scatter ghost points to local vector,using the 2-step process
165: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
166: By placing code between these two statements, computations can be
167: done while messages are in transition.
168: */
169: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
170: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
172: /*
173: Get pointers to vector data
174: */
175: DAVecGetArray(da,localX,&x);
176: DAVecGetArray(da,F,&f);
178: /*
179: Get local grid boundaries
180: */
181: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
183: /*
184: Compute function over the locally owned part of the grid
185: */
186: for (j=ys; j<ys+ym; j++) {
187: for (i=xs; i<xs+xm; i++) {
188: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
189: f[j][i] = x[j][i];
190: continue;
191: }
192: u = x[j][i];
193: uxx = (two*u - x[j][i-1] - x[j][i+1])*sx;
194: uyy = (two*u - x[j-1][i] - x[j+1][i])*sy;
195: /* f[j][i] = -(uxx + uyy); */
196: 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 +
197: (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
198: }
199: }
201: /*
202: Restore vectors
203: */
204: DAVecRestoreArray(da,localX,&x);
205: DAVecRestoreArray(da,F,&f);
206: DARestoreLocalVector(da,&localX);
207: PetscLogFlops(11.0*ym*xm);
208: return(0);
209: }
211: /* ------------------------------------------------------------------- */
214: PetscErrorCode FormInitialSolution(DA da,Vec U)
215: {
217: PetscInt i,j,xs,ys,xm,ym,Mx,My;
218: PetscScalar **u;
219: PetscReal hx,hy,x,y,r;
222: DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
223: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
225: hx = 1.0/(PetscReal)(Mx-1);
226: hy = 1.0/(PetscReal)(My-1);
228: /*
229: Get pointers to vector data
230: */
231: DAVecGetArray(da,U,&u);
233: /*
234: Get local grid boundaries
235: */
236: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
238: /*
239: Compute function over the locally owned part of the grid
240: */
241: for (j=ys; j<ys+ym; j++) {
242: y = j*hy;
243: for (i=xs; i<xs+xm; i++) {
244: x = i*hx;
245: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
246: if (r < .125) {
247: u[j][i] = PetscExpScalar(-30.0*r*r*r);
248: } else {
249: u[j][i] = 0.0;
250: }
251: }
252: }
254: /*
255: Restore vectors
256: */
257: DAVecRestoreArray(da,U,&u);
258: return(0);
259: }
263: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
264: {
266: PetscReal norm;
267: MPI_Comm comm;
270: VecNorm(v,NORM_2,&norm);
271: PetscObjectGetComm((PetscObject)ts,&comm);
272: PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
273: return(0);
274: }
278: /*
279: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
280: Input Parameters:
281: snes - the SNES context
282: its - iteration number
283: fnorm - 2-norm function value (may be estimated)
284: ctx - optional user-defined context for private data for the
285: monitor routine, as set by SNESMonitorSet()
286: */
287: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
288: {
290:
292: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
293: return(0);
294: }