diff --git a/.hgtags b/.hgtags index 63323ff40..313515843 100644 --- a/.hgtags +++ b/.hgtags @@ -2,3 +2,4 @@ bbf0a2ffbd22c4fbecf946c9c96e6c4fac5cbdae woc_pre_fld_merge 48b4e9d9d6b90f703e48e621b488136be2a0e9cf woc_fld_merge b86d8ba026d6a0ec30f15d8134add1e55fae2958 Wise10_GalaxyBirth 2d90aa38e06f00a531db45a43225cde1faf093f2 enzo-2.2 +a444b00827d3d4351dceea8197f1e446b9ada0da enzogold0001 diff --git a/doc/manual/source/parameters/problemtypes.rst b/doc/manual/source/parameters/problemtypes.rst index b29ee6cd9..37c852d08 100644 --- a/doc/manual/source/parameters/problemtypes.rst +++ b/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. diff --git a/run/Hydro/Hydro-1D/ShockInABox/ShockInABox.enzo b/run/Hydro/Hydro-1D/ShockInABox/ShockInABox.enzo new file mode 100644 index 000000000..b3434abe2 --- /dev/null +++ b/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 + diff --git a/run/Hydro/Hydro-1D/ShockInABox/ShockInABox.enzotest b/run/Hydro/Hydro-1D/ShockInABox/ShockInABox.enzotest new file mode 100644 index 000000000..d572732dd --- /dev/null +++ b/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 diff --git a/run/Hydro/Hydro-1D/ShockInABox/ShockInABox__test_shockinabox.py b/run/Hydro/Hydro-1D/ShockInABox/ShockInABox__test_shockinabox.py new file mode 100644 index 000000000..c09d74bb7 --- /dev/null +++ b/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) + diff --git a/run/Hydro/Hydro-1D/ShockInABox/input_shock.enzo b/run/Hydro/Hydro-1D/ShockInABox/input_shock.enzo new file mode 100644 index 000000000..040fc6b5c --- /dev/null +++ b/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 + diff --git a/run/Hydro/Hydro-1D/ShockInABox/make_plots.py b/run/Hydro/Hydro-1D/ShockInABox/make_plots.py new file mode 100644 index 000000000..4504e41ea --- /dev/null +++ b/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() + diff --git a/run/Hydro/Hydro-1D/ShockInABox/make_shock.py b/run/Hydro/Hydro-1D/ShockInABox/make_shock.py new file mode 100755 index 000000000..f2dc9a08f --- /dev/null +++ b/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) + diff --git a/run/test_runner.py b/run/test_runner.py index 1fdb94164..eed4a4038 100755 --- a/run/test_runner.py +++ b/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 diff --git a/src/enzo/ClusterInitialize.C b/src/enzo/ClusterInitialize.C index 9a37b4c3d..07608654c 100644 --- a/src/enzo/ClusterInitialize.C +++ b/src/enzo/ClusterInitialize.C @@ -15,9 +15,7 @@ // This routine intializes a new simulation based on the parameter file. // -#include -#include -#include +#include "preincludes.h" #include "macros_and_parameters.h" #include "typedefs.h" #include "global_data.h" diff --git a/src/enzo/CollapseTestInitialize.C b/src/enzo/CollapseTestInitialize.C index 833512487..019238e4d 100644 --- a/src/enzo/CollapseTestInitialize.C +++ b/src/enzo/CollapseTestInitialize.C @@ -87,10 +87,13 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, CollapseTestSphereCutOff[MAX_SPHERES], CollapseTestSphereAng1[MAX_SPHERES], CollapseTestSphereAng2[MAX_SPHERES], - CollapseTestSphereMetallicity[MAX_SPHERES]; + CollapseTestSphereMetallicity[MAX_SPHERES], + CollapseTestSphereSmoothRadius[MAX_SPHERES]; int CollapseTestSphereNumShells[MAX_SPHERES], CollapseTestSphereInitialLevel[MAX_SPHERES], - CollapseTestSphereType[MAX_SPHERES]; + CollapseTestSphereType[MAX_SPHERES], + CollapseTestSphereConstantPressure[MAX_SPHERES], + CollapseTestSphereSmoothSurface[MAX_SPHERES]; FLOAT CollapseTestSphereRadius[MAX_SPHERES], CollapseTestSphereCoreRadius[MAX_SPHERES], CollapseTestSpherePosition[MAX_SPHERES][MAX_DIMENSION]; @@ -107,6 +110,7 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, CollapseTestSphereAng1[sphere] = 0; CollapseTestSphereAng2[sphere] = 0; CollapseTestSphereNumShells[sphere] = 1; + CollapseTestSphereSmoothRadius[sphere] = 1.2; CollapseTestSphereMetallicity[sphere] = 0; CollapseTestSphereInitialLevel[sphere] = 0; @@ -116,6 +120,8 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, CollapseTestSphereVelocity[sphere][dim] = 0; } CollapseTestSphereType[sphere] = 0; + CollapseTestSphereConstantPressure[sphere] = FALSE; + CollapseTestSphereSmoothSurface[sphere] = FALSE; } for (dim = 0; dim < MAX_DIMENSION; dim++) CollapseTestUniformVelocity[dim] = 0; @@ -199,6 +205,15 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, if (sscanf(line, "CollapseTestSphereInitialLevel[%"ISYM"]", &sphere) > 0) ret += sscanf(line, "CollapseTestSphereInitialLevel[%"ISYM"] = %"ISYM, &sphere, &CollapseTestSphereInitialLevel[sphere]); + if (sscanf(line, "CollapseTestSphereConstantPressure[%"ISYM"]", &sphere) > 0) + ret += sscanf(line, "CollapseTestSphereConstantPressure[%"ISYM"] = %"ISYM, &sphere, + &CollapseTestSphereConstantPressure[sphere]); + if (sscanf(line, "CollapseTestSphereSmoothSurface[%"ISYM"]", &sphere) > 0) + ret += sscanf(line, "CollapseTestSphereSmoothSurface[%"ISYM"] = %"ISYM, &sphere, + &CollapseTestSphereSmoothSurface[sphere]); + if (sscanf(line, "CollapseTestSphereSmoothRadius[%"ISYM"]", &sphere) > 0) + ret += sscanf(line, "CollapseTestSphereSmoothRadius[%"FSYM"] = %"FSYM, &sphere, + &CollapseTestSphereSmoothRadius[sphere]); /* if the line is suspicious, issue a warning */ @@ -219,8 +234,9 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, CollapseTestSphereDispersion, CollapseTestSphereCutOff, CollapseTestSphereAng1, CollapseTestSphereAng2, CollapseTestSphereNumShells, - CollapseTestSphereType, CollapseTestUseParticles, - CollapseTestParticleMeanDensity, + CollapseTestSphereType, CollapseTestSphereConstantPressure, + CollapseTestSphereSmoothSurface, CollapseTestSphereSmoothRadius, + CollapseTestUseParticles, CollapseTestParticleMeanDensity, CollapseTestUniformVelocity, CollapseTestUseColour, CollapseTestUseMetals, CollapseTestInitialTemperature, CollapseTestInitialDensity, @@ -322,8 +338,9 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, CollapseTestSphereDispersion, CollapseTestSphereCutOff, CollapseTestSphereAng1, CollapseTestSphereAng2, CollapseTestSphereNumShells, - CollapseTestSphereType, CollapseTestUseParticles, - CollapseTestParticleMeanDensity, + CollapseTestSphereType, CollapseTestSphereConstantPressure, + CollapseTestSphereSmoothSurface, CollapseTestSphereSmoothRadius, + CollapseTestUseParticles, CollapseTestParticleMeanDensity, CollapseTestUniformVelocity, CollapseTestUseColour, CollapseTestUseMetals, CollapseTestInitialTemperature, CollapseTestInitialDensity, @@ -365,8 +382,9 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, CollapseTestSphereDispersion, CollapseTestSphereCutOff, CollapseTestSphereAng1, CollapseTestSphereAng2, CollapseTestSphereNumShells, - CollapseTestSphereType, CollapseTestUseParticles, - CollapseTestParticleMeanDensity, + CollapseTestSphereType, CollapseTestSphereConstantPressure, + CollapseTestSphereSmoothSurface, CollapseTestSphereSmoothRadius, + CollapseTestUseParticles, CollapseTestParticleMeanDensity, CollapseTestUniformVelocity, CollapseTestUseColour, CollapseTestUseMetals, CollapseTestInitialTemperature, CollapseTestInitialDensity, @@ -480,6 +498,12 @@ int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, CollapseTestSphereAng2[sphere]); fprintf(Outfptr, "CollapseTestSphereNumShells[%"ISYM"] = %"ISYM"\n", sphere, CollapseTestSphereNumShells[sphere]); + fprintf(Outfptr, "CollapseTestSphereConstantPressure[%"ISYM"] = %"ISYM"\n", sphere, + CollapseTestSphereConstantPressure[sphere]); + fprintf(Outfptr, "CollapseTestSphereSmoothSurface[%"ISYM"] = %"ISYM"\n", sphere, + CollapseTestSphereSmoothSurface[sphere]); + fprintf(Outfptr, "CollapseTestSphereSmoothRadius[%"ISYM"] = %"GOUTSYM"\n", sphere, + CollapseTestSphereSmoothRadius[sphere]); } } diff --git a/src/enzo/CommunicationTranspose.C b/src/enzo/CommunicationTranspose.C index faed0fdcc..82e121d4a 100644 --- a/src/enzo/CommunicationTranspose.C +++ b/src/enzo/CommunicationTranspose.C @@ -14,6 +14,7 @@ #include "mpi.h" #endif /* USE_MPI */ +#include "EnzoTiming.h" #include #include "ErrorExceptions.h" #include "macros_and_parameters.h" @@ -73,6 +74,7 @@ int CommunicationTranspose(region *FromRegion, int NumberOfFromRegions, region *ToRegion, int NumberOfToRegions, int TransposeOrder) { + TIMER_START("CommunicationTranspose"); int retval; switch (UnigridTranspose) { case 0: @@ -93,6 +95,7 @@ int CommunicationTranspose(region *FromRegion, int NumberOfFromRegions, default: ENZO_VFAIL("Invalid value for UnigridTranspose = %d", UnigridTranspose); } // ENDSWITCH + TIMER_STOP("CommunicationTranspose"); return retval; } diff --git a/src/enzo/Grid.h b/src/enzo/Grid.h index f85b98a4f..c4f76df33 100644 --- a/src/enzo/Grid.h +++ b/src/enzo/Grid.h @@ -1922,6 +1922,9 @@ int zEulerSweep(int j, int NumberOfSubgrids, fluxes *SubgridFluxes[], float SphereAng2[MAX_SPHERES], int SphereNumShells[MAX_SPHERES], int SphereType[MAX_SPHERES], + int SphereConstantPressure[MAX_SPHERES], + int SphereSmoothSurface[MAX_SPHERES], + float SphereSmoothRadius[MAX_SPHERES], int SphereUseParticles, float ParticleMeanDensity, float UniformVelocity[MAX_DIMENSION], diff --git a/src/enzo/Grid_ClusterInitializeGrid.C b/src/enzo/Grid_ClusterInitializeGrid.C index 7dc8d0999..a93870761 100644 --- a/src/enzo/Grid_ClusterInitializeGrid.C +++ b/src/enzo/Grid_ClusterInitializeGrid.C @@ -12,9 +12,7 @@ / ************************************************************************/ -#include -#include -#include +#include "preincludes.h" #include "macros_and_parameters.h" #include "typedefs.h" #include "global_data.h" diff --git a/src/enzo/Grid_CollapseTestInitializeGrid.C b/src/enzo/Grid_CollapseTestInitializeGrid.C index 4d2fdf804..9531578fc 100644 --- a/src/enzo/Grid_CollapseTestInitializeGrid.C +++ b/src/enzo/Grid_CollapseTestInitializeGrid.C @@ -78,6 +78,9 @@ int grid::CollapseTestInitializeGrid(int NumberOfSpheres, float SphereAng2[MAX_SPHERES], int SphereNumShells[MAX_SPHERES], int SphereType[MAX_SPHERES], + int SphereConstantPressure[MAX_SPHERES], + int SphereSmoothSurface[MAX_SPHERES], + float SphereSmoothRadius[MAX_SPHERES], int SphereUseParticles, float ParticleMeanDensity, float UniformVelocity[MAX_DIMENSION], @@ -275,7 +278,8 @@ int grid::CollapseTestInitializeGrid(int NumberOfSpheres, float density, dens1, old_density, Velocity[MAX_DIMENSION], DiskVelocity[MAX_DIMENSION], temperature, temp1, sigma, sigma1, - colour, weight, a, DMVelocity[MAX_DIMENSION], metallicity; + colour, weight, a, DMVelocity[MAX_DIMENSION], metallicity, + outer_radius; FLOAT r, rcyl, x, y = 0, z = 0; int n = 0, ibin; @@ -477,7 +481,9 @@ int grid::CollapseTestInitializeGrid(int NumberOfSpheres, /* Compute Cartesian coordinates for rotational properties */ - if (r < SphereRadius[sphere]) { + outer_radius = (SphereSmoothSurface[sphere] == TRUE) ? + SphereSmoothRadius[sphere]*SphereRadius[sphere] : SphereRadius[sphere]; + if (r < outer_radius) { /* Compute spherical coordinate theta */ @@ -755,10 +761,15 @@ int grid::CollapseTestInitializeGrid(int NumberOfSpheres, // temp1 = InitialTemperature; if (SphereType[sphere] != 7 && SphereType[sphere] != 9) - if (temp1 == InitialTemperature) - temperature = SphereTemperature[sphere]; - else + if (temp1 == InitialTemperature) { + if (SphereConstantPressure[sphere] == TRUE) { + temperature = SphereTemperature[sphere] * (SphereDensity[sphere] / dens1); + } else { + temperature = SphereTemperature[sphere]; + } + } else { temperature = temp1; + } sigma = sigma1; if (SphereType[sphere] != 10) @@ -828,6 +839,18 @@ int grid::CollapseTestInitializeGrid(int NumberOfSpheres, } + if (SphereSmoothSurface[sphere] == TRUE && + r < SphereSmoothRadius[sphere]*SphereRadius[sphere] && + r > SphereRadius[sphere]) { + float ramp = 1.0 - 1.0 * tanh((3.0/(SphereSmoothRadius[sphere]-1.0))* + (r/SphereRadius[sphere] - 1.0)); + ramp = max(ramp, 1.0/density); + density *= ramp; + if (SphereConstantPressure[sphere] == TRUE) { + temperature /= ramp; + } + } // end: if (SmoothSurface) + } // end: loop over spheres /* Set density. */ diff --git a/src/enzo/Grid_CosmologySimulationInitializeGrid.C b/src/enzo/Grid_CosmologySimulationInitializeGrid.C index c31f18b1f..389a9a609 100644 --- a/src/enzo/Grid_CosmologySimulationInitializeGrid.C +++ b/src/enzo/Grid_CosmologySimulationInitializeGrid.C @@ -178,7 +178,7 @@ int grid::CosmologySimulationInitializeGrid( if ( (ParallelRootGridIO != TRUE) || (ParallelParticleIO != TRUE && !CosmologySimulationCalculatePositions) ) { - ENZO_FAIL("ParallelRootGridIO and ParallelParticleIO MUST be set for > 64 cpus!"); + fprintf(stderr, "WARNING: ParallelRootGridIO and ParallelParticleIO are recommended for > 64 cpus to shrink data read-in time."); } } diff --git a/src/enzo/Grid_HydroShockTubesInitializeGrid.C b/src/enzo/Grid_HydroShockTubesInitializeGrid.C index 87b052c8d..1203b4ee5 100644 --- a/src/enzo/Grid_HydroShockTubesInitializeGrid.C +++ b/src/enzo/Grid_HydroShockTubesInitializeGrid.C @@ -24,6 +24,8 @@ int grid::HydroShockTubesInitializeGrid(float x0, ) { + int MachNum, PSTempNum, PSDenNum; + NumberOfBaryonFields = 0; FieldType[NumberOfBaryonFields++] = Density; FieldType[NumberOfBaryonFields++] = Velocity1; @@ -34,19 +36,24 @@ int grid::HydroShockTubesInitializeGrid(float x0, FieldType[NumberOfBaryonFields++] = InternalEnergy; } + if(ShockMethod){ + FieldType[MachNum = NumberOfBaryonFields++] = Mach; + if(StorePreShockFields){ + FieldType[PSTempNum = NumberOfBaryonFields++] = PreShockTemperature; + FieldType[PSDenNum = NumberOfBaryonFields++] = PreShockDensity; + } + } + if (ProcessorNumber != MyProcessorNumber) { return SUCCESS; } - int size = 1, activesize = 1, dim; + int size = 1, index, dim; for (dim = 0; dim < GridRank; dim++) size *= GridDimension[dim]; - for (dim = 0; dim < GridRank; dim++) - activesize *= (GridDimension[dim] - 2*NumberOfGhostZones); - int field; for (field = 0; field < NumberOfBaryonFields; field++) if (BaryonField[field] == NULL) @@ -63,29 +70,44 @@ int grid::HydroShockTubesInitializeGrid(float x0, FLOAT x; int i; + for (int k = 0; k < GridDimension[2]; k++) { + for (int j = 0; j < GridDimension[1]; j++) { for (i = 0; i < GridDimension[0]; i++) { x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; + index = GRIDINDEX_NOGHOST(i,j,k); if (x <= x0) { - BaryonField[iden ][i] = rhol; - BaryonField[ivx ][i] = vxl; - BaryonField[ivy ][i] = vyl; - BaryonField[ivz ][i] = vzl; - BaryonField[ietot][i] = etotl; + BaryonField[iden ][index] = rhol; + BaryonField[ivx ][index] = vxl; + BaryonField[ivy ][index] = vyl; + BaryonField[ivz ][index] = vzl; + BaryonField[ietot][index] = etotl; if (DualEnergyFormalism) { - BaryonField[ieint][i] = etotl - 0.5*(vxl*vxl+vyl*vyl+vzl*vzl); + BaryonField[ieint][index] = etotl - 0.5*(vxl*vxl+vyl*vyl+vzl*vzl); } } else { - BaryonField[iden ][i] = rhor; - BaryonField[ivx ][i] = vxr; - BaryonField[ivy ][i] = vyr; - BaryonField[ivz ][i] = vzr; - BaryonField[ietot][i] = etotr; + BaryonField[iden ][index] = rhor; + BaryonField[ivx ][index] = vxr; + BaryonField[ivy ][index] = vyr; + BaryonField[ivz ][index] = vzr; + BaryonField[ietot][index] = etotr; if (DualEnergyFormalism) { - BaryonField[ieint][i] = etotr - 0.5*(vxr*vxr+vyr*vyr+vzr*vzr); + BaryonField[ieint][index] = etotr - 0.5*(vxr*vxr+vyr*vyr+vzr*vzr); } } + + //Shock + if (ShockMethod) { + BaryonField[MachNum][index] = tiny_number; + if (StorePreShockFields) { + BaryonField[PSTempNum][index] = tiny_number; + BaryonField[PSDenNum][index] = tiny_number; + } + } + + } + } } return SUCCESS; @@ -102,6 +124,8 @@ int grid::HydroShockTubesInitializeGrid(float x0, float x1, ) { + int MachNum, PSTempNum, PSDenNum; + NumberOfBaryonFields = 0; FieldType[NumberOfBaryonFields++] = Density; FieldType[NumberOfBaryonFields++] = Velocity1; @@ -112,24 +136,27 @@ int grid::HydroShockTubesInitializeGrid(float x0, float x1, FieldType[NumberOfBaryonFields++] = InternalEnergy; } + if(ShockMethod){ + FieldType[MachNum = NumberOfBaryonFields++] = Mach; + if(StorePreShockFields){ + FieldType[PSTempNum = NumberOfBaryonFields++] = PreShockTemperature; + FieldType[PSDenNum = NumberOfBaryonFields++] = PreShockDensity; + } + } if (ProcessorNumber != MyProcessorNumber) { return SUCCESS; } - int size = 1, activesize = 1, dim; + int size = 1, dim, index; for (dim = 0; dim < GridRank; dim++) size *= GridDimension[dim]; - for (dim = 0; dim < GridRank; dim++) - activesize *= (GridDimension[dim] - 2*NumberOfGhostZones); - int field; for (field = 0; field < NumberOfBaryonFields; field++) if (BaryonField[field] == NULL) BaryonField[field] = new float[size]; - /* transform pressure to total energy */ float etotl, etotr, etotc, v2; @@ -144,39 +171,53 @@ int grid::HydroShockTubesInitializeGrid(float x0, float x1, FLOAT x; int i; + for (int k = 0; k < GridDimension[2]; k++) { + for (int j = 0; j < GridDimension[1]; j++) { for (i = 0; i < GridDimension[0]; i++) { x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i]; + index = GRIDINDEX_NOGHOST(i,j,k); if (x <= x0) { - BaryonField[iden ][i] = rhol; - BaryonField[ivx ][i] = vxl; - BaryonField[ivy ][i] = vyl; - BaryonField[ivz ][i] = vzl; - BaryonField[ietot][i] = etotl; + BaryonField[iden ][index] = rhol; + BaryonField[ivx ][index] = vxl; + BaryonField[ivy ][index] = vyl; + BaryonField[ivz ][index] = vzl; + BaryonField[ietot][index] = etotl; if (DualEnergyFormalism) { - BaryonField[ieint][i] = etotl - 0.5*(vxl*vxl+vyl*vyl+vzl*vzl); + BaryonField[ieint][index] = etotl - 0.5*(vxl*vxl+vyl*vyl+vzl*vzl); } } else if (x <= x1) { - BaryonField[iden ][i] = rhoc; - BaryonField[ivx ][i] = vxc; - BaryonField[ivy ][i] = vyc; - BaryonField[ivz ][i] = vzc; - BaryonField[ietot][i] = etotc; + BaryonField[iden ][index] = rhoc; + BaryonField[ivx ][index] = vxc; + BaryonField[ivy ][index] = vyc; + BaryonField[ivz ][index] = vzc; + BaryonField[ietot][index] = etotc; if (DualEnergyFormalism) { - BaryonField[ieint][i] = etotc - 0.5*(vxc*vxc+vyc*vyc+vzc*vzc); + BaryonField[ieint][index] = etotc - 0.5*(vxc*vxc+vyc*vyc+vzc*vzc); } } else { - BaryonField[iden ][i] = rhor; - BaryonField[ivx ][i] = vxr; - BaryonField[ivy ][i] = vyr; - BaryonField[ivz ][i] = vzr; - BaryonField[ietot][i] = etotr; + BaryonField[iden ][index] = rhor; + BaryonField[ivx ][index] = vxr; + BaryonField[ivy ][index] = vyr; + BaryonField[ivz ][index] = vzr; + BaryonField[ietot][index] = etotr; if (DualEnergyFormalism) { - BaryonField[ieint][i] = etotr - 0.5*(vxr*vxr+vyr*vyr+vzr*vzr); + BaryonField[ieint][index] = etotr - 0.5*(vxr*vxr+vyr*vyr+vzr*vzr); } } + + //Shock + if (ShockMethod) { + BaryonField[MachNum][index] = tiny_number; + if (StorePreShockFields) { + BaryonField[PSTempNum][index] = tiny_number; + BaryonField[PSDenNum][index] = tiny_number; + } + } + } + } } return SUCCESS; diff --git a/src/enzo/Grid_NestedCosmologySimulationInitializeGrid.C b/src/enzo/Grid_NestedCosmologySimulationInitializeGrid.C index 888a83920..d75c27d35 100644 --- a/src/enzo/Grid_NestedCosmologySimulationInitializeGrid.C +++ b/src/enzo/Grid_NestedCosmologySimulationInitializeGrid.C @@ -181,7 +181,7 @@ int grid::NestedCosmologySimulationInitializeGrid( { if ( (ParallelRootGridIO != TRUE) && (ParallelParticleIO != TRUE) ) { - // ENZO_FAIL("ParallelRootGridIO (and possibly ParallelParticleIO) should be set for > 64 cpus!\n"); + fprintf(stderr, "WARNING: ParallelRootGridIO and ParallelParticleIO are recommended for > 64 cpus to shrink data read-in time."); } } diff --git a/src/enzo/Grid_PhotonTestInitializeGrid.C b/src/enzo/Grid_PhotonTestInitializeGrid.C index 03d502238..2be6de981 100644 --- a/src/enzo/Grid_PhotonTestInitializeGrid.C +++ b/src/enzo/Grid_PhotonTestInitializeGrid.C @@ -60,6 +60,9 @@ int grid::PhotonTestInitializeGrid(int NumberOfSpheres, float SphereAng2[MAX_SPHERES], int SphereNumShells[MAX_SPHERES], int SphereType[MAX_SPHERES], + int SphereConstantPressure[MAX_SPHERES], + int SphereSmoothSurface[MAX_SPHERES], + float SphereSmoothRadius[MAX_SPHERES], float SphereHII[MAX_SPHERES], float SphereHeII[MAX_SPHERES], float SphereHeIII[MAX_SPHERES], @@ -298,7 +301,7 @@ int grid::PhotonTestInitializeGrid(int NumberOfSpheres, /* Loop over the mesh. */ float density, dens1, Velocity[MAX_DIMENSION], - temperature, temp1, sigma, sigma1, colour; + temperature, temp1, sigma, sigma1, colour, outer_radius; float HII_Fraction, HeII_Fraction, HeIII_Fraction, H2I_Fraction; FLOAT r, x, y = 0, z = 0; int n = 0; @@ -415,7 +418,9 @@ int grid::PhotonTestInitializeGrid(int NumberOfSpheres, pow(fabs(z-SpherePosition[sphere][2]), 2) ); r = max(r, 0.1*CellWidth[0][0]); - if (r < SphereRadius[sphere]) { + outer_radius = (SphereSmoothSurface[sphere] == TRUE) ? + SphereSmoothRadius[sphere]*SphereRadius[sphere] : SphereRadius[sphere]; + if (r < outer_radius) { /* Compute Cartesian coordinates for rotational properties */ @@ -648,9 +653,16 @@ int grid::PhotonTestInitializeGrid(int NumberOfSpheres, if (dens1 > density) { density = dens1; - if (temp1 == InitialTemperature) - temp1 = SphereTemperature[sphere]; - temperature = temp1; + if (SphereType[sphere] != 7 && SphereType[sphere] != 9) + if (temp1 == InitialTemperature) { + if (SphereConstantPressure[sphere] == TRUE) { + temperature = SphereTemperature[sphere] * (SphereDensity[sphere] / dens1); + } else { + temperature = SphereTemperature[sphere]; + } + } else { + temperature = temp1; + } sigma = sigma1; if (SphereType[sphere] != 6 && SphereType[sphere] != 10 && @@ -666,6 +678,19 @@ int grid::PhotonTestInitializeGrid(int NumberOfSpheres, } } // end: if (r < SphereRadius) + + if (SphereSmoothSurface[sphere] == TRUE && + r < SphereSmoothRadius[sphere]*SphereRadius[sphere] && + r > SphereRadius[sphere]) { + float ramp = 1.0 - 1.0 * tanh((3.0/(SphereSmoothRadius[sphere]-1.0))* + (r/SphereRadius[sphere] - 1.0)); + ramp = max(ramp, 1.0/density); + density *= ramp; + if (SphereConstantPressure[sphere] == TRUE) { + temperature /= ramp; + } + } // end: if (SmoothSurface) + } // end: loop over spheres /* Set density. */ diff --git a/src/enzo/HydroShockTubesInitialize.C b/src/enzo/HydroShockTubesInitialize.C index dbc794ae2..f9de16232 100644 --- a/src/enzo/HydroShockTubesInitialize.C +++ b/src/enzo/HydroShockTubesInitialize.C @@ -1,6 +1,4 @@ -#include -#include -#include +#include "preincludes.h" #include "macros_and_parameters.h" #include "typedefs.h" #include "global_data.h" @@ -35,6 +33,9 @@ int HydroShockTubesInitialize(FILE *fptr, FILE *Outfptr, char *Vel2Name = "y-velocity"; char *Vel3Name = "z-velocity"; char *ColourName = "colour"; + char *MachName = "Mach"; + char *PSTempName = "PreShock_Temperature"; + char *PSDenName = "PreShock_Density"; /* declarations */ @@ -209,6 +210,14 @@ int HydroShockTubesInitialize(FILE *fptr, FILE *Outfptr, DataLabel[count++] = GEName; } + if (ShockMethod) { + DataLabel[count++] = MachName; + if(StorePreShockFields){ + DataLabel[count++] = PSTempName; + DataLabel[count++] = PSDenName; + } + } + for (i = 0; i < count; i++) DataUnits[i] = NULL; diff --git a/src/enzo/InitializePythonInterface.C b/src/enzo/InitializePythonInterface.C index 9e6d737ea..c34537021 100644 --- a/src/enzo/InitializePythonInterface.C +++ b/src/enzo/InitializePythonInterface.C @@ -264,6 +264,13 @@ void ExportParameterFile(TopGridData *MetaData, FLOAT CurrentTime, FLOAT OldTime PyDict_SetItemString(yt_parameter_file, "TopGridDimensions", tgd_tuple); Py_XDECREF(tgd_tuple); Py_XDECREF(tgd0); Py_XDECREF(tgd1); Py_XDECREF(tgd2); + tgd0 = PyLong_FromLong((long) MetaData->LeftFaceBoundaryCondition[0]); + tgd1 = PyLong_FromLong((long) MetaData->LeftFaceBoundaryCondition[1]); + tgd2 = PyLong_FromLong((long) MetaData->LeftFaceBoundaryCondition[2]); + tgd_tuple = PyTuple_Pack(3, tgd0, tgd1, tgd2); + PyDict_SetItemString(yt_parameter_file, "LeftFaceBoundaryCondition", tgd_tuple); + Py_XDECREF(tgd_tuple); Py_XDECREF(tgd0); Py_XDECREF(tgd1); Py_XDECREF(tgd2); + /* Do the conversion factors */ for (int dim = 0; dim < MAX_NUMBER_OF_BARYON_FIELDS; dim++) { if (DataLabel[dim]) { diff --git a/src/enzo/Make.mach.arizona b/src/enzo/Make.mach.arizona index cd61cc3a2..cfc2cdb4e 100644 --- a/src/enzo/Make.mach.arizona +++ b/src/enzo/Make.mach.arizona @@ -1,8 +1,8 @@ #======================================================================= # -# FILE: Make.mach.UofA +# FILE: Make.mach.arizona # -# DESCRIPTION: Makefile settings for University of Arizona HPC Cluster +# DESCRIPTION: Makefile settings for University of Arizona SMP Cluster # # AUTHOR: Cameron Hummels (chummels@email.arizona.edu) # @@ -10,32 +10,32 @@ # #======================================================================= -MACH_TEXT = University of Arizona's HPC Clusters +MACH_TEXT = University of Arizona SMP Cluster MACH_VALID = 1 MACH_FILE = Make.mach.arizona -MACHINE_NOTES = "MACHINE_NOTES for University of Arizona's HPC Cluster: \ +MACHINE_NOTES = " MACHINE_NOTES for University of Arizona HPC Cluster: \ The following modules are needed to compile enzo here with \ the default machine file. \ - hdf5/1.8.8 \ - intel/2012.0.032 \ - mpt/2.04 (for SGI MPI compilers) \ \ You can run this command, \ - $ module load hdf5 mpt intel" + $ module load intel mpt hdf5" #----------------------------------------------------------------------- # Install paths (local variables) #----------------------------------------------------------------------- -LOCAL_HDF5_INSTALL = /uaopt/hdf5/1.8.8 -LOCAL_COMPILER_DIR = /gsfs1/uaopt/intel/2012.0.032/composer_xe_2011_sp1.6.233/compiler +LOCAL_MPI_INSTALL = +LOCAL_HDF5_INSTALL = /uaopt/hdf5/1.8.8/ +#LOCAL_HDF5_INSTALL = $(YT_DEST) +#LOCAL_PYTHON_INSTALL = $(YT_DEST) +LOCAL_COMPILER_DIR = #----------------------------------------------------------------------- # Compiler settings #----------------------------------------------------------------------- -MACH_CPP = cpp # C preprocessor command +MACH_CPP = icpc # C preprocessor command # With MPI # Preferred compiler is the SGI library @@ -70,7 +70,7 @@ MACH_DEFINES = -DLINUX -DH5_USE_16_API MACH_CPPFLAGS = -P -traditional MACH_CFLAGS = -MACH_CXXFLAGS = +MACH_CXXFLAGS = -DMPICH_IGNORE_CXX_SEEK MACH_FFLAGS = MACH_F90FLAGS = MACH_LDFLAGS = @@ -88,13 +88,12 @@ MACH_OPT_AGGRESSIVE = -O3 # Includes #----------------------------------------------------------------------- -LOCAL_INCLUDES_MPI = # MPI includes +LOCAL_INCLUDES_MPI = -I$(LOCAL_MPI_INSTALL)/include # MPI includes LOCAL_INCLUDES_HDF5 = -I$(LOCAL_HDF5_INSTALL)/include # HDF5 includes LOCAL_INCLUDES_HYPRE = # hypre includes LOCAL_INCLUDES_PAPI = # PAPI includes MACH_INCLUDES = $(LOCAL_INCLUDES_HDF5) - MACH_INCLUDES_MPI = $(LOCAL_INCLUDES_MPI) MACH_INCLUDES_HYPRE = $(LOCAL_INCLUDES_HYPRE) MACH_INCLUDES_PAPI = $(LOCAL_INCLUDES_PAPI) @@ -103,12 +102,13 @@ MACH_INCLUDES_PAPI = $(LOCAL_INCLUDES_PAPI) # Libraries #----------------------------------------------------------------------- -LOCAL_LIBS_MPI = # MPI libraries +#LOCAL_LIBS_MPI = -L$(LOCAL_MPI_INSTALL)/lib -lmpi # MPI libraries +LOCAL_LIBS_MPI = -L$(LOCAL_MPI_INSTALL)/lib -lmpi++ # MPI libraries LOCAL_LIBS_HDF5 = -L$(LOCAL_HDF5_INSTALL)/lib -lhdf5 -lz # HDF5 libraries LOCAL_LIBS_HYPRE = # hypre libraries LOCAL_LIBS_PAPI = # PAPI libraries -LOCAL_LIBS_MACH = -L $(LOCAL_COMPILER_DIR)/lib/intel64 -lifcore -lifport -lsvml # Machine-dependent libraries +LOCAL_LIBS_MACH = -L $(LOCAL_COMPILER_DIR) -limf -lifcore -lifport -lsvml -lgfortran # Machine-dependent libraries MACH_LIBS = $(LOCAL_LIBS_HDF5) $(LOCAL_LIBS_MACH) MACH_LIBS_MPI = $(LOCAL_LIBS_MPI) diff --git a/src/enzo/PhotonGrid_Methods.h b/src/enzo/PhotonGrid_Methods.h index 38f4cd09c..73d8561de 100644 --- a/src/enzo/PhotonGrid_Methods.h +++ b/src/enzo/PhotonGrid_Methods.h @@ -412,6 +412,9 @@ int PhotonTestInitializeGrid(int NumberOfSpheres, float SphereAng2[MAX_SPHERES], int SphereNumShells[MAX_SPHERES], int SphereType[MAX_SPHERES], + int SphereConstantPressure[MAX_SPHERES], + int SphereSmoothSurface[MAX_SPHERES], + float SphereSmoothRadius[MAX_SPHERES], float SphereHII[MAX_SPHERES], float SphereHeII[MAX_SPHERES], float SphereHeIII[MAX_SPHERES], diff --git a/src/enzo/PhotonTestInitialize.C b/src/enzo/PhotonTestInitialize.C index 0f40e5cae..ed316e139 100644 --- a/src/enzo/PhotonTestInitialize.C +++ b/src/enzo/PhotonTestInitialize.C @@ -94,7 +94,9 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, int PhotonTestUseParticles = FALSE; int PhotonTestUseColour = FALSE; float PhotonTestInitialTemperature = 1000; - int PhotonTestSphereType[MAX_SPHERES]; + int PhotonTestSphereType[MAX_SPHERES], + PhotonTestSphereConstantPressure[MAX_SPHERES], + PhotonTestSphereSmoothSurface[MAX_SPHERES]; float PhotonTestSphereDensity[MAX_SPHERES], PhotonTestSphereTemperature[MAX_SPHERES], PhotonTestSphereVelocity[MAX_SPHERES][MAX_DIMENSION], @@ -104,6 +106,7 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, PhotonTestSphereCutOff[MAX_SPHERES], PhotonTestSphereAng1[MAX_SPHERES], PhotonTestSphereAng2[MAX_SPHERES], + PhotonTestSphereSmoothRadius[MAX_SPHERES], PhotonTestSphereRadius[MAX_SPHERES], PhotonTestSphereCoreRadius[MAX_SPHERES], PhotonTestSphereHIIFraction[MAX_SPHERES], @@ -137,6 +140,7 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, PhotonTestSphereAng1[sphere] = 0; PhotonTestSphereAng2[sphere] = 0; PhotonTestSphereNumShells[sphere] = 1; + PhotonTestSphereSmoothRadius[sphere] = 1.2; PhotonTestSphereHIIFraction[sphere] = PhotonTestInitialFractionHII; PhotonTestSphereHeIIFraction[sphere] = PhotonTestInitialFractionHeII; PhotonTestSphereHeIIIFraction[sphere] = PhotonTestInitialFractionHeIII; @@ -148,6 +152,8 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, PhotonTestSphereVelocity[sphere][dim] = 0; } PhotonTestSphereType[sphere] = 0; + PhotonTestSphereConstantPressure[sphere] = FALSE; + PhotonTestSphereSmoothSurface[sphere] = FALSE; } for (dim = 0; dim < MAX_DIMENSION; dim++) PhotonTestUniformVelocity[dim] = 0; @@ -178,6 +184,15 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, if (sscanf(line, "PhotonTestSphereType[%"ISYM"]", &sphere) > 0) ret += sscanf(line, "PhotonTestSphereType[%"ISYM"] = %"ISYM, &sphere, &PhotonTestSphereType[sphere]); + if (sscanf(line, "PhotonTestSphereConstantPressure[%"ISYM"]", &sphere) > 0) + ret += sscanf(line, "PhotonTestSphereConstantPressure[%"ISYM"] = %"ISYM, &sphere, + &PhotonTestSphereConstantPressure[sphere]); + if (sscanf(line, "PhotonTestSphereSmoothSurface[%"ISYM"]", &sphere) > 0) + ret += sscanf(line, "PhotonTestSphereSmoothSurface[%"ISYM"] = %"ISYM, &sphere, + &PhotonTestSphereSmoothSurface[sphere]); + if (sscanf(line, "PhotonTestSphereSmoothRadius[%"ISYM"]", &sphere) > 0) + ret += sscanf(line, "PhotonTestSphereSmoothRadius[%"FSYM"] = %"FSYM, &sphere, + &PhotonTestSphereSmoothRadius[sphere]); if (sscanf(line, "PhotonTestSphereRadius[%"ISYM"]", &sphere) > 0) ret += sscanf(line, "PhotonTestSphereRadius[%"ISYM"] = %"FSYM, &sphere, &PhotonTestSphereRadius[sphere]); @@ -291,7 +306,8 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, PhotonTestFracKeplerianRot, PhotonTestSphereTurbulence, PhotonTestSphereCutOff, PhotonTestSphereAng1, PhotonTestSphereAng2, PhotonTestSphereNumShells, - PhotonTestSphereType, + PhotonTestSphereType, PhotonTestSphereConstantPressure, + PhotonTestSphereSmoothSurface, PhotonTestSphereSmoothRadius, PhotonTestSphereHIIFraction, PhotonTestSphereHeIIFraction, PhotonTestSphereHeIIIFraction, PhotonTestSphereH2IFraction, PhotonTestUseParticles, @@ -358,7 +374,8 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, PhotonTestFracKeplerianRot, PhotonTestSphereTurbulence, PhotonTestSphereCutOff, PhotonTestSphereAng1, PhotonTestSphereAng2, PhotonTestSphereNumShells, - PhotonTestSphereType, + PhotonTestSphereType, PhotonTestSphereConstantPressure, + PhotonTestSphereSmoothSurface, PhotonTestSphereSmoothRadius, PhotonTestSphereHIIFraction, PhotonTestSphereHeIIFraction, PhotonTestSphereHeIIIFraction, PhotonTestSphereH2IFraction, PhotonTestUseParticles, @@ -411,7 +428,8 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, PhotonTestFracKeplerianRot, PhotonTestSphereTurbulence, PhotonTestSphereCutOff, PhotonTestSphereAng1, PhotonTestSphereAng2, PhotonTestSphereNumShells, - PhotonTestSphereType, + PhotonTestSphereType, PhotonTestSphereConstantPressure, + PhotonTestSphereSmoothSurface, PhotonTestSphereSmoothRadius, PhotonTestSphereHIIFraction, PhotonTestSphereHeIIFraction, PhotonTestSphereHeIIIFraction, PhotonTestSphereH2IFraction, PhotonTestUseParticles, @@ -521,6 +539,12 @@ int PhotonTestInitialize(FILE *fptr, FILE *Outfptr, for (sphere = 0; sphere < PhotonTestNumberOfSpheres; sphere++) { fprintf(Outfptr, "PhotonTestSphereType[%"ISYM"] = %"ISYM"\n", sphere, PhotonTestSphereType[sphere]); + fprintf(Outfptr, "PhotonTestSphereConstantPressure[%"ISYM"] = %"ISYM"\n", sphere, + PhotonTestSphereConstantPressure[sphere]); + fprintf(Outfptr, "PhotonTestSphereSmoothSurface[%"ISYM"] = %"ISYM"\n", sphere, + PhotonTestSphereSmoothSurface[sphere]); + fprintf(Outfptr, "PhotonTestSphereSmoothRadius[%"ISYM"] = %"GOUTSYM"\n", sphere, + PhotonTestSphereSmoothRadius[sphere]); fprintf(Outfptr, "PhotonTestSphereRadius[%"ISYM"] = %"GOUTSYM"\n", sphere, PhotonTestSphereRadius[sphere]); fprintf(Outfptr, "PhotonTestSphereCoreRadius[%"ISYM"] = %"GOUTSYM"\n", sphere, diff --git a/src/enzo/ShockInABoxInitialize.C b/src/enzo/ShockInABoxInitialize.C index 96358b040..32f389532 100644 --- a/src/enzo/ShockInABoxInitialize.C +++ b/src/enzo/ShockInABoxInitialize.C @@ -38,6 +38,9 @@ int ShockInABoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, char *Vel2Name = "y-velocity"; char *Vel3Name = "z-velocity"; + char *MachName = "Mach"; + char *PSTempName = "PreShock_Temperature"; + char *PSDenName = "PreShock_Density"; /* declarations */ char line[MAX_LINE_LENGTH]; @@ -81,18 +84,18 @@ int ShockInABoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, ret += sscanf(line, "ShockInABoxBoundary = %"PSYM, &ShockInABoxBoundary); - ret += sscanf(line, "ShockInABoxLeftDensity = %"FSYM, + ret += sscanf(line, "ShockInABoxLeftDensity = %"ESYM, &ShockInABoxDensity[0]); - ret += sscanf(line, "ShockInABoxLeftPressure = %"FSYM, + ret += sscanf(line, "ShockInABoxLeftPressure = %"ESYM, &ShockInABoxPressure[0]); - ret += sscanf(line, "ShockInABoxLeftVelocity = %"FSYM, + ret += sscanf(line, "ShockInABoxLeftVelocity = %"ESYM, &ShockInABoxVelocity[0]); - ret += sscanf(line, "ShockInABoxRightDensity = %"FSYM, + ret += sscanf(line, "ShockInABoxRightDensity = %"ESYM, &ShockInABoxDensity[1]); - ret += sscanf(line, "ShockInABoxRightPressure = %"FSYM, + ret += sscanf(line, "ShockInABoxRightPressure = %"ESYM, &ShockInABoxPressure[1]); - ret += sscanf(line, "ShockInABoxRightVelocity = %"FSYM, + ret += sscanf(line, "ShockInABoxRightVelocity = %"ESYM, &ShockInABoxVelocity[1]); ret += sscanf(line, "ShockInABoxSubgridLeft = %"PSYM, @@ -103,7 +106,7 @@ int ShockInABoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, /* if the line is suspicious, issue a warning */ if (ret == 0 && strstr(line, "=") && strstr(line, "ShockInABox")) - fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line); + fprintf(stderr, "ShockInABox warning: the following parameter line was not interpreted:\n%s\n", line); } @@ -112,14 +115,6 @@ int ShockInABoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, if( ShockInABoxDirection != 0 ) ENZO_FAIL("Only ShockInABoxDirection=0 supported at the moment!"); -// if (TopGrid.GridData->ShockTubeInitializeGrid(ShockInABoxDirection, -// ShockInABoxBoundary, -// ShockInABoxDensity, -// ShockInABoxPressure, -// ShockInABoxVelocity) == FAIL) { -// ENZO_FAIL("Error in ShockTubeInitializeGrid.\n"); -// } - if (TopGrid.GridData-> HydroShockTubesInitializeGrid(ShockInABoxBoundary, ShockInABoxDensity[0], ShockInABoxDensity[1], @@ -174,17 +169,46 @@ int ShockInABoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, } } + + /* set up field names and units */ + + DataLabel[0] = DensName; + DataLabel[1] = Vel1Name; + DataLabel[2] = Vel2Name; + DataLabel[3] = Vel3Name; + DataLabel[4] = TEName; + + DataUnits[0] = NULL; + DataUnits[1] = NULL; + DataUnits[2] = NULL; + DataUnits[3] = NULL; + DataUnits[4] = NULL; + + int field_counter = 5; + + if (ShockMethod) { + DataLabel[field_counter++] = MachName; + if(StorePreShockFields){ + DataLabel[field_counter++] = PSTempName; + DataLabel[field_counter++] = PSDenName; + } + } + for (int j = 0; j < field_counter; j++) + DataUnits[j] = NULL; + /* Initialize the exterior. */ Exterior.Prepare(TopGrid.GridData); - float InflowValue[5], Dummy[5]; + float InflowValue[field_counter], Dummy[field_counter]; + for (int j = 0; j < field_counter; j++){ + InflowValue[j] = 0.0; + Dummy[j] = 0.0; + } InflowValue[0] = ShockInABoxDensity[0]; - InflowValue[1] = ShockInABoxPressure[0]/(Gamma-1.0)/ShockInABoxDensity[0] + InflowValue[1] = ShockInABoxVelocity[0]; + InflowValue[4] = ShockInABoxPressure[0]/(Gamma-1.0)/ShockInABoxDensity[0] + 0.5*POW(ShockInABoxVelocity[0], 2); - InflowValue[2] = ShockInABoxVelocity[0]; - InflowValue[3] = 0.0; - InflowValue[4] = 0.0; if (Exterior.InitializeExternalBoundaryFace(0, inflow, outflow, InflowValue, Dummy) == FAIL) { @@ -198,21 +222,6 @@ int ShockInABoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, Exterior.InitializeExternalBoundaryFace(2, reflecting, reflecting, Dummy, Dummy); - - /* set up field names and units */ - - DataLabel[0] = DensName; - DataLabel[1] = TEName; - DataLabel[2] = Vel1Name; - DataLabel[3] = Vel2Name; - DataLabel[4] = Vel3Name; - - DataUnits[0] = NULL; - DataUnits[1] = NULL; - DataUnits[2] = NULL; - DataUnits[3] = NULL; - DataUnits[4] = NULL; - /* Write parameters to parameter output file */ if (MyProcessorNumber == ROOT_PROCESSOR) { @@ -221,16 +230,16 @@ int ShockInABoxInitialize(FILE *fptr, FILE *Outfptr, HierarchyEntry &TopGrid, fprintf(Outfptr, "ShockInABoxBoundary = %"GOUTSYM"\n\n", ShockInABoxBoundary); - fprintf(Outfptr, "ShockInABoxLeftDensity = %"FSYM"\n", ShockInABoxDensity[0]); - fprintf(Outfptr, "ShockInABoxLeftPressure = %"FSYM"\n", + fprintf(Outfptr, "ShockInABoxLeftDensity = %"ESYM"\n", ShockInABoxDensity[0]); + fprintf(Outfptr, "ShockInABoxLeftPressure = %"ESYM"\n", ShockInABoxPressure[0]); - fprintf(Outfptr, "ShockInABoxLeftVelocity = %"FSYM"\n\n", + fprintf(Outfptr, "ShockInABoxLeftVelocity = %"ESYM"\n\n", ShockInABoxVelocity[0]); - fprintf(Outfptr, "ShockInABoxRightDensity = %"FSYM"\n", ShockInABoxDensity[1]); - fprintf(Outfptr, "ShockInABoxRightPressure = %"FSYM"\n", + fprintf(Outfptr, "ShockInABoxRightDensity = %"ESYM"\n", ShockInABoxDensity[1]); + fprintf(Outfptr, "ShockInABoxRightPressure = %"ESYM"\n", ShockInABoxPressure[1]); - fprintf(Outfptr, "ShockInABoxRightVelocity = %"FSYM"\n\n", + fprintf(Outfptr, "ShockInABoxRightVelocity = %"ESYM"\n\n", ShockInABoxVelocity[1]); } diff --git a/src/enzo/hydro_rk/EvolveLevel_RK2.C b/src/enzo/hydro_rk/EvolveLevel_RK2.C index b6d76b98c..bf737ae15 100644 --- a/src/enzo/hydro_rk/EvolveLevel_RK2.C +++ b/src/enzo/hydro_rk/EvolveLevel_RK2.C @@ -484,7 +484,6 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], } // ENDIF UseHydro - /* Solve the cooling and species rate equations. */ Grids[grid1]->GridData->MultiSpeciesHandler(); @@ -508,10 +507,12 @@ int EvolveLevel_RK2(TopGridData *MetaData, LevelHierarchyEntry *LevelArray[], /* Gravity: clean up AccelerationField. */ +#ifndef SAB if ((level != MaximumGravityRefinementLevel || MaximumGravityRefinementLevel == MaximumRefinementLevel)) Grids[grid1]->GridData->DeleteAccelerationField(); +#endif //!SAB Grids[grid1]->GridData->DeleteParticleAcceleration(); if (UseFloor) diff --git a/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C b/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C index eaebe8370..b8126c112 100644 --- a/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C +++ b/src/enzo/hydro_rk/Grid_MHDRK2_2ndStep.C @@ -122,6 +122,15 @@ int grid::MHDRK2_2ndStep(fluxes *SubgridFluxes[], delete [] dU[field]; } + /* If we're supposed to be outputting on Density, we need to update + the current maximum value of that Density. */ + + if(OutputOnDensity == 1){ + int DensNum = FindField(Density, FieldType, NumberOfBaryonFields); + for(int i = 0; i < size; i++) + CurrentMaximumDensity = max(BaryonField[DensNum][i], CurrentMaximumDensity); + } + return SUCCESS; } diff --git a/src/enzo/hydro_rk/Grid_RungeKutta2_2ndStep.C b/src/enzo/hydro_rk/Grid_RungeKutta2_2ndStep.C index ebb262370..204a2e2ad 100644 --- a/src/enzo/hydro_rk/Grid_RungeKutta2_2ndStep.C +++ b/src/enzo/hydro_rk/Grid_RungeKutta2_2ndStep.C @@ -163,6 +163,15 @@ int grid::RungeKutta2_2ndStep(fluxes *SubgridFluxes[], } // PerformanceTimers[1] += ReturnWallTime() - time1; + + /* If we're supposed to be outputting on Density, we need to update + the current maximum value of that Density. */ + + if(OutputOnDensity == 1){ + int DensNum = FindField(Density, FieldType, NumberOfBaryonFields); + for(int i = 0; i < size; i++) + CurrentMaximumDensity = max(BaryonField[DensNum][i], CurrentMaximumDensity); + } return SUCCESS; diff --git a/src/performance_tools/performance_tools.py b/src/performance_tools/performance_tools.py index a1e6fd58f..cb1b15520 100755 --- a/src/performance_tools/performance_tools.py +++ b/src/performance_tools/performance_tools.py @@ -59,14 +59,24 @@ def preserve_extrema(extrema, xdata, ydata): if extrema[4] == 0: minx = np.min(xdata) maxx = np.max(xdata) - miny = np.min(ydata[np.nonzero(ydata)]) - maxy = np.max(ydata[np.nonzero(ydata)]) + nonzero = ydata[np.nonzero(ydata)] + if len(nonzero) == 0: + miny = 1e-1 ### Reasonable defaults given that no data to plot + maxy = 1 + else: + miny = np.min(nonzero) + maxy = np.max(nonzero) ### Otherwise, preserve the existing extrema else: minx = np.min([extrema[0], np.min(xdata)]) maxx = np.max([extrema[1], np.max(xdata)]) - miny = np.min([extrema[2], np.min(ydata[np.nonzero(ydata)])]) - maxy = np.max([extrema[3], np.max(ydata[np.nonzero(ydata)])]) + nonzero = ydata[np.nonzero(ydata)] + if len(nonzero) == 0: + miny = extrema[2] + maxy = extrema[3] + else: + miny = np.min([extrema[2], np.min(nonzero)]) + maxy = np.max([extrema[3], np.max(nonzero)]) return [minx,maxx,miny,maxy,1] def smooth(x, window_len=11, window='hanning'):