Actual source code: biharmonic2.c


  2: static char help[] = "Solves biharmonic equation in 1d.\n";

  4: /*
  5:   Solves the equation biharmonic equation in split form

  7:     w = -kappa \Delta u
  8:     u_t =  \Delta w
  9:     -1  <= u <= 1
 10:     Periodic boundary conditions

 12: Evolve the biharmonic heat equation with bounds:  (same as biharmonic)
 13: ---------------
 14: ./biharmonic2 -ts_monitor -snes_monitor -ts_monitor_draw_solution  -pc_type lu  -draw_pause .1 -snes_converged_reason  -ts_type beuler  -da_refine 5 -draw_fields 1 -ts_dt 9.53674e-9

 16:     w = -kappa \Delta u  + u^3  - u
 17:     u_t =  \Delta w
 18:     -1  <= u <= 1
 19:     Periodic boundary conditions

 21: Evolve the Cahn-Hillard equations: (this fails after a few timesteps 12/17/2017)
 22: ---------------
 23: ./biharmonic2 -ts_monitor -snes_monitor -ts_monitor_draw_solution  -pc_type lu  -draw_pause .1 -snes_converged_reason   -ts_type beuler    -da_refine 6  -draw_fields 1  -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard

 25: */
 26: #include <petscdm.h>
 27: #include <petscdmda.h>
 28: #include <petscts.h>
 29: #include <petscdraw.h>

 31: /*
 32:    User-defined routines
 33: */
 34: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, Vec, void *), FormInitialSolution(DM, Vec, PetscReal);
 35: typedef struct {
 36:   PetscBool cahnhillard;
 37:   PetscReal kappa;
 38:   PetscInt  energy;
 39:   PetscReal tol;
 40:   PetscReal theta;
 41:   PetscReal theta_c;
 42: } UserCtx;

 44: int main(int argc, char **argv)
 45: {
 46:   TS            ts;   /* nonlinear solver */
 47:   Vec           x, r; /* solution, residual vectors */
 48:   Mat           J;    /* Jacobian matrix */
 49:   PetscInt      steps, Mx;
 50:   DM            da;
 51:   MatFDColoring matfdcoloring;
 52:   ISColoring    iscoloring;
 53:   PetscReal     dt;
 54:   PetscReal     vbounds[] = {-100000, 100000, -1.1, 1.1};
 55:   SNES          snes;
 56:   UserCtx       ctx;

 58:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 59:      Initialize program
 60:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 62:   PetscInitialize(&argc, &argv, (char *)0, help);
 63:   ctx.kappa = 1.0;
 64:   PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL);
 65:   ctx.cahnhillard = PETSC_FALSE;

 67:   PetscOptionsGetBool(NULL, NULL, "-cahn-hillard", &ctx.cahnhillard, NULL);
 68:   PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 2, vbounds);
 69:   PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 600, 600);
 70:   ctx.energy = 1;
 71:   /*PetscOptionsGetInt(NULL,NULL,"-energy",&ctx.energy,NULL);*/
 72:   PetscOptionsGetInt(NULL, NULL, "-energy", &ctx.energy, NULL);
 73:   ctx.tol = 1.0e-8;
 74:   PetscOptionsGetReal(NULL, NULL, "-tol", &ctx.tol, NULL);
 75:   ctx.theta   = .001;
 76:   ctx.theta_c = 1.0;
 77:   PetscOptionsGetReal(NULL, NULL, "-theta", &ctx.theta, NULL);
 78:   PetscOptionsGetReal(NULL, NULL, "-theta_c", &ctx.theta_c, NULL);

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:      Create distributed array (DMDA) to manage parallel grid and vectors
 82:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 83:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 2, 2, NULL, &da);
 84:   DMSetFromOptions(da);
 85:   DMSetUp(da);
 86:   DMDASetFieldName(da, 0, "Biharmonic heat equation: w = -kappa*u_xx");
 87:   DMDASetFieldName(da, 1, "Biharmonic heat equation: u");
 88:   DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
 89:   dt = 1.0 / (10. * ctx.kappa * Mx * Mx * Mx * Mx);

 91:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:      Extract global vectors from DMDA; then duplicate for remaining
 93:      vectors that are the same types
 94:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 95:   DMCreateGlobalVector(da, &x);
 96:   VecDuplicate(x, &r);

 98:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:      Create timestepping solver context
100:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101:   TSCreate(PETSC_COMM_WORLD, &ts);
102:   TSSetDM(ts, da);
103:   TSSetProblemType(ts, TS_NONLINEAR);
104:   TSSetIFunction(ts, NULL, FormFunction, &ctx);
105:   TSSetMaxTime(ts, .02);
106:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE);

108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109:      Create matrix data structure; set Jacobian evaluation routine

111: <     Set Jacobian matrix data structure and default Jacobian evaluation
112:      routine. User can override with:
113:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
114:                 (unless user explicitly sets preconditioner)
115:      -snes_mf_operator : form preconditioning matrix as set by the user,
116:                          but use matrix-free approx for Jacobian-vector
117:                          products within Newton-Krylov method

119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   TSGetSNES(ts, &snes);
121:   DMCreateColoring(da, IS_COLORING_GLOBAL, &iscoloring);
122:   DMSetMatType(da, MATAIJ);
123:   DMCreateMatrix(da, &J);
124:   MatFDColoringCreate(J, iscoloring, &matfdcoloring);
125:   MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, ts);
126:   MatFDColoringSetFromOptions(matfdcoloring);
127:   MatFDColoringSetUp(J, iscoloring, matfdcoloring);
128:   ISColoringDestroy(&iscoloring);
129:   SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring);

131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132:      Customize nonlinear solver
133:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134:   TSSetType(ts, TSBEULER);

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:      Set initial conditions
138:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139:   FormInitialSolution(da, x, ctx.kappa);
140:   TSSetTimeStep(ts, dt);
141:   TSSetSolution(ts, x);

143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:      Set runtime options
145:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146:   TSSetFromOptions(ts);

148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149:      Solve nonlinear system
150:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151:   TSSolve(ts, x);
152:   TSGetStepNumber(ts, &steps);

154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:      Free work space.  All PETSc objects should be destroyed when they
156:      are no longer needed.
157:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158:   MatDestroy(&J);
159:   MatFDColoringDestroy(&matfdcoloring);
160:   VecDestroy(&x);
161:   VecDestroy(&r);
162:   TSDestroy(&ts);
163:   DMDestroy(&da);

165:   PetscFinalize();
166:   return 0;
167: }

169: typedef struct {
170:   PetscScalar w, u;
171: } Field;
172: /* ------------------------------------------------------------------- */
173: /*
174:    FormFunction - Evaluates nonlinear function, F(x).

176:    Input Parameters:
177: .  ts - the TS context
178: .  X - input vector
179: .  ptr - optional user-defined context, as set by SNESSetFunction()

181:    Output Parameter:
182: .  F - function vector
183:  */
184: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec Xdot, Vec F, void *ptr)
185: {
186:   DM        da;
187:   PetscInt  i, Mx, xs, xm;
188:   PetscReal hx, sx;
189:   Field    *x, *xdot, *f;
190:   Vec       localX, localXdot;
191:   UserCtx  *ctx = (UserCtx *)ptr;

193:   TSGetDM(ts, &da);
194:   DMGetLocalVector(da, &localX);
195:   DMGetLocalVector(da, &localXdot);
196:   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);

198:   hx = 1.0 / (PetscReal)Mx;
199:   sx = 1.0 / (hx * hx);

201:   /*
202:      Scatter ghost points to local vector,using the 2-step process
203:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
204:      By placing code between these two statements, computations can be
205:      done while messages are in transition.
206:   */
207:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
208:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
209:   DMGlobalToLocalBegin(da, Xdot, INSERT_VALUES, localXdot);
210:   DMGlobalToLocalEnd(da, Xdot, INSERT_VALUES, localXdot);

212:   /*
213:      Get pointers to vector data
214:   */
215:   DMDAVecGetArrayRead(da, localX, &x);
216:   DMDAVecGetArrayRead(da, localXdot, &xdot);
217:   DMDAVecGetArray(da, F, &f);

219:   /*
220:      Get local grid boundaries
221:   */
222:   DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);

224:   /*
225:      Compute function over the locally owned part of the grid
226:   */
227:   for (i = xs; i < xs + xm; i++) {
228:     f[i].w = x[i].w + ctx->kappa * (x[i - 1].u + x[i + 1].u - 2.0 * x[i].u) * sx;
229:     if (ctx->cahnhillard) {
230:       switch (ctx->energy) {
231:       case 1: /* double well */
232:         f[i].w += -x[i].u * x[i].u * x[i].u + x[i].u;
233:         break;
234:       case 2: /* double obstacle */
235:         f[i].w += x[i].u;
236:         break;
237:       case 3: /* logarithmic */
238:         if (PetscRealPart(x[i].u) < -1.0 + 2.0 * ctx->tol) f[i].w += .5 * ctx->theta * (-PetscLogReal(ctx->tol) + PetscLogScalar((1.0 - x[i].u) / 2.0)) + ctx->theta_c * x[i].u;
239:         else if (PetscRealPart(x[i].u) > 1.0 - 2.0 * ctx->tol) f[i].w += .5 * ctx->theta * (-PetscLogScalar((1.0 + x[i].u) / 2.0) + PetscLogReal(ctx->tol)) + ctx->theta_c * x[i].u;
240:         else f[i].w += .5 * ctx->theta * (-PetscLogScalar((1.0 + x[i].u) / 2.0) + PetscLogScalar((1.0 - x[i].u) / 2.0)) + ctx->theta_c * x[i].u;
241:         break;
242:       }
243:     }
244:     f[i].u = xdot[i].u - (x[i - 1].w + x[i + 1].w - 2.0 * x[i].w) * sx;
245:   }

247:   /*
248:      Restore vectors
249:   */
250:   DMDAVecRestoreArrayRead(da, localXdot, &xdot);
251:   DMDAVecRestoreArrayRead(da, localX, &x);
252:   DMDAVecRestoreArray(da, F, &f);
253:   DMRestoreLocalVector(da, &localX);
254:   DMRestoreLocalVector(da, &localXdot);
255:   return 0;
256: }

258: /* ------------------------------------------------------------------- */
259: PetscErrorCode FormInitialSolution(DM da, Vec X, PetscReal kappa)
260: {
261:   PetscInt  i, xs, xm, Mx, xgs, xgm;
262:   Field    *x;
263:   PetscReal hx, xx, r, sx;
264:   Vec       Xg;

266:   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);

268:   hx = 1.0 / (PetscReal)Mx;
269:   sx = 1.0 / (hx * hx);

271:   /*
272:      Get pointers to vector data
273:   */
274:   DMCreateLocalVector(da, &Xg);
275:   DMDAVecGetArray(da, Xg, &x);

277:   /*
278:      Get local grid boundaries
279:   */
280:   DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
281:   DMDAGetGhostCorners(da, &xgs, NULL, NULL, &xgm, NULL, NULL);

283:   /*
284:      Compute u function over the locally owned part of the grid including ghost points
285:   */
286:   for (i = xgs; i < xgs + xgm; i++) {
287:     xx = i * hx;
288:     r  = PetscSqrtReal((xx - .5) * (xx - .5));
289:     if (r < .125) x[i].u = 1.0;
290:     else x[i].u = -.50;
291:     /* fill in x[i].w so that valgrind doesn't detect use of uninitialized memory */
292:     x[i].w = 0;
293:   }
294:   for (i = xs; i < xs + xm; i++) x[i].w = -kappa * (x[i - 1].u + x[i + 1].u - 2.0 * x[i].u) * sx;

296:   /*
297:      Restore vectors
298:   */
299:   DMDAVecRestoreArray(da, Xg, &x);

301:   /* Grab only the global part of the vector */
302:   VecSet(X, 0);
303:   DMLocalToGlobalBegin(da, Xg, ADD_VALUES, X);
304:   DMLocalToGlobalEnd(da, Xg, ADD_VALUES, X);
305:   VecDestroy(&Xg);
306:   return 0;
307: }

309: /*TEST

311:    build:
312:      requires: !complex !single

314:    test:
315:      args: -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason  -ts_type beuler  -da_refine 5 -ts_dt 9.53674e-9 -ts_max_steps 50
316:      requires: x

318: TEST*/