1    | /***************************************
2    |   $Header: /cvsroot/petscgraphics/chts.c,v 1.28 2004/06/30 15:38:17 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, PETSC_FALSE,
203  | 				PETSC_FALSE, PETSC_FALSE); CHKERRQ (ierr);
204  | 	  DPRINTF ("Done, sending to Geomview\n",0);
205  | 	  ierr = GeomviewDisplayTriangulation
206  | 	    (user->comm, minmax, "Cahn-Hilliard", PETSC_TRUE); CHKERRQ (ierr);
207  | 	  DPRINTF ("Done.\n",0);
208  | 	}
209  |       else
210  | 	{
211  | 	  ierr = VecView (X,user->theviewer); CHKERRQ (ierr);
212  | 	}
213  |     }
214  | 
215  |   if (user->save_data)
216  |     {
217  |       PetscReal deltat;
218  |       field_plot_type chtype=FIELD_SCALAR;
219  |       char filename [100], *paramdata [4],
220  | 	*paramnames [4] = { "time", "timestep", "gamma", "kappa" };
221  | 
222  |       for (ierr=0; ierr<4; ierr++)
223  | 	paramdata[ierr] = (char *) malloc (50*sizeof(char));
224  |       snprintf (filename, 99, "chtsout.time%.3d", stepno);
225  |       TSGetTimeStep (thets, &deltat);
226  |       snprintf (paramdata[0], 49, "%g", time);
227  |       snprintf (paramdata[1], 49, "%g", deltat);
228  |       snprintf (paramdata[2], 49, "%g", user->gamma);
229  |       snprintf (paramdata[3], 49, "%g", user->kappa);
230  |       DPRINTF ("Storing data using IlluMultiSave()\n",0);
231  |       IlluMultiSave (user->da, X, filename, 1.,1.,1., &chtype, 4, paramnames,
232  | 		     paramdata, COMPRESS_INT_NONE | COMPRESS_GZIP_FAST);
233  |     }
234  | 
235  |   DPRINTF ("Completed timestep monitor tasks.\n",0);
236  | 
237  |   return 0;
238  | }
239  | 
240  | 
241  | #undef __FUNCT__
242  | #define __FUNCT__ "main"
243  | 
244  | /*++++++++++++++++++++++++++++++++++++++
245  |   The usual main function.
246  | 
247  |   int main Returns 0 or error.
248  | 
249  |   int argc Number of args.
250  | 
251  |   char **argv The args.
252  |   ++++++++++++++++++++++++++++++++++++++*/
253  | 
254  | int main (int argc, char **argv)
255  | {
256  |   TS            thets;               /* the timestepping solver */
257  |   Vec           x;                   /* solution vector */
258  |   AppCtx        user;                /* user-defined work context */
259  |   int           dim, ierr;
260  |   PetscDraw     draw;
261  |   PetscTruth    matfree;
262  |   PetscReal     xmin, xmax;
263  |   int           temp;
264  |   PetscScalar   ftime, time_total_max = 100.0; /* default max total time */
265  |   int           steps = 0, time_steps_max = 5; /* default max timesteps */
266  | 
267  |   PetscInitialize (&argc, &argv, (char *)0, help);
268  |   ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &user.rank); CHKERRQ (ierr);
269  |   ierr = MPI_Comm_size (PETSC_COMM_WORLD, &user.size); CHKERRQ (ierr);
270  |   user.comm = PETSC_COMM_WORLD;
271  |   user.mc = 1; user.chvar = 0; /* Sets number of variables, which is conc */
272  | 
273  |   /* Create user context, set problem data, create vector data structures.
274  |      Also, set the initial condition */
275  | 
276  |   DPRINTF ("About to initialize problem\n",0);
277  |   ierr = InitializeProblem (&user, &x); CHKERRQ (ierr);
278  |   if (user.load_data > -1)
279  |     steps = user.load_data;
280  |   ierr = VecDuplicate (user.localX, &user.localF); CHKERRQ (ierr);
281  | 
282  |   /* Set up displays to show graph of the solution */
283  | 
284  |   if (!user.no_contours)
285  |     {
286  |       if (user.threedee)
287  | 	{
288  | 	  ierr = GeomviewBegin (PETSC_COMM_WORLD); CHKERRQ (ierr);
289  | 	}
290  |       else
291  | 	{
292  | 	  ierr = PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", PETSC_DECIDE,
293  | 				      PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
294  | 				      &user.theviewer); CHKERRQ (ierr);
295  | 	  ierr = PetscViewerDrawGetDraw (user.theviewer, 0, &draw);
296  | 	  CHKERRQ (ierr);
297  | 	  ierr = PetscDrawSetDoubleBuffer (draw); CHKERRQ (ierr);
298  | 	}
299  |     }
300  | 
301  |   /* Create timestepping solver context */
302  |   DPRINTF ("Creating timestepping context\n",0);
303  |   ierr = TSCreate (PETSC_COMM_WORLD, &thets); CHKERRQ (ierr);
304  |   ierr = TSSetProblemType (thets, TS_NONLINEAR); CHKERRQ (ierr);
305  |   ierr = VecGetSize (x, &dim); CHKERRQ (ierr);
306  |   ierr = PetscPrintf (user.comm, "global size = %d, kappa = %g, epsilon = %g, "
307  | 		      "gamma = %g, mparam = %g\n", dim, user.kappa,
308  | 		      user.epsilon, user.gamma, user.mparam); CHKERRQ (ierr);
309  | 
310  |   /* Set function evaluation routine and monitor */
311  |   DPRINTF ("Setting RHS function\n",0);
312  |   if (user.threedee)
313  |     ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_3d, &user);
314  |   else
315  |     ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_2d, &user);
316  |   CHKERRQ(ierr);
317  |   ierr = TSSetMonitor (thets, ch_ts_monitor, &user, PETSC_NULL); CHKERRQ(ierr);
318  | 
319  |   /* This condition prevents the construction of the Jacobian if we're
320  |      running matrix-free. */
321  |   ierr = PetscOptionsHasName (PETSC_NULL, "-snes_mf", &matfree); CHKERRQ(ierr);
322  | 
323  |   if (!matfree)
324  |     {
325  |       /* Set up the Jacobian using cahnhill.c subroutines */
326  |       DPRINTF ("Using analytical Jacobian from cahnhill.c\n",0);
327  |       if (user.threedee) {
328  | 	ierr = ch_jacobian_alpha_3d (&user); CHKERRQ (ierr);
329  | 	ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_3d,
330  | 				 &user); CHKERRQ (ierr); }
331  |       else {
332  | 	ierr = ch_jacobian_alpha_2d (&user); CHKERRQ (ierr);
333  | 	ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_2d,
334  | 				 &user); CHKERRQ (ierr);}
335  |       /*}*/
336  |     }
337  | 
338  |   /* Set solution vector and initial timestep (currently a fraction of the
339  |      explicit diffusion stability criterion */
340  |   ierr = TSSetInitialTimeStep (thets, 0.0, 0.1/(user.mx-1)/(user.mx-1));
341  |   CHKERRQ (ierr);
342  |   ierr = TSSetSolution (thets, x); CHKERRQ (ierr);
343  | 
344  |   /* Customize timestepping solver, default to Crank-Nicholson */
345  |   ierr = TSSetDuration (thets, time_steps_max, time_total_max); CHKERRQ (ierr);
346  |   ierr = TSSetType (thets, TS_CRANK_NICHOLSON); CHKERRQ (ierr);
347  |   ierr = TSSetFromOptions (thets); CHKERRQ (ierr);
348  | 
349  |   /* Run the solver */
350  |   DPRINTF ("Running the solver\n",0);
351  |   ierr = TSStep (thets, &steps, &ftime); CHKERRQ (ierr);
352  | 
353  |   /* Final clean-up */
354  |   DPRINTF ("Cleaning up\n",0);
355  |   if (!user.no_contours)
356  |     {
357  |       if (user.threedee)
358  | 	{
359  | 	  ierr = GeomviewEnd (PETSC_COMM_WORLD); CHKERRQ (ierr);
360  | 	}
361  |       else
362  | 	{
363  | 	  ierr = PetscViewerDestroy (user.theviewer); CHKERRQ (ierr);
364  | 	}
365  |     }
366  |   ierr = VecDestroy (user.localX); CHKERRQ (ierr);
367  |   ierr = VecDestroy (x); CHKERRQ (ierr);
368  |   ierr = VecDestroy (user.localF); CHKERRQ (ierr);
369  |   ierr = TSDestroy (thets); CHKERRQ (ierr);  
370  |   ierr = PetscFree (user.label); CHKERRQ (ierr);
371  |   ierr = DADestroy (user.da); CHKERRQ (ierr);
372  | 
373  |   PetscFinalize ();
374  |   printf ("Game over man!\n");
375  |   return 0;
376  | }
377  | 
378  | 
379  | #undef __FUNCT__
380  | #define __FUNCT__ "FormInitialCondition"
381  | 
382  | /*++++++++++++++++++++++++++++++++++++++
383  |   Like it says, put together the initial condition.
384  | 
385  |   int FormInitialCondition Returns zero or error.
386  | 
387  |   AppCtx *user The user context structure.
388  | 
389  |   Vec X Vector in which to place the initial condition.
390  |   ++++++++++++++++++++++++++++++++++++++*/
391  | 
392  | int FormInitialCondition (AppCtx *user, Vec X)
393  | {
394  |   int     i,j,k, row, mc, chvar, mx,my,mz, ierr, xs,ys,zs, xm,ym,zm;
395  |   int     gxm,gym,gzm, gxs,gys,gzs;
396  |   PetscScalar  mparam;
397  |   PetscScalar  *x;
398  |   Vec     localX = user->localX;
399  |   PetscRandom rand;
400  | 
401  |   mc = user->mc;
402  |   chvar = user->chvar;
403  |   mx = user->mx; my = user->my; mz = user->mz;
404  |   mparam = user->mparam;
405  | 
406  |   /* Get a pointer to vector data.
407  |        - For default PETSc vectors, VecGetArray() returns a pointer to
408  |          the data array.  Otherwise, the routine is implementation dependent.
409  |        - You MUST call VecRestoreArray() when you no longer need access to
410  |          the array. */
411  |   if (user->random)
412  |     {
413  |       DPRINTF ("Setting initial condition as random numbers\n",0);
414  |       ierr = PetscRandomCreate (user->comm, RANDOM_DEFAULT, &rand);
415  |       CHKERRQ (ierr);
416  |       ierr = PetscRandomSetInterval (rand, 0.35, 0.65); CHKERRQ (ierr);
417  |       ierr = VecSetRandom (rand, X); CHKERRQ (ierr);
418  |       ierr = PetscRandomDestroy (rand); CHKERRQ (ierr);
419  |     }
420  |   else if (user->load_data > -1)
421  |     {
422  |       int usermetacount;
423  |       char basename [1000], **usermetanames, **usermetadata;
424  |       sprintf (basename, "chtsout.time%.3d", user->load_data);
425  |       DPRINTF ("Loading data for timestep %d, basename %s\n", user->load_data,
426  | 	       basename);
427  |       IlluMultiRead (user->da, X, basename, &usermetacount, &usermetanames,
428  | 		     &usermetadata);
429  |       /* TODO: free these */
430  |       for (i=0; i<usermetacount; i++)
431  | 	PetscPrintf (PETSC_COMM_WORLD, "Parameter \"%s\" = \"%s\"\n",
432  | 		     usermetanames [i], usermetadata [i]);
433  |     }
434  |   else
435  |     {
436  |       DPRINTF ("Getting local array\n",0);
437  |       ierr = VecGetArray(localX,&x); CHKERRQ (ierr);
438  | 
439  |       /* Get local grid boundaries (for 2-dimensional DA):
440  | 	 xs, ys   - starting grid indices (no ghost points)
441  | 	 xm, ym   - widths of local grid (no ghost points)
442  | 	 gxs, gys - starting grid indices (including ghost points)
443  | 	 gxm, gym - widths of local grid (including ghost points) */
444  |       DPRINTF ("Getting corners and ghost corners\n",0);
445  |       ierr = DAGetCorners (user->da, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr);
446  |       ierr = DAGetGhostCorners (user->da, &gxs,&gys,&gzs, &gxm,&gym,&gzm);
447  |       CHKERRQ (ierr);
448  | 
449  |       /* Set up 2-D so it works */
450  |       if (!user->threedee) { zs = gzs = 0; zm = 1; mz = 10; }
451  | 
452  |       /* Compute initial condition over the locally owned part of the grid
453  | 	 Initial condition is motionless fluid and equilibrium temperature */
454  |       DPRINTF ("Looping to set initial condition\n",0);
455  |       for (k=zs; k<zs+zm; k++)
456  | 	for (j=ys; j<ys+ym; j++)
457  | 	  for (i=xs; i<xs+xm; i++)
458  | 	    {
459  | 	      row = i - gxs + (j - gys)*gxm + (k-gzs)*gxm*gym;
460  | 	      /* x[C(row)] = (i>=mx*0.43921) ? 1.0 : 0.0; */
461  | 	      x[C(row)]     = ((i<mx/3 || i>2*mx/3) && (j<my/3 || j>2*my/3) &&
462  | 			       (k<mz/3 || k>2*mz/3)) ? 1.0 : 0.0;
463  | 	      /* x[C(row)] = (i<mx/2 && j<my/2 && k<mz/2) ? 1.0 : 0.0; */
464  | 	      /* x[V(row)]     = 0.0;
465  | 		 x[Omega(row)] = 0.0;
466  | 		 x[Temp(row)]  = (mparam>0)*(PetscScalar)(i)/(PetscScalar)(mx-1); */
467  | 	    }
468  | 
469  |       /* Restore vector */
470  |       DPRINTF ("Restoring array to local vector\n",0);
471  |       ierr = VecRestoreArray (localX,&x); CHKERRQ (ierr);
472  | 
473  |       /* Insert values into global vector */
474  |       DPRINTF ("Inserting local vector values into global vector\n",0);
475  |       ierr = DALocalToGlobal (user->da,localX,INSERT_VALUES,X); CHKERRQ (ierr);
476  |     }
477  | 
478  |   return 0;
479  | }
480  | 
481  | 
482  | #undef __FUNCT__
483  | #define __FUNCT__ "InitializeProblem"
484  | 
485  | /*++++++++++++++++++++++++++++++++++++++
486  |   This takes the gory details of initialization out of the way (importing
487  |   parameters into the user context, etc.).
488  | 
489  |   int InitializeProblem Returns zero or error.
490  | 
491  |   AppCtx *user The user context to fill.
492  | 
493  |   Vec *xvec Vector into which to put the initial condition.
494  |   ++++++++++++++++++++++++++++++++++++++*/
495  | 
496  | int InitializeProblem (AppCtx *user, Vec *xvec)
497  | {
498  |   int        Nx,Ny,Nz;  /* number of processors in x-, y- and z- directions */
499  |   int        xs,xm,ys,ym,zs,zm, Nlocal,ierr;
500  |   Vec        xv;
501  |   PetscTruth twodee;
502  | 
503  |   /* Initialize problem parameters */
504  |   DPRINTF ("Initializing user->parameters\n",0);
505  |   user->mx = 20;
506  |   user->my = 20; 
507  |   user->mz = 20; 
508  |   ierr = PetscOptionsGetInt(PETSC_NULL, "-mx", &user->mx, PETSC_NULL);
509  |   CHKERRQ (ierr);
510  |   ierr = PetscOptionsGetInt(PETSC_NULL, "-my", &user->my, PETSC_NULL);
511  |   CHKERRQ (ierr);
512  |   ierr = PetscOptionsGetInt(PETSC_NULL, "-mz", &user->mz, PETSC_NULL);
513  |   CHKERRQ (ierr);
514  |   /* No. of components in the unknown vector and auxiliary vector */
515  |   user->mc = 1;
516  |   /* Problem parameters (kappa, epsilon, gamma, and mparam) */
517  |   user->kappa = 1.0;
518  |   ierr = PetscOptionsGetScalar(PETSC_NULL, "-kappa", &user->kappa, PETSC_NULL);
519  |   CHKERRQ (ierr);
520  |   user->epsilon = 1.0/(user->mx-1);
521  |   ierr = PetscOptionsGetScalar (PETSC_NULL, "-epsilon", &user->epsilon,
522  | 				PETSC_NULL); CHKERRQ (ierr);
523  |   user->gamma = 1.0;
524  |   ierr = PetscOptionsGetScalar(PETSC_NULL, "-gamma", &user->gamma, PETSC_NULL);
525  |   CHKERRQ (ierr);
526  |   user->mparam = 0.0;
527  |   ierr = PetscOptionsGetScalar (PETSC_NULL, "-mparam", &user->mparam,
528  | 				PETSC_NULL); CHKERRQ (ierr);
529  |   ierr = PetscOptionsHasName (PETSC_NULL, "-twodee", &twodee);
530  |   user->threedee = !twodee;
531  |   CHKERRQ (ierr);
532  |   ierr = PetscOptionsHasName (PETSC_NULL, "-printv", &user->print_vecs);
533  |   CHKERRQ (ierr);
534  |   ierr = PetscOptionsHasName (PETSC_NULL, "-printg", &user->print_grid);
535  |   CHKERRQ (ierr);
536  |   ierr = PetscOptionsHasName (PETSC_NULL, "-no_contours", &user->no_contours);
537  |   CHKERRQ (ierr);
538  |   ierr = PetscOptionsHasName (PETSC_NULL, "-random", &user->random);
539  |   CHKERRQ (ierr);
540  |   ierr = PetscOptionsHasName (PETSC_NULL, "-save_data", &user->save_data);
541  |   CHKERRQ (ierr);
542  |   user->load_data = -1;
543  |   ierr = PetscOptionsGetInt (PETSC_NULL, "-load_data", &user->load_data,
544  | 			     PETSC_NULL); CHKERRQ (ierr);
545  | 
546  |   /* Create distributed array (DA) to manage parallel grid and vectors
547  |      for principal unknowns (x) and governing residuals (f)
548  |      Note the stencil width of 2 for this 4th-order equation. */
549  |   DPRINTF ("Creating distributed array\n",0);
550  |   Nx = PETSC_DECIDE;
551  |   Ny = PETSC_DECIDE;
552  |   Nz = PETSC_DECIDE;
553  |   ierr = PetscOptionsGetInt (PETSC_NULL, "-Nx", &Nx, PETSC_NULL);CHKERRQ(ierr);
554  |   ierr = PetscOptionsGetInt (PETSC_NULL, "-Ny", &Ny, PETSC_NULL);CHKERRQ(ierr);
555  |   ierr = PetscOptionsGetInt (PETSC_NULL, "-Nz", &Nz, PETSC_NULL);CHKERRQ(ierr);
556  |   if (user->threedee)
557  |     {
558  |       DPRINTF ("Three Dee!\n",0);
559  |       user->period = DA_XYZPERIODIC;
560  |       ierr = DACreate3d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX,
561  | 			 user->mx, user->my, user->mz, Nx, Ny, Nz, user->mc, 2,
562  | 			 PETSC_NULL, PETSC_NULL, PETSC_NULL, &user->da);
563  |       CHKERRQ (ierr);
564  |     }
565  |   else
566  |     {
567  |       user->period = DA_XYPERIODIC;
568  |       ierr = DACreate2d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX,
569  | 			 user->mx, user->my, Nx, Ny, user->mc, 2,
570  | 			 PETSC_NULL, PETSC_NULL, &user->da); CHKERRQ (ierr);
571  |       user->mz = Nz = 1;
572  |     }
573  | 
574  |   /* Extract global vector from DA */
575  |   DPRINTF ("Extracting global vector from DA...\n",0);
576  |   ierr = DACreateGlobalVector(user->da,&xv); CHKERRQ (ierr);
577  | 
578  |   /* Label PDE components */
579  |   DPRINTF ("Labeling PDE components\n",0);
580  |   ierr = PetscMalloc (user->mc * sizeof(char*), &(user->label));CHKERRQ (ierr);
581  |   user->label[0] = "Concentration (C)";
582  |   /* user->label[1] = "Velocity (V)";
583  |      user->label[2] = "Omega";
584  |      user->label[3] = "Temperature"; */
585  |   ierr = DASetFieldName (user->da, 0, user->label[0]); CHKERRQ(ierr);
586  | 
587  |   user->x_old = 0;
588  | 
589  |   /* Get local vector */
590  |   DPRINTF ("Getting local vector\n",0);
591  |   ierr = DACreateLocalVector (user->da, &user->localX); CHKERRQ (ierr);
592  | 
593  |   /* Print grid info */
594  |   if (user->print_grid)
595  |     {
596  |       DPRINTF ("Printing grid information\n",0);
597  |       ierr = DAView(user->da,PETSC_VIEWER_STDOUT_SELF); CHKERRQ (ierr);
598  |       ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr);
599  |       if (!user->threedee) {
600  | 	zs = 0; zm = 1; }
601  |       ierr = PetscPrintf(PETSC_COMM_WORLD,
602  | 			 "global grid: %d X %d X %d with %d components per"
603  | 			 " node ==> global vector dimension %d\n",
604  | 			 user->mx, user->my, user->mz, user->mc,
605  | 			 user->mc*user->mx*user->my*user->mz);
606  |       CHKERRQ (ierr);
607  |       fflush(stdout);
608  |       ierr = VecGetLocalSize (xv, &Nlocal); CHKERRQ (ierr);
609  |       ierr = PetscSynchronizedPrintf
610  | 	(PETSC_COMM_WORLD,"[%d] local grid %d X %d X %d with %d components"
611  | 	 " per node ==> local vector dimension %d\n",
612  | 	 user->rank, xm, ym, zm, user->mc, Nlocal); CHKERRQ (ierr);
613  |       ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ (ierr);
614  |   }  
615  | 
616  |   /* Compute initial condition */
617  |   DPRINTF ("Forming inital condition\n",0);
618  |   ierr = FormInitialCondition (user, xv); CHKERRQ (ierr);
619  | 
620  |   *xvec = xv;
621  |   return 0;
622  | }