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