1 | /***************************************
2 | $Header$
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 = TSMonitorSet (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, TSCRANK_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 | }