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  | }