Actual source code: ex12.c
2: /* Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options] */
4: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
6: /*
7: Solves the equation
9: u_tt - \Delta u = 0
11: which we split into two first-order systems
13: u_t - v = 0 so that F(u,v).u = -v
14: v_t - \Delta u = 0 so that F(u,v).v = -Delta u
16: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
17: Include "petscts.h" so that we can use SNES solvers. Note that this
18: file automatically includes:
19: petscsys.h - base PETSc routines petscvec.h - vectors
20: petscmat.h - matrices
21: petscis.h - index sets petscksp.h - Krylov subspace methods
22: petscviewer.h - viewers petscpc.h - preconditioners
23: petscksp.h - linear solvers
24: */
25: #include <petscdmda.h>
26: #include <petscts.h>
29: /*
30: User-defined routines
31: */
38: int main(int argc,char **argv)
39: {
40: TS ts; /* nonlinear solver */
41: Vec x,r; /* solution, residual vectors */
42: Mat J; /* Jacobian matrix */
43: PetscInt steps,maxsteps = 100; /* iterations for convergence */
44: PetscErrorCode ierr;
45: DM da;
46: MatFDColoring matfdcoloring;
47: ISColoring iscoloring;
48: PetscReal ftime;
49: SNES ts_snes;
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Initialize program
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: PetscInitialize(&argc,&argv,(char *)0,help);
55: PetscOptionsGetInt(PETSC_NULL,"-max_steps",&maxsteps,PETSC_NULL);
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create distributed array (DMDA) to manage parallel grid and vectors
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
61: 2,1,PETSC_NULL,PETSC_NULL,&da);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Extract global vectors from DMDA; then duplicate for remaining
65: vectors that are the same types
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: DMCreateGlobalVector(da,&x);
68: VecDuplicate(x,&r);
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Create timestepping solver context
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: TSCreate(PETSC_COMM_WORLD,&ts);
74: TSSetProblemType(ts,TS_NONLINEAR);
75: TSSetRHSFunction(ts,PETSC_NULL,FormFunction,da);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Create matrix data structure; set Jacobian evaluation routine
80: Set Jacobian matrix data structure and default Jacobian evaluation
81: routine. User can override with:
82: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
83: (unless user explicitly sets preconditioner)
84: -snes_mf_operator : form preconditioning matrix as set by the user,
85: but use matrix-free approx for Jacobian-vector
86: products within Newton-Krylov method
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: DMGetColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
90: DMGetMatrix(da,MATAIJ,&J);
91: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
92: ISColoringDestroy(&iscoloring);
93: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,da);
94: MatFDColoringSetFromOptions(matfdcoloring);
95: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
97: TSSetDuration(ts,maxsteps,1.0);
98: TSMonitorSet(ts,MyTSMonitor,0,0);
99:
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Customize nonlinear solver
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: TSSetType(ts,TSBEULER);
104: TSGetSNES(ts,&ts_snes);
105: SNESMonitorSet(ts_snes,MySNESMonitor,PETSC_NULL,PETSC_NULL);
106:
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Set initial conditions
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: FormInitialSolution(da,x);
111: TSSetInitialTimeStep(ts,0.0,.0001);
112: TSSetSolution(ts,x);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Set runtime options
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: TSSetFromOptions(ts);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Solve nonlinear system
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: TSSolve(ts,x,&ftime);
123: TSGetTimeStepNumber(ts,&steps);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Free work space. All PETSc objects should be destroyed when they
127: are no longer needed.
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: MatDestroy(&J);
130: MatFDColoringDestroy(&matfdcoloring);
131: VecDestroy(&x);
132: VecDestroy(&r);
133: TSDestroy(&ts);
134: DMDestroy(&da);
136: PetscFinalize();
137: return(0);
138: }
139: /* ------------------------------------------------------------------- */
142: /*
143: FormFunction - Evaluates nonlinear function, F(x).
145: Input Parameters:
146: . ts - the TS context
147: . X - input vector
148: . ptr - optional user-defined context, as set by SNESSetFunction()
150: Output Parameter:
151: . F - function vector
152: */
153: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
154: {
155: DM da = (DM)ptr;
157: PetscInt i,j,Mx,My,xs,ys,xm,ym;
158: PetscReal two = 2.0,hx,hy,/*hxdhy,hydhx,*/sx,sy;
159: PetscScalar u,uxx,uyy,v,***x,***f;
160: Vec localX;
163: DMGetLocalVector(da,&localX);
164: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
165: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
167: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
168: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
169: /*hxdhy = hx/hy;*/
170: /*hydhx = hy/hx;*/
172: /*
173: Scatter ghost points to local vector,using the 2-step process
174: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
175: By placing code between these two statements, computations can be
176: done while messages are in transition.
177: */
178: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
179: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
181: /*
182: Get pointers to vector data
183: */
184: DMDAVecGetArrayDOF(da,localX,&x);
185: DMDAVecGetArrayDOF(da,F,&f);
187: /*
188: Get local grid boundaries
189: */
190: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
192: /*
193: Compute function over the locally owned part of the grid
194: */
195: for (j=ys; j<ys+ym; j++) {
196: for (i=xs; i<xs+xm; i++) {
197: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
198: f[j][i][0] = x[j][i][0];
199: f[j][i][1] = x[j][i][1];
200: continue;
201: }
202: u = x[j][i][0];
203: v = x[j][i][1];
204: uxx = (two*u - x[j][i-1][0] - x[j][i+1][0])*sx;
205: uyy = (two*u - x[j-1][i][0] - x[j+1][i][0])*sy;
206: f[j][i][0] = -v;
207: f[j][i][1] = -(uxx + uyy);
208: }
209: }
211: /*
212: Restore vectors
213: */
214: DMDAVecRestoreArrayDOF(da,localX,&x);
215: DMDAVecRestoreArrayDOF(da,F,&f);
216: DMRestoreLocalVector(da,&localX);
217: PetscLogFlops(11.0*ym*xm);
218: return(0);
219: }
221: /* ------------------------------------------------------------------- */
224: PetscErrorCode FormInitialSolution(DM da,Vec U)
225: {
227: PetscInt i,j,xs,ys,xm,ym,Mx,My;
228: PetscScalar ***u;
229: PetscReal hx,hy,x,y,r;
232: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
233: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
235: hx = 1.0/(PetscReal)(Mx-1);
236: hy = 1.0/(PetscReal)(My-1);
238: /*
239: Get pointers to vector data
240: */
241: DMDAVecGetArrayDOF(da,U,&u);
243: /*
244: Get local grid boundaries
245: */
246: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
248: /*
249: Compute function over the locally owned part of the grid
250: */
251: for (j=ys; j<ys+ym; j++) {
252: y = j*hy;
253: for (i=xs; i<xs+xm; i++) {
254: x = i*hx;
255: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
256: if (r < .125) {
257: u[j][i][0] = PetscExpScalar(-30.0*r*r*r);
258: u[j][i][1] = 0.0;
259: } else {
260: u[j][i][0] = 0.0;
261: u[j][i][1] = 0.0;
262: }
263: }
264: }
266: /*
267: Restore vectors
268: */
269: DMDAVecRestoreArrayDOF(da,U,&u);
270: return(0);
271: }
275: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
276: {
278: PetscReal norm;
279: MPI_Comm comm;
282: VecNorm(v,NORM_2,&norm);
283: PetscObjectGetComm((PetscObject)ts,&comm);
284: PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
285: return(0);
286: }
290: /*
291: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
292: Input Parameters:
293: snes - the SNES context
294: its - iteration number
295: fnorm - 2-norm function value (may be estimated)
296: ctx - optional user-defined context for private data for the
297: monitor routine, as set by SNESMonitorSet()
298: */
299: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
300: {
302:
304: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
305: return(0);
306: }