Actual source code: ex12.c
petsc-3.7.5 2017-01-01
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,PetscViewerAndFormat*);
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 */
42: PetscErrorCode ierr;
43: DM da;
44: PetscReal ftime;
45: SNES ts_snes;
46: PetscViewerAndFormat *vf;
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Initialize program
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: PetscInitialize(&argc,&argv,(char*)0,help);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create distributed array (DMDA) to manage parallel grid and vectors
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
57: 2,1,NULL,NULL,&da);
58: DMDASetFieldName(da,0,"u");
59: DMDASetFieldName(da,1,"v");
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Extract global vectors from DMDA; then duplicate for remaining
63: vectors that are the same types
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: DMCreateGlobalVector(da,&x);
66: VecDuplicate(x,&r);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create timestepping solver context
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: TSCreate(PETSC_COMM_WORLD,&ts);
72: TSSetDM(ts,da);
73: TSSetProblemType(ts,TS_NONLINEAR);
74: TSSetRHSFunction(ts,NULL,FormFunction,da);
76: TSSetDuration(ts,maxsteps,1.0);
77: TSMonitorSet(ts,MyTSMonitor,0,0);
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Customize nonlinear solver
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: TSSetType(ts,TSBEULER);
83: TSGetSNES(ts,&ts_snes);
84: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
85: SNESMonitorSet(ts_snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void *))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Set initial conditions
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: FormInitialSolution(da,x);
91: TSSetInitialTimeStep(ts,0.0,.0001);
92: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
93: TSSetSolution(ts,x);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Set runtime options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: TSSetFromOptions(ts);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Solve nonlinear system
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: TSSolve(ts,x);
104: TSGetSolveTime(ts,&ftime);
105: TSGetTimeStepNumber(ts,&steps);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Free work space. All PETSc objects should be destroyed when they
109: are no longer needed.
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: VecDestroy(&x);
112: VecDestroy(&r);
113: TSDestroy(&ts);
114: DMDestroy(&da);
116: PetscFinalize();
117: return(0);
118: }
119: /* ------------------------------------------------------------------- */
122: /*
123: FormFunction - Evaluates nonlinear function, F(x).
125: Input Parameters:
126: . ts - the TS context
127: . X - input vector
128: . ptr - optional user-defined context, as set by SNESSetFunction()
130: Output Parameter:
131: . F - function vector
132: */
133: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
134: {
135: DM da = (DM)ptr;
137: PetscInt i,j,Mx,My,xs,ys,xm,ym;
138: PetscReal hx,hy,/*hxdhy,hydhx,*/ sx,sy;
139: PetscScalar u,uxx,uyy,v,***x,***f;
140: Vec localX;
143: DMGetLocalVector(da,&localX);
144: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
146: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
147: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
148: /*hxdhy = hx/hy;*/
149: /*hydhx = hy/hx;*/
151: /*
152: Scatter ghost points to local vector,using the 2-step process
153: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
154: By placing code between these two statements, computations can be
155: done while messages are in transition.
156: */
157: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
158: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
160: /*
161: Get pointers to vector data
162: */
163: DMDAVecGetArrayDOF(da,localX,&x);
164: DMDAVecGetArrayDOF(da,F,&f);
166: /*
167: Get local grid boundaries
168: */
169: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
171: /*
172: Compute function over the locally owned part of the grid
173: */
174: for (j=ys; j<ys+ym; j++) {
175: for (i=xs; i<xs+xm; i++) {
176: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
177: f[j][i][0] = x[j][i][0];
178: f[j][i][1] = x[j][i][1];
179: continue;
180: }
181: u = x[j][i][0];
182: v = x[j][i][1];
183: uxx = (-2.0*u + x[j][i-1][0] + x[j][i+1][0])*sx;
184: uyy = (-2.0*u + x[j-1][i][0] + x[j+1][i][0])*sy;
185: f[j][i][0] = v;
186: f[j][i][1] = uxx + uyy;
187: }
188: }
190: /*
191: Restore vectors
192: */
193: DMDAVecRestoreArrayDOF(da,localX,&x);
194: DMDAVecRestoreArrayDOF(da,F,&f);
195: DMRestoreLocalVector(da,&localX);
196: PetscLogFlops(11.0*ym*xm);
197: return(0);
198: }
200: /* ------------------------------------------------------------------- */
203: PetscErrorCode FormInitialSolution(DM da,Vec U)
204: {
206: PetscInt i,j,xs,ys,xm,ym,Mx,My;
207: PetscScalar ***u;
208: PetscReal hx,hy,x,y,r;
211: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
212: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
214: hx = 1.0/(PetscReal)(Mx-1);
215: hy = 1.0/(PetscReal)(My-1);
217: /*
218: Get pointers to vector data
219: */
220: DMDAVecGetArrayDOF(da,U,&u);
222: /*
223: Get local grid boundaries
224: */
225: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
227: /*
228: Compute function over the locally owned part of the grid
229: */
230: for (j=ys; j<ys+ym; j++) {
231: y = j*hy;
232: for (i=xs; i<xs+xm; i++) {
233: x = i*hx;
234: r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5));
235: if (r < .125) {
236: u[j][i][0] = PetscExpReal(-30.0*r*r*r);
237: u[j][i][1] = 0.0;
238: } else {
239: u[j][i][0] = 0.0;
240: u[j][i][1] = 0.0;
241: }
242: }
243: }
245: /*
246: Restore vectors
247: */
248: DMDAVecRestoreArrayDOF(da,U,&u);
249: return(0);
250: }
254: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
255: {
257: PetscReal norm;
258: MPI_Comm comm;
261: VecNorm(v,NORM_2,&norm);
262: PetscObjectGetComm((PetscObject)ts,&comm);
263: if (step > -1) { /* -1 is used to indicate an interpolated value */
264: PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
265: }
266: return(0);
267: }
271: /*
272: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
273: Input Parameters:
274: snes - the SNES context
275: its - iteration number
276: fnorm - 2-norm function value (may be estimated)
277: ctx - optional user-defined context for private data for the
278: monitor routine, as set by SNESMonitorSet()
279: */
280: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf)
281: {
285: SNESMonitorDefaultShort(snes,its,fnorm,vf);
286: return(0);
287: }