Actual source code: ex12.c
petsc-3.6.1 2015-08-06
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 <petscdm.h>
24: #include <petscdmda.h>
25: #include <petscts.h>
28: /*
29: User-defined routines
30: */
31: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
32: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
33: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,void*);
37: int main(int argc,char **argv)
38: {
39: TS ts; /* nonlinear solver */
40: Vec x,r; /* solution, residual vectors */
41: PetscInt steps,maxsteps = 100; /* iterations for convergence */
43: DM da;
44: PetscReal ftime;
45: SNES ts_snes;
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Initialize program
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: PetscInitialize(&argc,&argv,(char*)0,help);
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Create distributed array (DMDA) to manage parallel grid and vectors
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
56: 2,1,NULL,NULL,&da);
57: DMDASetFieldName(da,0,"u");
58: DMDASetFieldName(da,1,"v");
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Extract global vectors from DMDA; then duplicate for remaining
62: vectors that are the same types
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: DMCreateGlobalVector(da,&x);
65: VecDuplicate(x,&r);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create timestepping solver context
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: TSCreate(PETSC_COMM_WORLD,&ts);
71: TSSetDM(ts,da);
72: TSSetProblemType(ts,TS_NONLINEAR);
73: TSSetRHSFunction(ts,NULL,FormFunction,da);
75: TSSetDuration(ts,maxsteps,1.0);
76: TSMonitorSet(ts,MyTSMonitor,0,0);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Customize nonlinear solver
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: TSSetType(ts,TSBEULER);
82: TSGetSNES(ts,&ts_snes);
83: SNESMonitorSet(ts_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);
101: TSGetSolveTime(ts,&ftime);
102: TSGetTimeStepNumber(ts,&steps);
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Free work space. All PETSc objects should be destroyed when they
106: are no longer needed.
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: VecDestroy(&x);
109: VecDestroy(&r);
110: TSDestroy(&ts);
111: DMDestroy(&da);
113: PetscFinalize();
114: return(0);
115: }
116: /* ------------------------------------------------------------------- */
119: /*
120: FormFunction - Evaluates nonlinear function, F(x).
122: Input Parameters:
123: . ts - the TS context
124: . X - input vector
125: . ptr - optional user-defined context, as set by SNESSetFunction()
127: Output Parameter:
128: . F - function vector
129: */
130: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
131: {
132: DM da = (DM)ptr;
134: PetscInt i,j,Mx,My,xs,ys,xm,ym;
135: PetscReal hx,hy,/*hxdhy,hydhx,*/ sx,sy;
136: PetscScalar u,uxx,uyy,v,***x,***f;
137: Vec localX;
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);
146: /*hxdhy = hx/hy;*/
147: /*hydhx = hy/hx;*/
149: /*
150: Scatter ghost points to local vector,using the 2-step process
151: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
152: By placing code between these two statements, computations can be
153: done while messages are in transition.
154: */
155: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
156: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
158: /*
159: Get pointers to vector data
160: */
161: DMDAVecGetArrayDOF(da,localX,&x);
162: DMDAVecGetArrayDOF(da,F,&f);
164: /*
165: Get local grid boundaries
166: */
167: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
169: /*
170: Compute function over the locally owned part of the grid
171: */
172: for (j=ys; j<ys+ym; j++) {
173: for (i=xs; i<xs+xm; i++) {
174: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
175: f[j][i][0] = x[j][i][0];
176: f[j][i][1] = x[j][i][1];
177: continue;
178: }
179: u = x[j][i][0];
180: v = x[j][i][1];
181: uxx = (-2.0*u + x[j][i-1][0] + x[j][i+1][0])*sx;
182: uyy = (-2.0*u + x[j-1][i][0] + x[j+1][i][0])*sy;
183: f[j][i][0] = v;
184: f[j][i][1] = uxx + uyy;
185: }
186: }
188: /*
189: Restore vectors
190: */
191: DMDAVecRestoreArrayDOF(da,localX,&x);
192: DMDAVecRestoreArrayDOF(da,F,&f);
193: DMRestoreLocalVector(da,&localX);
194: PetscLogFlops(11.0*ym*xm);
195: return(0);
196: }
198: /* ------------------------------------------------------------------- */
201: PetscErrorCode FormInitialSolution(DM da,Vec U)
202: {
204: PetscInt i,j,xs,ys,xm,ym,Mx,My;
205: PetscScalar ***u;
206: PetscReal hx,hy,x,y,r;
209: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
210: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
212: hx = 1.0/(PetscReal)(Mx-1);
213: hy = 1.0/(PetscReal)(My-1);
215: /*
216: Get pointers to vector data
217: */
218: DMDAVecGetArrayDOF(da,U,&u);
220: /*
221: Get local grid boundaries
222: */
223: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
225: /*
226: Compute function over the locally owned part of the grid
227: */
228: for (j=ys; j<ys+ym; j++) {
229: y = j*hy;
230: for (i=xs; i<xs+xm; i++) {
231: x = i*hx;
232: r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5));
233: if (r < .125) {
234: u[j][i][0] = PetscExpReal(-30.0*r*r*r);
235: u[j][i][1] = 0.0;
236: } else {
237: u[j][i][0] = 0.0;
238: u[j][i][1] = 0.0;
239: }
240: }
241: }
243: /*
244: Restore vectors
245: */
246: DMDAVecRestoreArrayDOF(da,U,&u);
247: return(0);
248: }
252: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
253: {
255: PetscReal norm;
256: MPI_Comm comm;
259: VecNorm(v,NORM_2,&norm);
260: PetscObjectGetComm((PetscObject)ts,&comm);
261: PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
262: return(0);
263: }
267: /*
268: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
269: Input Parameters:
270: snes - the SNES context
271: its - iteration number
272: fnorm - 2-norm function value (may be estimated)
273: ctx - optional user-defined context for private data for the
274: monitor routine, as set by SNESMonitorSet()
275: */
276: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
277: {
281: SNESMonitorDefaultShort(snes,its,fnorm,ctx);
282: return(0);
283: }