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: PetscErrorCode ierr;
38: DM da;
39: PetscReal ftime;
40: SNES ts_snes;
41: PetscBool usemonitor = PETSC_TRUE;
42: PetscViewerAndFormat *vf;
44: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: Initialize program
46: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
48: PetscOptionsGetBool(NULL,NULL,"-usemonitor",&usemonitor,NULL);
50: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: Create distributed array (DMDA) to manage parallel grid and vectors
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
54: DMSetFromOptions(da);
55: DMSetUp(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: TSSetMaxTime(ts,1.0);
75: if (usemonitor) {
76: TSMonitorSet(ts,MyTSMonitor,0,0);
77: }
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Customize nonlinear solver
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: TSSetType(ts,TSBEULER);
83: TSGetSNES(ts,&ts_snes);
84: if (usemonitor) {
85: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
86: SNESMonitorSet(ts_snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void *))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
87: }
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Set initial conditions
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: FormInitialSolution(da,x);
92: TSSetTimeStep(ts,.0001);
93: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
94: TSSetSolution(ts,x);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Set runtime options
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: TSSetFromOptions(ts);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Solve nonlinear system
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: TSSolve(ts,x);
105: TSGetSolveTime(ts,&ftime);
106: TSGetStepNumber(ts,&steps);
107: VecViewFromOptions(x,NULL,"-final_sol");
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Free work space. All PETSc objects should be destroyed when they
111: are no longer needed.
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: VecDestroy(&x);
114: VecDestroy(&r);
115: TSDestroy(&ts);
116: DMDestroy(&da);
118: PetscFinalize();
119: return ierr;
120: }
121: /* ------------------------------------------------------------------- */
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: /* ------------------------------------------------------------------- */
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,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 = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5));
232: if (r < .125) {
233: u[j][i][0] = PetscExpReal(-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: }
249: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
250: {
252: PetscReal norm;
253: MPI_Comm comm;
256: VecNorm(v,NORM_2,&norm);
257: PetscObjectGetComm((PetscObject)ts,&comm);
258: if (step > -1) { /* -1 is used to indicate an interpolated value */
259: PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
260: }
261: return(0);
262: }
264: /*
265: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
266: Input Parameters:
267: snes - the SNES context
268: its - iteration number
269: fnorm - 2-norm function value (may be estimated)
270: ctx - optional user-defined context for private data for the
271: monitor routine, as set by SNESMonitorSet()
272: */
273: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf)
274: {
278: SNESMonitorDefaultShort(snes,its,fnorm,vf);
279: return(0);
280: }
281: /*TEST
283: test:
284: args: -da_grid_x 20 -ts_max_time 3 -ts_dt 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -ksp_monitor_short
285: requires: !single
287: test:
288: suffix: 2
289: args: -da_grid_x 20 -ts_max_time 0.11 -ts_dt 1e-1 -ts_type glle -ts_monitor -ksp_monitor_short
290: requires: !single
292: test:
293: suffix: glvis_da_2d_vect
294: 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
295: requires: !single
297: test:
298: suffix: glvis_da_2d_vect_ll
299: 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
300: requires: !single
302: TEST*/