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: PetscErrorCode ierr;
33: PetscReal ftime;
34: AppCtx user; /* user-defined work context */
36: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
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 %D, 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 ierr;
103: }
105: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
106: {
107: PetscErrorCode ierr;
108: AppCtx *user = (AppCtx*)ptr;
109: PetscScalar *f;
110: const PetscScalar *x;
111: PetscInt i,mx;
112: PetscReal hx,eps;
115: mx = user->mx;
116: eps = user->param;
117: hx = (user->xright-user->xleft)/(mx-1);
118: VecGetArrayRead(X,&x);
119: VecGetArray(F,&f);
120: f[0] = 2.*eps*(x[1]-x[0])/(hx*hx); /*boundary*/
121: for (i=1;i<mx-1;i++) {
122: f[i]= eps*(x[i+1]-2.*x[i]+x[i-1])/(hx*hx);
123: }
124: f[mx-1] = 2.*eps*(x[mx-2]- x[mx-1])/(hx*hx); /*boundary*/
125: VecRestoreArrayRead(X,&x);
126: VecRestoreArray(F,&f);
127: return(0);
128: }
130: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
131: {
132: PetscErrorCode ierr;
133: AppCtx *user = (AppCtx*)ptr;
134: PetscScalar *f;
135: const PetscScalar *x,*xdot;
136: PetscInt i,mx;
139: mx = user->mx;
140: VecGetArrayRead(X,&x);
141: VecGetArrayRead(Xdot,&xdot);
142: VecGetArray(F,&f);
144: for (i=0;i<mx;i++) {
145: f[i]= xdot[i] - x[i]*(1.-x[i]*x[i]);
146: }
148: VecRestoreArrayRead(X,&x);
149: VecRestoreArrayRead(Xdot,&xdot);
150: VecRestoreArray(F,&f);
151: return(0);
152: }
154: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U, Vec Udot, PetscReal a, Mat J,Mat Jpre,void *ctx)
155: {
156: PetscErrorCode ierr;
157: AppCtx *user = (AppCtx *)ctx;
158: PetscScalar v;
159: const PetscScalar *x;
160: PetscInt i,col;
163: VecGetArrayRead(U,&x);
164: for (i=0; i < user->mx; i++) {
165: v = a - 1. + 3.*x[i]*x[i];
166: col = i;
167: MatSetValues(J,1,&i,1,&col,&v,INSERT_VALUES);
168: }
169: VecRestoreArrayRead(U,&x);
171: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
172: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
173: if (J != Jpre) {
174: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
175: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
176: }
177: /* MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/
178: return(0);
179: }
181: static PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ctx)
182: {
183: AppCtx *user = (AppCtx*)ctx;
184: PetscInt i;
185: PetscScalar *x;
186: PetscReal hx,x_map;
190: hx = (user->xright-user->xleft)/(PetscReal)(user->mx-1);
191: VecGetArray(U,&x);
192: for (i=0;i<user->mx;i++) {
193: x_map = user->xleft + i*hx;
194: if (x_map >= 0.7065) {
195: x[i] = PetscTanhReal((x_map-0.8)/(2.*PetscSqrtReal(user->param)));
196: } else if (x_map >= 0.4865) {
197: x[i] = PetscTanhReal((0.613-x_map)/(2.*PetscSqrtReal(user->param)));
198: } else if (x_map >= 0.28) {
199: x[i] = PetscTanhReal((x_map-0.36)/(2.*PetscSqrtReal(user->param)));
200: } else if (x_map >= -0.7) {
201: x[i] = PetscTanhReal((0.2-x_map)/(2.*PetscSqrtReal(user->param)));
202: } else if (x_map >= -1) {
203: x[i] = PetscTanhReal((x_map+0.9)/(2.*PetscSqrtReal(user->param)));
204: }
205: }
206: VecRestoreArray(U,&x);
207: return(0);
208: }
210: /*TEST
212: test:
213: 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
214: requires: x
215: timeoutfactor: 3
217: TEST*/