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: }