Actual source code: ex25.c
petsc-3.7.3 2016-08-01
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).
8: Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and multiple solves.
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
15: Initial conditions:
16: u(x,0) = 1 + sin(2 pi x)
17: v(x,0) = 3
19: Boundary conditions:
20: u(0,t) = u(1,t) = 1
21: v(0,t) = v(1,t) = 3
22: */
24: #include <petscdm.h>
25: #include <petscdmda.h>
26: #include <petscts.h>
28: typedef struct {
29: PetscScalar u,v;
30: } Field;
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: };
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);
48: int main(int argc,char **argv)
49: {
50: PetscInt cycle;
53: MPI_Init(&argc,&argv);
54: for (cycle=0; cycle<4; cycle++) {
55: Brusselator(argc,argv,cycle);
56: if (ierr) return 1;
57: }
58: MPI_Finalize();
59: return 0;
60: }
64: PetscErrorCode Brusselator(int argc,char **argv,PetscInt cycle)
65: {
66: TS ts; /* nonlinear solver */
67: Vec X; /* solution, residual vectors */
68: Mat J; /* Jacobian matrix */
69: PetscInt steps,maxsteps,mx;
70: PetscErrorCode ierr;
71: DM da;
72: PetscReal ftime,hx,dt,xmax,xmin;
73: struct _User user; /* user-defined work context */
74: TSConvergedReason reason;
76: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Create distributed array (DMDA) to manage parallel grid and vectors
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-11,2,2,NULL,&da);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Extract global vectors from DMDA;
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: DMCreateGlobalVector(da,&X);
88: /* Initialize user application context */
89: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Advection-reaction options","");
90: {
91: user.A = 1;
92: user.B = 3;
93: user.alpha = 0.1;
94: user.uleft = 1;
95: user.uright = 1;
96: user.vleft = 3;
97: user.vright = 3;
98: PetscOptionsReal("-A","Reaction rate","",user.A,&user.A,NULL);
99: PetscOptionsReal("-B","Reaction rate","",user.B,&user.B,NULL);
100: PetscOptionsReal("-alpha","Diffusion coefficient","",user.alpha,&user.alpha,NULL);
101: PetscOptionsReal("-uleft","Dirichlet boundary condition","",user.uleft,&user.uleft,NULL);
102: PetscOptionsReal("-uright","Dirichlet boundary condition","",user.uright,&user.uright,NULL);
103: PetscOptionsReal("-vleft","Dirichlet boundary condition","",user.vleft,&user.vleft,NULL);
104: PetscOptionsReal("-vright","Dirichlet boundary condition","",user.vright,&user.vright,NULL);
105: }
106: PetscOptionsEnd();
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Create timestepping solver context
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: TSCreate(PETSC_COMM_WORLD,&ts);
112: TSSetDM(ts,da);
113: TSSetType(ts,TSARKIMEX);
114: TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
115: TSSetIFunction(ts,NULL,FormIFunction,&user);
116: DMSetMatType(da,MATAIJ);
117: DMCreateMatrix(da,&J);
118: TSSetIJacobian(ts,J,J,FormIJacobian,&user);
120: ftime = 1.0;
121: maxsteps = 10000;
122: TSSetDuration(ts,maxsteps,ftime);
123: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Set initial conditions
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: FormInitialSolution(ts,X,&user);
129: TSSetSolution(ts,X);
130: VecGetSize(X,&mx);
131: hx = 1.0/(PetscReal)(mx/2-1);
132: dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
133: dt *= PetscPowRealInt(0.2,cycle); /* Shrink the time step in convergence study. */
134: TSSetInitialTimeStep(ts,0.0,dt);
135: TSSetTolerances(ts,1e-3*PetscPowRealInt(0.5,cycle),NULL,1e-3*PetscPowRealInt(0.5,cycle),NULL);
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Set runtime options
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: TSSetFromOptions(ts);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Solve nonlinear system
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: TSSolve(ts,X);
146: TSGetSolveTime(ts,&ftime);
147: TSGetTimeStepNumber(ts,&steps);
148: TSGetConvergedReason(ts,&reason);
149: VecMin(X,NULL,&xmin);
150: VecMax(X,NULL,&xmax);
151: 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);
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Free work space.
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: MatDestroy(&J);
157: VecDestroy(&X);
158: TSDestroy(&ts);
159: DMDestroy(&da);
160: PetscFinalize();
161: return ierr;
162: }
166: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
167: {
168: User user = (User)ptr;
169: DM da;
170: DMDALocalInfo info;
171: PetscInt i;
172: Field *x,*xdot,*f;
173: PetscReal hx;
174: Vec Xloc;
178: TSGetDM(ts,&da);
179: DMDAGetLocalInfo(da,&info);
180: hx = 1.0/(PetscReal)(info.mx-1);
182: /*
183: Scatter ghost points to local vector,using the 2-step process
184: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
185: By placing code between these two statements, computations can be
186: done while messages are in transition.
187: */
188: DMGetLocalVector(da,&Xloc);
189: DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
190: DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);
192: /* Get pointers to vector data */
193: DMDAVecGetArrayRead(da,Xloc,&x);
194: DMDAVecGetArrayRead(da,Xdot,&xdot);
195: DMDAVecGetArray(da,F,&f);
197: /* Compute function over the locally owned part of the grid */
198: for (i=info.xs; i<info.xs+info.xm; i++) {
199: if (i == 0) {
200: f[i].u = hx * (x[i].u - user->uleft);
201: f[i].v = hx * (x[i].v - user->vleft);
202: } else if (i == info.mx-1) {
203: f[i].u = hx * (x[i].u - user->uright);
204: f[i].v = hx * (x[i].v - user->vright);
205: } else {
206: f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx;
207: f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx;
208: }
209: }
211: /* Restore vectors */
212: DMDAVecRestoreArrayRead(da,Xloc,&x);
213: DMDAVecRestoreArrayRead(da,Xdot,&xdot);
214: DMDAVecRestoreArray(da,F,&f);
215: DMRestoreLocalVector(da,&Xloc);
216: return(0);
217: }
221: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
222: {
223: User user = (User)ptr;
224: DM da;
225: DMDALocalInfo info;
226: PetscInt i;
227: PetscReal hx;
228: Field *x,*f;
232: TSGetDM(ts,&da);
233: DMDAGetLocalInfo(da,&info);
234: hx = 1.0/(PetscReal)(info.mx-1);
236: /* Get pointers to vector data */
237: DMDAVecGetArrayRead(da,X,&x);
238: DMDAVecGetArray(da,F,&f);
240: /* Compute function over the locally owned part of the grid */
241: for (i=info.xs; i<info.xs+info.xm; i++) {
242: PetscScalar u = x[i].u,v = x[i].v;
243: f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
244: f[i].v = hx*(user->B*u - u*u*v);
245: }
247: /* Restore vectors */
248: DMDAVecRestoreArrayRead(da,X,&x);
249: DMDAVecRestoreArray(da,F,&f);
250: return(0);
251: }
253: /* --------------------------------------------------------------------- */
254: /*
255: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
256: */
259: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr)
260: {
261: User user = (User)ptr;
263: DMDALocalInfo info;
264: PetscInt i;
265: PetscReal hx;
266: DM da;
267: Field *x,*xdot;
270: TSGetDM(ts,&da);
271: DMDAGetLocalInfo(da,&info);
272: hx = 1.0/(PetscReal)(info.mx-1);
274: /* Get pointers to vector data */
275: DMDAVecGetArrayRead(da,X,&x);
276: DMDAVecGetArrayRead(da,Xdot,&xdot);
278: /* Compute function over the locally owned part of the grid */
279: for (i=info.xs; i<info.xs+info.xm; i++) {
280: if (i == 0 || i == info.mx-1) {
281: const PetscInt row = i,col = i;
282: const PetscScalar vals[2][2] = {{hx,0},{0,hx}};
283: MatSetValuesBlocked(Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);
284: } else {
285: const PetscInt row = i,col[] = {i-1,i,i+1};
286: const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;
287: const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}},
288: {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
289: MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);
290: }
291: }
293: /* Restore vectors */
294: DMDAVecRestoreArrayRead(da,X,&x);
295: DMDAVecRestoreArrayRead(da,Xdot,&xdot);
297: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
298: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
299: if (J != Jpre) {
300: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
301: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
302: }
303: return(0);
304: }
308: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
309: {
310: User user = (User)ctx;
311: DM da;
312: PetscInt i;
313: DMDALocalInfo info;
314: Field *x;
315: PetscReal hx;
319: TSGetDM(ts,&da);
320: DMDAGetLocalInfo(da,&info);
321: hx = 1.0/(PetscReal)(info.mx-1);
323: /* Get pointers to vector data */
324: DMDAVecGetArray(da,X,&x);
326: /* Compute function over the locally owned part of the grid */
327: for (i=info.xs; i<info.xs+info.xm; i++) {
328: PetscReal xi = i*hx;
329: x[i].u = user->uleft*(1.-xi) + user->uright*xi + PetscSinReal(2.*PETSC_PI*xi);
330: x[i].v = user->vleft*(1.-xi) + user->vright*xi;
331: }
332: DMDAVecRestoreArray(da,X,&x);
333: return(0);
334: }