Actual source code: ex17.c
petsc-3.3-p7 2013-05-11
1: static const char help[] = "Time-dependent PDE in 1d. Simplified from ex15.c for illustrating how to solve DAEs. \n";
2: /*
3: u_t = uxx
4: 0 < x < 1;
5: At t=0: u(x) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5)) < .125
6: u(x) = 0.0 if r >= .125
9: Boundary conditions:
10: Dirichlet BC:
11: At x=0, x=1, u = 0.0
13: Neumann BC:
14: At x=0, x=1: du(x,t)/dx = 0
16: Program usage:
17: mpiexec -n <procs> ./ex17 [-help] [all PETSc options]
18: e.g., mpiexec -n 2 ./ex17 -da_grid_x 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
19: ./ex17 -da_grid_x 40 -drawcontours -draw_pause .1
20: ./ex17 -da_grid_x 100 -drawcontours -draw_pause .1 -ts_type theta -ts_theta_theta 0.5 # Midpoint is not L-stable
21: ./ex17 -jac_type fd_coloring -drawcontours -draw_pause .1 -da_grid_x 500 -boundary 1
22: ./ex17 -da_grid_x 100 -drawcontours -draw_pause 1 -ts_type gl -ts_adapt_type none -ts_max_steps 2
23: */
25: #include <petscdmda.h>
26: #include <petscts.h>
28: typedef enum {JACOBIAN_ANALYTIC,JACOBIAN_FD_COLORING,JACOBIAN_FD_FULL} JacobianType;
29: static const char *const JacobianTypes[] = {"analytic","fd_coloring","fd_full","JacobianType","fd_",0};
31: /*
32: User-defined data structures and routines
33: */
34: typedef struct {
35: PetscBool drawcontours;
36: } MonitorCtx;
38: typedef struct {
39: PetscReal c;
40: PetscInt boundary; /* Type of boundary condition */
41: PetscBool viewJacobian;
42: } AppCtx;
44: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
45: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);
46: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
47: static PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
51: int main(int argc,char **argv)
52: {
53: TS ts; /* nonlinear solver */
54: Vec u; /* solution, residual vectors */
55: Mat J; /* Jacobian matrix */
56: PetscInt steps,maxsteps = 1000; /* iterations for convergence */
58: DM da;
59: MatFDColoring matfdcoloring = PETSC_NULL;
60: PetscReal ftime,dt;
61: MonitorCtx usermonitor; /* user-defined monitor context */
62: AppCtx user; /* user-defined work context */
63: JacobianType jacType;
65: PetscInitialize(&argc,&argv,(char *)0,help);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create distributed array (DMDA) to manage parallel grid and vectors
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-11,1,1,PETSC_NULL,&da);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Extract global vectors from DMDA;
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: DMCreateGlobalVector(da,&u);
77: /* Initialize user application context */
78: user.c = -30.0;
79: user.boundary = 0; /* 0: Dirichlet BC; 1: Neumann BC */
80: user.viewJacobian = PETSC_FALSE;
81: PetscOptionsGetInt(PETSC_NULL,"-boundary",&user.boundary,PETSC_NULL);
82: PetscOptionsHasName(PETSC_NULL,"-viewJacobian",&user.viewJacobian);
84: usermonitor.drawcontours = PETSC_FALSE;
85: PetscOptionsHasName(PETSC_NULL,"-drawcontours",&usermonitor.drawcontours);
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Create timestepping solver context
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: TSCreate(PETSC_COMM_WORLD,&ts);
91: TSSetProblemType(ts,TS_NONLINEAR);
92: TSSetType(ts,TSTHETA);
93: TSThetaSetTheta(ts,1.0); /* Make the Theta method behave like backward Euler */
94: TSSetIFunction(ts,PETSC_NULL,FormIFunction,&user);
96: DMCreateMatrix(da,MATAIJ,&J);
97: jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */
98: TSSetIJacobian(ts,J,J,FormIJacobian,&user);
100: TSSetDM(ts,da); /* Use TSGetDM() to access. Setting here allows easy use of geometric multigrid. */
102: ftime = 1.0;
103: TSSetDuration(ts,maxsteps,ftime);
104: TSMonitorSet(ts,MyTSMonitor,&usermonitor,PETSC_NULL);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Set initial conditions
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: FormInitialSolution(ts,u,&user);
110: TSSetSolution(ts,u);
111: dt = .01;
112: TSSetInitialTimeStep(ts,0.0,dt);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Set runtime options
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: TSSetFromOptions(ts);
119: /* Use slow fd Jacobian or fast fd Jacobian with colorings.
120: Note: this requirs snes which is not created until TSSetUp()/TSSetFromOptions() is called */
121: PetscOptionsBegin(((PetscObject)da)->comm,PETSC_NULL,"Options for Jacobian evaluation",PETSC_NULL);
122: PetscOptionsEnum("-jac_type","Type of Jacobian","",JacobianTypes,(PetscEnum)jacType,(PetscEnum*)&jacType,0);
123: PetscOptionsEnd();
124: if (jacType == JACOBIAN_FD_COLORING) {
125: SNES snes;
126: ISColoring iscoloring;
127: DMCreateColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
128: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
129: MatFDColoringSetFromOptions(matfdcoloring);
130: ISColoringDestroy(&iscoloring);
131: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))SNESTSFormFunction,ts);
132: TSGetSNES(ts,&snes);
133: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
134: } else if (jacType == JACOBIAN_FD_FULL){
135: SNES snes;
136: TSGetSNES(ts,&snes);
137: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobian,&user);
138: }
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Solve nonlinear system
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: TSSolve(ts,u,&ftime);
144: TSGetTimeStepNumber(ts,&steps);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Free work space.
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: MatDestroy(&J);
150: if (matfdcoloring) {MatFDColoringDestroy(&matfdcoloring);}
151: VecDestroy(&u);
152: TSDestroy(&ts);
153: DMDestroy(&da);
155: PetscFinalize();
156: return(0);
157: }
158: /* ------------------------------------------------------------------- */
161: static PetscErrorCode FormIFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
162: {
163: AppCtx *user=(AppCtx*)ptr;
164: DM da;
166: PetscInt i,Mx,xs,xm;
167: PetscReal hx,sx;
168: PetscScalar *u,*udot,*f;
169: Vec localU;
172: TSGetDM(ts,&da);
173: DMGetLocalVector(da,&localU);
174: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
175: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
177: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
179: /*
180: Scatter ghost points to local vector,using the 2-step process
181: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
182: By placing code between these two statements, computations can be
183: done while messages are in transition.
184: */
185: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
186: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
188: /* Get pointers to vector data */
189: DMDAVecGetArray(da,localU,&u);
190: DMDAVecGetArray(da,Udot,&udot);
191: DMDAVecGetArray(da,F,&f);
193: /* Get local grid boundaries */
194: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
196: /* Compute function over the locally owned part of the grid */
197: for (i=xs; i<xs+xm; i++) {
198: if (user->boundary == 0) { /* Dirichlet BC */
199: if (i == 0 || i == Mx-1) {
200: f[i] = u[i]; /* F = U */
201: } else {
202: f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
203: }
204: } else { /* Neumann BC */
205: if (i == 0) {
206: f[i] = u[0] - u[1];
207: } else if (i == Mx-1) {
208: f[i] = u[i] - u[i-1];
209: } else {
210: f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
211: }
212: }
213: }
215: /* Restore vectors */
216: DMDAVecRestoreArray(da,localU,&u);
217: DMDAVecRestoreArray(da,Udot,&udot);
218: DMDAVecRestoreArray(da,F,&f);
219: DMRestoreLocalVector(da,&localU);
220: return(0);
221: }
223: /* --------------------------------------------------------------------- */
224: /*
225: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
226: */
229: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *str,void *ctx)
230: {
232: PetscInt i,rstart,rend,Mx;
233: PetscReal hx,sx;
234: AppCtx *user = (AppCtx*)ctx;
235: DM da;
236: MatStencil col[3],row;
237: PetscInt nc;
238: PetscScalar vals[3];
241: TSGetDM(ts,&da);
242: MatGetOwnershipRange(*Jpre,&rstart,&rend);
243: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
244: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
245: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
246: for (i=rstart; i<rend; i++) {
247: nc = 0;
248: row.i = i;
249: if (user->boundary == 0 && (i == 0 || i == Mx-1)) {
250: col[nc].i = i; vals[nc++] = 1.0;
251: } else if (user->boundary > 0 && i == 0) { /* Left Neumann */
252: col[nc].i = i; vals[nc++] = 1.0;
253: col[nc].i = i+1; vals[nc++] = -1.0;
254: } else if (user->boundary > 0 && i == Mx-1) { /* Right Neumann */
255: col[nc].i = i-1; vals[nc++] = -1.0;
256: col[nc].i = i; vals[nc++] = 1.0;
257: } else { /* Interior */
258: col[nc].i = i-1; vals[nc++] = -1.0*sx;
259: col[nc].i = i; vals[nc++] = 2.0*sx + a;
260: col[nc].i = i+1; vals[nc++] = -1.0*sx;
261: }
262: MatSetValuesStencil(*Jpre,1,&row,nc,col,vals,INSERT_VALUES);
263: }
265: MatAssemblyBegin(*Jpre,MAT_FINAL_ASSEMBLY);
266: MatAssemblyEnd(*Jpre,MAT_FINAL_ASSEMBLY);
267: if (*J != *Jpre) {
268: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
269: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
270: }
271: if (user->viewJacobian){
272: PetscPrintf(((PetscObject)*Jpre)->comm,"Jpre:\n");
273: MatView(*Jpre,PETSC_VIEWER_STDOUT_WORLD);
274: }
275: return(0);
276: }
278: /* ------------------------------------------------------------------- */
281: PetscErrorCode FormInitialSolution(TS ts,Vec U,void* ptr)
282: {
283: AppCtx *user=(AppCtx*)ptr;
284: PetscReal c=user->c;
285: DM da;
287: PetscInt i,xs,xm,Mx;
288: PetscScalar *u;
289: PetscReal hx,x,r;
292: TSGetDM(ts,&da);
293: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
294: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
296: hx = 1.0/(PetscReal)(Mx-1);
298: /* Get pointers to vector data */
299: DMDAVecGetArray(da,U,&u);
301: /* Get local grid boundaries */
302: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
304: /* Compute function over the locally owned part of the grid */
305: for (i=xs; i<xs+xm; i++) {
306: x = i*hx;
307: r = PetscSqrtScalar((x-.5)*(x-.5));
308: if (r < .125) {
309: u[i] = PetscExpScalar(c*r*r*r);
310: } else {
311: u[i] = 0.0;
312: }
313: }
315: /* Restore vectors */
316: DMDAVecRestoreArray(da,U,&u);
317: return(0);
318: }
322: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ptr)
323: {
325: PetscReal norm,vmax,vmin;
326: MPI_Comm comm;
327: MonitorCtx *user = (MonitorCtx*)ptr;
330: VecNorm(v,NORM_1,&norm);
331: VecMax(v,PETSC_NULL,&vmax);
332: VecMin(v,PETSC_NULL,&vmin);
333: PetscObjectGetComm((PetscObject)ts,&comm);
334: PetscPrintf(comm,"timestep %D: time %G, solution norm %G, max %G, min %G\n",step,ptime,norm,vmax,vmin);
336: if (user->drawcontours){
337: VecView(v,PETSC_VIEWER_DRAW_WORLD);
338: }
339: return(0);
340: }