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=PetscSqrtReal((x-.5)*(x-.5)) < .125
  6:            u(x) = 0.0           if r >= .125

  8:    Boundary conditions:
  9:    Dirichlet BC:
 10:    At x=0, x=1, u = 0.0

 12:    Neumann BC:
 13:    At x=0, x=1: du(x,t)/dx = 0

 15:    mpiexec -n 2 ./ex17 -da_grid_x 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
 16:          ./ex17 -da_grid_x 40 -monitor_solution
 17:          ./ex17 -da_grid_x 100  -ts_type theta -ts_theta_theta 0.5 # Midpoint is not L-stable
 18:          ./ex17 -jac_type fd_coloring  -da_grid_x 500 -boundary 1
 19:          ./ex17 -da_grid_x 100  -ts_type gl -ts_adapt_type none -ts_max_steps 2
 20: */

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

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

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

 42: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
 43: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
 44: static PetscErrorCode FormInitialSolution(TS, Vec, void *);

 46: int main(int argc, char **argv)
 47: {
 48:   TS           ts; /* nonlinear solver */
 49:   Vec          u;  /* solution, residual vectors */
 50:   Mat          J;  /* Jacobian matrix */
 51:   PetscInt     nsteps;
 52:   PetscReal    vmin, vmax, norm;
 53:   DM           da;
 54:   PetscReal    ftime, dt;
 55:   AppCtx       user; /* user-defined work context */
 56:   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);
 65:   DMSetFromOptions(da);
 66:   DMSetUp(da);

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

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

 78:   PetscOptionsGetInt(NULL, NULL, "-boundary", &user.boundary, NULL);
 79:   PetscOptionsHasName(NULL, NULL, "-viewJacobian", &user.viewJacobian);

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

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

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

 96:   ftime = 1.0;
 97:   TSSetMaxTime(ts, ftime);
 98:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Set initial conditions
102:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   FormInitialSolution(ts, u, &user);
104:   TSSetSolution(ts, u);
105:   dt = .01;
106:   TSSetTimeStep(ts, dt);

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

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

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

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

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

166:   TSGetDM(ts, &da);
167:   DMGetLocalVector(da, &localU);
168:   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);

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

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

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

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

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

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

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

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

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

270: /* ------------------------------------------------------------------- */
271: PetscErrorCode FormInitialSolution(TS ts, Vec U, void *ptr)
272: {
273:   AppCtx      *user = (AppCtx *)ptr;
274:   PetscReal    c    = user->c;
275:   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: }

305: /*TEST

307:     test:
308:       requires: !single
309:       args: -da_grid_x 40 -ts_max_steps 2 -snes_monitor_short -ksp_monitor_short -ts_monitor

311:     test:
312:       suffix: 2
313:       requires: !single
314:       args: -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5

316: TEST*/