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