Actual source code: allen_cahn.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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*);


 28: int main(int argc, char **argv)
 29: {
 30:   TS                ts;
 31:   Vec               x; /*solution vector*/
 32:   Mat               A; /*Jacobian*/
 33:   PetscInt          steps,mx;
 34:   PetscErrorCode    ierr;
 35:   PetscReal         ftime;
 36:   AppCtx            user;       /* user-defined work context */

 38:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;

 40:   /* Initialize user application context */
 41:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Allen-Cahn equation","");
 42:   user.param = 9e-4;
 43:   user.xleft = -1.;
 44:   user.xright = 2.;
 45:   user.mx = 400;
 46:   PetscOptionsReal("-eps","parameter","",user.param,&user.param,NULL);
 47:   PetscOptionsEnd();

 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:    Set runtime options
 51:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 52:   /*
 53:    * PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
 54:    */
 55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:    Create necessary matrix and vectors, solve same ODE on every process
 57:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 58:   MatCreate(PETSC_COMM_WORLD,&A);
 59:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,user.mx,user.mx);
 60:   MatSetFromOptions(A);
 61:   MatSetUp(A);
 62:   MatCreateVecs(A,&x,NULL);

 64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:    Create time stepping solver context
 66:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 67:   TSCreate(PETSC_COMM_WORLD,&ts);
 68:   TSSetType(ts,TSEIMEX);
 69:   TSSetRHSFunction(ts,NULL,RHSFunction,&user);
 70:   TSSetIFunction(ts,NULL,FormIFunction,&user);
 71:   TSSetIJacobian(ts,A,A,FormIJacobian,&user);
 72:   ftime = 22;
 73:   TSSetMaxTime(ts,ftime);
 74:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

 76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:    Set initial conditions
 78:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 79:   FormInitialSolution(ts,x,&user);
 80:   TSSetSolution(ts,x);
 81:   VecGetSize(x,&mx);

 83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84:    Set runtime options
 85:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 86:   TSSetFromOptions(ts);

 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:    Solve nonlinear system
 90:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 91:   TSSolve(ts,x);
 92:   TSGetTime(ts,&ftime);
 93:   TSGetStepNumber(ts,&steps);
 94:   PetscPrintf(PETSC_COMM_WORLD,"eps %g, steps %D, ftime %g\n",(double)user.param,steps,(double)ftime);
 95:   /*   VecView(x,PETSC_VIEWER_STDOUT_WORLD);*/

 97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 98:    Free work space.
 99:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100:   MatDestroy(&A);
101:   VecDestroy(&x);
102:   TSDestroy(&ts);
103:   PetscFinalize();
104:   return ierr;
105: }


108: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
109: {
110:   PetscErrorCode    ierr;
111:   AppCtx            *user = (AppCtx*)ptr;
112:   PetscScalar       *f;
113:   const PetscScalar *x;
114:   PetscInt          i,mx;
115:   PetscReal         hx,eps;

118:   mx = user->mx;
119:   eps = user->param;
120:   hx = (user->xright-user->xleft)/(mx-1);
121:   VecGetArrayRead(X,&x);
122:   VecGetArray(F,&f);
123:   f[0] = 2.*eps*(x[1]-x[0])/(hx*hx); /*boundary*/
124:   for(i=1;i<mx-1;i++) {
125:     f[i]= eps*(x[i+1]-2.*x[i]+x[i-1])/(hx*hx);
126:   }
127:   f[mx-1] = 2.*eps*(x[mx-2]- x[mx-1])/(hx*hx); /*boundary*/
128:   VecRestoreArrayRead(X,&x);
129:   VecRestoreArray(F,&f);
130:   return(0);
131: }


134: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
135: {
136:   PetscErrorCode    ierr;
137:   AppCtx            *user = (AppCtx*)ptr;
138:   PetscScalar       *f;
139:   const PetscScalar *x,*xdot;
140:   PetscInt          i,mx;

143:   mx = user->mx;
144:   VecGetArrayRead(X,&x);
145:   VecGetArrayRead(Xdot,&xdot);
146:   VecGetArray(F,&f);

148:   for(i=0;i<mx;i++) {
149:     f[i]= xdot[i] - x[i]*(1.-x[i]*x[i]);
150:   }

152:   VecRestoreArrayRead(X,&x);
153:   VecRestoreArrayRead(Xdot,&xdot);
154:   VecRestoreArray(F,&f);
155:   return(0);
156: }


159: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U, Vec Udot, PetscReal a, Mat J,Mat Jpre,void *ctx)
160: {
161:   PetscErrorCode    ierr;
162:   AppCtx            *user = (AppCtx *)ctx;
163:   PetscScalar       v;
164:   const PetscScalar *x;
165:   PetscInt          i,col;

168:   VecGetArrayRead(U,&x);
169:   for(i=0; i < user->mx; i++) {
170:     v = a - 1. + 3.*x[i]*x[i];
171:     col = i;
172:     MatSetValues(J,1,&i,1,&col,&v,INSERT_VALUES);
173:   }
174:   VecRestoreArrayRead(U,&x);

176:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
177:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
178:   if (J != Jpre) {
179:     MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
180:     MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
181:   }
182:   /*  MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/
183:   return(0);
184: }


187: static PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ctx)
188: {
189:   AppCtx *user = (AppCtx*)ctx;
190:   PetscInt       i;
191:   PetscScalar    *x;
192:   PetscReal      hx,x_map;

196:   hx = (user->xright-user->xleft)/(PetscReal)(user->mx-1);
197:   VecGetArray(U,&x);
198:   for(i=0;i<user->mx;i++) {
199:     x_map = user->xleft + i*hx;
200:     if(x_map >= 0.7065) {
201:       x[i] = PetscTanhReal((x_map-0.8)/(2.*PetscSqrtReal(user->param)));
202:     } else if(x_map >= 0.4865) {
203:       x[i] = PetscTanhReal((0.613-x_map)/(2.*PetscSqrtReal(user->param)));
204:     } else if(x_map >= 0.28) {
205:       x[i] = PetscTanhReal((x_map-0.36)/(2.*PetscSqrtReal(user->param)));
206:     } else if(x_map >= -0.7) {
207:       x[i] = PetscTanhReal((0.2-x_map)/(2.*PetscSqrtReal(user->param)));
208:     } else if(x_map >= -1) {
209:       x[i] = PetscTanhReal((x_map+0.9)/(2.*PetscSqrtReal(user->param)));
210:     }
211:   }
212:   VecRestoreArray(U,&x);
213:   return(0);
214: }

216: /*TEST

218:      test:
219:        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
220:        requires: x
221:        timeoutfactor: 3

223: TEST*/