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
main(int argc,char ** argv)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
Brusselator(int argc,char ** argv,PetscInt cycle)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
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void * ptr)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
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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 */
FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void * ptr)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
FormInitialSolution(TS ts,Vec X,void * ctx)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