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 PetscReal ftime; 33 AppCtx user; /* user-defined work context */ 34 35 PetscFunctionBeginUser; 36 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 37 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 PetscCall(PetscOptionsReal("-eps","parameter","",user.param,&user.param,NULL)); 45 PetscOptionsEnd(); 46 47 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 48 Set runtime options 49 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 50 /* 51 * PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 52 */ 53 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 54 Create necessary matrix and vectors, solve same ODE on every process 55 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 56 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 57 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,user.mx,user.mx)); 58 PetscCall(MatSetFromOptions(A)); 59 PetscCall(MatSetUp(A)); 60 PetscCall(MatCreateVecs(A,&x,NULL)); 61 62 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 63 Create time stepping solver context 64 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 65 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 66 PetscCall(TSSetType(ts,TSEIMEX)); 67 PetscCall(TSSetRHSFunction(ts,NULL,RHSFunction,&user)); 68 PetscCall(TSSetIFunction(ts,NULL,FormIFunction,&user)); 69 PetscCall(TSSetIJacobian(ts,A,A,FormIJacobian,&user)); 70 ftime = 22; 71 PetscCall(TSSetMaxTime(ts,ftime)); 72 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 73 74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 75 Set initial conditions 76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 77 PetscCall(FormInitialSolution(ts,x,&user)); 78 PetscCall(TSSetSolution(ts,x)); 79 PetscCall(VecGetSize(x,&mx)); 80 81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 Set runtime options 83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 84 PetscCall(TSSetFromOptions(ts)); 85 86 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 87 Solve nonlinear system 88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 89 PetscCall(TSSolve(ts,x)); 90 PetscCall(TSGetTime(ts,&ftime)); 91 PetscCall(TSGetStepNumber(ts,&steps)); 92 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"eps %g, steps %" PetscInt_FMT ", ftime %g\n",(double)user.param,steps,(double)ftime)); 93 /* PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));*/ 94 95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96 Free work space. 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 PetscCall(MatDestroy(&A)); 99 PetscCall(VecDestroy(&x)); 100 PetscCall(TSDestroy(&ts)); 101 PetscCall(PetscFinalize()); 102 return 0; 103 } 104 105 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 106 { 107 AppCtx *user = (AppCtx*)ptr; 108 PetscScalar *f; 109 const PetscScalar *x; 110 PetscInt i,mx; 111 PetscReal hx,eps; 112 113 PetscFunctionBegin; 114 mx = user->mx; 115 eps = user->param; 116 hx = (user->xright-user->xleft)/(mx-1); 117 PetscCall(VecGetArrayRead(X,&x)); 118 PetscCall(VecGetArray(F,&f)); 119 f[0] = 2.*eps*(x[1]-x[0])/(hx*hx); /*boundary*/ 120 for (i=1;i<mx-1;i++) { 121 f[i]= eps*(x[i+1]-2.*x[i]+x[i-1])/(hx*hx); 122 } 123 f[mx-1] = 2.*eps*(x[mx-2]- x[mx-1])/(hx*hx); /*boundary*/ 124 PetscCall(VecRestoreArrayRead(X,&x)); 125 PetscCall(VecRestoreArray(F,&f)); 126 PetscFunctionReturn(0); 127 } 128 129 static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr) 130 { 131 AppCtx *user = (AppCtx*)ptr; 132 PetscScalar *f; 133 const PetscScalar *x,*xdot; 134 PetscInt i,mx; 135 136 PetscFunctionBegin; 137 mx = user->mx; 138 PetscCall(VecGetArrayRead(X,&x)); 139 PetscCall(VecGetArrayRead(Xdot,&xdot)); 140 PetscCall(VecGetArray(F,&f)); 141 142 for (i=0;i<mx;i++) { 143 f[i]= xdot[i] - x[i]*(1.-x[i]*x[i]); 144 } 145 146 PetscCall(VecRestoreArrayRead(X,&x)); 147 PetscCall(VecRestoreArrayRead(Xdot,&xdot)); 148 PetscCall(VecRestoreArray(F,&f)); 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec U, Vec Udot, PetscReal a, Mat J,Mat Jpre,void *ctx) 153 { 154 AppCtx *user = (AppCtx *)ctx; 155 PetscScalar v; 156 const PetscScalar *x; 157 PetscInt i,col; 158 159 PetscFunctionBegin; 160 PetscCall(VecGetArrayRead(U,&x)); 161 for (i=0; i < user->mx; i++) { 162 v = a - 1. + 3.*x[i]*x[i]; 163 col = i; 164 PetscCall(MatSetValues(J,1,&i,1,&col,&v,INSERT_VALUES)); 165 } 166 PetscCall(VecRestoreArrayRead(U,&x)); 167 168 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 169 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 170 if (J != Jpre) { 171 PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY)); 172 PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY)); 173 } 174 /* MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/ 175 PetscFunctionReturn(0); 176 } 177 178 static PetscErrorCode FormInitialSolution(TS ts,Vec U,void *ctx) 179 { 180 AppCtx *user = (AppCtx*)ctx; 181 PetscInt i; 182 PetscScalar *x; 183 PetscReal hx,x_map; 184 185 PetscFunctionBegin; 186 hx = (user->xright-user->xleft)/(PetscReal)(user->mx-1); 187 PetscCall(VecGetArray(U,&x)); 188 for (i=0;i<user->mx;i++) { 189 x_map = user->xleft + i*hx; 190 if (x_map >= 0.7065) { 191 x[i] = PetscTanhReal((x_map-0.8)/(2.*PetscSqrtReal(user->param))); 192 } else if (x_map >= 0.4865) { 193 x[i] = PetscTanhReal((0.613-x_map)/(2.*PetscSqrtReal(user->param))); 194 } else if (x_map >= 0.28) { 195 x[i] = PetscTanhReal((x_map-0.36)/(2.*PetscSqrtReal(user->param))); 196 } else if (x_map >= -0.7) { 197 x[i] = PetscTanhReal((0.2-x_map)/(2.*PetscSqrtReal(user->param))); 198 } else if (x_map >= -1) { 199 x[i] = PetscTanhReal((x_map+0.9)/(2.*PetscSqrtReal(user->param))); 200 } 201 } 202 PetscCall(VecRestoreArray(U,&x)); 203 PetscFunctionReturn(0); 204 } 205 206 /*TEST 207 208 test: 209 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 210 requires: x 211 timeoutfactor: 3 212 213 TEST*/ 214