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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    217:   hx = 1.0 / (PetscReal)Mx;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    403: PetscErrorCode MyDestroy(void **ptr)
    404: {
    405:   UserCtx *ctx = *(UserCtx **)ptr;

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

    412: /*TEST

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

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

    422: TEST*/