Actual source code: ex12.c
1: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
2: /*
3: Solves the equation
5: u_tt - \Delta u = 0
7: which we split into two first-order systems
9: u_t - v = 0 so that F(u,v).u = v
10: v_t - \Delta u = 0 so that F(u,v).v = Delta u
12: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
13: Include "petscts.h" so that we can use SNES solvers. Note that this
14: file automatically includes:
15: petscsys.h - base PETSc routines petscvec.h - vectors
16: petscmat.h - matrices
17: petscis.h - index sets petscksp.h - Krylov subspace methods
18: petscviewer.h - viewers petscpc.h - preconditioners
19: petscksp.h - linear solvers
20: */
21: #include <petscdm.h>
22: #include <petscdmda.h>
23: #include <petscts.h>
25: /*
26: User-defined routines
27: */
28: extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec);
29: extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
30: extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,PetscViewerAndFormat*);
32: int main(int argc,char **argv)
33: {
34: TS ts; /* nonlinear solver */
35: Vec x,r; /* solution, residual vectors */
36: PetscInt steps; /* iterations for convergence */
37: DM da;
38: PetscReal ftime;
39: SNES ts_snes;
40: PetscBool usemonitor = PETSC_TRUE;
41: PetscViewerAndFormat *vf;
43: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Initialize program
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46: PetscInitialize(&argc,&argv,(char*)0,help);
47: PetscOptionsGetBool(NULL,NULL,"-usemonitor",&usemonitor,NULL);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create distributed array (DMDA) to manage parallel grid and vectors
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
53: DMSetFromOptions(da);
54: DMSetUp(da);
55: DMDASetFieldName(da,0,"u");
56: DMDASetFieldName(da,1,"v");
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Extract global vectors from DMDA; then duplicate for remaining
60: vectors that are the same types
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62: DMCreateGlobalVector(da,&x);
63: VecDuplicate(x,&r);
65: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: Create timestepping solver context
67: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: TSCreate(PETSC_COMM_WORLD,&ts);
69: TSSetDM(ts,da);
70: TSSetProblemType(ts,TS_NONLINEAR);
71: TSSetRHSFunction(ts,NULL,FormFunction,da);
73: TSSetMaxTime(ts,1.0);
74: if (usemonitor) {
75: TSMonitorSet(ts,MyTSMonitor,0,0);
76: }
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Customize nonlinear solver
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: TSSetType(ts,TSBEULER);
82: TSGetSNES(ts,&ts_snes);
83: if (usemonitor) {
84: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
85: SNESMonitorSet(ts_snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void *))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
86: }
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Set initial conditions
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: FormInitialSolution(da,x);
91: TSSetTimeStep(ts,.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: TSGetStepNumber(ts,&steps);
106: VecViewFromOptions(x,NULL,"-final_sol");
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Free work space. All PETSc objects should be destroyed when they
110: are no longer needed.
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: VecDestroy(&x);
113: VecDestroy(&r);
114: TSDestroy(&ts);
115: DMDestroy(&da);
117: PetscFinalize();
118: return 0;
119: }
120: /* ------------------------------------------------------------------- */
121: /*
122: FormFunction - Evaluates nonlinear function, F(x).
124: Input Parameters:
125: . ts - the TS context
126: . X - input vector
127: . ptr - optional user-defined context, as set by SNESSetFunction()
129: Output Parameter:
130: . F - function vector
131: */
132: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
133: {
134: DM da = (DM)ptr;
135: PetscInt i,j,Mx,My,xs,ys,xm,ym;
136: PetscReal hx,hy,/*hxdhy,hydhx,*/ sx,sy;
137: PetscScalar u,uxx,uyy,v,***x,***f;
138: Vec localX;
141: DMGetLocalVector(da,&localX);
142: 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);
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: /* ------------------------------------------------------------------- */
199: PetscErrorCode FormInitialSolution(DM da,Vec U)
200: {
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,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: DMDAVecGetArrayDOF(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 = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5));
229: if (r < .125) {
230: u[j][i][0] = PetscExpReal(-30.0*r*r*r);
231: u[j][i][1] = 0.0;
232: } else {
233: u[j][i][0] = 0.0;
234: u[j][i][1] = 0.0;
235: }
236: }
237: }
239: /*
240: Restore vectors
241: */
242: DMDAVecRestoreArrayDOF(da,U,&u);
243: return 0;
244: }
246: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
247: {
248: PetscReal norm;
249: MPI_Comm comm;
252: VecNorm(v,NORM_2,&norm);
253: PetscObjectGetComm((PetscObject)ts,&comm);
254: if (step > -1) { /* -1 is used to indicate an interpolated value */
255: PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
256: }
257: return 0;
258: }
260: /*
261: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
262: Input Parameters:
263: snes - the SNES context
264: its - iteration number
265: fnorm - 2-norm function value (may be estimated)
266: ctx - optional user-defined context for private data for the
267: monitor routine, as set by SNESMonitorSet()
268: */
269: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf)
270: {
272: SNESMonitorDefaultShort(snes,its,fnorm,vf);
273: return 0;
274: }
275: /*TEST
277: test:
278: args: -da_grid_x 20 -ts_max_time 3 -ts_dt 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -ksp_monitor_short
279: requires: !single
281: test:
282: suffix: 2
283: args: -da_grid_x 20 -ts_max_time 0.11 -ts_dt 1e-1 -ts_type glle -ts_monitor -ksp_monitor_short
284: requires: !single
286: test:
287: suffix: glvis_da_2d_vect
288: args: -usemonitor 0 -da_grid_x 20 -ts_max_time 0.3 -ts_dt 1e-1 -ts_type glle -final_sol glvis: -viewer_glvis_dm_da_bs 2,0
289: requires: !single
291: test:
292: suffix: glvis_da_2d_vect_ll
293: args: -usemonitor 0 -da_grid_x 20 -ts_max_time 0.3 -ts_dt 1e-1 -ts_type glle -final_sol glvis: -viewer_glvis_dm_da_bs 2,0 -viewer_glvis_dm_da_ll
294: requires: !single
296: TEST*/