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:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 61:   PetscFunctionBeginUser;
 62:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 63:   ctx.kappa = 1.0;
 64:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL));
 65:   ctx.cahnhillard = PETSC_FALSE;

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

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:      Create distributed array (DMDA) to manage parallel grid and vectors
 82:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 83:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 2, 2, NULL, &da));
 84:   PetscCall(DMSetFromOptions(da));
 85:   PetscCall(DMSetUp(da));
 86:   PetscCall(DMDASetFieldName(da, 0, "Biharmonic heat equation: w = -kappa*u_xx"));
 87:   PetscCall(DMDASetFieldName(da, 1, "Biharmonic heat equation: u"));
 88:   PetscCall(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:   PetscCall(DMCreateGlobalVector(da, &x));
 96:   PetscCall(VecDuplicate(x, &r));

 98:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:      Create timestepping solver context
100:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
102:   PetscCall(TSSetDM(ts, da));
103:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
104:   PetscCall(TSSetIFunction(ts, NULL, FormFunction, &ctx));
105:   PetscCall(TSSetMaxTime(ts, .02));
106:   PetscCall(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:   PetscCall(TSGetSNES(ts, &snes));
121:   PetscCall(DMCreateColoring(da, IS_COLORING_GLOBAL, &iscoloring));
122:   PetscCall(DMSetMatType(da, MATAIJ));
123:   PetscCall(DMCreateMatrix(da, &J));
124:   PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring));
125:   PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, ts));
126:   PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
127:   PetscCall(MatFDColoringSetUp(J, iscoloring, matfdcoloring));
128:   PetscCall(ISColoringDestroy(&iscoloring));
129:   PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring));

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

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

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

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

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

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

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

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

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

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

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

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

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

267:   PetscFunctionBegin;
268:   PetscCall(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));

270:   hx = 1.0 / (PetscReal)Mx;
271:   sx = 1.0 / (hx * hx);

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

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

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

298:   /*
299:      Restore vectors
300:   */
301:   PetscCall(DMDAVecRestoreArray(da, Xg, &x));

303:   /* Grab only the global part of the vector */
304:   PetscCall(VecSet(X, 0));
305:   PetscCall(DMLocalToGlobalBegin(da, Xg, ADD_VALUES, X));
306:   PetscCall(DMLocalToGlobalEnd(da, Xg, ADD_VALUES, X));
307:   PetscCall(VecDestroy(&Xg));
308:   PetscFunctionReturn(PETSC_SUCCESS);
309: }

311: /*TEST

313:    build:
314:      requires: !complex !single

316:    test:
317:      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
318:      requires: x

320: TEST*/