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