1 static char help[] = "Solves the time dependent Allen-Cahn equation with IMEX methods";
2
3 /*
4 * allen_cahn.c
5 *
6 * Created on: Jun 8, 2012
7 * Author: Hong Zhang
8 */
9
10 #include <petscts.h>
11
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;
20
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, PetscCtx ctx);
24 static PetscErrorCode FormInitialSolution(TS, Vec, void *);
25
main(int argc,char ** argv)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 */
34
35 PetscFunctionBeginUser;
36 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37
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();
46
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));
61
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));
73
74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75 Set initial conditions
76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77 PetscCall(FormInitialSolution(ts, x, &user));
78 PetscCall(TSSetSolution(ts, x));
79 PetscCall(VecGetSize(x, &mx));
80
81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 Set runtime options
83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84 PetscCall(TSSetFromOptions(ts));
85
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));*/
94
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 }
104
RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
112
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 }
126
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void * ptr)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;
133
134 PetscFunctionBegin;
135 mx = user->mx;
136 PetscCall(VecGetArrayRead(X, &x));
137 PetscCall(VecGetArrayRead(Xdot, &xdot));
138 PetscCall(VecGetArray(F, &f));
139
140 for (i = 0; i < mx; i++) f[i] = xdot[i] - x[i] * (1. - x[i] * x[i]);
141
142 PetscCall(VecRestoreArrayRead(X, &x));
143 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
144 PetscCall(VecRestoreArray(F, &f));
145 PetscFunctionReturn(PETSC_SUCCESS);
146 }
147
FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,PetscCtx ctx)148 static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, PetscCtx ctx)
149 {
150 AppCtx *user = (AppCtx *)ctx;
151 PetscScalar v;
152 const PetscScalar *x;
153 PetscInt i, col;
154
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));
163
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 }
173
FormInitialSolution(TS ts,Vec U,PetscCtx ctx)174 static PetscErrorCode FormInitialSolution(TS ts, Vec U, PetscCtx ctx)
175 {
176 AppCtx *user = (AppCtx *)ctx;
177 PetscInt i;
178 PetscScalar *x;
179 PetscReal hx, x_map;
180
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 }
201
202 /*TEST
203
204 test:
205 args: -ts_rtol 1e-04 -ts_time_step 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
208
209 TEST*/
210