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