Actual source code: ex7.c
2: /* Program usage: mpirun -np <procs> ex5 [-help] [all PETSc options] */
4: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
7: /*
8: Include "petscda.h" so that we can use distributed arrays (DAs).
9: Include "petscts.h" so that we can use SNES solvers. Note that this
10: file automatically includes:
11: petsc.h - base PETSc routines petscvec.h - vectors
12: petscsys.h - system routines petscmat.h - matrices
13: petscis.h - index sets petscksp.h - Krylov subspace methods
14: petscviewer.h - viewers petscpc.h - preconditioners
15: petscksp.h - linear solvers
16: */
17: #include petscda.h
18: #include petscts.h
21: /*
22: User-defined routines
23: */
29: int main(int argc,char **argv)
30: {
31: TS ts; /* nonlinear solver */
32: Vec x,r; /* solution, residual vectors */
33: Mat J; /* Jacobian matrix */
34: PetscInt steps,maxsteps = 100; /* iterations for convergence */
35: PetscErrorCode ierr;
36: DA da;
37: MatFDColoring matfdcoloring;
38: ISColoring iscoloring;
39: PetscReal ftime;
41: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: Initialize program
43: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
45: PetscInitialize(&argc,&argv,(char *)0,help);
47: PetscOptionsGetInt(PETSC_NULL,"-max_steps",&maxsteps,PETSC_NULL);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create distributed array (DA) to manage parallel grid and vectors
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,
53: 1,1,PETSC_NULL,PETSC_NULL,&da);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Extract global vectors from DA; then duplicate for remaining
57: vectors that are the same types
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: DACreateGlobalVector(da,&x);
60: VecDuplicate(x,&r);
62: TSCreate(PETSC_COMM_WORLD,&ts);
63: TSSetProblemType(ts,TS_NONLINEAR);
64: TSSetRHSFunction(ts,FormFunction,da);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create matrix data structure; set Jacobian evaluation routine
70: Set Jacobian matrix data structure and default Jacobian evaluation
71: routine. User can override with:
72: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
73: (unless user explicitly sets preconditioner)
74: -snes_mf_operator : form preconditioning matrix as set by the user,
75: but use matrix-free approx for Jacobian-vector
76: products within Newton-Krylov method
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: DAGetColoring(da,IS_COLORING_GLOBAL,&iscoloring);
80: DAGetMatrix(da,MATAIJ,&J);
81: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
82: ISColoringDestroy(iscoloring);
83: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,da);
84: MatFDColoringSetFromOptions(matfdcoloring);
85: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Customize nonlinear solver; set runtime options
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: TSSetInitialTimeStep(ts,0.0,.0001);
90: TSSetType(ts,TS_BEULER);
91: TSSetDuration(ts,maxsteps,1.0);
92: TSSetFromOptions(ts);
93: TSMonitorSet(ts,Monitor,0,0);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Set initial conditions
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: FormInitialSolution(da,x);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Solve nonlinear system
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: TSSetSolution(ts,x);
104: TSStep(ts,&steps,&ftime);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Free work space. All PETSc objects should be destroyed when they
109: are no longer needed.
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: MatDestroy(J);
113: MatFDColoringDestroy(matfdcoloring);
114: VecDestroy(x);
115: VecDestroy(r);
116: TSDestroy(ts);
117: DADestroy(da);
118: PetscFinalize();
120: return(0);
121: }
122: /* ------------------------------------------------------------------- */
125: /*
126: FormFunction - Evaluates nonlinear function, F(x).
128: Input Parameters:
129: . ts - the TS context
130: . X - input vector
131: . ptr - optional user-defined context, as set by SNESSetFunction()
133: Output Parameter:
134: . F - function vector
135: */
136: PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr)
137: {
138: DA da = (DA)ptr;
140: PetscInt i,j,Mx,My,xs,ys,xm,ym;
141: PetscReal two = 2.0,hx,hy,hxdhy,hydhx,sx,sy;
142: PetscScalar u,uxx,uyy,**x,**f;
143: Vec localX;
146: DAGetLocalVector(da,&localX);
147: DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
148: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
150: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
151: hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
152: hxdhy = hx/hy;
153: hydhx = hy/hx;
155: /*
156: Scatter ghost points to local vector,using the 2-step process
157: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
158: By placing code between these two statements, computations can be
159: done while messages are in transition.
160: */
161: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
162: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
164: /*
165: Get pointers to vector data
166: */
167: DAVecGetArray(da,localX,&x);
168: DAVecGetArray(da,F,&f);
170: /*
171: Get local grid boundaries
172: */
173: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
175: /*
176: Compute function over the locally owned part of the grid
177: */
178: for (j=ys; j<ys+ym; j++) {
179: for (i=xs; i<xs+xm; i++) {
180: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
181: f[j][i] = x[j][i];
182: continue;
183: }
184: u = x[j][i];
185: uxx = (two*u - x[j][i-1] - x[j][i+1])*sx;
186: uyy = (two*u - x[j-1][i] - x[j+1][i])*sy;
187: /* f[j][i] = -(uxx + uyy); */
188: f[j][i] = -u*(uxx + uyy) - (4.0 - 1.0)*((x[j][i+1] - x[j][i-1])*(x[j][i+1] - x[j][i-1])*.25*sx +
189: (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
190: }
191: }
193: /*
194: Restore vectors
195: */
196: DAVecRestoreArray(da,localX,&x);
197: DAVecRestoreArray(da,F,&f);
198: DARestoreLocalVector(da,&localX);
199: PetscLogFlops(11*ym*xm);
200: return(0);
201: }
203: /* ------------------------------------------------------------------- */
206: PetscErrorCode FormInitialSolution(DA da,Vec U)
207: {
209: PetscInt i,j,xs,ys,xm,ym,Mx,My;
210: PetscScalar **u;
211: PetscReal hx,hy,x,y,r;
214: DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
215: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
217: hx = 1.0/(PetscReal)(Mx-1);
218: hy = 1.0/(PetscReal)(My-1);
220: /*
221: Get pointers to vector data
222: */
223: DAVecGetArray(da,U,&u);
225: /*
226: Get local grid boundaries
227: */
228: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
230: /*
231: Compute function over the locally owned part of the grid
232: */
233: for (j=ys; j<ys+ym; j++) {
234: y = j*hy;
235: for (i=xs; i<xs+xm; i++) {
236: x = i*hx;
237: r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
238: if (r < .125) {
239: u[j][i] = PetscExpScalar(-30.0*r*r*r);
240: } else {
241: u[j][i] = 0.0;
242: }
243: }
244: }
246: /*
247: Restore vectors
248: */
249: DAVecRestoreArray(da,U,&u);
250: return(0);
251: }
255: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx)
256: {
258: PetscReal norm;
259: MPI_Comm comm;
262: VecNorm(v,NORM_2,&norm);
263: PetscObjectGetComm((PetscObject)ts,&comm);
264: PetscPrintf(comm,"timestep %D time %G norm %G\n",step,ptime,norm);
265: return(0);
266: }