Actual source code: ex7.c
petsc-3.3-p7 2013-05-11
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 "petscdmda.h" so that we can use distributed arrays (DMDAs).
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 <petscdmda.h>
18: #include <petscts.h>
21: /*
22: User-defined routines
23: */
24: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
25: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
26: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,void *);
30: int main(int argc,char **argv)
31: {
32: TS ts; /* time integrator */
33: SNES snes; /* nonlinear solver */
34: Vec x,r; /* solution, residual vectors */
35: Mat J; /* Jacobian matrix */
36: PetscInt steps,maxsteps = 100; /* iterations for convergence */
37: PetscErrorCode ierr;
38: DM da;
39: MatFDColoring matfdcoloring;
40: ISColoring iscoloring;
41: PetscReal ftime;
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 (DMDA) to manage parallel grid and vectors
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
53: 1,1,PETSC_NULL,PETSC_NULL,&da);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Extract global vectors from DMDA; then duplicate for remaining
57: vectors that are the same types
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: DMCreateGlobalVector(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,PETSC_NULL,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: DMCreateColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
82: DMCreateMatrix(da,MATAIJ,&J);
83: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
84: ISColoringDestroy(&iscoloring);
85: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))SNESTSFormFunction,ts);
86: MatFDColoringSetFromOptions(matfdcoloring);
87: TSGetSNES(ts,&snes);
88: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
90: TSSetDuration(ts,maxsteps,1.0);
91: TSMonitorSet(ts,MyTSMonitor,0,0);
92:
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Customize nonlinear solver
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: TSSetType(ts,TSBEULER);
97: SNESMonitorSet(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: TSSolve(ts,x,&ftime);
115: TSGetTimeStepNumber(ts,&steps);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Free work space. All PETSc objects should be destroyed when they
119: are no longer needed.
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: MatDestroy(&J);
122: MatFDColoringDestroy(&matfdcoloring);
123: VecDestroy(&x);
124: VecDestroy(&r);
125: TSDestroy(&ts);
126: DMDestroy(&da);
128: PetscFinalize();
129: return(0);
130: }
131: /* ------------------------------------------------------------------- */
134: /*
135: FormFunction - Evaluates nonlinear function, F(x).
137: Input Parameters:
138: . ts - the TS context
139: . X - input vector
140: . ptr - optional user-defined context, as set by SNESSetFunction()
142: Output Parameter:
143: . F - function vector
144: */
145: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
146: {
147: DM da = (DM)ptr;
149: PetscInt i,j,Mx,My,xs,ys,xm,ym;
150: PetscReal two = 2.0,hx,hy,sx,sy;
151: PetscScalar u,uxx,uyy,**x,**f;
152: Vec localX;
155: DMGetLocalVector(da,&localX);
156: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
157: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
159: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
160: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
162: /*
163: Scatter ghost points to local vector,using the 2-step process
164: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
165: By placing code between these two statements, computations can be
166: done while messages are in transition.
167: */
168: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
169: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
171: /*
172: Get pointers to vector data
173: */
174: DMDAVecGetArray(da,localX,&x);
175: DMDAVecGetArray(da,F,&f);
177: /*
178: Get local grid boundaries
179: */
180: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
182: /*
183: Compute function over the locally owned part of the grid
184: */
185: for (j=ys; j<ys+ym; j++) {
186: for (i=xs; i<xs+xm; i++) {
187: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
188: f[j][i] = x[j][i];
189: continue;
190: }
191: u = x[j][i];
192: uxx = (two*u - x[j][i-1] - x[j][i+1])*sx;
193: uyy = (two*u - x[j-1][i] - x[j+1][i])*sy;
194: /* f[j][i] = -(uxx + uyy); */
195: 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 +
196: (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
197: }
198: }
200: /*
201: Restore vectors
202: */
203: DMDAVecRestoreArray(da,localX,&x);
204: DMDAVecRestoreArray(da,F,&f);
205: DMRestoreLocalVector(da,&localX);
206: PetscLogFlops(11.0*ym*xm);
207: return(0);
208: }
210: /* ------------------------------------------------------------------- */
213: PetscErrorCode FormInitialSolution(DM da,Vec U)
214: {
216: PetscInt i,j,xs,ys,xm,ym,Mx,My;
217: PetscScalar **u;
218: PetscReal hx,hy,x,y,r;
221: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
222: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
224: hx = 1.0/(PetscReal)(Mx-1);
225: hy = 1.0/(PetscReal)(My-1);
227: /*
228: Get pointers to vector data
229: */
230: DMDAVecGetArray(da,U,&u);
232: /*
233: Get local grid boundaries
234: */
235: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
237: /*
238: Compute function over the locally owned part of the grid
239: */
240: for (j=ys; j<ys+ym; j++) {
241: y = j*hy;
242: for (i=xs; i<xs+xm; i++) {
243: x = i*hx;
244: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
245: if (r < .125) {
246: u[j][i] = PetscExpScalar(-30.0*r*r*r);
247: } else {
248: u[j][i] = 0.0;
249: }
250: }
251: }
253: /*
254: Restore vectors
255: */
256: DMDAVecRestoreArray(da,U,&u);
257: return(0);
258: }
262: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
263: {
265: PetscReal norm;
266: MPI_Comm comm;
269: VecNorm(v,NORM_2,&norm);
270: PetscObjectGetComm((PetscObject)ts,&comm);
271: PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
272: return(0);
273: }
277: /*
278: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
279: Input Parameters:
280: snes - the SNES context
281: its - iteration number
282: fnorm - 2-norm function value (may be estimated)
283: ctx - optional user-defined context for private data for the
284: monitor routine, as set by SNESMonitorSet()
285: */
286: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
287: {
289:
291: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
292: return(0);
293: }