root/trunk/matml/src/castabox/castabox.c

Revision 230, 22.2 kB (checked in by powell, 3 years ago)

Fixes for new PETSc 2.3.0, illuminator 0.10.0.
  • Property svn:keywords set to Author Date Id Revision
Line 
1/* Cast-a-Box software, Copyright 2003 Adam Powell
2 *
3 * This program is free software; you can redistribute it and/or modify it
4 * under the terms of the GNU General Public License as published by the Free
5 * Software Foundation; either version 2 of the License, or (at your option)
6 * any later version.
7 *
8 * This program is distributed in the hope that it will be useful, but WITHOUT
9 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
11 * more details.
12 *
13 * The license can be obtained from the web at http://www.gnu.org/licenses/ .
14 *
15 * If obtained via MIT OpenCourseWare (http://ocw.mit.edu/), you may use it
16 * under the terms of the MIT Open Courseware license instead; this can be
17 * obtained from the web at http://ocw.mit.edu/global/terms-of-use.html .
18 *
19 * $Header$
20 */
21
22#include <stdlib.h> /* For atoi() */
23#include <string.h>
24/* This also #includes petscda.h, petscvec.h and several other dependencies. */
25#include <illuminator.h>
26
27/* Compile with -DDEBUG to print debugging information */
28#undef DPRINTF
29#ifdef DEBUG
30#define DPRINTF(fmt, args...) ierr = PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args); CHKERRQ (ierr)
31#else
32#define DPRINTF(fmt, args...)
33#endif
34
35/* This is printed (along with a TON of other stuff) when you run:
36   castabox -help */
37static char help[] =
38"Cast-a-box version 0.5 by Adam Powell September 25, 2005\n\
39\n\
40This very simple program simulates 3-D heat conduction and solidification of\n\
41a substance in one quarter of a box mold.  It uses the enthalpy method for\n\
42solidification front tracking (which also permits modification to model a\n\
43freezing range) and explicit finite differencing with a uniform grid to\n\
44discretize the equations.  There is no fluid flow in this model.\n\
45\n\
46Internally, it uses PETSc distributed arrays to store the temperature and\n\
47enthalpy fields, so it runs efficiently in parallel (where available).\n\
48For visualization, it uses the Illuminator library.  For info, see URLs:\n\
49  http://www.mcs.anl.gov/petsc/\n\
50  http://lyre.mit.edu/~powell/illuminator.html\n\
51\n\
52You can set parameters at the command line, such as grid and 1/4 box sizes:\n\
53  -grid_x ## [default 5]            -width_x ## [default 0.5]\n\
54  -grid_y ## [default 5]            -width_y ## [default 0.5]\n\
55  -grid_z ## [default 10]           -width_z ## [default 1.0]\n\
56\n\
57Material properties can be set using:\n\
58  -liquid_k ##   [default 1.0]      -solid_k ##     [default 1.0]\n\
59  -liquid_rho ## [default 1.0]      -solid_rho ##   [default 1.0]\n\
60  -liquid_cp ##  [default 1.0]      -solid_cp ##    [default 1.0]\n\
61  -melt_temp ##  [default 1000.0]   -latent_heat ## [default 100.0]\n\
62\n\
63Finally, initial and boundary conditions are set by:\n\
64  -init_temp ##   [default 1200.0]  -init_enthalpy ## [default 300.0]\n\
65  -xmax_tenv ##   [default 0.0]     -xmax_h ##        [default 1.0]\n\
66  -ymax_tenv ##   [default 0.0]     -ymax_h ##        [default 1.0]\n\
67  -bottom_tenv ## [default 0.0]     -bottom_h ##      [default 1.0]\n\
68  -top_tenv ##    [default 0.0]     -top_h ##         [default 0.5]\n\
69  -top_epssig ##  [default 5.67e-11]\n\
70[The top surface has a mixed convection-radiation boundary condition:\n\
71q = h (T-Tenv) + epssig (T^4 - Tenv^4)]\n\
72\n\
73The maximum possible explicit timestep size is calculated automatically from\n\
74the larger of the liquid and solid thermal diffusivities and grid spacing.\n\
75\n";
76
77static char help_commands[] =
78"Cast-a-box commands:\n\
79  h or ?       Print this helpful information\n\
80  g nnn        Go nnn timesteps\n\
81  p [T1 T2 T3] Plot specified temperature contours (none to print settings)\n\
82  o            Turn plotting off\n\
83  q or x       Exit\n\
84";
85
86/* This just makes it easier to pass all of the parameters between functions */
87typedef struct {
88  int grid_x, grid_y, grid_z;                  /* Grid dimensions */
89  PetscScalar width_x, width_y, width_z;       /* (1/4) box dimensions */
90  PetscScalar deltax_m1, deltay_m1, deltaz_m1; /* Inverse square spacings */
91  PetscScalar liquid_k, liquid_rho, liquid_cp; /* Liquid properties */
92  PetscScalar solid_k, solid_rho, solid_cp;    /* Solid properties */
93  PetscScalar melt_temp, latent_heat;          /* Melting behavior */
94  PetscScalar xmax_tenv, ymax_tenv, bottom_tenv, top_tenv; /* Boundary temps */
95  PetscScalar xmax_h, ymax_h, bottom_h, top_h, top_epssig; /* Boundary hs */
96} parameter_values;
97
98/* Tell which field variable in the distributed array is which */
99typedef struct {
100  PetscScalar T, H;
101} Field;
102#define TEMPERATURE 0
103#define ENTHALPY 1
104
105
106/* This is the T(H) function, the inverse of the H(T) function */
107static inline PetscScalar temperature(PetscScalar enthalpy, parameter_values p)
108{
109  /* This uses reference enthalpy=0 for solid at the melting point */
110  if (enthalpy < 0)
111    return p.melt_temp + enthalpy/p.solid_rho/p.solid_cp;
112
113  /* This will have to be modified for a freezing range */
114  if (enthalpy > p.latent_heat)
115    return p.melt_temp + (enthalpy-p.latent_heat)/p.liquid_rho/p.liquid_cp;
116
117  /* This too */
118  return p.melt_temp;
119}
120
121
122/* This gives the liquid conductivity in the liquid, solid in the
123   solid, and a linear interpolation in between */
124static inline PetscScalar conductivity
125(PetscScalar enthalpy, parameter_values p)
126{
127  if (enthalpy < 0)
128    return p.solid_k;
129
130  if (enthalpy > p.latent_heat)
131    return p.liquid_k;
132
133  return p.liquid_k * (enthalpy)/p.latent_heat +
134    p.solid_k * (1.-(enthalpy)/p.latent_heat);
135}
136
137
138/* The do_timestep() function does a single timestep */
139#undef __FUNCT__
140#define __FUNCT__ "do_timestep"
141int do_timestep (DA box_da, parameter_values p, PetscScalar deltat)
142{
143  Vec global, localold;
144  static Vec localnew=(Vec)NULL;
145  Field ***localold_array, ***localnew_array;
146  int xs,ys,zs, xm,ym,zm, gxs,gys,gzs, gxm,gym,gzm, ix,iy,iz, ierr;
147
148  DPRINTF ("Getting the vectors, sending global vector to local.\n",0);
149  ierr = DAGetGlobalVector (box_da, &global); CHKERRQ (ierr);
150  ierr = DAGetLocalVector (box_da, &localold); CHKERRQ (ierr);
151  ierr = DAGlobalToLocalBegin (box_da, global, INSERT_VALUES, localold);
152  CHKERRQ (ierr);
153  ierr = DAGlobalToLocalEnd (box_da, global, INSERT_VALUES, localold);
154  CHKERRQ (ierr);
155  if (localnew == NULL)
156    {
157      DPRINTF ("Creating the static vector to hold the new timestep\n",0);
158      ierr = VecDuplicate (localold, &localnew); CHKERRQ (ierr);
159    }
160
161  DPRINTF ("Getting corner configuration and creating local arrays.\n",0);
162  ierr = DAGetCorners (box_da, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr);
163  ierr = DAGetGhostCorners (box_da, &gxs,&gys,&gzs, &gxm,&gym,&gzm);
164  CHKERRQ (ierr);
165  ierr = DAVecGetArray (box_da, localold, (void **)&localold_array);
166  CHKERRQ (ierr);
167  ierr = DAVecGetArray (box_da, localnew, (void **)&localnew_array);
168  CHKERRQ (ierr);
169
170  DPRINTF ("Doing a timestep of size %g.\n", deltat);
171  for (iz=zs; iz<zs+zm; iz++)
172    for (iy=ys; iy<ys+ym; iy++)
173      for (ix=xs; ix<xs+xm; ix++)
174        {
175          PetscScalar flux_zp1,flux_zm1, flux_yp1,flux_ym1, flux_xp1,flux_xm1;
176
177          /* Calculate flux on x-1 and x+1 faces using average conductivity */
178          if (ix == 0) /* Symmetry boundary condition at x=0 */
179            flux_xm1 = 0.;
180          else
181            flux_xm1 = 0.5 * (conductivity (localold_array[iz][iy][ix-1].H, p)+
182                              conductivity (localold_array[iz][iy][ix].H, p)) *
183              (localold_array[iz][iy][ix-1].T - localold_array[iz][iy][ix].T) *
184              p.deltax_m1;
185
186          if (ix == p.grid_x-1) /* Outer x boundary condition */
187            flux_xp1 = p.xmax_h * (localold_array[iz][iy][ix].T - p.xmax_tenv);
188          else
189            flux_xp1 = 0.5 * (conductivity (localold_array[iz][iy][ix].H, p) +
190                              conductivity (localold_array[iz][iy][ix+1].H,p))*
191              (localold_array[iz][iy][ix].T - localold_array[iz][iy][ix+1].T) *
192              p.deltax_m1;
193
194          /* Calculate flux on y-1 and y+1 faces using average conductivity */
195          if (iy == 0) /* Symmetry boundary condition at y=0 */
196            flux_ym1 = 0.;
197          else
198            flux_ym1 = 0.5 * (conductivity (localold_array[iz][iy-1][ix].H, p)+
199                              conductivity (localold_array[iz][iy][ix].H, p)) *
200              (localold_array[iz][iy-1][ix].T - localold_array[iz][iy][ix].T) *
201              p.deltay_m1;
202
203          if (iy == p.grid_y-1) /* Outer y boundary condition */
204            flux_yp1 = p.ymax_h * (localold_array[iz][iy][ix].T - p.ymax_tenv);
205          else
206            flux_yp1 = 0.5 * (conductivity (localold_array[iz][iy][ix].H, p) +
207                              conductivity (localold_array[iz][iy+1][ix].H,p))*
208              (localold_array[iz][iy][ix].T - localold_array[iz][iy+1][ix].T) *
209              p.deltay_m1;
210
211          /* Calculate flux on z-1 and z+1 faces using average conductivity */
212          if (iz == 0) /* Bottom boundary condition */
213            flux_zm1 = p.bottom_h*(p.bottom_tenv-localold_array[iz][iy][ix].T);
214          else
215            flux_zm1 = 0.5 * (conductivity (localold_array[iz-1][iy][ix].H, p)+
216                              conductivity (localold_array[iz][iy][ix].H, p)) *
217              (localold_array[iz-1][iy][ix].T - localold_array[iz][iy][ix].T) *
218              p.deltaz_m1;
219
220          if (iz == p.grid_z-1) /* Top boundary condition: h, radiation */
221            {
222              flux_zp1 = p.top_h * (localold_array[iz][iy][ix].T - p.top_tenv)+
223                p.top_epssig *
224                (localold_array[iz][iy][ix].T * localold_array[iz][iy][ix].T *
225                 localold_array[iz][iy][ix].T * localold_array[iz][iy][ix].T -
226                 p.top_tenv * p.top_tenv * p.top_tenv * p.top_tenv);
227              /* DPRINTF ("Top: ix=%d, iy=%d, T=%g, flux=%g\n", ix, iy,
228                 localold_array[iz][iy][ix].T, flux_zp1); */
229            }
230          else
231            flux_zp1 = 0.5 * (conductivity (localold_array[iz][iy][ix].H, p) +
232                              conductivity (localold_array[iz+1][iy][ix].H,p))*
233              (localold_array[iz][iy][ix].T - localold_array[iz+1][iy][ix].T) *
234              p.deltaz_m1;
235
236          /* Calculate new enthalpy */
237          localnew_array[iz][iy][ix].H = localold_array[iz][iy][ix].H + deltat*
238            ((flux_xm1-flux_xp1)*p.deltax_m1 +
239             (flux_ym1-flux_yp1)*p.deltay_m1 +
240             (flux_zm1-flux_zp1)*p.deltaz_m1);
241
242          /* Calculate new temperature */
243          localnew_array[iz][iy][ix].T =
244            temperature (localnew_array[iz][iy][ix].H, p);
245        }
246
247  DPRINTF ("Restoring arrays, copying new values into the old vector.\n",0);
248  ierr = DAVecRestoreArray (box_da, localnew, (void **)&localnew_array);
249  CHKERRQ (ierr);
250  ierr = DAVecRestoreArray (box_da, localold, (void **)&localold_array);
251  CHKERRQ (ierr);
252  ierr = VecCopy (localnew, localold);
253
254  DPRINTF ("Sending local back to global, restoring vectors.\n",0);
255  ierr = DALocalToGlobal (box_da, localold, INSERT_VALUES, global);
256  CHKERRQ(ierr);
257  ierr = DARestoreLocalVector (box_da, &localold); CHKERRQ (ierr);
258  ierr = DARestoreGlobalVector (box_da, &global); CHKERRQ (ierr);
259  return 0;
260}
261
262
263/* The timestep_loop() function handles user input and does the timstepping. */
264#undef __FUNCT__
265#define __FUNCT__ "timestep_loop"
266int timestep_loop (DA box_da, parameter_values p)
267{
268  Vec global;
269  PetscScalar deltat, current_time=0., Tmax, Tmin,
270    plot_temps[6] = { 1195., 1150., 1100., 1050., 999.99, 950. },
271    plot_colors[24] = { 1.,0.,0.,.5, 1.,1.,0.,.5, 0.,1.,0.,.5, 0.,1.,1.,.5,
272                    0.,0.,1.,.5, 1.,0.,1.,.5 };
273  int ts=0, plots=6, keepon=1, ierr;
274  char outstring [101];
275  ISurface Surf;
276  IDisplay Disp;
277
278  DPRINTF ("Starting the Geomview display using Illuminator.\n",0);
279  ierr = SurfCreate (&Surf); CHKERRQ (ierr);
280  ierr = GeomviewBegin (PETSC_COMM_WORLD, &Disp); CHKERRQ (ierr);
281
282  DPRINTF ("Calculating maximum stable explicit timestep size.\n",0);
283  {
284    double liquid_alpha, solid_alpha;
285    liquid_alpha = p.liquid_k/p.liquid_rho/p.liquid_cp;
286    solid_alpha = p.solid_k/p.solid_rho/p.solid_cp;
287
288    deltat = 1./2. / PetscMax (liquid_alpha, solid_alpha) /
289      (p.deltax_m1*p.deltax_m1 + p.deltay_m1*p.deltay_m1 +
290       p.deltaz_m1*p.deltaz_m1);
291  }
292
293  ierr = DAGetGlobalVector (box_da, &global); CHKERRQ (ierr);
294  ierr = VecStrideMax (global, TEMPERATURE, PETSC_NULL, &Tmax);
295  CHKERRQ (ierr);
296  ierr = VecStrideMin (global, TEMPERATURE, PETSC_NULL, &Tmin);
297  CHKERRQ (ierr);
298  ierr = DARestoreGlobalVector (box_da, &global); CHKERRQ (ierr);
299
300  snprintf (outstring,100,"Timestep %d (time %g): deltat %g, Tmax %g, Tmin %g",
301            ts, current_time, deltat, Tmax, Tmin);
302#ifdef DEBUG
303  ierr = PetscPrintf (PETSC_COMM_WORLD, "%s\n", outstring); CHKERRQ (ierr);
304#else
305  strncat (outstring, "            ", 79-strlen(outstring));
306  ierr = PetscPrintf (PETSC_COMM_WORLD, "\r%s\n", outstring); CHKERRQ (ierr);
307#endif
308
309  while (keepon)
310    {
311      char instring[100];
312
313      ierr = PetscPrintf (PETSC_COMM_WORLD,
314                          "Cast-a-box command (h for help): ");
315      CHKERRQ (ierr);
316      ierr = PetscSynchronizedFGets (PETSC_COMM_WORLD, stdin, 99, instring);
317      CHKERRQ (ierr);
318
319      switch (*instring)
320        {
321        case 'h':
322        case 'H':
323        case '?':
324          ierr = PetscPrintf (PETSC_COMM_WORLD, help_commands); CHKERRQ (ierr);
325          break;
326
327        case 'q':
328        case 'Q':
329        case 'x':
330        case 'X':
331          keepon=0;
332          break;
333
334        case 'p':
335        case 'P':
336          {
337            int count=0, newplots=0;
338
339            while (newplots<6 && instring[count] != 0)
340              {
341                while ((instring[count] < '0' || instring[count] > '9') &&
342                       instring[count] != '-' && instring[count] != '.' &&
343                       instring[count] != '\0')
344                  count++;
345                if (instring[count])
346                  {
347                    sscanf (instring+count, "%lf\n", plot_temps+newplots);
348                    newplots++;
349                    while ((instring[count] >= '0' && instring[count] <= '9')||
350                           instring[count] == '-' || instring[count] == '.')
351                      count++;
352                  }
353              }
354
355            if (newplots)
356              {
357                PetscScalar minmax[6] =
358                  { 0., p.width_x, 0., p.width_y, 0., p.width_z };
359
360                ierr = DATriangulate (Surf, box_da,global, TEMPERATURE, minmax,
361                                      plots=newplots, plot_temps, plot_colors,
362                                      PETSC_FALSE, PETSC_FALSE, PETSC_FALSE);
363                CHKERRQ (ierr);
364                /* Work-around for illuminator non-periodic bug */
365                minmax[1] = p.width_x * (p.grid_x-1) / p.grid_x;
366                minmax[3] = p.width_y * (p.grid_y-1) / p.grid_y;
367                minmax[5] = p.width_z * (p.grid_z-1) / p.grid_z;
368                ierr = GeomviewDisplayTriangulation
369                  (PETSC_COMM_WORLD, Surf, Disp, minmax, "Temperature",
370                   PETSC_FALSE); CHKERRQ(ierr);
371              }
372            else
373              {
374                ierr = PetscPrintf (PETSC_COMM_WORLD,
375                                    "Current plot temperatures:");
376                CHKERRQ (ierr);
377                for (count=0; count<plots; count++)
378                  {
379                    ierr = PetscPrintf (PETSC_COMM_WORLD, " %g",
380                                        plot_temps[count]); CHKERRQ (ierr);
381                  }
382                ierr = PetscPrintf (PETSC_COMM_WORLD, "\n"); CHKERRQ (ierr);
383              }
384
385            break;
386          }
387
388        case 'g':
389        case 'G':
390        case '\0':
391          {
392            int count=0, i;
393            while ((instring[count] < '0' || instring[count] > '9') &&
394                   instring[count] != '\0')
395              count++;
396
397            DPRINTF ("count = %d\n", count);
398
399            if (instring[count] == '\0')
400              count = 1;
401            else
402              count = atoi (instring+count);
403
404            DPRINTF ("Doing %d timesteps\n", count);
405
406            for (i=0; i<count; i++)
407              {
408                ierr = do_timestep (box_da, p, deltat); CHKERRQ (ierr);
409
410                ierr = DAGetGlobalVector (box_da, &global); CHKERRQ (ierr);
411                ierr = VecStrideMax (global, TEMPERATURE, PETSC_NULL, &Tmax);
412                CHKERRQ (ierr);
413                ierr = VecStrideMin (global, TEMPERATURE, PETSC_NULL, &Tmin);
414                CHKERRQ (ierr);
415                ierr = DARestoreGlobalVector (box_da, &global); CHKERRQ (ierr);
416
417                current_time += deltat;
418                ts++;
419
420                snprintf (outstring,100,"Timestep %d (time %g): deltat %g, Tmax %g, Tmin %g",
421                          ts, current_time, deltat, Tmax, Tmin);
422#ifdef DEBUG
423                ierr = PetscPrintf (PETSC_COMM_WORLD, "%s\n", outstring); CHKERRQ (ierr);
424#else
425                strncat (outstring, "            ", 79-strlen(outstring));
426                ierr = PetscPrintf (PETSC_COMM_WORLD, "\r%s", outstring); CHKERRQ (ierr);
427#endif
428
429                if (plots)
430                  {
431                    PetscScalar minmax[6] =
432                      { 0., p.width_x, 0., p.width_y, 0., p.width_z };
433
434                    ierr = DATriangulate (Surf, box_da, global, TEMPERATURE,
435                                          minmax, plots,plot_temps,plot_colors,
436                                          PETSC_FALSE,PETSC_FALSE,PETSC_FALSE);
437                    CHKERRQ (ierr);
438                    /* Work-around for illuminator non-periodic bug */
439                    minmax[1] = p.width_x * (p.grid_x-1) / p.grid_x;
440                    minmax[3] = p.width_y * (p.grid_y-1) / p.grid_y;
441                    minmax[5] = p.width_z * (p.grid_z-1) / p.grid_z;
442                    ierr = GeomviewDisplayTriangulation
443                      (PETSC_COMM_WORLD, Surf, Disp, minmax, "Temperature",
444                       PETSC_FALSE); CHKERRQ (ierr);
445                    ierr = SurfClear (Surf); CHKERRQ (ierr);
446                  }
447              }
448#ifndef DEBUG
449            if (count)
450              {
451                ierr = PetscPrintf (PETSC_COMM_WORLD, "\n"); CHKERRQ (ierr);
452              }
453#endif
454            break;
455          }
456
457        default:
458          ierr = PetscPrintf (PETSC_COMM_WORLD, "Unrecognized: %s\n",instring);
459          CHKERRQ (ierr);
460        }
461    }
462
463  ierr = GeomviewEnd (PETSC_COMM_WORLD, Disp); CHKERRQ (ierr);
464  return 0;
465}
466
467
468/* The main() function reads in the parameters, calculates the timestep size,
469 * creates the distributed array, and starts Geomview for visualization, then
470 * calls the timestep loop to do the simulation. */
471#undef __FUNCT__
472#define __FUNCT__ "main"
473int main (int argc, char *argv[])
474{
475  DA box_da;
476  parameter_values p;
477  int ierr;
478
479  /* Start with PETSc initialization */
480  ierr = PetscInitialize (&argc, &argv, NULL, help); CHKERRQ (ierr);
481
482  DPRINTF ("Getting the parameters from the command line into p.\n",0);
483  p.grid_x = p.grid_y = 5;
484  p.grid_z = 10;
485  ierr = PetscOptionsGetInt (PETSC_NULL, "-grid_x", &p.grid_x, PETSC_NULL);
486  CHKERRQ (ierr);
487  ierr = PetscOptionsGetInt (PETSC_NULL, "-grid_y", &p.grid_y, PETSC_NULL);
488  CHKERRQ (ierr);
489  ierr = PetscOptionsGetInt (PETSC_NULL, "-grid_z", &p.grid_z, PETSC_NULL);
490  CHKERRQ (ierr);
491
492  p.width_x = p.width_y = 0.5;
493  p.width_z = 1.;
494  ierr = PetscOptionsGetScalar (PETSC_NULL, "-width_x", &p.width_x,
495                                PETSC_NULL); CHKERRQ (ierr);
496  ierr = PetscOptionsGetScalar (PETSC_NULL, "-width_y", &p.width_y,
497                                PETSC_NULL); CHKERRQ (ierr);
498  ierr = PetscOptionsGetScalar (PETSC_NULL, "-width_z", &p.width_z,
499                                PETSC_NULL); CHKERRQ (ierr);
500  p.deltax_m1 = (PetscScalar) p.grid_x / p.width_x;
501  p.deltay_m1 = (PetscScalar) p.grid_y / p.width_y;
502  p.deltaz_m1 = (PetscScalar) p.grid_z / p.width_z;
503
504  p.liquid_k = p.liquid_rho = p.liquid_cp = 1.;
505  ierr = PetscOptionsGetScalar (PETSC_NULL, "-liquid_k", &p.liquid_k,
506                                PETSC_NULL); CHKERRQ (ierr);
507  ierr = PetscOptionsGetScalar (PETSC_NULL, "-liquid_rho", &p.liquid_rho,
508                                PETSC_NULL); CHKERRQ (ierr);
509  ierr = PetscOptionsGetScalar (PETSC_NULL, "-liquid_cp", &p.liquid_cp,
510                                PETSC_NULL); CHKERRQ (ierr);
511
512  p.solid_k = p.solid_rho = p.solid_cp = 1.;
513  ierr = PetscOptionsGetScalar (PETSC_NULL, "-solid_k", &p.solid_k,
514                                PETSC_NULL); CHKERRQ (ierr);
515  ierr = PetscOptionsGetScalar (PETSC_NULL, "-solid_rho", &p.solid_rho,
516                                PETSC_NULL); CHKERRQ (ierr);
517  ierr = PetscOptionsGetScalar (PETSC_NULL, "-solid_cp", &p.solid_cp,
518                                PETSC_NULL); CHKERRQ (ierr);
519
520  p.melt_temp = 1000.;
521  p.latent_heat = 100.;
522  ierr = PetscOptionsGetScalar (PETSC_NULL, "-melt_temp", &p.melt_temp,
523                                PETSC_NULL); CHKERRQ (ierr);
524  ierr = PetscOptionsGetScalar (PETSC_NULL, "-latent_heat", &p.latent_heat,
525                                PETSC_NULL); CHKERRQ (ierr);
526
527  p.xmax_tenv = p.ymax_tenv = p.bottom_tenv = p.top_tenv = 0.;
528  ierr = PetscOptionsGetScalar (PETSC_NULL, "-xmax_tenv", &p.xmax_tenv,
529                                PETSC_NULL); CHKERRQ (ierr);
530  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ymax_tenv", &p.ymax_tenv,
531                                PETSC_NULL); CHKERRQ (ierr);
532  ierr = PetscOptionsGetScalar (PETSC_NULL, "-bottom_tenv", &p.bottom_tenv,
533                                PETSC_NULL); CHKERRQ (ierr);
534  ierr = PetscOptionsGetScalar (PETSC_NULL, "-top_tenv", &p.top_tenv,
535                                PETSC_NULL); CHKERRQ (ierr);
536
537  p.xmax_h = p.ymax_h = p.bottom_h = 1.;
538  p.top_h = 0.5;
539  p.top_epssig = 5.67e-11;
540  ierr = PetscOptionsGetScalar (PETSC_NULL, "-xmax_h", &p.xmax_h, PETSC_NULL);
541  CHKERRQ (ierr);
542  ierr = PetscOptionsGetScalar (PETSC_NULL, "-ymax_h", &p.ymax_h, PETSC_NULL);
543  CHKERRQ (ierr);
544  ierr = PetscOptionsGetScalar (PETSC_NULL, "-bottom_h", &p.bottom_h,
545                                PETSC_NULL); CHKERRQ (ierr);
546  ierr = PetscOptionsGetScalar (PETSC_NULL, "-top_h", &p.top_h,
547                                PETSC_NULL); CHKERRQ (ierr);
548  ierr = PetscOptionsGetScalar (PETSC_NULL, "-top_epssig", &p.top_epssig,
549                                PETSC_NULL); CHKERRQ (ierr);
550
551  DPRINTF ("Simulation parameters:\n",0);
552  DPRINTF ("Grid dimensions: %dx%dx%d, Physical dimensions: %gx%gx%g\n",
553           p.grid_x, p.grid_y, p.grid_z, p.width_x, p.width_y, p.width_z);
554  DPRINTF ("Liquid properties: k=%g, rho=%g, cp=%g, melt temp=%g\n",
555           p.liquid_k, p.liquid_rho, p.liquid_cp, p.melt_temp);
556  DPRINTF ("Solid properties:  k=%g, rho=%g, cp=%g, latent heat=%g\n",
557           p.solid_k, p.solid_rho, p.solid_cp, p.latent_heat);
558  DPRINTF ("Env. temperatures: xmax=%g, ymax=%g, top=%g, bottom=%g\n",
559           p.xmax_tenv, p.ymax_tenv, p.top_tenv, p.bottom_tenv);
560  DPRINTF ("h values: xmax=%g, ymax=%g, top=%g, bottom=%g, epssig=%g\n",
561           p.xmax_h, p.ymax_h, p.bottom_h, p.top_h, p.top_epssig);
562
563  DPRINTF ("Creating the distributed array with box stencil of width 1,\n"
564           "  2 field variables (temperature and enthalpy).\n",0);
565  ierr = DACreate3d (PETSC_COMM_WORLD, DA_NONPERIODIC, DA_STENCIL_BOX,
566                     p.grid_x, p.grid_y, p.grid_z,
567                     PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
568                     2, 1, PETSC_NULL, PETSC_NULL, PETSC_NULL, &box_da);
569  CHKERRQ (ierr);
570  ierr = DASetFieldName (box_da, TEMPERATURE, "Temperature"); CHKERRQ (ierr);
571  ierr = DASetFieldName (box_da, ENTHALPY, "Enthalpy"); CHKERRQ (ierr);
572
573  DPRINTF ("Setting the initial conditions.\n",0);
574  {
575    Vec global;
576    Field ***global_array, initial;
577    int xs,ys,zs, xm,ym,zm, ix,iy,iz;
578
579    /* Default: 200 degree initial superheat */
580    initial.H = p.latent_heat + 200.*p.liquid_rho*p.liquid_cp;
581    ierr = PetscOptionsGetScalar (PETSC_NULL, "-init_enthalpy", &initial.H,
582                                  PETSC_NULL); CHKERRQ (ierr);
583    /* Provided initial temperature overrides enthalpy */
584    initial.T = temperature (initial.H, p);
585    ierr = PetscOptionsGetScalar (PETSC_NULL, "-init_temp", &initial.T,
586                                  PETSC_NULL); CHKERRQ (ierr);
587    if (initial.T != temperature (initial.H, p))
588      initial.H = ((initial.T<p.melt_temp) ?
589                   p.solid_rho*p.solid_cp*(initial.T-p.melt_temp) :
590                   p.liquid_rho*p.liquid_cp*(initial.T-p.melt_temp) +
591                   p.latent_heat);
592
593    ierr = DAGetGlobalVector (box_da, &global); CHKERRQ (ierr);
594    ierr = DAVecGetArray (box_da, global, (void **)&global_array); CHKERRQ (ierr);
595    ierr = DAGetCorners(box_da, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr);
596
597    for (iz=zs; iz<zs+zm; iz++)
598      for (iy=ys; iy<ys+ym; iy++)
599        for (ix=xs; ix<xs+xm; ix++)
600          global_array[iz][iy][ix] = initial;
601
602    ierr = DAVecRestoreArray (box_da, global, (void **)&global_array); CHKERRQ (ierr);
603    ierr = DARestoreGlobalVector (box_da, &global); CHKERRQ (ierr);
604  }
605
606  DPRINTF ("Starting the timestep loop...\n",0);
607  ierr = timestep_loop (box_da, p);
608
609  DPRINTF ("Shutting down.\n",0);
610  ierr = DADestroy (box_da); CHKERRQ (ierr);
611  ierr = PetscFinalize (); CHKERRQ (ierr);
612  return 0;
613}
Note: See TracBrowser for help on using the browser.