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