Actual source code: allen_cahn.c

      1: static char help[] = "Solves the time dependent Allen-Cahn equation with IMEX methods";

      3: /*
      4:  * allen_cahn.c
      5:  *
      6:  *  Created on: Jun 8, 2012
      7:  *      Author: Hong Zhang
      8:  */

     10: #include <petscts.h>

     12: /*
     13:  * application context
     14:  */
     15: typedef struct {
     16:   PetscReal param;         /* parameter */
     17:   PetscReal xleft, xright; /* range in x-direction */
     18:   PetscInt  mx;            /* Discretization in x-direction */
     19: } AppCtx;

     21: static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
     22: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
     23: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *ctx);
     24: static PetscErrorCode FormInitialSolution(TS, Vec, void *);

     26: int main(int argc, char **argv)
     27: {
     28:   TS        ts;
     29:   Vec       x; /* solution vector */
     30:   Mat       A; /* Jacobian */
     31:   PetscInt  steps, mx;
     32:   PetscReal ftime;
     33:   AppCtx    user; /* user-defined work context */

     35:   PetscFunctionBeginUser;
     36:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

     38:   /* Initialize user application context */
     39:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Allen-Cahn equation", "");
     40:   user.param  = 9e-4;
     41:   user.xleft  = -1.;
     42:   user.xright = 2.;
     43:   user.mx     = 400;
     44:   PetscCall(PetscOptionsReal("-eps", "parameter", "", user.param, &user.param, NULL));
     45:   PetscOptionsEnd();

     47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     48:    Set runtime options
     49:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     50:   /*
     51:    * PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL));
     52:    */
     53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     54:    Create necessary matrix and vectors, solve same ODE on every process
     55:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     56:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     57:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, user.mx, user.mx));
     58:   PetscCall(MatSetFromOptions(A));
     59:   PetscCall(MatSetUp(A));
     60:   PetscCall(MatCreateVecs(A, &x, NULL));

     62:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     63:    Create time stepping solver context
     64:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     65:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     66:   PetscCall(TSSetType(ts, TSEIMEX));
     67:   PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user));
     68:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
     69:   PetscCall(TSSetIJacobian(ts, A, A, FormIJacobian, &user));
     70:   ftime = 22;
     71:   PetscCall(TSSetMaxTime(ts, ftime));
     72:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

     74:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     75:    Set initial conditions
     76:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     77:   PetscCall(FormInitialSolution(ts, x, &user));
     78:   PetscCall(TSSetSolution(ts, x));
     79:   PetscCall(VecGetSize(x, &mx));

     81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     82:    Set runtime options
     83:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     84:   PetscCall(TSSetFromOptions(ts));

     86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     87:    Solve nonlinear system
     88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     89:   PetscCall(TSSolve(ts, x));
     90:   PetscCall(TSGetTime(ts, &ftime));
     91:   PetscCall(TSGetStepNumber(ts, &steps));
     92:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "eps %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.param, steps, (double)ftime));
     93:   /*   PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));*/

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:    Free work space.
     97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     98:   PetscCall(MatDestroy(&A));
     99:   PetscCall(VecDestroy(&x));
    100:   PetscCall(TSDestroy(&ts));
    101:   PetscCall(PetscFinalize());
    102:   return 0;
    103: }

    105: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
    106: {
    107:   AppCtx            *user = (AppCtx *)ptr;
    108:   PetscScalar       *f;
    109:   const PetscScalar *x;
    110:   PetscInt           i, mx;
    111:   PetscReal          hx, eps;

    113:   PetscFunctionBegin;
    114:   mx  = user->mx;
    115:   eps = user->param;
    116:   hx  = (user->xright - user->xleft) / (mx - 1);
    117:   PetscCall(VecGetArrayRead(X, &x));
    118:   PetscCall(VecGetArray(F, &f));
    119:   f[0] = 2. * eps * (x[1] - x[0]) / (hx * hx); /*boundary*/
    120:   for (i = 1; i < mx - 1; i++) f[i] = eps * (x[i + 1] - 2. * x[i] + x[i - 1]) / (hx * hx);
    121:   f[mx - 1] = 2. * eps * (x[mx - 2] - x[mx - 1]) / (hx * hx); /*boundary*/
    122:   PetscCall(VecRestoreArrayRead(X, &x));
    123:   PetscCall(VecRestoreArray(F, &f));
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
    128: {
    129:   AppCtx            *user = (AppCtx *)ptr;
    130:   PetscScalar       *f;
    131:   const PetscScalar *x, *xdot;
    132:   PetscInt           i, mx;

    134:   PetscFunctionBegin;
    135:   mx = user->mx;
    136:   PetscCall(VecGetArrayRead(X, &x));
    137:   PetscCall(VecGetArrayRead(Xdot, &xdot));
    138:   PetscCall(VecGetArray(F, &f));

    140:   for (i = 0; i < mx; i++) f[i] = xdot[i] - x[i] * (1. - x[i] * x[i]);

    142:   PetscCall(VecRestoreArrayRead(X, &x));
    143:   PetscCall(VecRestoreArrayRead(Xdot, &xdot));
    144:   PetscCall(VecRestoreArray(F, &f));
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }

    148: static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, void *ctx)
    149: {
    150:   AppCtx            *user = (AppCtx *)ctx;
    151:   PetscScalar        v;
    152:   const PetscScalar *x;
    153:   PetscInt           i, col;

    155:   PetscFunctionBegin;
    156:   PetscCall(VecGetArrayRead(U, &x));
    157:   for (i = 0; i < user->mx; i++) {
    158:     v   = a - 1. + 3. * x[i] * x[i];
    159:     col = i;
    160:     PetscCall(MatSetValues(J, 1, &i, 1, &col, &v, INSERT_VALUES));
    161:   }
    162:   PetscCall(VecRestoreArrayRead(U, &x));

    164:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    165:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    166:   if (J != Jpre) {
    167:     PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    168:     PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    169:   }
    170:   /*  MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/
    171:   PetscFunctionReturn(PETSC_SUCCESS);
    172: }

    174: static PetscErrorCode FormInitialSolution(TS ts, Vec U, void *ctx)
    175: {
    176:   AppCtx      *user = (AppCtx *)ctx;
    177:   PetscInt     i;
    178:   PetscScalar *x;
    179:   PetscReal    hx, x_map;

    181:   PetscFunctionBegin;
    182:   hx = (user->xright - user->xleft) / (PetscReal)(user->mx - 1);
    183:   PetscCall(VecGetArray(U, &x));
    184:   for (i = 0; i < user->mx; i++) {
    185:     x_map = user->xleft + i * hx;
    186:     if (x_map >= 0.7065) {
    187:       x[i] = PetscTanhReal((x_map - 0.8) / (2. * PetscSqrtReal(user->param)));
    188:     } else if (x_map >= 0.4865) {
    189:       x[i] = PetscTanhReal((0.613 - x_map) / (2. * PetscSqrtReal(user->param)));
    190:     } else if (x_map >= 0.28) {
    191:       x[i] = PetscTanhReal((x_map - 0.36) / (2. * PetscSqrtReal(user->param)));
    192:     } else if (x_map >= -0.7) {
    193:       x[i] = PetscTanhReal((0.2 - x_map) / (2. * PetscSqrtReal(user->param)));
    194:     } else if (x_map >= -1) {
    195:       x[i] = PetscTanhReal((x_map + 0.9) / (2. * PetscSqrtReal(user->param)));
    196:     }
    197:   }
    198:   PetscCall(VecRestoreArray(U, &x));
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: /*TEST

    204:      test:
    205:        args: -ts_rtol 1e-04 -ts_dt 0.025 -pc_type lu -ksp_error_if_not_converged TRUE -ts_type eimex -ts_adapt_type none -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_draw_solution
    206:        requires: x
    207:        timeoutfactor: 3

    209: TEST*/