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