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