Actual source code: ex17.c
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=sqrt((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: DM da;
40: PetscReal c;
41: PetscInt boundary; /* Type of boundary condition */
42: PetscBool viewJacobian;
43: } AppCtx;
45: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
46: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);
47: static PetscErrorCode FormInitialSolution(Vec,void*);
48: static PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*);
52: int main(int argc,char **argv)
53: {
54: TS ts; /* nonlinear solver */
55: Vec u; /* solution, residual vectors */
56: Mat J; /* Jacobian matrix */
57: PetscInt steps,maxsteps = 1000; /* iterations for convergence */
59: DM da;
60: MatFDColoring matfdcoloring = PETSC_NULL;
61: PetscReal ftime,dt;
62: MonitorCtx usermonitor; /* user-defined monitor context */
63: AppCtx user; /* user-defined work context */
64: JacobianType jacType;
66: PetscInitialize(&argc,&argv,(char *)0,help);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create distributed array (DMDA) to manage parallel grid and vectors
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-11,1,1,PETSC_NULL,&da);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Extract global vectors from DMDA;
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: DMCreateGlobalVector(da,&u);
78: /* Initialize user application context */
79: user.da = da;
80: user.c = -30.0;
81: user.boundary = 0; /* 0: Dirichlet BC; 1: Neumann BC */
82: user.viewJacobian = PETSC_FALSE;
83: PetscOptionsGetInt(PETSC_NULL,"-boundary",&user.boundary,PETSC_NULL);
84: PetscOptionsHasName(PETSC_NULL,"-viewJacobian",&user.viewJacobian);
86: usermonitor.drawcontours = PETSC_FALSE;
87: PetscOptionsHasName(PETSC_NULL,"-drawcontours",&usermonitor.drawcontours);
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Create timestepping solver context
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: TSCreate(PETSC_COMM_WORLD,&ts);
93: TSSetProblemType(ts,TS_NONLINEAR);
94: TSSetType(ts,TSTHETA);
95: TSThetaSetTheta(ts,1.0); /* Make the Theta method behave like backward Euler */
96: TSSetIFunction(ts,PETSC_NULL,FormIFunction,&user);
98: DMGetMatrix(da,MATAIJ,&J);
99: jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */
100: TSSetIJacobian(ts,J,J,FormIJacobian,&user);
102: ftime = 1.0;
103: TSSetDuration(ts,maxsteps,ftime);
104: TSMonitorSet(ts,MyTSMonitor,&usermonitor,PETSC_NULL);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Set initial conditions
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: FormInitialSolution(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: DMGetColoring(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: //TSSetFromOptions(ts);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Solve nonlinear system
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: TSSolve(ts,u,&ftime);
146: TSGetTimeStepNumber(ts,&steps);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Free work space.
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: MatDestroy(&J);
152: if (matfdcoloring) {MatFDColoringDestroy(&matfdcoloring);}
153: VecDestroy(&u);
154: TSDestroy(&ts);
155: DMDestroy(&da);
157: PetscFinalize();
158: return(0);
159: }
160: /* ------------------------------------------------------------------- */
163: static PetscErrorCode FormIFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
164: {
165: AppCtx *user=(AppCtx*)ptr;
166: DM da = (DM)user->da;
168: PetscInt i,Mx,xs,xm;
169: PetscReal hx,sx;
170: PetscScalar *u,*udot,*f;
171: Vec localU;
174: DMGetLocalVector(da,&localU);
175: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
176: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
178: hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
180: /*
181: Scatter ghost points to local vector,using the 2-step process
182: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
183: By placing code between these two statements, computations can be
184: done while messages are in transition.
185: */
186: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
187: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
189: /* Get pointers to vector data */
190: DMDAVecGetArray(da,localU,&u);
191: DMDAVecGetArray(da,Udot,&udot);
192: DMDAVecGetArray(da,F,&f);
194: /* Get local grid boundaries */
195: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
197: /* Compute function over the locally owned part of the grid */
198: for (i=xs; i<xs+xm; i++) {
199: if (user->boundary == 0) { /* Dirichlet BC */
200: if (i == 0 || i == Mx-1) {
201: f[i] = u[i]; /* F = U */
202: } else {
203: f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
204: }
205: } else { /* Neumann BC */
206: if (i == 0) {
207: f[i] = u[0] - u[1];
208: } else if (i == Mx-1) {
209: f[i] = u[i] - u[i-1];
210: } else {
211: f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
212: }
213: }
214: }
216: /* Restore vectors */
217: DMDAVecRestoreArray(da,localU,&u);
218: DMDAVecRestoreArray(da,Udot,&udot);
219: DMDAVecRestoreArray(da,F,&f);
220: DMRestoreLocalVector(da,&localU);
221: return(0);
222: }
224: /* --------------------------------------------------------------------- */
225: /*
226: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
227: */
230: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *str,void *ctx)
231: {
233: PetscInt i,rstart,rend,Mx;
234: PetscReal hx,sx;
235: AppCtx *user = (AppCtx*)ctx;
236: DM da = (DM)user->da;
237: MatStencil col[3],row;
238: PetscInt nc;
239: PetscScalar vals[3];
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(Vec U,void* ptr)
282: {
283: AppCtx *user=(AppCtx*)ptr;
284: DM da=user->da;
285: PetscReal c=user->c;
287: PetscInt i,xs,xm,Mx;
288: PetscScalar *u;
289: PetscReal hx,x,r;
292: DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
293: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
295: hx = 1.0/(PetscReal)(Mx-1);
297: /* Get pointers to vector data */
298: DMDAVecGetArray(da,U,&u);
300: /* Get local grid boundaries */
301: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
303: /* Compute function over the locally owned part of the grid */
304: for (i=xs; i<xs+xm; i++) {
305: x = i*hx;
306: r = PetscSqrtScalar((x-.5)*(x-.5));
307: if (r < .125) {
308: u[i] = PetscExpScalar(c*r*r*r);
309: } else {
310: u[i] = 0.0;
311: }
312: }
314: /* Restore vectors */
315: DMDAVecRestoreArray(da,U,&u);
316: return(0);
317: }
321: PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ptr)
322: {
324: PetscReal norm,vmax,vmin;
325: MPI_Comm comm;
326: MonitorCtx *user = (MonitorCtx*)ptr;
329: VecNorm(v,NORM_1,&norm);
330: VecMax(v,PETSC_NULL,&vmax);
331: VecMin(v,PETSC_NULL,&vmin);
332: PetscObjectGetComm((PetscObject)ts,&comm);
333: PetscPrintf(comm,"timestep %D: time %G, solution norm %G, max %G, min %G\n",step,ptime,norm,vmax,vmin);
335: if (user->drawcontours){
336: VecView(v,PETSC_VIEWER_DRAW_WORLD);
337: }
338: return(0);
339: }