/* Cast-a-Box software, Copyright 2003 Adam Powell * * This program is free software; you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by the Free * Software Foundation; either version 2 of the License, or (at your option) * any later version. * * This program is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for * more details. * * The license can be obtained from the web at http://www.gnu.org/licenses/ . * * If obtained via MIT OpenCourseWare (http://ocw.mit.edu/), you may use it * under the terms of the MIT Open Courseware license instead; this can be * obtained from the web at http://ocw.mit.edu/global/terms-of-use.html . * * $Header$ */ #include /* For atoi() */ #include /* This also #includes petscda.h, petscvec.h and several other dependencies. */ #include /* Compile with -DDEBUG to print debugging information */ #undef DPRINTF #ifdef DEBUG #define DPRINTF(fmt, args...) ierr = PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args); CHKERRQ (ierr) #else #define DPRINTF(fmt, args...) #endif /* This is printed (along with a TON of other stuff) when you run: castabox -help */ static char help[] = "Cast-a-box version 0.5 by Adam Powell September 25, 2005\n\ \n\ This very simple program simulates 3-D heat conduction and solidification of\n\ a substance in one quarter of a box mold. It uses the enthalpy method for\n\ solidification front tracking (which also permits modification to model a\n\ freezing range) and explicit finite differencing with a uniform grid to\n\ discretize the equations. There is no fluid flow in this model.\n\ \n\ Internally, it uses PETSc distributed arrays to store the temperature and\n\ enthalpy fields, so it runs efficiently in parallel (where available).\n\ For visualization, it uses the Illuminator library. For info, see URLs:\n\ http://www.mcs.anl.gov/petsc/\n\ http://lyre.mit.edu/~powell/illuminator.html\n\ \n\ You can set parameters at the command line, such as grid and 1/4 box sizes:\n\ -grid_x ## [default 5] -width_x ## [default 0.5]\n\ -grid_y ## [default 5] -width_y ## [default 0.5]\n\ -grid_z ## [default 10] -width_z ## [default 1.0]\n\ \n\ Material properties can be set using:\n\ -liquid_k ## [default 1.0] -solid_k ## [default 1.0]\n\ -liquid_rho ## [default 1.0] -solid_rho ## [default 1.0]\n\ -liquid_cp ## [default 1.0] -solid_cp ## [default 1.0]\n\ -melt_temp ## [default 1000.0] -latent_heat ## [default 100.0]\n\ \n\ Finally, initial and boundary conditions are set by:\n\ -init_temp ## [default 1200.0] -init_enthalpy ## [default 300.0]\n\ -xmax_tenv ## [default 0.0] -xmax_h ## [default 1.0]\n\ -ymax_tenv ## [default 0.0] -ymax_h ## [default 1.0]\n\ -bottom_tenv ## [default 0.0] -bottom_h ## [default 1.0]\n\ -top_tenv ## [default 0.0] -top_h ## [default 0.5]\n\ -top_epssig ## [default 5.67e-11]\n\ [The top surface has a mixed convection-radiation boundary condition:\n\ q = h (T-Tenv) + epssig (T^4 - Tenv^4)]\n\ \n\ The maximum possible explicit timestep size is calculated automatically from\n\ the larger of the liquid and solid thermal diffusivities and grid spacing.\n\ \n"; static char help_commands[] = "Cast-a-box commands:\n\ h or ? Print this helpful information\n\ g nnn Go nnn timesteps\n\ p [T1 T2 T3] Plot specified temperature contours (none to print settings)\n\ o Turn plotting off\n\ q or x Exit\n\ "; /* This just makes it easier to pass all of the parameters between functions */ typedef struct { int grid_x, grid_y, grid_z; /* Grid dimensions */ PetscScalar width_x, width_y, width_z; /* (1/4) box dimensions */ PetscScalar deltax_m1, deltay_m1, deltaz_m1; /* Inverse square spacings */ PetscScalar liquid_k, liquid_rho, liquid_cp; /* Liquid properties */ PetscScalar solid_k, solid_rho, solid_cp; /* Solid properties */ PetscScalar melt_temp, latent_heat; /* Melting behavior */ PetscScalar xmax_tenv, ymax_tenv, bottom_tenv, top_tenv; /* Boundary temps */ PetscScalar xmax_h, ymax_h, bottom_h, top_h, top_epssig; /* Boundary hs */ } parameter_values; /* Tell which field variable in the distributed array is which */ typedef struct { PetscScalar T, H; } Field; #define TEMPERATURE 0 #define ENTHALPY 1 /* This is the T(H) function, the inverse of the H(T) function */ static inline PetscScalar temperature(PetscScalar enthalpy, parameter_values p) { /* This uses reference enthalpy=0 for solid at the melting point */ if (enthalpy < 0) return p.melt_temp + enthalpy/p.solid_rho/p.solid_cp; /* This will have to be modified for a freezing range */ if (enthalpy > p.latent_heat) return p.melt_temp + (enthalpy-p.latent_heat)/p.liquid_rho/p.liquid_cp; /* This too */ return p.melt_temp; } /* This gives the liquid conductivity in the liquid, solid in the solid, and a linear interpolation in between */ static inline PetscScalar conductivity (PetscScalar enthalpy, parameter_values p) { if (enthalpy < 0) return p.solid_k; if (enthalpy > p.latent_heat) return p.liquid_k; return p.liquid_k * (enthalpy)/p.latent_heat + p.solid_k * (1.-(enthalpy)/p.latent_heat); } /* The do_timestep() function does a single timestep */ #undef __FUNCT__ #define __FUNCT__ "do_timestep" int do_timestep (DA box_da, parameter_values p, PetscScalar deltat) { Vec global, localold; static Vec localnew=(Vec)NULL; Field ***localold_array, ***localnew_array; int xs,ys,zs, xm,ym,zm, gxs,gys,gzs, gxm,gym,gzm, ix,iy,iz, ierr; DPRINTF ("Getting the vectors, sending global vector to local.\n",0); ierr = DAGetGlobalVector (box_da, &global); CHKERRQ (ierr); ierr = DAGetLocalVector (box_da, &localold); CHKERRQ (ierr); ierr = DAGlobalToLocalBegin (box_da, global, INSERT_VALUES, localold); CHKERRQ (ierr); ierr = DAGlobalToLocalEnd (box_da, global, INSERT_VALUES, localold); CHKERRQ (ierr); if (localnew == NULL) { DPRINTF ("Creating the static vector to hold the new timestep\n",0); ierr = VecDuplicate (localold, &localnew); CHKERRQ (ierr); } DPRINTF ("Getting corner configuration and creating local arrays.\n",0); ierr = DAGetCorners (box_da, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr); ierr = DAGetGhostCorners (box_da, &gxs,&gys,&gzs, &gxm,&gym,&gzm); CHKERRQ (ierr); ierr = DAVecGetArray (box_da, localold, (void **)&localold_array); CHKERRQ (ierr); ierr = DAVecGetArray (box_da, localnew, (void **)&localnew_array); CHKERRQ (ierr); DPRINTF ("Doing a timestep of size %g.\n", deltat); for (iz=zs; iz '9') && instring[count] != '-' && instring[count] != '.' && instring[count] != '\0') count++; if (instring[count]) { sscanf (instring+count, "%lf\n", plot_temps+newplots); newplots++; while ((instring[count] >= '0' && instring[count] <= '9')|| instring[count] == '-' || instring[count] == '.') count++; } } if (newplots) { PetscScalar minmax[6] = { 0., p.width_x, 0., p.width_y, 0., p.width_z }; ierr = DATriangulate (Surf, box_da,global, TEMPERATURE, minmax, plots=newplots, plot_temps, plot_colors, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE); CHKERRQ (ierr); /* Work-around for illuminator non-periodic bug */ minmax[1] = p.width_x * (p.grid_x-1) / p.grid_x; minmax[3] = p.width_y * (p.grid_y-1) / p.grid_y; minmax[5] = p.width_z * (p.grid_z-1) / p.grid_z; ierr = GeomviewDisplayTriangulation (PETSC_COMM_WORLD, Surf, Disp, minmax, "Temperature", PETSC_FALSE); CHKERRQ(ierr); } else { ierr = PetscPrintf (PETSC_COMM_WORLD, "Current plot temperatures:"); CHKERRQ (ierr); for (count=0; count '9') && instring[count] != '\0') count++; DPRINTF ("count = %d\n", count); if (instring[count] == '\0') count = 1; else count = atoi (instring+count); DPRINTF ("Doing %d timesteps\n", count); for (i=0; i