Actual source code: ex12.c
petsc-3.3-p7 2013-05-11
2: /* Program usage: mpiexec -n <procs> ex12 [-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: */
32: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
33: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
34: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,void *);
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);
62: DMDASetFieldName(da,0,"u");
63: DMDASetFieldName(da,1,"v");
65: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: Extract global vectors from DMDA; then duplicate for remaining
67: vectors that are the same types
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69: DMCreateGlobalVector(da,&x);
70: VecDuplicate(x,&r);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create timestepping solver context
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: TSCreate(PETSC_COMM_WORLD,&ts);
76: TSSetDM(ts,da);
77: TSSetProblemType(ts,TS_NONLINEAR);
78: TSSetRHSFunction(ts,PETSC_NULL,FormFunction,da);
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Create matrix data structure; set Jacobian evaluation routine
83: Set Jacobian matrix data structure and default Jacobian evaluation
84: routine. User can override with:
85: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
86: (unless user explicitly sets preconditioner)
87: -snes_mf_operator : form preconditioning matrix as set by the user,
88: but use matrix-free approx for Jacobian-vector
89: products within Newton-Krylov method
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: DMCreateColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
93: DMCreateMatrix(da,MATAIJ,&J);
94: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
95: ISColoringDestroy(&iscoloring);
96: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
97: MatFDColoringSetFromOptions(matfdcoloring);
98: TSGetSNES(ts,&ts_snes);
99: SNESSetJacobian(ts_snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
101: TSSetDuration(ts,maxsteps,1.0);
102: TSMonitorSet(ts,MyTSMonitor,0,0);
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Customize nonlinear solver
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: TSSetType(ts,TSBEULER);
108: TSGetSNES(ts,&ts_snes);
109: SNESMonitorSet(ts_snes,MySNESMonitor,PETSC_NULL,PETSC_NULL);
110:
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Set initial conditions
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: FormInitialSolution(da,x);
115: TSSetInitialTimeStep(ts,0.0,.0001);
116: TSSetSolution(ts,x);
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Set runtime options
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: TSSetFromOptions(ts);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Solve nonlinear system
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: TSSolve(ts,x,&ftime);
127: TSGetTimeStepNumber(ts,&steps);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Free work space. All PETSc objects should be destroyed when they
131: are no longer needed.
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: MatDestroy(&J);
134: MatFDColoringDestroy(&matfdcoloring);
135: VecDestroy(&x);
136: VecDestroy(&r);
137: TSDestroy(&ts);
138: DMDestroy(&da);
140: PetscFinalize();
141: return(0);
142: }
143: /* ------------------------------------------------------------------- */
146: /*
147: FormFunction - Evaluates nonlinear function, F(x).
149: Input Parameters:
150: . ts - the TS context
151: . X - input vector
152: . ptr - optional user-defined context, as set by SNESSetFunction()
154: Output Parameter:
155: . F - function vector
156: */
157: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
158: {
159: DM da = (DM)ptr;
161: PetscInt i,j,Mx,My,xs,ys,xm,ym;
162: PetscReal hx,hy,/*hxdhy,hydhx,*/sx,sy;
163: PetscScalar u,uxx,uyy,v,***x,***f;
164: Vec localX;
167: DMGetLocalVector(da,&localX);
168: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
169: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
171: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
172: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
173: /*hxdhy = hx/hy;*/
174: /*hydhx = hy/hx;*/
176: /*
177: Scatter ghost points to local vector,using the 2-step process
178: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
179: By placing code between these two statements, computations can be
180: done while messages are in transition.
181: */
182: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
183: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
185: /*
186: Get pointers to vector data
187: */
188: DMDAVecGetArrayDOF(da,localX,&x);
189: DMDAVecGetArrayDOF(da,F,&f);
191: /*
192: Get local grid boundaries
193: */
194: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
196: /*
197: Compute function over the locally owned part of the grid
198: */
199: for (j=ys; j<ys+ym; j++) {
200: for (i=xs; i<xs+xm; i++) {
201: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
202: f[j][i][0] = x[j][i][0];
203: f[j][i][1] = x[j][i][1];
204: continue;
205: }
206: u = x[j][i][0];
207: v = x[j][i][1];
208: uxx = (-2.0*u + x[j][i-1][0] + x[j][i+1][0])*sx;
209: uyy = (-2.0*u + x[j-1][i][0] + x[j+1][i][0])*sy;
210: f[j][i][0] = v;
211: f[j][i][1] = uxx + uyy;
212: }
213: }
215: /*
216: Restore vectors
217: */
218: DMDAVecRestoreArrayDOF(da,localX,&x);
219: DMDAVecRestoreArrayDOF(da,F,&f);
220: DMRestoreLocalVector(da,&localX);
221: PetscLogFlops(11.0*ym*xm);
222: return(0);
223: }
225: /* ------------------------------------------------------------------- */
228: PetscErrorCode FormInitialSolution(DM da,Vec U)
229: {
231: PetscInt i,j,xs,ys,xm,ym,Mx,My;
232: PetscScalar ***u;
233: PetscReal hx,hy,x,y,r;
236: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
237: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
239: hx = 1.0/(PetscReal)(Mx-1);
240: hy = 1.0/(PetscReal)(My-1);
242: /*
243: Get pointers to vector data
244: */
245: DMDAVecGetArrayDOF(da,U,&u);
247: /*
248: Get local grid boundaries
249: */
250: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
252: /*
253: Compute function over the locally owned part of the grid
254: */
255: for (j=ys; j<ys+ym; j++) {
256: y = j*hy;
257: for (i=xs; i<xs+xm; i++) {
258: x = i*hx;
259: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
260: if (r < .125) {
261: u[j][i][0] = PetscExpScalar(-30.0*r*r*r);
262: u[j][i][1] = 0.0;
263: } else {
264: u[j][i][0] = 0.0;
265: u[j][i][1] = 0.0;
266: }
267: }
268: }
270: /*
271: Restore vectors
272: */
273: DMDAVecRestoreArrayDOF(da,U,&u);
274: return(0);
275: }
279: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
280: {
282: PetscReal norm;
283: MPI_Comm comm;
286: VecNorm(v,NORM_2,&norm);
287: PetscObjectGetComm((PetscObject)ts,&comm);
288: PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
289: return(0);
290: }
294: /*
295: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
296: Input Parameters:
297: snes - the SNES context
298: its - iteration number
299: fnorm - 2-norm function value (may be estimated)
300: ctx - optional user-defined context for private data for the
301: monitor routine, as set by SNESMonitorSet()
302: */
303: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
304: {
306:
308: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
309: return(0);
310: }