Actual source code: ex7.c
petsc-3.5.4 2015-05-23
2: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
5: /*
6: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
7: Include "petscts.h" so that we can use SNES solvers. Note that this
8: file automatically includes:
9: petscsys.h - base PETSc routines petscvec.h - vectors
10: petscmat.h - matrices
11: petscis.h - index sets petscksp.h - Krylov subspace methods
12: petscviewer.h - viewers petscpc.h - preconditioners
13: petscksp.h - linear solvers
14: */
15: #include <petscdm.h>
16: #include <petscdmda.h>
17: #include <petscts.h>
20: /*
21: User-defined routines
22: */
23: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
24: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
25: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,void*);
29: int main(int argc,char **argv)
30: {
31: TS ts; /* time integrator */
32: SNES snes;
33: Vec x,r; /* solution, residual vectors */
34: PetscInt maxsteps = 100; /* iterations for convergence */
36: DM da;
38: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39: Initialize program
40: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41: PetscInitialize(&argc,&argv,(char*)0,help);
43: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Create distributed array (DMDA) to manage parallel grid and vectors
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
47: 1,1,NULL,NULL,&da);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Extract global vectors from DMDA; then duplicate for remaining
51: vectors that are the same types
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: DMCreateGlobalVector(da,&x);
54: VecDuplicate(x,&r);
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Create timestepping solver context
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: TSCreate(PETSC_COMM_WORLD,&ts);
60: TSSetProblemType(ts,TS_NONLINEAR);
61: TSSetRHSFunction(ts,NULL,FormFunction,da);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Create matrix data structure; set Jacobian evaluation routine
66: Set Jacobian matrix data structure and default Jacobian evaluation
67: routine. User can override with:
68: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
69: (unless user explicitly sets preconditioner)
70: -snes_mf_operator : form preconditioning matrix as set by the user,
71: but use matrix-free approx for Jacobian-vector
72: products within Newton-Krylov method
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: TSSetDuration(ts,maxsteps,1.0);
77: TSMonitorSet(ts,MyTSMonitor,0,0);
78: TSSetDM(ts,da);
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Customize nonlinear solver
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: TSSetType(ts,TSBEULER);
83: TSGetSNES(ts,&snes);
84: SNESMonitorSet(snes,MySNESMonitor,NULL,NULL);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Set initial conditions
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: FormInitialSolution(da,x);
90: TSSetInitialTimeStep(ts,0.0,.0001);
91: TSSetSolution(ts,x);
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Set runtime options
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: TSSetFromOptions(ts);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Solve nonlinear system
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: TSSolve(ts,x);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Free work space. All PETSc objects should be destroyed when they
105: are no longer needed.
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: VecDestroy(&x);
108: VecDestroy(&r);
109: TSDestroy(&ts);
110: DMDestroy(&da);
112: PetscFinalize();
113: return(0);
114: }
115: /* ------------------------------------------------------------------- */
118: /*
119: FormFunction - Evaluates nonlinear function, F(x).
121: Input Parameters:
122: . ts - the TS context
123: . X - input vector
124: . ptr - optional user-defined context, as set by SNESSetFunction()
126: Output Parameter:
127: . F - function vector
128: */
129: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
130: {
131: DM da;
133: PetscInt i,j,Mx,My,xs,ys,xm,ym;
134: PetscReal two = 2.0,hx,hy,sx,sy;
135: PetscScalar u,uxx,uyy,**x,**f;
136: Vec localX;
139: TSGetDM(ts,&da);
140: DMGetLocalVector(da,&localX);
141: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
142: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
144: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
145: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
147: /*
148: Scatter ghost points to local vector,using the 2-step process
149: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
150: By placing code between these two statements, computations can be
151: done while messages are in transition.
152: */
153: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
154: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
156: /*
157: Get pointers to vector data
158: */
159: DMDAVecGetArray(da,localX,&x);
160: DMDAVecGetArray(da,F,&f);
162: /*
163: Get local grid boundaries
164: */
165: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
167: /*
168: Compute function over the locally owned part of the grid
169: */
170: for (j=ys; j<ys+ym; j++) {
171: for (i=xs; i<xs+xm; i++) {
172: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
173: f[j][i] = x[j][i];
174: continue;
175: }
176: u = x[j][i];
177: uxx = (two*u - x[j][i-1] - x[j][i+1])*sx;
178: uyy = (two*u - x[j-1][i] - x[j+1][i])*sy;
179: /* f[j][i] = -(uxx + uyy); */
180: 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 +
181: (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
182: }
183: }
185: /*
186: Restore vectors
187: */
188: DMDAVecRestoreArray(da,localX,&x);
189: DMDAVecRestoreArray(da,F,&f);
190: DMRestoreLocalVector(da,&localX);
191: PetscLogFlops(11.0*ym*xm);
192: return(0);
193: }
195: /* ------------------------------------------------------------------- */
198: PetscErrorCode FormInitialSolution(DM da,Vec U)
199: {
201: PetscInt i,j,xs,ys,xm,ym,Mx,My;
202: PetscScalar **u;
203: PetscReal hx,hy,x,y,r;
206: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
207: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
209: hx = 1.0/(PetscReal)(Mx-1);
210: hy = 1.0/(PetscReal)(My-1);
212: /*
213: Get pointers to vector data
214: */
215: DMDAVecGetArray(da,U,&u);
217: /*
218: Get local grid boundaries
219: */
220: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
222: /*
223: Compute function over the locally owned part of the grid
224: */
225: for (j=ys; j<ys+ym; j++) {
226: y = j*hy;
227: for (i=xs; i<xs+xm; i++) {
228: x = i*hx;
229: r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5));
230: if (r < .125) u[j][i] = PetscExpReal(-30.0*r*r*r);
231: else u[j][i] = 0.0;
232: }
233: }
235: /*
236: Restore vectors
237: */
238: DMDAVecRestoreArray(da,U,&u);
239: return(0);
240: }
244: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
245: {
247: PetscReal norm;
248: MPI_Comm comm;
251: VecNorm(v,NORM_2,&norm);
252: PetscObjectGetComm((PetscObject)ts,&comm);
253: PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
254: return(0);
255: }
259: /*
260: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
261: Input Parameters:
262: snes - the SNES context
263: its - iteration number
264: fnorm - 2-norm function value (may be estimated)
265: ctx - optional user-defined context for private data for the
266: monitor routine, as set by SNESMonitorSet()
267: */
268: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
269: {
273: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
274: return(0);
275: }