Merged enzo/enzo-dev into week-of-code

--HG--
branch : week-of-code

.hgtags

@@ -2,3 +2,4 @@ bbf0a2ffbd22c4fbecf946c9c96e6c4fac5cbdae woc_pre_fld_merge
 48b4e9d9d6b90f703e48e621b488136be2a0e9cf woc_fld_merge
 b86d8ba026d6a0ec30f15d8134add1e55fae2958 Wise10_GalaxyBirth
 2d90aa38e06f00a531db45a43225cde1faf093f2 enzo-2.2
+a444b00827d3d4351dceea8197f1e446b9ada0da enzogold0001

doc/manual/source/parameters/problemtypes.rst

@@ -684,6 +684,15 @@ Collapse Test (27)
     ``CollapseTestSphereAng2`` are set, the rotational axis linearly
     changes with radius between ``CollapseTestSphereAng1`` and
     ``CollapseTestSphereAng2``.  Units in radians. Default: 0.
+``CollapseTestSphereConstantPressure`` (external)
+    Constant pressure inside the sphere that is equal to the pressure
+    at the outer radius.  Default: 0
+``CollapseTestSphereSmoothSurface`` (external)
+    The density interface between the ambient and sphere medium is
+    smoothed with a hyperbolic tangent.  Default: 0
+``CollapseTestSmoothRadius`` (external)
+    The outer radius of the smoothed interface.  This parameter is in
+    units of the sphere radius.  Default: 1.2
 ``CollapseTestSphereInitialLevel`` (external)
     Failed experiment to try to force refinement to a specified level.
     Not working. Default: 0.

run/Hydro/Hydro-1D/ShockInABox/ShockInABox.enzo

@@ -0,0 +1,71 @@
+# 
+# ShockInABox Problem.
+# Shock Mach Number: 3.000000
+# PreShock Temperature: 2.727273e+06
+# PostShock Temperature: 1.000000e+07
+# PreShock Density: 1.000000e+00
+# 
+
+#
+# AMR PROBLEM DEFINITION FILE: 1D Shock Propogation test
+#
+#  define problem
+#
+ProblemType            = 5        // Shock In A Box
+TopGridRank            = 1
+TopGridDimensions      = 32
+
+ShockInABoxBoundary      = 0.5
+
+LeftFaceBoundaryCondition = 2 0 0  // set left faces to inflow 
+RightFaceBoundaryCondition = 1 0 0  // set right faces to outflow
+
+#
+#  set I/O and stop/start parameters
+#
+StopTime               = 1.0e0
+dtDataDump             = 1.0e-1
+
+#
+#  set Hydro parameters
+#
+HydroMethod            = 0
+Gamma                  = 1.66667
+CourantSafetyNumber    = 0.8
+PPMDiffusionParameter  = 1       // diffusion on 
+PPMFlatteningParameter = 1       // flattening on
+PPMSteepeningParameter = 1       // steepening on
+
+#
+#  set grid refinement parameters
+#
+StaticHierarchy           = 0    // static hierarchy
+MaximumRefinementLevel    = 2
+RefineBy                  = 2    // refinement factor
+CellFlaggingMethod        = 1
+MinimumEfficiency         = 0.5
+
+#
+#  set some misc global parameters
+#
+OutputTemperature = 1
+
+# 
+#  Turn on Shock Finding
+#
+ShockTemperatureFloor=1.0e0
+StorePreShockFields = 1
+FindShocksOnlyOnOutput = 1
+ShockMethod = 1
+
+
+DensityUnits = 1.67453400e-24
+LengthUnits = 3.08567758e+24
+TimeUnits = 3.15576000e+16
+ShockInABoxLeftDensity = 1.00000000e+00
+ShockInABoxLeftVelocity = 6.90120768e-01
+ShockInABoxLeftPressure = 3.91989033e-02
+ShockInABoxRightDensity = 3.00000000e+00
+ShockInABoxRightVelocity = 1.78920199e-01
+ShockInABoxRightPressure = 4.31187936e-01
+

run/Hydro/Hydro-1D/ShockInABox/ShockInABox.enzotest

@@ -0,0 +1,11 @@
+name = 'ShockInABox'
+answer_testing_script = 'ShockInABox__test_shockinabox.py' 
+nprocs = 1
+runtime = 'short'
+hydro = True
+gravity = False
+dimensionality = 1
+max_time_minutes = 1
+fullsuite = True
+pushsuite = True
+quicksuite = True

run/Hydro/Hydro-1D/ShockInABox/ShockInABox__test_shockinabox.py

@@ -0,0 +1,16 @@
+from yt.mods import *
+from yt.funcs import *
+from yt.testing import *
+from yt.frontends.enzo.answer_testing_support import \
+    requires_outputlog
+from yt.utilities.answer_testing.framework import AssertWrapper
+import os
+
+# Verifies that OutputLog exists
+@requires_outputlog(os.path.dirname(__file__), "ShockInABox.enzo")
+def test_shockinabox():
+    pf = load('DD0010/data0010')
+    mach = pf.h.find_max('Mach')[0]
+    myname = 'ShockInABox_Mach'
+    yield AssertWrapper(myname, assert_allclose, mach, 3.0, 1.0e-2)
+

run/Hydro/Hydro-1D/ShockInABox/input_shock.enzo

@@ -0,0 +1,51 @@
+#
+# AMR PROBLEM DEFINITION FILE: 1D Shock Propogation test
+#
+#  define problem
+#
+ProblemType            = 5        // Shock In A Box
+TopGridRank            = 1
+TopGridDimensions      = 32
+
+ShockInABoxBoundary      = 0.5
+
+LeftFaceBoundaryCondition = 2 0 0  // set left faces to inflow 
+RightFaceBoundaryCondition = 1 0 0  // set right faces to outflow
+
+#
+#  set I/O and stop/start parameters
+#
+StopTime               = 1.0e0
+dtDataDump             = 1.0e-1
+
+#
+#  set Hydro parameters
+#
+HydroMethod            = 0
+CourantSafetyNumber    = 0.8
+PPMDiffusionParameter  = 1       // diffusion on 
+PPMFlatteningParameter = 1       // flattening on
+PPMSteepeningParameter = 1       // steepening on
+
+#
+#  set grid refinement parameters
+#
+StaticHierarchy           = 0    // static hierarchy
+MaximumRefinementLevel    = 2
+RefineBy                  = 2    // refinement factor
+CellFlaggingMethod        = 1
+MinimumEfficiency         = 0.5
+
+#
+#  set some misc global parameters
+#
+OutputTemperature = 1
+
+# 
+#  Turn on Shock Finding
+#
+ShockTemperatureFloor=1.0e0
+StorePreShockFields = 1
+FindShocksOnlyOnOutput = 1
+ShockMethod = 1
+

run/Hydro/Hydro-1D/ShockInABox/make_plots.py

@@ -0,0 +1,47 @@
+from yt.mods import *
+import pylab
+
+### define simulation output directory and filename base
+output_dir_base = 'DD'
+datafile_base = 'data'
+
+### load data
+ts = TimeSeriesData.from_filenames("*/*.hierarchy")
+for pf in ts:
+    pylab.clf()
+    print pf.current_time
+
+    ### extract an ortho_ray (1D solution vector)
+    ray = pf.h.ortho_ray(0, [0.5, 0.5])
+
+    pylab.figure(1, figsize=(10,8))
+
+    # Density Plot
+    pylab.subplot(2,2,1)
+    pylab.semilogy(ray['x'],ray['Density'], 'k')
+    pylab.xlabel('Position')
+    pylab.ylabel('Density')
+
+    # Temperature Plot
+    pylab.subplot(2,2,2)
+    pylab.semilogy(ray['x'],ray['Temperature'], 'b')
+    pylab.xlabel('Position')
+    pylab.ylabel('Temperature')
+
+    # Mach Plot
+    pylab.subplot(2,2,3)
+    pylab.plot(ray['x'],ray['Mach'], 'k')
+    pylab.xlabel('x')
+    pylab.ylabel('Mach')
+
+    # Mach Plot
+    pylab.subplot(2,2,4)
+    pylab.plot(ray['x'],ray['VelocityMagnitude'], 'k')
+    pylab.xlabel('x')
+    pylab.ylabel('|v|')
+
+    ### Save plot
+    pylab.savefig('%s_thermal.png' % pf)
+
+pylab.clf()
+

run/Hydro/Hydro-1D/ShockInABox/make_shock.py

@@ -0,0 +1,100 @@
+#!/usr/bin/python
+import sys
+import os
+
+kboltz = 1.3806504e-16  # erg K^-1
+sec_per_Gyr  = 31.5576e15
+mpc_per_cm    = 3.24077929e-25
+cm_per_mpc    = 1.0 / mpc_per_cm
+mh = 1.674534e-24  # g
+mu = 0.6
+
+def setup_shockbox(length_units, time_units, d1, T1, m, gamma=None):
+    if gamma is None: gamma = 5./3.
+
+    dens_units = mh
+    temp_units = dens_units*(length_units/time_units)**2/kboltz
+    p1 = (T1/temp_units)*d1/mu
+    d2 = d1*((gamma+1.)*m*m)/((gamma-1.)*m*m + 2.);
+    p2 = p1*(2.0*gamma*m*m - (gamma-1.))/(gamma+1.);
+    c1 = (gamma*p1/(d1))**0.5;
+    v1 = 0.0
+    v2 = m*c1*(1.-d1/d2);
+    shockspeed = 1.0 * c1 * m;
+
+    vel_units = length_units/time_units
+
+    lines = []
+    lines.append('\n')
+    lines.append('DensityUnits = %0.8e\n' % dens_units)
+    lines.append('LengthUnits = %0.8e\n' % length_units)
+    lines.append('TimeUnits = %0.8e\n' % time_units)
+    lines.append('Gamma = %f\n' % gamma)
+    lines.append('ShockInABoxLeftDensity = %0.8e\n' % (d1)) 
+    lines.append('ShockInABoxLeftVelocity = %0.8e\n' % ((shockspeed - v1))) 
+    lines.append('ShockInABoxLeftPressure = %0.8e\n' % (p1) ) 
+
+    lines.append('ShockInABoxRightDensity = %0.8e\n' % (d2)) 
+    lines.append('ShockInABoxRightVelocity = %0.8e\n' % ((shockspeed - v2))) 
+    lines.append('ShockInABoxRightPressure = %0.8e\n' % (p2) ) 
+    lines.append('\n')
+
+    return lines
+
+def get_header(d1, T1, T2, M):
+    lines = []
+    lines.append('# \n')
+    lines.append('# Custom ShockInABox Problem.\n')
+    lines.append('# Shock Mach Number: %f\n' % M)
+    lines.append('# PreShock Temperature: %e\n' % T1)
+    lines.append('# PostShock Temperature: %e\n' % T2)
+    lines.append('# PreShock Density: %e\n' % d1)
+    lines.append('# \n\n')
+    return lines
+
+def write_lines(outfile, lines):
+    f = file(outfile,'w')
+    f.writelines(lines)
+    f.close()
+
+def get_lines(infile):
+    f = file(infile,'r')
+    lines = f.readlines()
+    f.close()
+    return lines
+
+def add_lines(infile, outfile, lines):
+    orig_lines = get_lines(infile)
+    of = file(outfile,'w')
+    of.writelines(orig_lines)
+    of.writelines(lines)
+    of.close()
+
+def TempJump(mach, Gamma=None):
+    if Gamma is None:
+        Gamma = 5.0/3.0
+    M2 = mach*mach;
+    TJ = (2.0*Gamma*M2-(Gamma-1.0))*((Gamma-1.0)*M2+2.0)/(M2*(Gamma+1.0)**2);
+    return TJ
+
+simtime = 1.0*sec_per_Gyr # 1 Gyr
+boxsize = 1.0*cm_per_mpc # 1 Mpc box
+pre_shock_den = 1.0e0 # part/cc
+postT = 1.0e7 # K
+mach = 3.0 
+gas_gamma = 5./3.
+
+infile = 'input_shock.enzo'
+myname = 'CustomShockBox.enzo'
+
+# No modification is needed below here.
+
+preT = postT/TempJump(mach, gas_gamma)
+
+header = get_header(pre_shock_den, preT, postT, mach)
+inlines = get_lines(infile)
+customlines = setup_shockbox(boxsize, simtime, pre_shock_den, preT, mach,
+                             gamma=gas_gamma)
+
+write_lines(myname, header+inlines+customlines)
+

run/test_runner.py

@@ -47,7 +47,7 @@
 
 # Set the filename for the latest version of the gold standard
 # and for the default local standard output
-ytcfg["yt", "gold_standard_filename"] = str("enzogold0000")
+ytcfg["yt", "gold_standard_filename"] = str("enzogold0001")
 ytcfg["yt", "local_standard_filename"] = str("enzolocaldev")
 from yt.utilities.answer_testing.framework import \
     AnswerTesting

src/enzo/ClusterInitialize.C

@@ -15,9 +15,7 @@
 // This routine intializes a new simulation based on the parameter file.
 //
 
-#include <string.h>
-#include <stdio.h>
-#include <math.h>
+#include "preincludes.h"
 #include "macros_and_parameters.h"
 #include "typedefs.h"
 #include "global_data.h"