Actual source code: ex12.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
4: /*
5: Solves the equation
7: u_tt - \Delta u = 0
9: which we split into two first-order systems
11: u_t - v = 0 so that F(u,v).u = v
12: v_t - \Delta u = 0 so that F(u,v).v = Delta u
14: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
15: Include "petscts.h" so that we can use SNES solvers. Note that this
16: file automatically includes:
17: petscsys.h - base PETSc routines petscvec.h - vectors
18: petscmat.h - matrices
19: petscis.h - index sets petscksp.h - Krylov subspace methods
20: petscviewer.h - viewers petscpc.h - preconditioners
21: petscksp.h - linear solvers
22: */
23: #include <petscdmda.h>
24: #include <petscts.h>
27: /*
28: User-defined routines
29: */
30: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
31: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
32: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,void*);
36: int main(int argc,char **argv)
37: {
38: TS ts; /* nonlinear solver */
39: Vec x,r; /* solution, residual vectors */
40: PetscInt steps,maxsteps = 100; /* iterations for convergence */
42: DM da;
43: PetscReal ftime;
44: SNES ts_snes;
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Initialize program
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: PetscInitialize(&argc,&argv,(char*)0,help);
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create distributed array (DMDA) to manage parallel grid and vectors
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
55: 2,1,NULL,NULL,&da);
56: DMDASetFieldName(da,0,"u");
57: DMDASetFieldName(da,1,"v");
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Extract global vectors from DMDA; then duplicate for remaining
61: vectors that are the same types
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: DMCreateGlobalVector(da,&x);
64: VecDuplicate(x,&r);
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create timestepping solver context
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69: TSCreate(PETSC_COMM_WORLD,&ts);
70: TSSetDM(ts,da);
71: TSSetProblemType(ts,TS_NONLINEAR);
72: TSSetRHSFunction(ts,NULL,FormFunction,da);
74: TSSetDuration(ts,maxsteps,1.0);
75: TSMonitorSet(ts,MyTSMonitor,0,0);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Customize nonlinear solver
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: TSSetType(ts,TSBEULER);
81: TSGetSNES(ts,&ts_snes);
82: SNESMonitorSet(ts_snes,MySNESMonitor,NULL,NULL);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Set initial conditions
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: FormInitialSolution(da,x);
88: TSSetInitialTimeStep(ts,0.0,.0001);
89: TSSetSolution(ts,x);
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Set runtime options
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: TSSetFromOptions(ts);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Solve nonlinear system
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: TSSolve(ts,x);
100: TSGetSolveTime(ts,&ftime);
101: TSGetTimeStepNumber(ts,&steps);
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 = (DM)ptr;
133: PetscInt i,j,Mx,My,xs,ys,xm,ym;
134: PetscReal hx,hy,/*hxdhy,hydhx,*/ sx,sy;
135: PetscScalar u,uxx,uyy,v,***x,***f;
136: Vec localX;
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);
145: /*hxdhy = hx/hy;*/
146: /*hydhx = hy/hx;*/
148: /*
149: Scatter ghost points to local vector,using the 2-step process
150: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
151: By placing code between these two statements, computations can be
152: done while messages are in transition.
153: */
154: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
155: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
157: /*
158: Get pointers to vector data
159: */
160: DMDAVecGetArrayDOF(da,localX,&x);
161: DMDAVecGetArrayDOF(da,F,&f);
163: /*
164: Get local grid boundaries
165: */
166: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
168: /*
169: Compute function over the locally owned part of the grid
170: */
171: for (j=ys; j<ys+ym; j++) {
172: for (i=xs; i<xs+xm; i++) {
173: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
174: f[j][i][0] = x[j][i][0];
175: f[j][i][1] = x[j][i][1];
176: continue;
177: }
178: u = x[j][i][0];
179: v = x[j][i][1];
180: uxx = (-2.0*u + x[j][i-1][0] + x[j][i+1][0])*sx;
181: uyy = (-2.0*u + x[j-1][i][0] + x[j+1][i][0])*sy;
182: f[j][i][0] = v;
183: f[j][i][1] = uxx + uyy;
184: }
185: }
187: /*
188: Restore vectors
189: */
190: DMDAVecRestoreArrayDOF(da,localX,&x);
191: DMDAVecRestoreArrayDOF(da,F,&f);
192: DMRestoreLocalVector(da,&localX);
193: PetscLogFlops(11.0*ym*xm);
194: return(0);
195: }
197: /* ------------------------------------------------------------------- */
200: PetscErrorCode FormInitialSolution(DM da,Vec U)
201: {
203: PetscInt i,j,xs,ys,xm,ym,Mx,My;
204: PetscScalar ***u;
205: PetscReal hx,hy,x,y,r;
208: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
209: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
211: hx = 1.0/(PetscReal)(Mx-1);
212: hy = 1.0/(PetscReal)(My-1);
214: /*
215: Get pointers to vector data
216: */
217: DMDAVecGetArrayDOF(da,U,&u);
219: /*
220: Get local grid boundaries
221: */
222: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
224: /*
225: Compute function over the locally owned part of the grid
226: */
227: for (j=ys; j<ys+ym; j++) {
228: y = j*hy;
229: for (i=xs; i<xs+xm; i++) {
230: x = i*hx;
231: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
232: if (r < .125) {
233: u[j][i][0] = PetscExpScalar(-30.0*r*r*r);
234: u[j][i][1] = 0.0;
235: } else {
236: u[j][i][0] = 0.0;
237: u[j][i][1] = 0.0;
238: }
239: }
240: }
242: /*
243: Restore vectors
244: */
245: DMDAVecRestoreArrayDOF(da,U,&u);
246: return(0);
247: }
251: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
252: {
254: PetscReal norm;
255: MPI_Comm comm;
258: VecNorm(v,NORM_2,&norm);
259: PetscObjectGetComm((PetscObject)ts,&comm);
260: PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
261: return(0);
262: }
266: /*
267: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
268: Input Parameters:
269: snes - the SNES context
270: its - iteration number
271: fnorm - 2-norm function value (may be estimated)
272: ctx - optional user-defined context for private data for the
273: monitor routine, as set by SNESMonitorSet()
274: */
275: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
276: {
280: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
281: return(0);
282: }