Actual source code: ex7.c
petsc-3.4.5 2014-06-29
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 <petscdmda.h>
16: #include <petscts.h>
19: /*
20: User-defined routines
21: */
22: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
23: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
24: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,void*);
28: int main(int argc,char **argv)
29: {
30: TS ts; /* time integrator */
31: SNES snes;
32: Vec x,r; /* solution, residual vectors */
33: PetscInt maxsteps = 100; /* iterations for convergence */
35: DM da;
37: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: Initialize program
39: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40: PetscInitialize(&argc,&argv,(char*)0,help);
42: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: Create distributed array (DMDA) to manage parallel grid and vectors
44: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
45: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
46: 1,1,NULL,NULL,&da);
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Extract global vectors from DMDA; then duplicate for remaining
50: vectors that are the same types
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: DMCreateGlobalVector(da,&x);
53: VecDuplicate(x,&r);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Create timestepping solver context
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: TSCreate(PETSC_COMM_WORLD,&ts);
59: TSSetProblemType(ts,TS_NONLINEAR);
60: TSSetRHSFunction(ts,NULL,FormFunction,da);
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create matrix data structure; set Jacobian evaluation routine
65: Set Jacobian matrix data structure and default Jacobian evaluation
66: routine. User can override with:
67: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
68: (unless user explicitly sets preconditioner)
69: -snes_mf_operator : form preconditioning matrix as set by the user,
70: but use matrix-free approx for Jacobian-vector
71: products within Newton-Krylov method
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: TSSetDuration(ts,maxsteps,1.0);
76: TSMonitorSet(ts,MyTSMonitor,0,0);
77: TSSetDM(ts,da);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Customize nonlinear solver
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: TSSetType(ts,TSBEULER);
82: TSGetSNES(ts,&snes);
83: SNESMonitorSet(snes,MySNESMonitor,NULL,NULL);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Set initial conditions
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: FormInitialSolution(da,x);
89: TSSetInitialTimeStep(ts,0.0,.0001);
90: TSSetSolution(ts,x);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Set runtime options
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: TSSetFromOptions(ts);
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Solve nonlinear system
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: TSSolve(ts,x);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Free work space. All PETSc objects should be destroyed when they
104: are no longer needed.
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: VecDestroy(&x);
107: VecDestroy(&r);
108: TSDestroy(&ts);
109: DMDestroy(&da);
111: PetscFinalize();
112: return(0);
113: }
114: /* ------------------------------------------------------------------- */
117: /*
118: FormFunction - Evaluates nonlinear function, F(x).
120: Input Parameters:
121: . ts - the TS context
122: . X - input vector
123: . ptr - optional user-defined context, as set by SNESSetFunction()
125: Output Parameter:
126: . F - function vector
127: */
128: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
129: {
130: DM da;
132: PetscInt i,j,Mx,My,xs,ys,xm,ym;
133: PetscReal two = 2.0,hx,hy,sx,sy;
134: PetscScalar u,uxx,uyy,**x,**f;
135: Vec localX;
138: TSGetDM(ts,&da);
139: DMGetLocalVector(da,&localX);
140: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
141: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
143: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
144: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
146: /*
147: Scatter ghost points to local vector,using the 2-step process
148: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
149: By placing code between these two statements, computations can be
150: done while messages are in transition.
151: */
152: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
153: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
155: /*
156: Get pointers to vector data
157: */
158: DMDAVecGetArray(da,localX,&x);
159: DMDAVecGetArray(da,F,&f);
161: /*
162: Get local grid boundaries
163: */
164: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
166: /*
167: Compute function over the locally owned part of the grid
168: */
169: for (j=ys; j<ys+ym; j++) {
170: for (i=xs; i<xs+xm; i++) {
171: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
172: f[j][i] = x[j][i];
173: continue;
174: }
175: u = x[j][i];
176: uxx = (two*u - x[j][i-1] - x[j][i+1])*sx;
177: uyy = (two*u - x[j-1][i] - x[j+1][i])*sy;
178: /* f[j][i] = -(uxx + uyy); */
179: 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 +
180: (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
181: }
182: }
184: /*
185: Restore vectors
186: */
187: DMDAVecRestoreArray(da,localX,&x);
188: DMDAVecRestoreArray(da,F,&f);
189: DMRestoreLocalVector(da,&localX);
190: PetscLogFlops(11.0*ym*xm);
191: return(0);
192: }
194: /* ------------------------------------------------------------------- */
197: PetscErrorCode FormInitialSolution(DM da,Vec U)
198: {
200: PetscInt i,j,xs,ys,xm,ym,Mx,My;
201: PetscScalar **u;
202: PetscReal hx,hy,x,y,r;
205: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
206: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
208: hx = 1.0/(PetscReal)(Mx-1);
209: hy = 1.0/(PetscReal)(My-1);
211: /*
212: Get pointers to vector data
213: */
214: DMDAVecGetArray(da,U,&u);
216: /*
217: Get local grid boundaries
218: */
219: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
221: /*
222: Compute function over the locally owned part of the grid
223: */
224: for (j=ys; j<ys+ym; j++) {
225: y = j*hy;
226: for (i=xs; i<xs+xm; i++) {
227: x = i*hx;
228: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
229: if (r < .125) u[j][i] = PetscExpScalar(-30.0*r*r*r);
230: else u[j][i] = 0.0;
231: }
232: }
234: /*
235: Restore vectors
236: */
237: DMDAVecRestoreArray(da,U,&u);
238: return(0);
239: }
243: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
244: {
246: PetscReal norm;
247: MPI_Comm comm;
250: VecNorm(v,NORM_2,&norm);
251: PetscObjectGetComm((PetscObject)ts,&comm);
252: PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
253: return(0);
254: }
258: /*
259: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
260: Input Parameters:
261: snes - the SNES context
262: its - iteration number
263: fnorm - 2-norm function value (may be estimated)
264: ctx - optional user-defined context for private data for the
265: monitor routine, as set by SNESMonitorSet()
266: */
267: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
268: {
272: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
273: return(0);
274: }