Actual source code: allen_cahn.c
petsc-3.12.0 2019-09-29
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*/