1 | /***************************************
2 | $Header: /cvsroot/petscgraphics/chts.c,v 1.26 2003/07/25 15:41:51 hazelsct Exp $
3 |
4 | This is the Cahn Hilliard timestepping code. It is provided here as an
5 | example usage of the Illuminator Distributed Visualization Library.
6 | ***************************************/
7 |
8 | static char help[] = "Solves a nonlinear system in parallel using PETSc timestepping routines.\n\
9 | \n\
10 | The 2D or 3D transient Cahn-Hilliard phase field equation is solved on a 1x1\n\
11 | square or 1x1x1 cube.\n\
12 | The model is governed by the following parameters:\n\
13 | -twodee : obvious (default is 3-D)\n\
14 | -random : random initial condition (default is a box)\n\
15 | -kappa <kap>, where <kap> = diffusivity\n\
16 | -epsilon <eps>, where <eps> = interface thickness (default 1/mx)\n\
17 | -gamma <gam>, where <gam> = interfacial tension (default 1)\n\
18 | -mparam <m>, where <m> = free energy parameter, bet 0 & 1/2 for 0 stable\n\
19 | Model parameters alpha and beta are set from epsilon and gamma according to:\n\
20 | alpha = gam*eps beta=gam/eps\n\
21 | Low-level options:\n\
22 | -mx <xg>, where <xg> = number of grid points in the x-direction\n\
23 | -my <yg>, where <yg> = number of grid points in the y-direction\n\
24 | -mz <zg>, where <zg> = number of grid points in the z-direction\n\
25 | -printg : print grid information\n\
26 | Graphics of the contours of C are drawn by default at each timestep:\n\
27 | -no_contours : do not draw contour plots of solution\n\
28 | Parallelism can be invoked based on the DA construct:\n\
29 | -Nx <npx>, where <npx> = number of processors in the x-direction\n\
30 | -Ny <npy>, where <npy> = number of processors in the y-direction\n\
31 | -Nz <npz>, where <npz> = number of processors in the z-direction\n\n";
32 |
33 |
34 | /* Including cahnhill.h includes the necessary PETSc header files */
35 | #include "cahnhill.h"
36 | #include "illuminator.h"
37 |
38 |
39 | /* Set #if to 1 to debug, 0 otherwise */
40 | #undef DPRINTF
41 | #ifdef DEBUG
42 | #define DPRINTF(fmt, args...) ierr = PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args); CHKERRQ(ierr)
43 | #else
44 | #define DPRINTF(fmt, args...)
45 | #endif
46 |
47 |
48 | /* Routines given below in this file */
49 | int FormInitialCondition(AppCtx*,Vec);
50 | int UserLevelEnd(AppCtx*,Vec);
51 | int InitializeProblem(AppCtx*,Vec*);
52 |
53 |
54 | /*++++++++++++++++++++++++++++++++++++++
55 | Wrapper for
56 | +latex+{\tt ch\_residual\_vector\_2d()} in {\tt cahnhill.c}.
57 | +html+ <tt>ch_residual_vector_2d()</tt> in <tt>cahnhill.c</tt>.
58 |
59 | int ch_ts_residual_vector_2d Usual return: zero or error.
60 |
61 | TS thets Timestepping context, ignored here.
62 |
63 | PetscScalar time Current time, also ignored.
64 |
65 | Vec X Current solution vector.
66 |
67 | Vec F Function vector to return.
68 |
69 | void *ptr User data pointer.
70 | ++++++++++++++++++++++++++++++++++++++*/
71 |
72 | int ch_ts_residual_vector_2d
73 | (TS thets, PetscScalar time, Vec X, Vec F, void *ptr)
74 | { return ch_residual_vector_2d (X, F, ptr); }
75 |
76 |
77 | /*++++++++++++++++++++++++++++++++++++++
78 | Wrapper for
79 | +latex+{\tt ch\_residual\_vector\_3d()} in {\tt cahnhill.c}.
80 | +html+ <tt>ch_residual_vector_3d()</tt> in <tt>cahnhill.c</tt>.
81 |
82 | int ch_ts_residual_vector_3d Usual return: zero or error.
83 |
84 | TS thets Timestepping context, ignored here.
85 |
86 | PetscScalar time Current time, also ignored.
87 |
88 | Vec X Current solution vector.
89 |
90 | Vec F Function vector to return.
91 |
92 | void *ptr User data pointer.
93 | ++++++++++++++++++++++++++++++++++++++*/
94 |
95 | int ch_ts_residual_vector_3d
96 | (TS thets, PetscScalar time, Vec X, Vec F, void *ptr)
97 | { return ch_residual_vector_3d (X, F, ptr); }
98 |
99 |
100 | /*++++++++++++++++++++++++++++++++++++++
101 | Wrapper for
102 | +latex+{\tt ch\_jacobian\_2d()} in {\tt cahnhill.c}.
103 | +html+ <tt>ch_jacobian_2d()</tt> in <tt>cahnhill.c</tt>.
104 |
105 | int ch_ts_jacobian_2d Usual return: zero or error.
106 |
107 | TS thets Timestepping context, ignored here.
108 |
109 | PetscScalar time Current time, also ignored.
110 |
111 | Vec X Current solution vector.
112 |
113 | Mat *A Place to put the new Jacobian.
114 |
115 | Mat *B Place to put the new conditioning matrix.
116 |
117 | MatStructure *flag Flag describing the volatility of the structure.
118 |
119 | void *ptr User data pointer.
120 | ++++++++++++++++++++++++++++++++++++++*/
121 |
122 | int ch_ts_jacobian_2d (TS thets, PetscScalar time, Vec X, Mat *A, Mat *B,
123 | MatStructure *flag, void *ptr) {
124 | return ch_jacobian_2d (X, A, B, flag, ptr); }
125 |
126 |
127 | /*++++++++++++++++++++++++++++++++++++++
128 | Wrapper for
129 | +latex+{\tt ch\_jacobian\_3d()} in {\tt cahnhill.c}.
130 | +html+ <tt>ch_jacobian_3d()</tt> in <tt>cahnhill.c</tt>.
131 |
132 | int ch_ts_jacobian_3d Usual return: zero or error.
133 |
134 | TS thets Timestepping context, ignored here.
135 |
136 | PetscScalar time Current time, also ignored.
137 |
138 | Vec X Current solution vector.
139 |
140 | Mat *A Place to put the new Jacobian.
141 |
142 | Mat *B Place to put the new conditioning matrix.
143 |
144 | MatStructure *flag Flag describing the volatility of the structure.
145 |
146 | void *ptr User data pointer.
147 | ++++++++++++++++++++++++++++++++++++++*/
148 |
149 | int ch_ts_jacobian_3d (TS thets, PetscScalar time, Vec X, Mat *A, Mat *B,
150 | MatStructure *flag, void *ptr) {
151 | return ch_jacobian_3d (X, A, B, flag, ptr); }
152 |
153 |
154 | #undef __FUNCT__
155 | #define __FUNCT__ "ch_ts_monitor"
156 |
157 | /*++++++++++++++++++++++++++++++++++++++
158 | Monitor routine which displays the current state, using
159 | +latex+{\tt Illuminator}'s {\tt geomview}
160 | +html+ <tt>Illuminator</tt>'s <tt>geomview</tt>
161 | front-end (unless -no_contours is used); and also saves it using
162 | +latex+{\tt IlluMultiSave()}
163 | +html+ <tt>IlluMultiSave()</tt>
164 | if -save_data is specified.
165 |
166 | int ch_ts_monitor Usual return: zero or error.
167 |
168 | TS thets Timestepping context, ignored here.
169 |
170 | int stepno Current time step number.
171 |
172 | PetscScalar time Current time.
173 |
174 | Vec X Vector of current solved field values.
175 |
176 | void *ptr User data pointer.
177 | ++++++++++++++++++++++++++++++++++++++*/
178 |
179 | int ch_ts_monitor (TS thets, int stepno, PetscScalar time, Vec X, void *ptr)
180 | {
181 | AppCtx *user;
182 | int temp, ierr;
183 | PetscReal xmin, xmax;
184 | PetscScalar minmax[6] = { 0.,1., 0.,1., 0.,1. };
185 | /* Example colors and isoquant values to pass to DATriangulate() */
186 | /* PetscScalar colors[16] = { 1,0,0,.5, 1,1,0,.5, 0,1,0,.5, 0,0,1,.5 };
187 | PetscScalar isovals[4] = { .2, .4, .6, .8 }; */
188 |
189 | user = (AppCtx *)ptr;
190 |
191 | ierr = VecMin (X, &temp, &xmin); CHKERRQ (ierr);
192 | ierr = VecMax (X, &temp, &xmax); CHKERRQ (ierr);
193 | PetscPrintf (user->comm,"Step %d, Time %g, C values from %g to %g\n",
194 | stepno, time, xmin, xmax);
195 |
196 | if (!user->no_contours)
197 | {
198 | if (user->threedee)
199 | {
200 | DPRINTF ("Starting triangulation\n",0);
201 | ierr = DATriangulate (user->da, X, user->chvar, minmax, PETSC_DECIDE,
202 | PETSC_NULL, PETSC_NULL); CHKERRQ (ierr);
203 | DPRINTF ("Done, sending to Geomview\n",0);
204 | ierr = GeomviewDisplayTriangulation
205 | (user->comm, minmax, "Cahn-Hilliard", PETSC_TRUE); CHKERRQ (ierr);
206 | DPRINTF ("Done.\n",0);
207 | }
208 | else
209 | {
210 | ierr = VecView (X,user->theviewer); CHKERRQ (ierr);
211 | }
212 | }
213 |
214 | if (user->save_data)
215 | {
216 | PetscReal deltat;
217 | field_plot_type chtype=FIELD_SCALAR_CONTOURS;
218 | char filename [100], *paramdata [4],
219 | *paramnames [4] = { "time", "timestep", "gamma", "kappa" };
220 |
221 | for (ierr=0; ierr<4; ierr++)
222 | paramdata[ierr] = (char *) malloc (50*sizeof(char));
223 | snprintf (filename, 99, "chtsout.time%.3d", stepno);
224 | TSGetTimeStep (thets, &deltat);
225 | snprintf (paramdata[0], 49, "%g", time);
226 | snprintf (paramdata[1], 49, "%g", deltat);
227 | snprintf (paramdata[2], 49, "%g", user->gamma);
228 | snprintf (paramdata[3], 49, "%g", user->kappa);
229 | DPRINTF ("Storing data using IlluMultiSave()\n",0);
230 | IlluMultiSave (user->da, X, filename, 1.,1.,1., &chtype, 4, paramnames,
231 | paramdata, COMPRESS_INT_NONE | COMPRESS_GZIP_FAST);
232 | }
233 |
234 | DPRINTF ("Completed timestep monitor tasks.\n",0);
235 |
236 | return 0;
237 | }
238 |
239 |
240 | #undef __FUNCT__
241 | #define __FUNCT__ "main"
242 |
243 | /*++++++++++++++++++++++++++++++++++++++
244 | The usual main function.
245 |
246 | int main Returns 0 or error.
247 |
248 | int argc Number of args.
249 |
250 | char **argv The args.
251 | ++++++++++++++++++++++++++++++++++++++*/
252 |
253 | int main (int argc, char **argv)
254 | {
255 | TS thets; /* the timestepping solver */
256 | Vec x; /* solution vector */
257 | AppCtx user; /* user-defined work context */
258 | int dim, ierr;
259 | PetscDraw draw;
260 | PetscTruth matfree;
261 | PetscReal xmin, xmax;
262 | int temp;
263 | PetscScalar ftime, time_total_max = 100.0; /* default max total time */
264 | int steps = 0, time_steps_max = 5; /* default max timesteps */
265 |
266 | PetscInitialize (&argc, &argv, (char *)0, help);
267 | ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &user.rank); CHKERRQ (ierr);
268 | ierr = MPI_Comm_size (PETSC_COMM_WORLD, &user.size); CHKERRQ (ierr);
269 | user.comm = PETSC_COMM_WORLD;
270 | user.mc = 1; user.chvar = 0; /* Sets number of variables, which is conc */
271 |
272 | /* Create user context, set problem data, create vector data structures.
273 | Also, set the initial condition */
274 |
275 | DPRINTF ("About to initialize problem\n",0);
276 | ierr = InitializeProblem (&user, &x); CHKERRQ (ierr);
277 | if (user.load_data > -1)
278 | steps = user.load_data;
279 | ierr = VecDuplicate (user.localX, &user.localF); CHKERRQ (ierr);
280 |
281 | /* Set up displays to show graph of the solution */
282 |
283 | if (!user.no_contours)
284 | {
285 | if (user.threedee)
286 | {
287 | ierr = GeomviewBegin (PETSC_COMM_WORLD); CHKERRQ (ierr);
288 | }
289 | else
290 | {
291 | ierr = PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", PETSC_DECIDE,
292 | PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
293 | &user.theviewer); CHKERRQ (ierr);
294 | ierr = PetscViewerDrawGetDraw (user.theviewer, 0, &draw);
295 | CHKERRQ (ierr);
296 | ierr = PetscDrawSetDoubleBuffer (draw); CHKERRQ (ierr);
297 | }
298 | }
299 |
300 | /* Create timestepping solver context */
301 | DPRINTF ("Creating timestepping context\n",0);
302 | ierr = TSCreate (PETSC_COMM_WORLD, &thets); CHKERRQ (ierr);
303 | ierr = TSSetProblemType (thets, TS_NONLINEAR); CHKERRQ (ierr);
304 | ierr = VecGetSize (x, &dim); CHKERRQ (ierr);
305 | ierr = PetscPrintf (user.comm, "global size = %d, kappa = %g, epsilon = %g, "
306 | "gamma = %g, mparam = %g\n", dim, user.kappa,
307 | user.epsilon, user.gamma, user.mparam); CHKERRQ (ierr);
308 |
309 | /* Set function evaluation routine and monitor */
310 | DPRINTF ("Setting RHS function\n",0);
311 | if (user.threedee)
312 | ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_3d, &user);
313 | else
314 | ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_2d, &user);
315 | CHKERRQ(ierr);
316 | ierr = TSSetMonitor (thets, ch_ts_monitor, &user, PETSC_NULL); CHKERRQ(ierr);
317 |
318 | /* This condition prevents the construction of the Jacobian if we're
319 | running matrix-free. */
320 | ierr = PetscOptionsHasName (PETSC_NULL, "-snes_mf", &matfree); CHKERRQ(ierr);
321 |
322 | if (!matfree)
323 | {
324 | /* Set up the Jacobian using cahnhill.c subroutines */
325 | DPRINTF ("Using analytical Jacobian from cahnhill.c\n",0);
326 | if (user.threedee) {
327 | ierr = ch_jacobian_alpha_3d (&user); CHKERRQ (ierr);
328 | ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_3d,
329 | &user); CHKERRQ (ierr); }
330 | else {
331 | ierr = ch_jacobian_alpha_2d (&user); CHKERRQ (ierr);
332 | ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_2d,
333 | &user); CHKERRQ (ierr);}
334 | /*}*/
335 | }
336 |
337 | /* Set solution vector and initial timestep (currently a fraction of the
338 | explicit diffusion stability criterion */
339 | ierr = TSSetInitialTimeStep (thets, 0.0, 0.1/(user.mx-1)/(user.mx-1));
340 | CHKERRQ (ierr);
341 | ierr = TSSetSolution (thets, x); CHKERRQ (ierr);
342 |
343 | /* Customize timestepping solver, default to Crank-Nicholson */
344 | ierr = TSSetDuration (thets, time_steps_max, time_total_max); CHKERRQ (ierr);
345 | ierr = TSSetType (thets, TS_CRANK_NICHOLSON); CHKERRQ (ierr);
346 | ierr = TSSetFromOptions (thets); CHKERRQ (ierr);
347 |
348 | /* Run the solver */
349 | DPRINTF ("Running the solver\n",0);
350 | ierr = TSStep (thets, &steps, &ftime); CHKERRQ (ierr);
351 |
352 | /* Final clean-up */
353 | DPRINTF ("Cleaning up\n",0);
354 | if (!user.no_contours)
355 | {
356 | if (user.threedee)
357 | {
358 | ierr = GeomviewEnd (PETSC_COMM_WORLD); CHKERRQ (ierr);
359 | }
360 | else
361 | {
362 | ierr = PetscViewerDestroy (user.theviewer); CHKERRQ (ierr);
363 | }
364 | }
365 | ierr = VecDestroy (user.localX); CHKERRQ (ierr);
366 | ierr = VecDestroy (x); CHKERRQ (ierr);
367 | ierr = VecDestroy (user.localF); CHKERRQ (ierr);
368 | ierr = TSDestroy (thets); CHKERRQ (ierr);
369 | ierr = PetscFree (user.label); CHKERRQ (ierr);
370 | ierr = DADestroy (user.da); CHKERRQ (ierr);
371 |
372 | PetscFinalize ();
373 | printf ("Game over man!\n");
374 | return 0;
375 | }
376 |
377 |
378 | #undef __FUNCT__
379 | #define __FUNCT__ "FormInitialCondition"
380 |
381 | /*++++++++++++++++++++++++++++++++++++++
382 | Like it says, put together the initial condition.
383 |
384 | int FormInitialCondition Returns zero or error.
385 |
386 | AppCtx *user The user context structure.
387 |
388 | Vec X Vector in which to place the initial condition.
389 | ++++++++++++++++++++++++++++++++++++++*/
390 |
391 | int FormInitialCondition (AppCtx *user, Vec X)
392 | {
393 | int i,j,k, row, mc, chvar, mx,my,mz, ierr, xs,ys,zs, xm,ym,zm;
394 | int gxm,gym,gzm, gxs,gys,gzs;
395 | PetscScalar mparam;
396 | PetscScalar *x;
397 | Vec localX = user->localX;
398 | PetscRandom rand;
399 |
400 | mc = user->mc;
401 | chvar = user->chvar;
402 | mx = user->mx; my = user->my; mz = user->mz;
403 | mparam = user->mparam;
404 |
405 | /* Get a pointer to vector data.
406 | - For default PETSc vectors, VecGetArray() returns a pointer to
407 | the data array. Otherwise, the routine is implementation dependent.
408 | - You MUST call VecRestoreArray() when you no longer need access to
409 | the array. */
410 | if (user->random)
411 | {
412 | DPRINTF ("Setting initial condition as random numbers\n",0);
413 | ierr = PetscRandomCreate (user->comm, RANDOM_DEFAULT, &rand);
414 | CHKERRQ (ierr);
415 | ierr = PetscRandomSetInterval (rand, 0.35, 0.65); CHKERRQ (ierr);
416 | ierr = VecSetRandom (rand, X); CHKERRQ (ierr);
417 | ierr = PetscRandomDestroy (rand); CHKERRQ (ierr);
418 | }
419 | else if (user->load_data > -1)
420 | {
421 | int usermetacount;
422 | char basename [1000], **usermetanames, **usermetadata;
423 | sprintf (basename, "chtsout.time%.3d", user->load_data);
424 | DPRINTF ("Loading data for timestep %d, basename %s\n", user->load_data,
425 | basename);
426 | IlluMultiRead (user->da, X, basename, &usermetacount, &usermetanames,
427 | &usermetadata);
428 | /* TODO: free these */
429 | for (i=0; i<usermetacount; i++)
430 | PetscPrintf (PETSC_COMM_WORLD, "Parameter \"%s\" = \"%s\"\n",
431 | usermetanames [i], usermetadata [i]);
432 | }
433 | else
434 | {
435 | DPRINTF ("Getting local array\n",0);
436 | ierr = VecGetArray(localX,&x); CHKERRQ (ierr);
437 |
438 | /* Get local grid boundaries (for 2-dimensional DA):
439 | xs, ys - starting grid indices (no ghost points)
440 | xm, ym - widths of local grid (no ghost points)
441 | gxs, gys - starting grid indices (including ghost points)
442 | gxm, gym - widths of local grid (including ghost points) */
443 | DPRINTF ("Getting corners and ghost corners\n",0);
444 | ierr = DAGetCorners (user->da, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr);
445 | ierr = DAGetGhostCorners (user->da, &gxs,&gys,&gzs, &gxm,&gym,&gzm);
446 | CHKERRQ (ierr);
447 |
448 | /* Set up 2-D so it works */
449 | if (!user->threedee) { zs = gzs = 0; zm = 1; mz = 10; }
450 |
451 | /* Compute initial condition over the locally owned part of the grid
452 | Initial condition is motionless fluid and equilibrium temperature */
453 | DPRINTF ("Looping to set initial condition\n",0);
454 | for (k=zs; k<zs+zm; k++)
455 | for (j=ys; j<ys+ym; j++)
456 | for (i=xs; i<xs+xm; i++)
457 | {
458 | row = i - gxs + (j - gys)*gxm + (k-gzs)*gxm*gym;
459 | /* x[C(row)] = (i>=mx*0.43921) ? 1.0 : 0.0; */
460 | x[C(row)] = ((i<mx/3 || i>2*mx/3) && (j<my/3 || j>2*my/3) &&
461 | (k<mz/3 || k>2*mz/3)) ? 1.0 : 0.0;
462 | /* x[C(row)] = (i<mx/2 && j<my/2 && k<mz/2) ? 1.0 : 0.0; */
463 | /* x[V(row)] = 0.0;
464 | x[Omega(row)] = 0.0;
465 | x[Temp(row)] = (mparam>0)*(PetscScalar)(i)/(PetscScalar)(mx-1); */
466 | }
467 |
468 | /* Restore vector */
469 | DPRINTF ("Restoring array to local vector\n",0);
470 | ierr = VecRestoreArray (localX,&x); CHKERRQ (ierr);
471 |
472 | /* Insert values into global vector */
473 | DPRINTF ("Inserting local vector values into global vector\n",0);
474 | ierr = DALocalToGlobal (user->da,localX,INSERT_VALUES,X); CHKERRQ (ierr);
475 | }
476 |
477 | return 0;
478 | }
479 |
480 |
481 | #undef __FUNCT__
482 | #define __FUNCT__ "InitializeProblem"
483 |
484 | /*++++++++++++++++++++++++++++++++++++++
485 | This takes the gory details of initialization out of the way (importing
486 | parameters into the user context, etc.).
487 |
488 | int InitializeProblem Returns zero or error.
489 |
490 | AppCtx *user The user context to fill.
491 |
492 | Vec *xvec Vector into which to put the initial condition.
493 | ++++++++++++++++++++++++++++++++++++++*/
494 |
495 | int InitializeProblem (AppCtx *user, Vec *xvec)
496 | {
497 | int Nx,Ny,Nz; /* number of processors in x-, y- and z- directions */
498 | int xs,xm,ys,ym,zs,zm, Nlocal,ierr;
499 | Vec xv;
500 | PetscTruth twodee;
501 |
502 | /* Initialize problem parameters */
503 | DPRINTF ("Initializing user->parameters\n",0);
504 | user->mx = 20;
505 | user->my = 20;
506 | user->mz = 20;
507 | ierr = PetscOptionsGetInt(PETSC_NULL, "-mx", &user->mx, PETSC_NULL);
508 | CHKERRQ (ierr);
509 | ierr = PetscOptionsGetInt(PETSC_NULL, "-my", &user->my, PETSC_NULL);
510 | CHKERRQ (ierr);
511 | ierr = PetscOptionsGetInt(PETSC_NULL, "-mz", &user->mz, PETSC_NULL);
512 | CHKERRQ (ierr);
513 | /* No. of components in the unknown vector and auxiliary vector */
514 | user->mc = 1;
515 | /* Problem parameters (kappa, epsilon, gamma, and mparam) */
516 | user->kappa = 1.0;
517 | ierr = PetscOptionsGetScalar(PETSC_NULL, "-kappa", &user->kappa, PETSC_NULL);
518 | CHKERRQ (ierr);
519 | user->epsilon = 1.0/(user->mx-1);
520 | ierr = PetscOptionsGetScalar (PETSC_NULL, "-epsilon", &user->epsilon,
521 | PETSC_NULL); CHKERRQ (ierr);
522 | user->gamma = 1.0;
523 | ierr = PetscOptionsGetScalar(PETSC_NULL, "-gamma", &user->gamma, PETSC_NULL);
524 | CHKERRQ (ierr);
525 | user->mparam = 0.0;
526 | ierr = PetscOptionsGetScalar (PETSC_NULL, "-mparam", &user->mparam,
527 | PETSC_NULL); CHKERRQ (ierr);
528 | ierr = PetscOptionsHasName (PETSC_NULL, "-twodee", &twodee);
529 | user->threedee = !twodee;
530 | CHKERRQ (ierr);
531 | ierr = PetscOptionsHasName (PETSC_NULL, "-printv", &user->print_vecs);
532 | CHKERRQ (ierr);
533 | ierr = PetscOptionsHasName (PETSC_NULL, "-printg", &user->print_grid);
534 | CHKERRQ (ierr);
535 | ierr = PetscOptionsHasName (PETSC_NULL, "-no_contours", &user->no_contours);
536 | CHKERRQ (ierr);
537 | ierr = PetscOptionsHasName (PETSC_NULL, "-random", &user->random);
538 | CHKERRQ (ierr);
539 | ierr = PetscOptionsHasName (PETSC_NULL, "-save_data", &user->save_data);
540 | CHKERRQ (ierr);
541 | user->load_data = -1;
542 | ierr = PetscOptionsGetInt (PETSC_NULL, "-load_data", &user->load_data,
543 | PETSC_NULL); CHKERRQ (ierr);
544 |
545 | /* Create distributed array (DA) to manage parallel grid and vectors
546 | for principal unknowns (x) and governing residuals (f)
547 | Note the stencil width of 2 for this 4th-order equation. */
548 | DPRINTF ("Creating distributed array\n",0);
549 | Nx = PETSC_DECIDE;
550 | Ny = PETSC_DECIDE;
551 | Nz = PETSC_DECIDE;
552 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Nx", &Nx, PETSC_NULL);CHKERRQ(ierr);
553 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Ny", &Ny, PETSC_NULL);CHKERRQ(ierr);
554 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Nz", &Nz, PETSC_NULL);CHKERRQ(ierr);
555 | if (user->threedee)
556 | {
557 | DPRINTF ("Three Dee!\n",0);
558 | user->period = DA_XYZPERIODIC;
559 | ierr = DACreate3d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX,
560 | user->mx, user->my, user->mz, Nx, Ny, Nz, user->mc, 2,
561 | PETSC_NULL, PETSC_NULL, PETSC_NULL, &user->da);
562 | CHKERRQ (ierr);
563 | }
564 | else
565 | {
566 | user->period = DA_XYPERIODIC;
567 | ierr = DACreate2d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX,
568 | user->mx, user->my, Nx, Ny, user->mc, 2,
569 | PETSC_NULL, PETSC_NULL, &user->da); CHKERRQ (ierr);
570 | user->mz = Nz = 1;
571 | }
572 |
573 | /* Extract global vector from DA */
574 | DPRINTF ("Extracting global vector from DA...\n",0);
575 | ierr = DACreateGlobalVector(user->da,&xv); CHKERRQ (ierr);
576 |
577 | /* Label PDE components */
578 | DPRINTF ("Labeling PDE components\n",0);
579 | ierr = PetscMalloc (user->mc * sizeof(char*), &(user->label));CHKERRQ (ierr);
580 | user->label[0] = "Concentration (C)";
581 | /* user->label[1] = "Velocity (V)";
582 | user->label[2] = "Omega";
583 | user->label[3] = "Temperature"; */
584 | ierr = DASetFieldName (user->da, 0, user->label[0]); CHKERRQ(ierr);
585 |
586 | user->x_old = 0;
587 |
588 | /* Get local vector */
589 | DPRINTF ("Getting local vector\n",0);
590 | ierr = DACreateLocalVector (user->da, &user->localX); CHKERRQ (ierr);
591 |
592 | /* Print grid info */
593 | if (user->print_grid)
594 | {
595 | DPRINTF ("Printing grid information\n",0);
596 | ierr = DAView(user->da,PETSC_VIEWER_STDOUT_SELF); CHKERRQ (ierr);
597 | ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr);
598 | if (!user->threedee) {
599 | zs = 0; zm = 1; }
600 | ierr = PetscPrintf(PETSC_COMM_WORLD,
601 | "global grid: %d X %d X %d with %d components per"
602 | " node ==> global vector dimension %d\n",
603 | user->mx, user->my, user->mz, user->mc,
604 | user->mc*user->mx*user->my*user->mz);
605 | CHKERRQ (ierr);
606 | fflush(stdout);
607 | ierr = VecGetLocalSize (xv, &Nlocal); CHKERRQ (ierr);
608 | ierr = PetscSynchronizedPrintf
609 | (PETSC_COMM_WORLD,"[%d] local grid %d X %d X %d with %d components"
610 | " per node ==> local vector dimension %d\n",
611 | user->rank, xm, ym, zm, user->mc, Nlocal); CHKERRQ (ierr);
612 | ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ (ierr);
613 | }
614 |
615 | /* Compute initial condition */
616 | DPRINTF ("Forming inital condition\n",0);
617 | ierr = FormInitialCondition (user, xv); CHKERRQ (ierr);
618 |
619 | *xvec = xv;
620 | return 0;
621 | }