Actual source code: heat.c

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

      3: /*
      4:   Solves the equation

      6:     u_t = kappa  \Delta u
      7:        Periodic boundary conditions

      9: Evolve the  heat equation:
     10: ---------------
     11: ./heat -ts_monitor -snes_monitor  -pc_type lu  -draw_pause .1 -snes_converged_reason  -ts_type cn  -da_refine 5 -mymonitor

     13: Evolve the  Allen-Cahn equation:
     14: ---------------
     15: ./heat -ts_monitor -snes_monitor  -pc_type lu  -draw_pause .1 -snes_converged_reason  -ts_type cn  -da_refine 5   -allen-cahn -kappa .001 -ts_max_time 5  -mymonitor

     17: Evolve the  Allen-Cahn equation: zoom in on part of the domain
     18: ---------------
     19: ./heat -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason     -ts_type cn  -da_refine 5   -allen-cahn -kappa .001 -ts_max_time 5  -zoom .25,.45  -mymonitor

     21: The option -square_initial indicates it should use a square wave initial condition otherwise it loads the file InitialSolution.heat as the initial solution. You should run with
     22: ./heat -square_initial -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason    -ts_type cn  -da_refine 9 -ts_max_time 1.e-4 -ts_time_step .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25  -ts_max_steps 15
     23: to generate InitialSolution.heat

     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, PetscCtx), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, PetscCtx);
     35: extern PetscCtxDestroyFn MyDestroy;
     36: typedef struct {
     37:   PetscReal           kappa;
     38:   PetscBool           allencahn;
     39:   PetscDrawViewPorts *ports;
     40: } AppCtx;

     42: int main(int argc, char **argv)
     43: {
     44:   TS          ts;   /* time integrator */
     45:   Vec         x, r; /* solution, residual vectors */
     46:   PetscInt    steps, Mx;
     47:   DM          da;
     48:   PetscReal   dt;
     49:   AppCtx      ctx;
     50:   PetscBool   mymonitor;
     51:   PetscViewer viewer;
     52:   PetscBool   flg;

     54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     55:      Initialize program
     56:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     57:   PetscFunctionBeginUser;
     58:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     59:   ctx.kappa = 1.0;
     60:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL));
     61:   ctx.allencahn = PETSC_FALSE;
     62:   PetscCall(PetscOptionsHasName(NULL, NULL, "-allen-cahn", &ctx.allencahn));
     63:   PetscCall(PetscOptionsHasName(NULL, NULL, "-mymonitor", &mymonitor));

     65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     66:      Create distributed array (DMDA) to manage parallel grid and vectors
     67:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     68:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da));
     69:   PetscCall(DMSetFromOptions(da));
     70:   PetscCall(DMSetUp(da));
     71:   PetscCall(DMDASetFieldName(da, 0, "Heat equation: u"));
     72:   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
     73:   dt = 1.0 / (ctx.kappa * Mx * Mx);

     75:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76:      Extract global vectors from DMDA; then duplicate for remaining
     77:      vectors that are the same types
     78:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     79:   PetscCall(DMCreateGlobalVector(da, &x));
     80:   PetscCall(VecDuplicate(x, &r));

     82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     83:      Create timestepping solver context
     84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     85:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     86:   PetscCall(TSSetDM(ts, da));
     87:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
     88:   PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &ctx));

     90:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     91:      Customize nonlinear solver
     92:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     93:   PetscCall(TSSetType(ts, TSCN));

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:      Set initial conditions
     97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     98:   PetscCall(FormInitialSolution(da, x));
     99:   PetscCall(TSSetTimeStep(ts, dt));
    100:   PetscCall(TSSetMaxTime(ts, .02));
    101:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE));
    102:   PetscCall(TSSetSolution(ts, x));

    104:   if (mymonitor) {
    105:     ctx.ports = NULL;
    106:     PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy));
    107:   }

    109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    110:      Set runtime options
    111:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    112:   PetscCall(TSSetFromOptions(ts));

    114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    115:      Solve nonlinear system
    116:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    117:   PetscCall(TSSolve(ts, x));
    118:   PetscCall(TSGetStepNumber(ts, &steps));
    119:   PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg));
    120:   if (flg) {
    121:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_WRITE, &viewer));
    122:     PetscCall(VecView(x, viewer));
    123:     PetscCall(PetscViewerDestroy(&viewer));
    124:   }

    126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    127:      Free work space.  All PETSc objects should be destroyed when they
    128:      are no longer needed.
    129:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    130:   PetscCall(VecDestroy(&x));
    131:   PetscCall(VecDestroy(&r));
    132:   PetscCall(TSDestroy(&ts));
    133:   PetscCall(DMDestroy(&da));

    135:   PetscCall(PetscFinalize());
    136:   return 0;
    137: }
    138: /* ------------------------------------------------------------------- */
    139: /*
    140:    FormFunction - Evaluates nonlinear function, F(x).

    142:    Input Parameters:
    143: .  ts - the TS context
    144: .  X - input vector
    145: .  ptr - optional user-defined context, as set by SNESSetFunction()

    147:    Output Parameter:
    148: .  F - function vector
    149:  */
    150: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, PetscCtx Ctx)
    151: {
    152:   DM           da;
    153:   PetscInt     i, Mx, xs, xm;
    154:   PetscReal    hx, sx;
    155:   PetscScalar *x, *f;
    156:   Vec          localX;
    157:   AppCtx      *ctx = (AppCtx *)Ctx;

    159:   PetscFunctionBegin;
    160:   PetscCall(TSGetDM(ts, &da));
    161:   PetscCall(DMGetLocalVector(da, &localX));
    162:   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));

    164:   hx = 1.0 / (PetscReal)Mx;
    165:   sx = 1.0 / (hx * hx);

    167:   /*
    168:      Scatter ghost points to local vector,using the 2-step process
    169:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    170:      By placing code between these two statements, computations can be
    171:      done while messages are in transition.
    172:   */
    173:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
    174:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

    176:   /*
    177:      Get pointers to vector data
    178:   */
    179:   PetscCall(DMDAVecGetArrayRead(da, localX, &x));
    180:   PetscCall(DMDAVecGetArray(da, F, &f));

    182:   /*
    183:      Get local grid boundaries
    184:   */
    185:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

    187:   /*
    188:      Compute function over the locally owned part of the grid
    189:   */
    190:   for (i = xs; i < xs + xm; i++) {
    191:     f[i] = ctx->kappa * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    192:     if (ctx->allencahn) f[i] += (x[i] - x[i] * x[i] * x[i]);
    193:   }

    195:   /*
    196:      Restore vectors
    197:   */
    198:   PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
    199:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    200:   PetscCall(DMRestoreLocalVector(da, &localX));
    201:   PetscFunctionReturn(PETSC_SUCCESS);
    202: }

    204: /* ------------------------------------------------------------------- */
    205: PetscErrorCode FormInitialSolution(DM da, Vec U)
    206: {
    207:   PetscInt           i, xs, xm, Mx, scale = 1, N;
    208:   PetscScalar       *u;
    209:   const PetscScalar *f;
    210:   PetscReal          hx, x, r;
    211:   Vec                finesolution;
    212:   PetscViewer        viewer;
    213:   PetscBool          flg;

    215:   PetscFunctionBegin;
    216:   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));

    218:   hx = 1.0 / (PetscReal)Mx;

    220:   /*
    221:      Get pointers to vector data
    222:   */
    223:   PetscCall(DMDAVecGetArray(da, U, &u));

    225:   /*
    226:      Get local grid boundaries
    227:   */
    228:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

    230:   /*  InitialSolution is obtained with
    231:       ./heat -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason    -ts_type cn  -da_refine 9 -ts_max_time 1.e-4 -ts_time_step .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25  -ts_max_steps 15
    232:   */
    233:   PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg));
    234:   if (!flg) {
    235:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer));
    236:     PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution));
    237:     PetscCall(VecLoad(finesolution, viewer));
    238:     PetscCall(PetscViewerDestroy(&viewer));
    239:     PetscCall(VecGetSize(finesolution, &N));
    240:     scale = N / Mx;
    241:     PetscCall(VecGetArrayRead(finesolution, &f));
    242:   }

    244:   /*
    245:      Compute function over the locally owned part of the grid
    246:   */
    247:   for (i = xs; i < xs + xm; i++) {
    248:     x = i * hx;
    249:     r = PetscSqrtReal((x - .5) * (x - .5));
    250:     if (r < .125) u[i] = 1.0;
    251:     else u[i] = -.5;

    253:     /* With the initial condition above the method is first order in space */
    254:     /* this is a smooth initial condition so the method becomes second order in space */
    255:     /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
    256:     /*  u[i] = f[scale*i];*/
    257:     if (!flg) u[i] = f[scale * i];
    258:   }
    259:   if (!flg) {
    260:     PetscCall(VecRestoreArrayRead(finesolution, &f));
    261:     PetscCall(VecDestroy(&finesolution));
    262:   }

    264:   /*
    265:      Restore vectors
    266:   */
    267:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    268:   PetscFunctionReturn(PETSC_SUCCESS);
    269: }

    271: /*
    272:     This routine is not parallel
    273: */
    274: PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, PetscCtx Ctx)
    275: {
    276:   AppCtx             *ctx = (AppCtx *)Ctx;
    277:   PetscDrawLG         lg;
    278:   PetscScalar        *u;
    279:   PetscInt            Mx, i, xs, xm, cnt;
    280:   PetscReal           x, y, hx, pause, sx, len, max, xx[2], yy[2];
    281:   PetscDraw           draw;
    282:   Vec                 localU;
    283:   DM                  da;
    284:   int                 colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE};
    285:   const char *const   legend[] = {"-kappa (\\grad u,\\grad u)", "(1 - u^2)^2"};
    286:   PetscDrawAxis       axis;
    287:   PetscDrawViewPorts *ports;
    288:   PetscReal           vbounds[] = {-1.1, 1.1};

    290:   PetscFunctionBegin;
    291:   PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds));
    292:   PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1200, 800));
    293:   PetscCall(TSGetDM(ts, &da));
    294:   PetscCall(DMGetLocalVector(da, &localU));
    295:   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));
    296:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
    297:   hx = 1.0 / (PetscReal)Mx;
    298:   sx = 1.0 / (hx * hx);
    299:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
    300:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
    301:   PetscCall(DMDAVecGetArrayRead(da, localU, &u));

    303:   PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg));
    304:   PetscCall(PetscDrawLGGetDraw(lg, &draw));
    305:   PetscCall(PetscDrawCheckResizedWindow(draw));
    306:   if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports));
    307:   ports = ctx->ports;
    308:   PetscCall(PetscDrawLGGetAxis(lg, &axis));
    309:   PetscCall(PetscDrawLGReset(lg));

    311:   xx[0] = 0.0;
    312:   xx[1] = 1.0;
    313:   cnt   = 2;
    314:   PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL));
    315:   xs = xx[0] / hx;
    316:   xm = (xx[1] - xx[0]) / hx;

    318:   /*
    319:       Plot the  energies
    320:   */
    321:   PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->allencahn ? 1 : 0)));
    322:   PetscCall(PetscDrawLGSetColors(lg, colors + 1));
    323:   PetscCall(PetscDrawViewPortsSet(ports, 2));
    324:   x = hx * xs;
    325:   for (i = xs; i < xs + xm; i++) {
    326:     xx[0] = xx[1] = x;
    327:     yy[0]         = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);
    328:     if (ctx->allencahn) yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
    329:     PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    330:     x += hx;
    331:   }
    332:   PetscCall(PetscDrawGetPause(draw, &pause));
    333:   PetscCall(PetscDrawSetPause(draw, 0.0));
    334:   PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", ""));
    335:   PetscCall(PetscDrawLGSetLegend(lg, legend));
    336:   PetscCall(PetscDrawLGDraw(lg));

    338:   /*
    339:       Plot the  forces
    340:   */
    341:   PetscCall(PetscDrawViewPortsSet(ports, 1));
    342:   PetscCall(PetscDrawLGReset(lg));
    343:   x   = xs * hx;
    344:   max = 0.;
    345:   for (i = xs; i < xs + xm; i++) {
    346:     xx[0] = xx[1] = x;
    347:     yy[0]         = PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    348:     max           = PetscMax(max, PetscAbs(yy[0]));
    349:     if (ctx->allencahn) {
    350:       yy[1] = PetscRealPart(u[i] - u[i] * u[i] * u[i]);
    351:       max   = PetscMax(max, PetscAbs(yy[1]));
    352:     }
    353:     PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    354:     x += hx;
    355:   }
    356:   PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", ""));
    357:   PetscCall(PetscDrawLGSetLegend(lg, NULL));
    358:   PetscCall(PetscDrawLGDraw(lg));

    360:   /*
    361:         Plot the solution
    362:   */
    363:   PetscCall(PetscDrawLGSetDimension(lg, 1));
    364:   PetscCall(PetscDrawViewPortsSet(ports, 0));
    365:   PetscCall(PetscDrawLGReset(lg));
    366:   x = hx * xs;
    367:   PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1));
    368:   PetscCall(PetscDrawLGSetColors(lg, colors));
    369:   for (i = xs; i < xs + xm; i++) {
    370:     xx[0] = x;
    371:     yy[0] = PetscRealPart(u[i]);
    372:     PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    373:     x += hx;
    374:   }
    375:   PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", ""));
    376:   PetscCall(PetscDrawLGDraw(lg));

    378:   /*
    379:       Print the  forces as arrows on the solution
    380:   */
    381:   x   = hx * xs;
    382:   cnt = xm / 60;
    383:   cnt = (!cnt) ? 1 : cnt;

    385:   for (i = xs; i < xs + xm; i += cnt) {
    386:     y   = PetscRealPart(u[i]);
    387:     len = .5 * PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
    388:     PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED));
    389:     if (ctx->allencahn) {
    390:       len = .5 * PetscRealPart(u[i] - u[i] * u[i] * u[i]) / max;
    391:       PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE));
    392:     }
    393:     x += cnt * hx;
    394:   }
    395:   PetscCall(DMDAVecRestoreArrayRead(da, localU, &x));
    396:   PetscCall(DMRestoreLocalVector(da, &localU));
    397:   PetscCall(PetscDrawStringSetSize(draw, .2, .2));
    398:   PetscCall(PetscDrawFlush(draw));
    399:   PetscCall(PetscDrawSetPause(draw, pause));
    400:   PetscCall(PetscDrawPause(draw));
    401:   PetscFunctionReturn(PETSC_SUCCESS);
    402: }

    404: PetscErrorCode MyDestroy(PetscCtxRt Ctx)
    405: {
    406:   AppCtx *ctx = *(AppCtx **)Ctx;

    408:   PetscFunctionBegin;
    409:   PetscCall(PetscDrawViewPortsDestroy(ctx->ports));
    410:   PetscFunctionReturn(PETSC_SUCCESS);
    411: }

    413: /*TEST

    415:    test:
    416:      args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial

    418:    test:
    419:      suffix: 2
    420:      args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001
    421:      requires: x

    423: TEST*/