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 */

 36:   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:   PetscOptionsReal("-eps", "parameter", "", user.param, &user.param, NULL);
 45:   PetscOptionsEnd();

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

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

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

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

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

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:    Free work space.
 97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   MatDestroy(&A);
 99:   VecDestroy(&x);
100:   TSDestroy(&ts);
101:   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:   mx  = user->mx;
114:   eps = user->param;
115:   hx  = (user->xright - user->xleft) / (mx - 1);
116:   VecGetArrayRead(X, &x);
117:   VecGetArray(F, &f);
118:   f[0] = 2. * eps * (x[1] - x[0]) / (hx * hx); /*boundary*/
119:   for (i = 1; i < mx - 1; i++) f[i] = eps * (x[i + 1] - 2. * x[i] + x[i - 1]) / (hx * hx);
120:   f[mx - 1] = 2. * eps * (x[mx - 2] - x[mx - 1]) / (hx * hx); /*boundary*/
121:   VecRestoreArrayRead(X, &x);
122:   VecRestoreArray(F, &f);
123:   return 0;
124: }

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

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

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

140:   VecRestoreArrayRead(X, &x);
141:   VecRestoreArrayRead(Xdot, &xdot);
142:   VecRestoreArray(F, &f);
143:   return 0;
144: }

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

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

161:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
162:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
163:   if (J != Jpre) {
164:     MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
165:     MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
166:   }
167:   /*  MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/
168:   return 0;
169: }

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

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

198: /*TEST

200:      test:
201:        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
202:        requires: x
203:        timeoutfactor: 3

205: TEST*/