Actual source code: ex12.c
petsc-3.8.4 2018-03-24
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*);
35: int main(int argc,char **argv)
36: {
37: TS ts; /* nonlinear solver */
38: Vec x,r; /* solution, residual vectors */
39: PetscInt steps; /* iterations for convergence */
40: PetscErrorCode ierr;
41: DM da;
42: PetscReal ftime;
43: SNES ts_snes;
44: PetscViewerAndFormat *vf;
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Initialize program
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create distributed array (DMDA) to manage parallel grid and vectors
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
55: DMSetFromOptions(da);
56: DMSetUp(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: TSSetMaxTime(ts,1.0);
76: TSMonitorSet(ts,MyTSMonitor,0,0);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Customize nonlinear solver
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: TSSetType(ts,TSBEULER);
82: TSGetSNES(ts,&ts_snes);
83: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
84: SNESMonitorSet(ts_snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void *))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Set initial conditions
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: FormInitialSolution(da,x);
90: TSSetTimeStep(ts,.0001);
91: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
92: TSSetSolution(ts,x);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Set runtime options
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: TSSetFromOptions(ts);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Solve nonlinear system
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: TSSolve(ts,x);
103: TSGetSolveTime(ts,&ftime);
104: TSGetStepNumber(ts,&steps);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Free work space. All PETSc objects should be destroyed when they
108: are no longer needed.
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: VecDestroy(&x);
111: VecDestroy(&r);
112: TSDestroy(&ts);
113: DMDestroy(&da);
115: PetscFinalize();
116: return ierr;
117: }
118: /* ------------------------------------------------------------------- */
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,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: /* ------------------------------------------------------------------- */
198: PetscErrorCode FormInitialSolution(DM da,Vec U)
199: {
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: {
249: PetscReal norm;
250: MPI_Comm comm;
253: VecNorm(v,NORM_2,&norm);
254: PetscObjectGetComm((PetscObject)ts,&comm);
255: if (step > -1) { /* -1 is used to indicate an interpolated value */
256: PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);
257: }
258: return(0);
259: }
261: /*
262: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
263: Input Parameters:
264: snes - the SNES context
265: its - iteration number
266: fnorm - 2-norm function value (may be estimated)
267: ctx - optional user-defined context for private data for the
268: monitor routine, as set by SNESMonitorSet()
269: */
270: PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf)
271: {
275: SNESMonitorDefaultShort(snes,its,fnorm,vf);
276: return(0);
277: }