1 static char help[] ="Solves the time dependent Allen-Cahn equation with IMEX methods"; 2 3 /* 4 * allen_cahn.c 5 * 6 * Created on: Jun 8, 2012 7 * Author: Hong Zhang 8 */ 9 10 #include <petscts.h> 11 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; 20 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*); 25 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 */ 35 36 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 37 38 /* Initialize user application context */ 39 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Allen-Cahn equation","");CHKERRQ(ierr); 40 user.param = 9e-4; 41 user.xleft = -1.; 42 user.xright = 2.; 43 user.mx = 400; 44 ierr = PetscOptionsReal("-eps","parameter","",user.param,&user.param,NULL);CHKERRQ(ierr); 45 ierr = PetscOptionsEnd();CHKERRQ(ierr); 46 47 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 48 Set runtime options 49 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 50 /* 51 * ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr); 52 */ 53 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 54 Create necessary matrix and vectors, solve same ODE on every process 55 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 56 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 57 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,user.mx,user.mx);CHKERRQ(ierr); 58 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 59 ierr = MatSetUp(A);CHKERRQ(ierr); 60 ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr); 61 62 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 63 Create time stepping solver context 64 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 65 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 66 ierr = TSSetType(ts,TSEIMEX);CHKERRQ(ierr); 67 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); 68 ierr = TSSetIFunction(ts,NULL,FormIFunction,&user);CHKERRQ(ierr); 69 ierr = TSSetIJacobian(ts,A,A,FormIJacobian,&user);CHKERRQ(ierr); 70 ftime = 22; 71 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 72 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 73 74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 75 Set initial conditions 76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 77 ierr = FormInitialSolution(ts,x,&user);CHKERRQ(ierr); 78 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 79 ierr = VecGetSize(x,&mx);CHKERRQ(ierr); 80 81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 Set runtime options 83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 84 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 85 86 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 87 Solve nonlinear system 88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 89 ierr = TSSolve(ts,x);CHKERRQ(ierr); 90 ierr = TSGetTime(ts,&ftime);CHKERRQ(ierr); 91 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 92 ierr = PetscPrintf(PETSC_COMM_WORLD,"eps %g, steps %D, ftime %g\n",(double)user.param,steps,(double)ftime);CHKERRQ(ierr); 93 /* ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);*/ 94 95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96 Free work space. 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 ierr = MatDestroy(&A);CHKERRQ(ierr); 99 ierr = VecDestroy(&x);CHKERRQ(ierr); 100 ierr = TSDestroy(&ts);CHKERRQ(ierr); 101 ierr = PetscFinalize(); 102 return ierr; 103 } 104 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; 113 114 PetscFunctionBegin; 115 mx = user->mx; 116 eps = user->param; 117 hx = (user->xright-user->xleft)/(mx-1); 118 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 119 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 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 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 126 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 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; 137 138 PetscFunctionBegin; 139 mx = user->mx; 140 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 141 ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr); 142 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 143 144 for (i=0;i<mx;i++) { 145 f[i]= xdot[i] - x[i]*(1.-x[i]*x[i]); 146 } 147 148 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 149 ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr); 150 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 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; 161 162 PetscFunctionBegin; 163 ierr = VecGetArrayRead(U,&x);CHKERRQ(ierr); 164 for (i=0; i < user->mx; i++) { 165 v = a - 1. + 3.*x[i]*x[i]; 166 col = i; 167 ierr = MatSetValues(J,1,&i,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 168 } 169 ierr = VecRestoreArrayRead(U,&x);CHKERRQ(ierr); 170 171 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 172 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 173 if (J != Jpre) { 174 ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 175 ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 176 } 177 /* MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/ 178 PetscFunctionReturn(0); 179 } 180 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; 187 PetscErrorCode ierr; 188 189 PetscFunctionBegin; 190 hx = (user->xright-user->xleft)/(PetscReal)(user->mx-1); 191 ierr = VecGetArray(U,&x);CHKERRQ(ierr); 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 ierr = VecRestoreArray(U,&x);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 210 /*TEST 211 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 216 217 TEST*/ 218