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