Actual source code: allen_cahn.c
petsc-3.5.4 2015-05-23
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;
22: static PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
23: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
24: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void *ctx);
25: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
30: int main(int argc, char **argv)
31: {
32: TS ts;
33: Vec x; /*solution vector*/
34: Mat A; /*Jacobian*/
35: PetscInt steps,maxsteps,mx;
36: PetscErrorCode ierr;
37: PetscReal ftime;
38: AppCtx user; /* user-defined work context */
40: PetscInitialize(&argc,&argv,NULL,help);
42: /* Initialize user application context */
43: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Allen-Cahn equation","");
44: user.param = 9e-4;
45: user.xleft = -1.;
46: user.xright = 2.;
47: user.mx = 400;
48: PetscOptionsReal("-eps","parameter","",user.param,&user.param,NULL);
49: PetscOptionsEnd();
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Set runtime options
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: /*
55: * PetscOptionsGetBool(NULL,"-monitor",&monitor,NULL);
56: */
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create necessary matrix and vectors, solve same ODE on every process
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: MatCreate(PETSC_COMM_WORLD,&A);
61: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,user.mx,user.mx);
62: MatSetFromOptions(A);
63: MatSetUp(A);
64: MatGetVecs(A,&x,NULL);
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create time stepping solver context
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69: TSCreate(PETSC_COMM_WORLD,&ts);
70: TSSetType(ts,TSEIMEX);
71: TSSetRHSFunction(ts,NULL,RHSFunction,&user);
72: TSSetIFunction(ts,NULL,FormIFunction,&user);
73: TSSetIJacobian(ts,A,A,FormIJacobian,&user);
74: ftime = 142;
75: maxsteps = 100000;
76: TSSetDuration(ts,maxsteps,ftime);
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Set initial conditions
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: FormInitialSolution(ts,x,&user);
82: TSSetSolution(ts,x);
83: VecGetSize(x,&mx);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Set runtime options
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: TSSetFromOptions(ts);
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Solve nonlinear system
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: TSSolve(ts,x);
94: TSGetTime(ts,&ftime);
95: TSGetTimeStepNumber(ts,&steps);
96: PetscPrintf(PETSC_COMM_WORLD,"eps %g, steps %D, ftime %g\n",(double)user.param,steps,(double)ftime);
97: /* VecView(x,PETSC_VIEWER_STDOUT_WORLD);*/
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Free work space.
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: MatDestroy(&A);
103: VecDestroy(&x);
104: TSDestroy(&ts);
105: PetscFinalize();
106: return(0);
107: }
112: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
113: {
115: AppCtx *user = (AppCtx*)ptr;
116: PetscScalar *x,*f;
117: PetscInt i,mx;
118: PetscReal hx,eps;
121: mx = user->mx;
122: eps = user->param;
123: hx = (user->xright-user->xleft)/(mx-1);
124: VecGetArray(X,&x);
125: VecGetArray(F,&f);
126: f[0] = 2.*eps*(x[1]-x[0])/(hx*hx); /*boundary*/
127: for(i=1;i<mx-1;i++) {
128: f[i]= eps*(x[i+1]-2.*x[i]+x[i-1])/(hx*hx);
129: }
130: f[mx-1] = 2.*eps*(x[mx-2]- x[mx-1])/(hx*hx); /*boundary*/
131: VecRestoreArray(X,&x);
132: VecRestoreArray(F,&f);
133: return(0);
134: }
139: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
140: {
142: AppCtx *user = (AppCtx*)ptr;
143: PetscScalar *x,*xdot,*f;
144: PetscInt i,mx;
147: mx = user->mx;
148: VecGetArray(X,&x);
149: VecGetArray(Xdot,&xdot);
150: VecGetArray(F,&f);
152: for(i=0;i<mx;i++) {
153: f[i]= xdot[i] - x[i]*(1.-x[i]*x[i]);
154: }
156: VecRestoreArray(X,&x);
157: VecRestoreArray(F,&f);
158: return(0);
159: }
164: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U, Vec Udot, PetscReal a, Mat J,Mat Jpre,void *ctx)
165: {
167: AppCtx *user = (AppCtx *)ctx;
168: PetscScalar *x,v;
169: PetscInt i,col;
171: VecGetArray(U,&x);
174: for(i=0; i < user->mx; i++) {
175: v = a - 1. + 3.*x[i]*x[i];
176: col = i;
177: MatSetValues(J,1,&i,1,&col,&v,INSERT_VALUES);
178: }
180: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
181: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
182: if (J != Jpre) {
183: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
184: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
185: }
186: /* MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/
187: return(0);
188: }
193: static PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ctx)
194: {
195: AppCtx *user = (AppCtx*)ctx;
196: PetscInt i;
197: PetscScalar *x;
198: PetscReal hx,x_map;
202: hx = (user->xright-user->xleft)/(PetscReal)(user->mx-1);
203: VecGetArray(U,&x);
204: for(i=0;i<user->mx;i++) {
205: x_map = user->xleft + i*hx;
206: if(x_map >= 0.7065) {
207: x[i] = tanh((x_map-0.8)/(2.*PetscSqrtReal(user->param)));
208: } else if(x_map >= 0.4865) {
209: x[i] = tanh((0.613-x_map)/(2.*PetscSqrtReal(user->param)));
210: } else if(x_map >= 0.28) {
211: x[i] = tanh((x_map-0.36)/(2.*PetscSqrtReal(user->param)));
212: } else if(x_map >= -0.7) {
213: x[i] = tanh((0.2-x_map)/(2.*PetscSqrtReal(user->param)));
214: } else if(x_map >= -1) {
215: x[i] = tanh((x_map+0.9)/(2.*PetscSqrtReal(user->param)));
216: }
217: }
218: VecRestoreArray(U,&x);
219: return(0);
220: }