1 static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods.\n"; 2 /* 3 u_t - alpha u_xx = A + u^2 v - (B+1) u 4 v_t - alpha v_xx = B u - u^2 v 5 0 < x < 1; 6 A = 1, B = 3, alpha = 1/50 7 8 Initial conditions: 9 u(x,0) = 1 + sin(2 pi x) 10 v(x,0) = 3 11 12 Boundary conditions: 13 u(0,t) = u(1,t) = 1 14 v(0,t) = v(1,t) = 3 15 */ 16 17 #include <petscdm.h> 18 #include <petscdmda.h> 19 #include <petscts.h> 20 21 typedef struct { 22 PetscScalar u,v; 23 } Field; 24 25 typedef struct _User *User; 26 struct _User { 27 PetscReal A,B; /* Reaction coefficients */ 28 PetscReal alpha; /* Diffusion coefficient */ 29 PetscReal uleft,uright; /* Dirichlet boundary conditions */ 30 PetscReal vleft,vright; /* Dirichlet boundary conditions */ 31 }; 32 33 static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*); 34 static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*); 35 static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 36 static PetscErrorCode FormInitialSolution(TS,Vec,void*); 37 38 int main(int argc,char **argv) 39 { 40 TS ts; /* nonlinear solver */ 41 Vec X; /* solution, residual vectors */ 42 Mat J; /* Jacobian matrix */ 43 PetscInt steps,mx; 44 PetscErrorCode ierr; 45 DM da; 46 PetscReal ftime,hx,dt; 47 struct _User user; /* user-defined work context */ 48 TSConvergedReason reason; 49 50 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 51 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 52 Create distributed array (DMDA) to manage parallel grid and vectors 53 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 54 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,2,2,NULL,&da);CHKERRQ(ierr); 55 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 56 ierr = DMSetUp(da);CHKERRQ(ierr); 57 58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 59 Extract global vectors from DMDA; 60 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 61 ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr); 62 63 /* Initialize user application context */ 64 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options",""); 65 { 66 user.A = 1; 67 user.B = 3; 68 user.alpha = 0.02; 69 user.uleft = 1; 70 user.uright = 1; 71 user.vleft = 3; 72 user.vright = 3; 73 ierr = PetscOptionsReal("-A","Reaction rate","",user.A,&user.A,NULL);CHKERRQ(ierr); 74 ierr = PetscOptionsReal("-B","Reaction rate","",user.B,&user.B,NULL);CHKERRQ(ierr); 75 ierr = PetscOptionsReal("-alpha","Diffusion coefficient","",user.alpha,&user.alpha,NULL);CHKERRQ(ierr); 76 ierr = PetscOptionsReal("-uleft","Dirichlet boundary condition","",user.uleft,&user.uleft,NULL);CHKERRQ(ierr); 77 ierr = PetscOptionsReal("-uright","Dirichlet boundary condition","",user.uright,&user.uright,NULL);CHKERRQ(ierr); 78 ierr = PetscOptionsReal("-vleft","Dirichlet boundary condition","",user.vleft,&user.vleft,NULL);CHKERRQ(ierr); 79 ierr = PetscOptionsReal("-vright","Dirichlet boundary condition","",user.vright,&user.vright,NULL);CHKERRQ(ierr); 80 } 81 ierr = PetscOptionsEnd();CHKERRQ(ierr); 82 83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 84 Create timestepping solver context 85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 86 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 87 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 88 ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr); 89 ierr = TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);CHKERRQ(ierr); 90 ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr); 91 ierr = TSSetIFunction(ts,NULL,FormIFunction,&user);CHKERRQ(ierr); 92 ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr); 93 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr); 94 ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr); 95 96 ftime = 10.0; 97 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 98 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 99 100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 Set initial conditions 102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 103 ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr); 104 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 105 ierr = VecGetSize(X,&mx);CHKERRQ(ierr); 106 hx = 1.0/(PetscReal)(mx/2-1); 107 dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */ 108 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 109 110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 111 Set runtime options 112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 113 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 114 115 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 116 Solve nonlinear system 117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 118 ierr = TSSolve(ts,X);CHKERRQ(ierr); 119 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 120 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 121 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 122 ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);CHKERRQ(ierr); 123 124 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 125 Free work space. 126 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 127 ierr = MatDestroy(&J);CHKERRQ(ierr); 128 ierr = VecDestroy(&X);CHKERRQ(ierr); 129 ierr = TSDestroy(&ts);CHKERRQ(ierr); 130 ierr = DMDestroy(&da);CHKERRQ(ierr); 131 ierr = PetscFinalize(); 132 return ierr; 133 } 134 135 static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr) 136 { 137 User user = (User)ptr; 138 DM da; 139 DMDALocalInfo info; 140 PetscInt i; 141 Field *x,*xdot,*f; 142 PetscReal hx; 143 Vec Xloc; 144 PetscErrorCode ierr; 145 146 PetscFunctionBeginUser; 147 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 148 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 149 hx = 1.0/(PetscReal)(info.mx-1); 150 151 /* 152 Scatter ghost points to local vector,using the 2-step process 153 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 154 By placing code between these two statements, computations can be 155 done while messages are in transition. 156 */ 157 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); 158 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 159 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 160 161 /* Get pointers to vector data */ 162 ierr = DMDAVecGetArrayRead(da,Xloc,&x);CHKERRQ(ierr); 163 ierr = DMDAVecGetArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); 164 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 165 166 /* Compute function over the locally owned part of the grid */ 167 for (i=info.xs; i<info.xs+info.xm; i++) { 168 if (i == 0) { 169 f[i].u = hx * (x[i].u - user->uleft); 170 f[i].v = hx * (x[i].v - user->vleft); 171 } else if (i == info.mx-1) { 172 f[i].u = hx * (x[i].u - user->uright); 173 f[i].v = hx * (x[i].v - user->vright); 174 } else { 175 f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx; 176 f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx; 177 } 178 } 179 180 /* Restore vectors */ 181 ierr = DMDAVecRestoreArrayRead(da,Xloc,&x);CHKERRQ(ierr); 182 ierr = DMDAVecRestoreArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); 183 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 184 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 185 PetscFunctionReturn(0); 186 } 187 188 static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 189 { 190 User user = (User)ptr; 191 DM da; 192 DMDALocalInfo info; 193 PetscInt i; 194 PetscReal hx; 195 Field *x,*f; 196 PetscErrorCode ierr; 197 198 PetscFunctionBeginUser; 199 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 200 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 201 hx = 1.0/(PetscReal)(info.mx-1); 202 203 /* Get pointers to vector data */ 204 ierr = DMDAVecGetArrayRead(da,X,&x);CHKERRQ(ierr); 205 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 206 207 /* Compute function over the locally owned part of the grid */ 208 for (i=info.xs; i<info.xs+info.xm; i++) { 209 PetscScalar u = x[i].u,v = x[i].v; 210 f[i].u = hx*(user->A + u*u*v - (user->B+1)*u); 211 f[i].v = hx*(user->B*u - u*u*v); 212 } 213 214 /* Restore vectors */ 215 ierr = DMDAVecRestoreArrayRead(da,X,&x);CHKERRQ(ierr); 216 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 /* --------------------------------------------------------------------- */ 221 /* 222 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot 223 */ 224 PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr) 225 { 226 User user = (User)ptr; 227 PetscErrorCode ierr; 228 DMDALocalInfo info; 229 PetscInt i; 230 PetscReal hx; 231 DM da; 232 Field *x,*xdot; 233 234 PetscFunctionBeginUser; 235 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 236 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 237 hx = 1.0/(PetscReal)(info.mx-1); 238 239 /* Get pointers to vector data */ 240 ierr = DMDAVecGetArrayRead(da,X,&x);CHKERRQ(ierr); 241 ierr = DMDAVecGetArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); 242 243 /* Compute function over the locally owned part of the grid */ 244 for (i=info.xs; i<info.xs+info.xm; i++) { 245 if (i == 0 || i == info.mx-1) { 246 const PetscInt row = i,col = i; 247 const PetscScalar vals[2][2] = {{hx,0},{0,hx}}; 248 ierr = MatSetValuesBlocked(Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);CHKERRQ(ierr); 249 } else { 250 const PetscInt row = i,col[] = {i-1,i,i+1}; 251 const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx; 252 const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}}, 253 {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}}; 254 ierr = MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);CHKERRQ(ierr); 255 } 256 } 257 258 /* Restore vectors */ 259 ierr = DMDAVecRestoreArrayRead(da,X,&x);CHKERRQ(ierr); 260 ierr = DMDAVecRestoreArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); 261 262 ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 263 ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 264 if (J != Jpre) { 265 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 266 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 267 } 268 PetscFunctionReturn(0); 269 } 270 271 PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx) 272 { 273 User user = (User)ctx; 274 DM da; 275 PetscInt i; 276 DMDALocalInfo info; 277 Field *x; 278 PetscReal hx; 279 PetscErrorCode ierr; 280 281 PetscFunctionBeginUser; 282 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 283 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 284 hx = 1.0/(PetscReal)(info.mx-1); 285 286 /* Get pointers to vector data */ 287 ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 288 289 /* Compute function over the locally owned part of the grid */ 290 for (i=info.xs; i<info.xs+info.xm; i++) { 291 PetscReal xi = i*hx; 292 x[i].u = user->uleft*(1.-xi) + user->uright*xi + PetscSinReal(2.*PETSC_PI*xi); 293 x[i].v = user->vleft*(1.-xi) + user->vright*xi; 294 } 295 ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 296 PetscFunctionReturn(0); 297 } 298 299 /*TEST 300 301 test: 302 args: -nox -da_grid_x 20 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none 303 304 TEST*/ 305