static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods.\n"; /* u_t - alpha u_xx = A + u^2 v - (B+1) u v_t - alpha v_xx = B u - u^2 v 0 < x < 1; A = 1, B = 3, alpha = 1/50 Initial conditions: u(x,0) = 1 + sin(2 pi x) v(x,0) = 3 Boundary conditions: u(0,t) = u(1,t) = 1 v(0,t) = v(1,t) = 3 */ #include #include #include typedef struct { PetscScalar u,v; } Field; typedef struct _User *User; struct _User { PetscReal A,B; /* Reaction coefficients */ PetscReal alpha; /* Diffusion coefficient */ PetscReal uleft,uright; /* Dirichlet boundary conditions */ PetscReal vleft,vright; /* Dirichlet boundary conditions */ }; static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*); static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*); static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); static PetscErrorCode FormInitialSolution(TS,Vec,void*); int main(int argc,char **argv) { TS ts; /* nonlinear solver */ Vec X; /* solution, residual vectors */ Mat J; /* Jacobian matrix */ PetscInt steps,mx; PetscErrorCode ierr; DM da; PetscReal ftime,hx,dt; struct _User user; /* user-defined work context */ TSConvergedReason reason; ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create distributed array (DMDA) to manage parallel grid and vectors - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,2,2,NULL,&da);CHKERRQ(ierr); ierr = DMSetFromOptions(da);CHKERRQ(ierr); ierr = DMSetUp(da);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Extract global vectors from DMDA; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr); /* Initialize user application context */ ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options",""); { user.A = 1; user.B = 3; user.alpha = 0.02; user.uleft = 1; user.uright = 1; user.vleft = 3; user.vright = 3; ierr = PetscOptionsReal("-A","Reaction rate","",user.A,&user.A,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-B","Reaction rate","",user.B,&user.B,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-alpha","Diffusion coefficient","",user.alpha,&user.alpha,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-uleft","Dirichlet boundary condition","",user.uleft,&user.uleft,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-uright","Dirichlet boundary condition","",user.uright,&user.uright,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-vleft","Dirichlet boundary condition","",user.vleft,&user.vleft,NULL);CHKERRQ(ierr); ierr = PetscOptionsReal("-vright","Dirichlet boundary condition","",user.vright,&user.vright,NULL);CHKERRQ(ierr); } ierr = PetscOptionsEnd();CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create timestepping solver context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); ierr = TSSetDM(ts,da);CHKERRQ(ierr); ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr); ierr = TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);CHKERRQ(ierr); ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr); ierr = TSSetIFunction(ts,NULL,FormIFunction,&user);CHKERRQ(ierr); ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr); ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr); ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr); ftime = 10.0; ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set initial conditions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr); ierr = TSSetSolution(ts,X);CHKERRQ(ierr); ierr = VecGetSize(X,&mx);CHKERRQ(ierr); hx = 1.0/(PetscReal)(mx/2-1); dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */ ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set runtime options - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve nonlinear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = TSSolve(ts,X);CHKERRQ(ierr); ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Free work space. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = VecDestroy(&X);CHKERRQ(ierr); ierr = TSDestroy(&ts);CHKERRQ(ierr); ierr = DMDestroy(&da);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; } static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr) { User user = (User)ptr; DM da; DMDALocalInfo info; PetscInt i; Field *x,*xdot,*f; PetscReal hx; Vec Xloc; PetscErrorCode ierr; PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr); ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); hx = 1.0/(PetscReal)(info.mx-1); /* Scatter ghost points to local vector,using the 2-step process DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). By placing code between these two statements, computations can be done while messages are in transition. */ ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); /* Get pointers to vector data */ ierr = DMDAVecGetArrayRead(da,Xloc,&x);CHKERRQ(ierr); ierr = DMDAVecGetArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); /* Compute function over the locally owned part of the grid */ for (i=info.xs; iuleft); f[i].v = hx * (x[i].v - user->vleft); } else if (i == info.mx-1) { f[i].u = hx * (x[i].u - user->uright); f[i].v = hx * (x[i].v - user->vright); } else { f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx; f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx; } } /* Restore vectors */ ierr = DMDAVecRestoreArrayRead(da,Xloc,&x);CHKERRQ(ierr); ierr = DMDAVecRestoreArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); PetscFunctionReturn(0); } static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) { User user = (User)ptr; DM da; DMDALocalInfo info; PetscInt i; PetscReal hx; Field *x,*f; PetscErrorCode ierr; PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr); ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); hx = 1.0/(PetscReal)(info.mx-1); /* Get pointers to vector data */ ierr = DMDAVecGetArrayRead(da,X,&x);CHKERRQ(ierr); ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); /* Compute function over the locally owned part of the grid */ for (i=info.xs; iA + u*u*v - (user->B+1)*u); f[i].v = hx*(user->B*u - u*u*v); } /* Restore vectors */ ierr = DMDAVecRestoreArrayRead(da,X,&x);CHKERRQ(ierr); ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); PetscFunctionReturn(0); } /* --------------------------------------------------------------------- */ /* IJacobian - Compute IJacobian = dF/dU + a dF/dUdot */ PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr) { User user = (User)ptr; PetscErrorCode ierr; DMDALocalInfo info; PetscInt i; PetscReal hx; DM da; Field *x,*xdot; PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr); ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); hx = 1.0/(PetscReal)(info.mx-1); /* Get pointers to vector data */ ierr = DMDAVecGetArrayRead(da,X,&x);CHKERRQ(ierr); ierr = DMDAVecGetArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); /* Compute function over the locally owned part of the grid */ for (i=info.xs; ialpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx; const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}}, {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}}; ierr = MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);CHKERRQ(ierr); } } /* Restore vectors */ ierr = DMDAVecRestoreArrayRead(da,X,&x);CHKERRQ(ierr); ierr = DMDAVecRestoreArrayRead(da,Xdot,&xdot);CHKERRQ(ierr); ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (J != Jpre) { ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); } PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx) { User user = (User)ctx; DM da; PetscInt i; DMDALocalInfo info; Field *x; PetscReal hx; PetscErrorCode ierr; PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr); ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); hx = 1.0/(PetscReal)(info.mx-1); /* Get pointers to vector data */ ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); /* Compute function over the locally owned part of the grid */ for (i=info.xs; iuleft*(1.-xi) + user->uright*xi + PetscSinReal(2.*PETSC_PI*xi); x[i].v = user->vleft*(1.-xi) + user->vright*xi; } ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); PetscFunctionReturn(0); } /*TEST test: 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 TEST*/