Actual source code: ex17.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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:    mpiexec -n 2 ./ex17 -da_grid_x 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
 17:          ./ex17 -da_grid_x 40 -monitor_solution
 18:          ./ex17 -da_grid_x 100  -ts_type theta -ts_theta_theta 0.5 # Midpoint is not L-stable
 19:          ./ex17 -jac_type fd_coloring  -da_grid_x 500 -boundary 1
 20:          ./ex17 -da_grid_x 100  -ts_type gl -ts_adapt_type none -ts_max_steps 2
 21: */

 23: #include <petscdm.h>
 24: #include <petscdmda.h>
 25: #include <petscts.h>

 27: typedef enum {JACOBIAN_ANALYTIC,JACOBIAN_FD_COLORING,JACOBIAN_FD_FULL} JacobianType;
 28: static const char *const JacobianTypes[] = {"analytic","fd_coloring","fd_full","JacobianType","fd_",0};

 30: /*
 31:    User-defined data structures and routines
 32: */
 33: typedef struct {
 34:   PetscReal c;
 35:   PetscInt  boundary;            /* Type of boundary condition */
 36:   PetscBool viewJacobian;
 37: } AppCtx;

 39: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
 40: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
 41: static PetscErrorCode FormInitialSolution(TS,Vec,void*);

 45: int main(int argc,char **argv)
 46: {
 47:   TS             ts;                   /* nonlinear solver */
 48:   Vec            u;                    /* solution, residual vectors */
 49:   Mat            J;                    /* Jacobian matrix */
 50:   PetscInt       maxsteps = 1000;     /* iterations for convergence */
 51:   PetscInt       nsteps;
 52:   PetscReal      vmin,vmax,norm;
 54:   DM             da;
 55:   PetscReal      ftime,dt;
 56:   AppCtx         user;              /* user-defined work context */
 57:   JacobianType   jacType;

 59:   PetscInitialize(&argc,&argv,(char*)0,help);

 61:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62:      Create distributed array (DMDA) to manage parallel grid and vectors
 63:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 64:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-11,1,1,NULL,&da);

 66:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:      Extract global vectors from DMDA;
 68:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 69:   DMCreateGlobalVector(da,&u);

 71:   /* Initialize user application context */
 72:   user.c            = -30.0;
 73:   user.boundary     = 0;  /* 0: Dirichlet BC; 1: Neumann BC */
 74:   user.viewJacobian = PETSC_FALSE;

 76:   PetscOptionsGetInt(NULL,NULL,"-boundary",&user.boundary,NULL);
 77:   PetscOptionsHasName(NULL,NULL,"-viewJacobian",&user.viewJacobian);

 79:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80:      Create timestepping solver context
 81:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 82:   TSCreate(PETSC_COMM_WORLD,&ts);
 83:   TSSetProblemType(ts,TS_NONLINEAR);
 84:   TSSetType(ts,TSTHETA);
 85:   TSThetaSetTheta(ts,1.0); /* Make the Theta method behave like backward Euler */
 86:   TSSetIFunction(ts,NULL,FormIFunction,&user);

 88:   DMSetMatType(da,MATAIJ);
 89:   DMCreateMatrix(da,&J);
 90:   jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */

 92:   TSSetDM(ts,da); /* Use TSGetDM() to access. Setting here allows easy use of geometric multigrid. */

 94:   ftime = 1.0;
 95:   TSSetDuration(ts,maxsteps,ftime);
 96:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

 98:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:      Set initial conditions
100:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101:   FormInitialSolution(ts,u,&user);
102:   TSSetSolution(ts,u);
103:   dt   = .01;
104:   TSSetInitialTimeStep(ts,0.0,dt);


107:   /* Use slow fd Jacobian or fast fd Jacobian with colorings.
108:      Note: this requirs snes which is not created until TSSetUp()/TSSetFromOptions() is called */
109:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Options for Jacobian evaluation",NULL);
110:   PetscOptionsEnum("-jac_type","Type of Jacobian","",JacobianTypes,(PetscEnum)jacType,(PetscEnum*)&jacType,0);
111:   PetscOptionsEnd();
112:   if (jacType == JACOBIAN_ANALYTIC) {
113:     TSSetIJacobian(ts,J,J,FormIJacobian,&user);
114:   } else if (jacType == JACOBIAN_FD_COLORING) {
115:     SNES snes;
116:     TSGetSNES(ts,&snes);
117:     SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,0);
118:   } else if (jacType == JACOBIAN_FD_FULL) {
119:     SNES snes;
120:     TSGetSNES(ts,&snes);
121:     SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,&user);
122:   }

124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125:      Set runtime options
126:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127:   TSSetFromOptions(ts);

129:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130:      Integrate ODE system
131:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132:   TSSolve(ts,u);

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:    Compute diagnostics of the solution
136:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137:   VecNorm(u,NORM_1,&norm);
138:   VecMax(u,NULL,&vmax);
139:   VecMin(u,NULL,&vmin);
140:   TSGetTimeStepNumber(ts,&nsteps);
141:   TSGetTime(ts,&ftime);
142:   PetscPrintf(PETSC_COMM_WORLD,"timestep %D: time %g, solution norm %g, max %g, min %g\n",nsteps,(double)ftime,(double)norm,(double)vmax,(double)vmin);

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:      Free work space.
146:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147:   MatDestroy(&J);
148:   VecDestroy(&u);
149:   TSDestroy(&ts);
150:   DMDestroy(&da);
151:   PetscFinalize();
152:   return(0);
153: }
154: /* ------------------------------------------------------------------- */
157: static PetscErrorCode FormIFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
158: {
159:   AppCtx         *user=(AppCtx*)ptr;
160:   DM             da;
162:   PetscInt       i,Mx,xs,xm;
163:   PetscReal      hx,sx;
164:   PetscScalar    *u,*udot,*f;
165:   Vec            localU;

168:   TSGetDM(ts,&da);
169:   DMGetLocalVector(da,&localU);
170:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
171:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

173:   hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);

175:   /*
176:      Scatter ghost points to local vector,using the 2-step process
177:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
178:      By placing code between these two statements, computations can be
179:      done while messages are in transition.
180:   */
181:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
182:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

184:   /* Get pointers to vector data */
185:   DMDAVecGetArrayRead(da,localU,&u);
186:   DMDAVecGetArrayRead(da,Udot,&udot);
187:   DMDAVecGetArray(da,F,&f);

189:   /* Get local grid boundaries */
190:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

192:   /* Compute function over the locally owned part of the grid */
193:   for (i=xs; i<xs+xm; i++) {
194:     if (user->boundary == 0) { /* Dirichlet BC */
195:       if (i == 0 || i == Mx-1) f[i] = u[i]; /* F = U */
196:       else                     f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
197:     } else { /* Neumann BC */
198:       if (i == 0)         f[i] = u[0] - u[1];
199:       else if (i == Mx-1) f[i] = u[i] - u[i-1];
200:       else                f[i] = udot[i] + (2.*u[i] - u[i-1] - u[i+1])*sx;
201:     }
202:   }

204:   /* Restore vectors */
205:   DMDAVecRestoreArrayRead(da,localU,&u);
206:   DMDAVecRestoreArrayRead(da,Udot,&udot);
207:   DMDAVecRestoreArray(da,F,&f);
208:   DMRestoreLocalVector(da,&localU);
209:   return(0);
210: }

212: /* --------------------------------------------------------------------- */
213: /*
214:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
215: */
218: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,void *ctx)
219: {
221:   PetscInt       i,rstart,rend,Mx;
222:   PetscReal      hx,sx;
223:   AppCtx         *user = (AppCtx*)ctx;
224:   DM             da;
225:   MatStencil     col[3],row;
226:   PetscInt       nc;
227:   PetscScalar    vals[3];

230:   TSGetDM(ts,&da);
231:   MatGetOwnershipRange(Jpre,&rstart,&rend);
232:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
233:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
234:   hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
235:   for (i=rstart; i<rend; i++) {
236:     nc    = 0;
237:     row.i = i;
238:     if (user->boundary == 0 && (i == 0 || i == Mx-1)) {
239:       col[nc].i = i; vals[nc++] = 1.0;
240:     } else if (user->boundary > 0 && i == 0) { /* Left Neumann */
241:       col[nc].i = i;   vals[nc++] = 1.0;
242:       col[nc].i = i+1; vals[nc++] = -1.0;
243:     } else if (user->boundary > 0 && i == Mx-1) { /* Right Neumann */
244:       col[nc].i = i-1; vals[nc++] = -1.0;
245:       col[nc].i = i;   vals[nc++] = 1.0;
246:     } else {                    /* Interior */
247:       col[nc].i = i-1; vals[nc++] = -1.0*sx;
248:       col[nc].i = i;   vals[nc++] = 2.0*sx + a;
249:       col[nc].i = i+1; vals[nc++] = -1.0*sx;
250:     }
251:     MatSetValuesStencil(Jpre,1,&row,nc,col,vals,INSERT_VALUES);
252:   }

254:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
255:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
256:   if (J != Jpre) {
257:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
258:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
259:   }
260:   if (user->viewJacobian) {
261:     PetscPrintf(PETSC_COMM_WORLD,"Jpre:\n");
262:     MatView(Jpre,PETSC_VIEWER_STDOUT_WORLD);
263:   }
264:   return(0);
265: }

267: /* ------------------------------------------------------------------- */
270: PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ptr)
271: {
272:   AppCtx         *user=(AppCtx*)ptr;
273:   PetscReal      c    =user->c;
274:   DM             da;
276:   PetscInt       i,xs,xm,Mx;
277:   PetscScalar    *u;
278:   PetscReal      hx,x,r;

281:   TSGetDM(ts,&da);
282:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

284:   hx = 1.0/(PetscReal)(Mx-1);

286:   /* Get pointers to vector data */
287:   DMDAVecGetArray(da,U,&u);

289:   /* Get local grid boundaries */
290:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

292:   /* Compute function over the locally owned part of the grid */
293:   for (i=xs; i<xs+xm; i++) {
294:     x = i*hx;
295:     r = PetscSqrtReal((x-.5)*(x-.5));
296:     if (r < .125) u[i] = PetscExpReal(c*r*r*r);
297:     else          u[i] = 0.0;
298:   }

300:   /* Restore vectors */
301:   DMDAVecRestoreArray(da,U,&u);
302:   return(0);
303: }