commit/enzo-dev: 8 new changesets

[commits...@bitbucket.org]

8 new commits in enzo-dev:

https://bitbucket.org/enzo/enzo-dev/commits/4ac1f78e7435/
Changeset: 4ac1f78e7435
Branch: week-of-code
User: jwise77
Date: 2019-02-26 19:24:33+00:00
Summary: Adding support for DM particles in the ZeldovichPancake problem
Affected #: 1 file

diff -r 49746a73ee3e -r 4ac1f78e7435 src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
--- a/src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
+++ b/src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
@@ -45,7 +45,7 @@
/* declarations */

float Amplitude, AmplitudeVel, kx, xLagrange, xEulerian, xEulerianOld;
- int dim, field, i, index;
+ int dim, field, i, j, k, index;
const float Pi = 3.14159;

/* error check */
@@ -53,9 +53,6 @@
if (ZeldovichPancakeDirection < 0 || ZeldovichPancakeDirection >= GridRank) {
ENZO_FAIL("ZeldovichPancakeDirection is improperly set.\n");
}
- if (ZeldovichPancakeOmegaCDMNow != 0) {
- ENZO_FAIL("Dark matter not yet supported.\n");
- }

/* create fields */

@@ -103,9 +100,11 @@

/* Determine the size of the fields. */

- int size = 1;
- for (dim = 0; dim < GridRank; dim++)
+ int size = 1, activesize = 1;
+ for (dim = 0; dim < GridRank; dim++) {
size *= GridDimension[dim];
+ activesize *= GridEndIndex[dim] - GridStartIndex[dim] + 1;
+ }

/* Allocate space for the fields. */

@@ -119,7 +118,14 @@

int NumberOfZones = GridEndIndex[ZeldovichPancakeDirection] -
GridStartIndex[ZeldovichPancakeDirection] + 1;
-
+
+ /* Set up particles */
+
+ if (ZeldovichPancakeOmegaCDMNow > 0) {
+ this->AllocateNewParticles(NumberOfZones);
+ this->NumberOfParticles = NumberOfZones;
+ }
+
/* Compute the amplitude of perturbations. */

Amplitude = (1+ZeldovichPancakeCollapseRedshift) / (1+InitialRedshift);
@@ -196,8 +202,8 @@
if (HydroMethod == MHD_RK) {
BaryonField[iPhi ][i] = 0.0;
}
- }
-
+
+ } // ENDFOR zones

/* set transverse velocities (i.e. erase any incorrectly set velocities). */

@@ -205,6 +211,38 @@
if (field != ZeldovichPancakeDirection+vel)
for (i = 0; i < size; i++)
BaryonField[field][i] = 0.0;
+
+ /* Set up particles if needed */
+
+ if (ZeldovichPancakeOmegaCDMNow > 0) {
+
+ int ipart = 0;
+ FLOAT pos0[MAX_DIMENSION], dr[MAX_DIMENSION];
+ for (dim = 0 ; dim < GridRank; dim++)
+ dr[dim] = (GridRank > dim) ? CellWidth[dim][0] :
+ DomainRightEdge[dim] - DomainLeftEdge[dim];
+
+ for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
+ pos0[2] = GridLeftEdge[2] + (0.5+k-GridStartIndex[2]) * dr[2];
+ for (j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
+ pos0[1] = GridLeftEdge[1] + (0.5+j-GridStartIndex[1]) * dr[1];
+ index = GRIDINDEX_NOGHOST(GridStartIndex[0], j, k);
+ for (i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, index++) {
+ pos0[0] = GridLeftEdge[0] + (0.5+i-GridStartIndex[0]) * dr[0];
+ ParticleMass[ipart] = ZeldovichPancakeOmegaCDMNow;
+ ParticleNumber[ipart] = ipart;
+ ParticleType[ipart] = PARTICLE_TYPE_DARK_MATTER;
+ for (dim = 0; dim < GridRank; dim++) {
+ ParticlePosition[dim][ipart] = pos0[dim] +
+ BaryonField[vel+dim][index] * InitialTimeInCodeUnits;
+ ParticleVelocity[dim][ipart] = BaryonField[vel+dim][index];
+ }
+ ipart++;
+ } // ENDFOR i
+ } // ENDFOR j
+ } // ENDFOR k
+
+ } // ENDIF particles

return SUCCESS;
}


https://bitbucket.org/enzo/enzo-dev/commits/ecced5aa9f43/
Changeset: ecced5aa9f43
Branch: week-of-code
User: jwise77
Date: 2019-02-27 16:01:35+00:00
Summary: Fixing bugs for particle movement in 1D and 2D simulations
Affected #: 1 file

diff -r 4ac1f78e7435 -r ecced5aa9f43 src/enzo/Grid_TransferSubgridParticles.C
--- a/src/enzo/Grid_TransferSubgridParticles.C
+++ b/src/enzo/Grid_TransferSubgridParticles.C
@@ -65,7 +65,7 @@

/* Set boundaries (with and without ghost zones) */

- int StartIndex[] = {1,1,1}, EndIndex[] = {1,1,1};
+ int StartIndex[] = {0,0,0}, EndIndex[] = {0,0,0};
if (IncludeGhostZones)
for (dim = 0; dim < GridRank; dim++) {
StartIndex[dim] = 0;
@@ -87,9 +87,9 @@
/* Compute index of particle position. */

i0 = int((ParticlePosition[0][i] - CellLeftEdge[0][0])/CellWidth[0][0]);
- if (GridRank > 0)
+ if (GridRank > 1)
j0 = int((ParticlePosition[1][i] - CellLeftEdge[1][0])/CellWidth[1][0]);
- if (GridRank > 1)
+ if (GridRank > 2)
k0 = int((ParticlePosition[2][i] - CellLeftEdge[2][0])/CellWidth[2][0]);

i0 = max(min(EndIndex[0], i0), StartIndex[0]);


https://bitbucket.org/enzo/enzo-dev/commits/48fbcd4ec204/
Changeset: 48fbcd4ec204
Branch: week-of-code
User: jwise77
Date: 2019-02-27 16:05:01+00:00
Summary: Merging
Affected #: 209 files

diff -r ecced5aa9f43 -r 48fbcd4ec204 README
--- a/README
+++ b/README
@@ -51,7 +51,7 @@
* Philipp Edelmann pede...@mpa-garching.mpg.de
* Andrew Emerick aemer...@gmail.com
* Nathan Goldbaum ngol...@ucsc.edu
- * Phillipp Grete pgr...@astro.physik.uni-goettingen.de
+ * Philipp Grete gr...@pa.msu.edu
* John Forbes jcfo...@ucsc.edu
* Oliver Hahn ha...@phys.ethz.ch
* Robert Harkness hark...@sdsc.edu

diff -r ecced5aa9f43 -r 48fbcd4ec204 doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -122,6 +122,8 @@

* `Conduction`_

+* `Subgrid-scale (SGS) turbulence model`_
+
* `Inline Analysis`_

* `Inline Halo Finding`_
@@ -3194,6 +3196,140 @@
``ConductionDynamicRebuildMinLevel`` (external)
The minimum level on which the dynamic hierarcy rebuild is performed. Default: 0.

+
+.. _sgs_parameters:
+
+Subgrid-scale (SGS) turbulence model
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The following parameter allow the use of an SGS turbulence model in
+Enzo, see test problem
+``run/MHD/3D/StochasticForcing/StochasticForcing.enzo``.
+
+It is recommended to not arbitrarily mix model terms, but either
+stick to one model family (nonlinear, dissipative, or scale-similarity)
+or conduct additional *a priori* test first.
+
+Best fit model coefficients based on *a priori* testing of compressible
+MHD turbulence for a wide range of sonic Mach numbers (0.2 to 20) can
+be found in Table II in `Grete et al. (2016) Physics of Plasmas 23 062317
+<https://dx.doi.org/10.1063/1.4954304>`_,
+where all models are presented in more detail.
+
+Overall, the nonlinear model (effectively parameter free)
+with an explicit 3-point stencil showed the
+best performance in decaying MHD test problem, see
+`Grete et al. (2017) Phys. Rev. E. 95 033206
+<https://dx.doi.org/10.1103/PhysRevE.95.033206>`_.
+
+
+``UseSGSModel`` (external)
+ This parameter generally turns the SGS machinery on (even though
+ no SGS term is added by default as every term needs a coefficient,
+ see below).
+ 1: Turn on. Default: 0
+
+Explicit filtering
+^^^^^^^^^^^^^^^^^^
+
+All SGS models rely on the notion that they are calculated based on
+filtered/resolved quantities.
+The spatial discretization itself acts as one filter.
+However, in shock-capturing schemes it is questionable how "resolved"
+quantities are at the grid-scale.
+The following three variables enable the explicit filtering of the grid-scale
+quantites as they are used in the SGS terms.
+
+See Table 1 in `Grete et al. (2017) Phys. Rev. E. 95 033206
+<https://dx.doi.org/10.1103/PhysRevE.95.033206>`_ for coefficients
+of a discrete box filter.
+The recommended values correspond to a discrete representation of a box filter
+on a 3-point stencil.
+
+
+``SGSFilterWidth`` (external)
+ Width (in units of cell widths) of the discrete filter.
+ Default: 0;
+ Recommended: 2.711;
+
+``SGSFilterStencil`` (external)
+ Discrete width of filter stencil in numbers of cells.
+ Default: 0;
+ Recommended: 3;
+ Maximum: 7;
+
+``SGSFilterWeights`` (external)
+ Symmetic filter weights that are used in the stencil. List of four floats.
+ First number corresponds to weight of central point X_i,
+ second number corresponds to weight of points X_i+1 and X_i-1, and so on.
+ Default: 0. 0. 0. 0.;
+ Recommended: 0.40150 0.29925 0.00000 0.0;
+
+Nonlinear model
+^^^^^^^^^^^^^^^
+
+``SGScoeffNLu`` (external)
+ Coefficient for nonlinear Reynolds stress model.
+ Default: 0;
+ Recommended: 1;
+
+``SGScoeffNLb`` (external)
+ Coefficient for nonlinear Maxwell stress (only MHD).
+ Default: 0;
+ Recommended: 1;
+
+``SGScoeffNLemfCompr`` (external)
+ Coefficient for nonlinear compressive EMF model (only MHD).
+ Default: 0;
+ Recommended: 1;
+
+Dissipative model
+^^^^^^^^^^^^^^^^^
+
+``SGScoeffEVStarEnS2Star`` (external)
+ Coefficient for traceless eddy-viscosity Reynolds stress model
+ (scaled by realizability condition in the kinetic SGS energy).
+ Default: 0;
+ Recommended: 0.01;
+
+``SGScoeffEnS2StarTrace`` (external)
+ Coefficient for the trace of the eddy-viscosity Reynolds stress model,
+ i.e., the kinetic SGS energy (derived from realizability condition).
+ Default: 0;
+ Recommended: 0.08;
+
+``SGScoeffERS2M2Star`` (external)
+ Coefficient for eddy-resistivity EMF model (only MHD; scaled by
+ realizable SGS energies)
+ Default: 0;
+ Recommended: 0.012;
+
+``SGScoeffERS2J2`` (external)
+ Coefficient for eddy-resistivity EMF model (only MHD; scaled by
+ Smagorinsky SGS energies)
+ Default: 0;
+
+Scale-similarity model
+^^^^^^^^^^^^^^^^^^^^^^
+
+``SGScoeffSSu`` (external)
+ Coefficient for scale-similarity Reynolds stress model.
+ Default: 0;
+ Recommended: 0.67;
+
+``SGScoeffSSb`` (external)
+ Coefficient for scale-similarity Maxwell stress (only MHD).
+ Default: 0;
+ Recommended: 0.9;
+
+``SGScoeffNLemfCompr`` (external)
+ Coefficient for scale-similarity EMF model (only MHD).
+ Default: 0;
+ Recommended: 0.89;
+
+
+
+
.. _other_parameters:

Other Parameters

diff -r ecced5aa9f43 -r 48fbcd4ec204 run/MHD/3D/StochasticForcing/StochasticForcing.enzo
--- a/run/MHD/3D/StochasticForcing/StochasticForcing.enzo
+++ b/run/MHD/3D/StochasticForcing/StochasticForcing.enzo
@@ -1,10 +1,13 @@
#
-# PROBLEM DEFINITION FILE: MHD/HD turbulence problem with stochastic forcing
+# PROBLEM DEFINITION FILE: MHD/HD turbulence problem with stochastic forcing
+# with subgrid-scale (SGS) turbulence model
# Philipp Grete 2014
#
# Typical "turbulence-in-a-box" problem with non-static driving field.
# For details on stochastic forcing, see Schmidt et al. 2009 A&A 494, 127-145
# http://dx.doi.org/10.1051/0004-6361:200809967
+# For details on the SGS model, see Grete et al. (2017) Phys. Rev. E. 95 033206
+# https://dx.doi.org/10.1103/PhysRevE.95.033206
#
# Works/properly tested only on 3D uniform grids with MUSCL type solvers and MHDCT at this point.
# For hydro use HydroMethod 3
@@ -25,7 +28,7 @@
#
# set I/O and stop/start parameters
#
-StopTime = 0.86
+StopTime = 0.516
dtDataDump = 0.0215
DataDumpName = data
ParallelRootGridIO = 1
@@ -64,4 +67,23 @@
DrivenFlowPressure = 1.0 # initial uniform pressure
DrivenFlowMagField = 3.1622777 # initial uniform field (x direction)

+UseSGSModel = 1 # use SGS model

+SGSFilterWidth = 2.7110
+SGSFilterStencil = 3
+SGSFilterWeights = 0.40150 0.29925 0.00000 0.0
+
+SGScoeffNLemfCompr = 1.0 # compr. unscaled nonliner model EMF
+SGScoeffNLu = 1.0 # compr. unscaled nonlinear model tauU
+SGScoeffNLb = 0.0 # compr. unscaled nonlinear model tauB
+#SGScoeffERS2J2 = 0.0 # eddy resistivity EMF model scaled by Smag.
+#energies
+##SGScoeffERS2M2Star = 0.012 # eddy resistivity EMF model scaled by
+#realiz. energies
+##SGScoeffEVStarEnS2Star = 0.01 # eddy viscosity tauUstar scaled by
+#realiz. energies
+##SGScoeffEnS2StarTrace = 0.08 # tauUtrace - SGS energy coeff from
+#realiz. energies
+##SGScoeffSSu = 0.67
+##SGScoeffSSb = 0.9
+##SGScoeffSSemf = 0.89

diff -r ecced5aa9f43 -r 48fbcd4ec204 run/MHD/3D/StochasticForcing/StochasticForcing.enzotest
--- a/run/MHD/3D/StochasticForcing/StochasticForcing.enzotest
+++ b/run/MHD/3D/StochasticForcing/StochasticForcing.enzotest
@@ -9,5 +9,5 @@
dimensionality = 3
max_time_minutes = 20
fullsuite = True
-pushsuite = False
+pushsuite = True
quicksuite = False

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/AMRH5writer.C
--- a/src/enzo/AMRH5writer.C
+++ b/src/enzo/AMRH5writer.C
@@ -24,6 +24,7 @@
#include "typedefs.h"
#include "AMRH5writer.h"
#include "StarParticleData.h"
+#include "phys_constants.h"

AMRHDF5Writer::AMRHDF5Writer() :
h5DataType(H5T_NATIVE_FLOAT),
@@ -151,14 +152,14 @@
error = true;
}

- // The 3.14159 header is to test the endian since this is a binary
+ // The pi header is to test the endian since this is a binary
// file. The root cell spacing is necessary to calculate the dx for
// all other grids.
- float myPi = 3.14159;
char nFields8bit = (char) NumberOfAllFields;
char stag8bit = (char) stag;
int field;
- fwrite(&myPi, sizeof(float), 1, index_file);
+ double test_pi = pi;
+ fwrite(&test_pi, sizeof(float), 1, index_file);
fwrite(&root_dx, sizeof(float), 1, index_file);
fwrite(&nFields8bit, sizeof(char), 1, index_file);
fwrite(&stag8bit, sizeof(char), 1, index_file);

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/AdiabaticExpansionInitialize.C
--- a/src/enzo/AdiabaticExpansionInitialize.C
+++ b/src/enzo/AdiabaticExpansionInitialize.C
@@ -29,6 +29,7 @@
#include "Hierarchy.h"
#include "TopGridData.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"

/* Set the mean molecular mass as in Grid_ComputeTemperatureField.C */

@@ -135,7 +136,7 @@
ENZO_FAIL("Error in GetUnits.");
}
PressureUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
CRUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);

for (int dim = 0; dim < MAX_DIMENSION; dim++)

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/AdjustMustRefineParticlesRefineToLevel.C
--- a/src/enzo/AdjustMustRefineParticlesRefineToLevel.C
+++ b/src/enzo/AdjustMustRefineParticlesRefineToLevel.C
@@ -41,6 +41,7 @@
#include "TopGridData.h"
#include "LevelHierarchy.h"
#include "Star.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -61,7 +62,7 @@
ENZO_FAIL("Your setup of assymmetrical topgrid is never tested.\n");
}

- const double pc = 3.086e18, ln2 = 0.69314718;
+ const double ln2 = 0.69314718;
int MustRefineParticlesRefineToLevel_prev = MustRefineParticlesRefineToLevel;

/* Compute Units. */
@@ -80,7 +81,7 @@

MustRefineParticlesRefineToLevel =
(int) floor( log( LengthUnits / MetaData->TopGridDims[0]
- / MustRefineParticlesRefineToLevelAutoAdjust / pc) / ln2);
+ / MustRefineParticlesRefineToLevelAutoAdjust / pc_cm) / ln2);

if (MyProcessorNumber == ROOT_PROCESSOR) {


diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/AnalysisBaseClass.C
--- a/src/enzo/AnalysisBaseClass.C
+++ b/src/enzo/AnalysisBaseClass.C
@@ -15,6 +15,7 @@
#include "TopGridData.h"
#include "CosmologyParameters.h"
#include "AnalysisBaseClass.h"
+#include "phys_constants.h"

std::map<HierarchyEntry *, int> OriginalGridID;
#ifdef USE_HDF5_GROUPS
@@ -31,16 +32,8 @@
HierarchyEntry *topgrid,
int level,
FLOAT *Left,
- FLOAT *Right):
-#ifdef USE_BAD_VALUES
- G(6.67e-8),
- MPC_CM(3.0856e+24),
- MSOLAR_G(1.989e+33)
-#else
- G(6.6742e-8),
- MPC_CM(3.0856776e+24),
- MSOLAR_G(1.989e+33)
-#endif
+ FLOAT *Right)
+
{

int i;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/AnalysisBaseClass.h
--- a/src/enzo/AnalysisBaseClass.h
+++ b/src/enzo/AnalysisBaseClass.h
@@ -33,10 +33,6 @@

protected:

- const float G;
- const float MPC_CM;
- const float MSOLAR_G;
-
TopGridData *MetaData;
HierarchyEntry *TopGrid;
int MaximumAnalysisLevel;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/CalculateSubtractionParameters.C
--- a/src/enzo/CalculateSubtractionParameters.C
+++ b/src/enzo/CalculateSubtractionParameters.C
@@ -29,6 +29,7 @@
#include "TopGridData.h"
#include "LevelHierarchy.h"
#include "CommunicationUtilities.h"
+#include "phys_constants.h"

int CommunicationBroadcastValue(int *Value, int BroadcastProcessor);

@@ -43,8 +44,7 @@
float &Radius, double &Subtraction)
{

- const double pc = 3.086e18, Msun = 1.989e33, Grav = 6.673e-8, yr = 3.1557e7, Myr = 3.1557e13,
- k_b = 1.38e-16, m_h = 1.673e-24, c = 3.0e10, sigma_T = 6.65e-25, h=0.70;
+ const double h=0.70;

float mdot, AccretedMass, SafetyFactor;
float MassEnclosed = 0, Metallicity = 0, ColdGasMass = 0, OneOverRSquaredSum, AvgVelocity[MAX_DIMENSION];
@@ -129,7 +129,7 @@
density = star_CurrentGrid->ReturnBaryonField(DensNum)[index];

if (MultiSpecies == 0) {
- number_density = density * DensityUnits / (Mu * m_h);
+ number_density = density * DensityUnits / (Mu * mh);
mu = Mu;
} else {
number_density =
@@ -157,7 +157,7 @@
because only integer can be broadcasted at the moment;
anyway, this shouldn't make a huge difference overall. */

- c_s = (int)(sqrt(Gamma * k_b * temperature[index] / (mu * m_h)));
+ c_s = (int)(sqrt(Gamma * kboltz * temperature[index] / (mu * mh)));

fprintf(stdout, "index = %d, temp = %g, mu = %g, density = %g, number_density = %g, c_s = %d\n",
index, temperature[index], mu, density, number_density, c_s);
@@ -182,7 +182,7 @@
/* When calculating Bondi radius, let's not bother to get accurate temperature and mu
just use the MBHAccretionFixedTemperature and default mu */

- c_s = (int)(sqrt(Gamma * k_b * MBHAccretionFixedTemperature / (Mu * m_h)));
+ c_s = (int)(sqrt(Gamma * kboltz * MBHAccretionFixedTemperature / (Mu * mh)));
#endif


@@ -192,14 +192,14 @@

old_mass = (float)(star_mass);
mdot = isnan(star_last_accretion_rate) ? 0.0 : star_last_accretion_rate;
- AccretedMass = mdot * dtForThisStar * TimeUnits; //in Msun
+ AccretedMass = mdot * dtForThisStar * TimeUnits; //in SolarMass


if (MBHAccretionRadius > 0)
- Radius = MBHAccretionRadius * pc / LengthUnits;
+ Radius = MBHAccretionRadius * pc_cm / LengthUnits;
else {
SafetyFactor = -MBHAccretionRadius;
- Radius = SafetyFactor * 2.0 * Grav * old_mass * Msun / (c_s * c_s) / LengthUnits;
+ Radius = SafetyFactor * 2.0 * GravConst * old_mass * SolarMass / (c_s * c_s) / LengthUnits;
}

Radius = min(max(Radius, 4*StarLevelCellWidth), 100*StarLevelCellWidth);
@@ -284,8 +284,8 @@
here Subtraction is the density that should be subtracted.
if one wants to use this, change Grid_SAMFS.C accordingly. */

- Subtraction = -AccretedMass * Msun /
- (4*M_PI/3.0 * pow(Radius*LengthUnits, 3)) / DensityUnits;
+ Subtraction = -AccretedMass * SolarMass /
+ (4*pi/3.0 * pow(Radius*LengthUnits, 3)) / DensityUnits;
#endif



diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/ClusterSMBHSumGasMass.C
--- a/src/enzo/ClusterSMBHSumGasMass.C
+++ b/src/enzo/ClusterSMBHSumGasMass.C
@@ -84,7 +84,7 @@
float ColdGasMassMsun=ClusterSMBHColdGasMass*MassUnits/SolarMass;
if (ClusterSMBHCalculateGasMass >1 ){ //2-calculate&remove,3-calculate&remove&re-orient,4-Bondi
if (ClusterSMBHCalculateGasMass == 3){
- ClusterSMBHJetDim = floor(Time*TimeUnits/(1.0e6*3.1557e7*ClusterSMBHJetPrecessionPeriod)); //ClusterSMBHJetPrecessionPeriod is now the Dim changing period.
+ ClusterSMBHJetDim = floor(Time*TimeUnits/(1.0e6*yr_s*ClusterSMBHJetPrecessionPeriod)); //ClusterSMBHJetPrecessionPeriod is now the Dim changing period.
}
if (ColdGasMassMsun > 0.000001) {
ClusterSMBHFeedbackSwitch = TRUE;
@@ -94,7 +94,7 @@
}
else {
ClusterSMBHJetMdot = (ColdGasMassMsun/(ClusterSMBHAccretionTime*1e6))/2.0; // AccretionTime from Myr to yr; reset Mdot, still in Msun/yr. Devide it by 2 because Mdot is for only one jet.
- ClusterSMBHJetEdot = (ClusterSMBHAccretionEpsilon*ClusterSMBHJetMdot * SolarMass/3.1557e7) * POW(clight,2)/1.0e44; //for one jet
+ ClusterSMBHJetEdot = (ClusterSMBHAccretionEpsilon*ClusterSMBHJetMdot * SolarMass/yr_s) * POW(clight,2)/1.0e44; //for one jet
}
}
else
@@ -108,7 +108,7 @@
ClusterSMBHFeedbackSwitch = (ColdGasMassMsun <= ClusterSMBHEnoughColdGas && LastClusterSMBHFeedbackSwitch == FALSE) ? FALSE : TRUE;

if (LastClusterSMBHFeedbackSwitch == FALSE && ClusterSMBHFeedbackSwitch == TRUE) {
- ClusterSMBHStartTime = Time + ClusterSMBHTramp*1.0e6*3.1557e7/TimeUnits;
+ ClusterSMBHStartTime = Time + ClusterSMBHTramp*1.0e6*yr_s/TimeUnits;
if (ClusterSMBHJetPrecessionPeriod < 0.00001 & ClusterSMBHJetPrecessionPeriod > -0.00001) //if precession off (set to 0), change the angle of the jets
ClusterSMBHJetAnglePhi += 0.5;
if (ClusterSMBHJetPrecessionPeriod < -0.00001) //if precession negative, change the jet dimension

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/CommunicationMergeStarParticle.C
--- a/src/enzo/CommunicationMergeStarParticle.C
+++ b/src/enzo/CommunicationMergeStarParticle.C
@@ -33,6 +33,7 @@
#include "TopGridData.h"
#include "Hierarchy.h"
#include "LevelHierarchy.h"
+#include "phys_constants.h"

/* function prototypes */

@@ -148,8 +149,7 @@
GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits,
&VelocityUnits, &MassUnits, 1.0);

- double Msun = 1.989e33;
- if(UsePhysicalUnit) MassUnits = DensityUnits*pow(LengthUnits,3)/Msun;
+ if(UsePhysicalUnit) MassUnits = DensityUnits*pow(LengthUnits,3)/SolarMass;

/*
if (debug) {

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/ConductionBubbleInitialize.C
--- a/src/enzo/ConductionBubbleInitialize.C
+++ b/src/enzo/ConductionBubbleInitialize.C
@@ -29,6 +29,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -108,7 +109,7 @@
total energy (that's not done automatically, since we don't know for sure
what the user is going to do). */
if (HydroMethod == MHD_RK){
- float MagneticUnits = sqrt(DensityUnits*4.0*M_PI)*VelocityUnits;
+ float MagneticUnits = sqrt(DensityUnits*4.0*pi)*VelocityUnits;
ConductionBubbleInitialUniformBField[0] /= MagneticUnits;
ConductionBubbleInitialUniformBField[1] /= MagneticUnits;
ConductionBubbleInitialUniformBField[2] /= MagneticUnits;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/ConductionCloudInitialize.C
--- a/src/enzo/ConductionCloudInitialize.C
+++ b/src/enzo/ConductionCloudInitialize.C
@@ -25,6 +25,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -123,11 +124,11 @@
return FAIL;
}

- float Boltzmann = 1.38e-16, mu = 0.6, mh=1.67e-24;
+ float mu = 0.6;

if(ConductionCloudTemperature > 1.0){

- ConductionCloudTotalEnergy = (Boltzmann*ConductionCloudTemperature)/((Gamma - 1.0)*mu*mh);
+ ConductionCloudTotalEnergy = (kboltz*ConductionCloudTemperature)/((Gamma - 1.0)*mu*mh);
ConductionCloudTotalEnergy /= (VelocityUnits*VelocityUnits);
printf("ConductionCloudTotalEnergy is %e and ConductionCloudTemperature is %e\n\n",ConductionCloudTotalEnergy, ConductionCloudTemperature);
fflush(stdout);

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/ConductionTestInitialize.C
--- a/src/enzo/ConductionTestInitialize.C
+++ b/src/enzo/ConductionTestInitialize.C
@@ -27,6 +27,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -105,12 +106,12 @@
return FAIL;
}

- float Boltzmann = 1.38e-16, mu = 0.6, mh=1.67e-24;
+ float mu = 0.6;

// User can set
if(ConductionTestTemperature > 1.0){

- ConductionTestTotalEnergy = (Boltzmann*ConductionTestTemperature)/((Gamma - 1.0)*mu*mh);
+ ConductionTestTotalEnergy = (kboltz*ConductionTestTemperature)/((Gamma - 1.0)*mu*mh);
ConductionTestTotalEnergy /= (VelocityUnits*VelocityUnits);
ConductionTestGasEnergy = ConductionTestTotalEnergy;
printf("ConductionTestTotalEnergy is %e and ConductionTestTemperature is %e\n\n",ConductionTestTotalEnergy, ConductionTestTemperature);
@@ -118,7 +119,7 @@
}

if (HydroMethod == MHD_RK){
- float MagneticUnits = sqrt(DensityUnits*4.0*M_PI)*VelocityUnits;
+ float MagneticUnits = sqrt(DensityUnits*4.0*pi)*VelocityUnits;
ConductionTestInitialUniformBField[0] /= MagneticUnits;
ConductionTestInitialUniformBField[1] /= MagneticUnits;
ConductionTestInitialUniformBField[2] /= MagneticUnits;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/CosmoIonizationInitialize.C
--- a/src/enzo/CosmoIonizationInitialize.C
+++ b/src/enzo/CosmoIonizationInitialize.C
@@ -29,6 +29,7 @@
#include "Hierarchy.h"
#include "LevelHierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"

// default constants
#define DEFAULT_MU 0.6 // mean molecular mass
@@ -134,8 +135,7 @@

// convert input temperature to internal energy
RadHydroTemperature = max(RadHydroTemperature,MIN_TEMP); // enforce minimum
- float mp = 1.67262171e-24; // proton mass [g]
- float kb = 1.3806504e-16; // boltzmann constant [erg/K]
+
float nH, HI, HII, nHe, HeI, HeII, HeIII, ne, num_dens, mu;
if (RadHydroChemistry == 1) {
HI = 1.0 - RadHydroInitialFractionHII;
@@ -157,7 +157,7 @@
mu = 1.0/num_dens;
}
// compute the internal energy
- float RadHydroIEnergy = kb*RadHydroTemperature/mu/mp/(Gamma-1.0);
+ float RadHydroIEnergy = kboltz*RadHydroTemperature/mu/mh/(Gamma-1.0);

// set up the grid(s) on this level
HierarchyEntry *Temp = &TopGrid;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/CosmologyGetUnits.C
--- a/src/enzo/CosmologyGetUnits.C
+++ b/src/enzo/CosmologyGetUnits.C
@@ -39,6 +39,7 @@
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"

/* function prototypes */

@@ -66,7 +67,7 @@
*DensityUnits = 1.8788e-29*OmegaMatterNow*POW(HubbleConstantNow,2)*
POW(1 + CurrentRedshift,3);

- *LengthUnits = 3.085678e24*ComovingBoxSize/HubbleConstantNow/
+ *LengthUnits = Mpc_cm*ComovingBoxSize/HubbleConstantNow/
(1 + CurrentRedshift);

*TemperatureUnits = 1.81723e6*POW(ComovingBoxSize,2)*OmegaMatterNow*

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/CosmologySimulationInitialize.C
--- a/src/enzo/CosmologySimulationInitialize.C
+++ b/src/enzo/CosmologySimulationInitialize.C
@@ -39,6 +39,7 @@
#include "CosmologyParameters.h"
#include "fortran.def"
#include "CommunicationUtilities.h"
+#include "phys_constants.h"

// Function prototypes

@@ -419,7 +420,7 @@
ENZO_FAIL("Error in GetUnits.");
}
PressureUnits = DensityUnits * VelocityUnits*VelocityUnits;
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);

for (int dim = 0; dim < MAX_DIMENSION; dim++) {
if (CosmologySimulationInitialUniformBField[dim] != 0.0 && HydroMethod != 4 && HydroMethod != 6)

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
--- a/src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
+++ b/src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
@@ -27,6 +27,7 @@
#include "Grid.h"
#include "MHD2DTestGlobalData.h"
#include "hydro_rk/EOS.h"
+#include "phys_constants.h"

/* function prototypes */

@@ -117,7 +118,7 @@
EOSSoundSpeed*EOSSoundSpeed*rho ;
EOS(pres, rho, eintl, h, cs, dpdrho, dpde, 0, 1); // compute eintl
ramp = 1./(1.+exp(-2/RampWidth*(y-0.5)));
- f = cos(2.*M_PI*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*M_PI*x*3);
+ f = cos(2.*pi*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*pi*x*3);
vx = f * (vxl + ramp*(vxu-vxl)) ;
vy = vyl + ramp*(vyu - vyl);
bx = (Bxl+ ramp*(Bxu-Bxl)) ;
@@ -152,7 +153,7 @@
EOSSoundSpeed*EOSSoundSpeed*rho ;
EOS(pres, rho, eintl, h, cs, dpdrho, dpde, 0, 1); // compute eintl
ramp = 1./(1.+exp(-2/RampWidth*(y-0.5)));
- f = cos(2.*M_PI*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*M_PI*x*3);
+ f = cos(2.*pi*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*pi*x*3);
vx = f * (vxl + ramp*(vxu-vxl)) ;
vy = vyl + ramp*(vyu - vyl);
bx = (Bxl+ ramp*(Bxu-Bxl)) ;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/FOF.C
--- a/src/enzo/FOF.C
+++ b/src/enzo/FOF.C
@@ -26,6 +26,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
+#include "phys_constants.h"

#include "Fluxes.h"
#include "GridList.h"
@@ -148,9 +149,9 @@

void set_units(FOFData &AllVars) {

- AllVars.UnitLength_in_cm = 3.085678e21;
- AllVars.UnitMass_in_g = 1.989e43;
- AllVars.UnitVelocity_in_cm_per_s = 1.0e5;
+ AllVars.UnitLength_in_cm = kpc_cm;
+ AllVars.UnitMass_in_g = SolarMass * 1.0E10;
+ AllVars.UnitVelocity_in_cm_per_s = km_cm;

AllVars.Theta = 0.8; /* opening angle for potential computation */
AllVars.DesDensityNgb = 32;
@@ -1395,7 +1396,7 @@
FILE *in, *out;
int i, Ngroups, *Npart_g;
float *CM, *Mass;
- double SolarMass = UnitMass_in_g / 1.989e33;
+ double Msun = UnitMass_in_g / SolarMass;
float RootBoxSize[3];

if (MyProcessorNumber == 0) {
@@ -1427,7 +1428,7 @@
"Halo", "Npart", "Mass [Msun]", "CM(x)", "CM(y)", "CM(z)");
for (i = 0; i < Ngroups; i++)
fprintf(out, "%8d %12d %12.6g %12.6f %12.6f %12.6f\n",
- i, Npart_g[i], Mass[i]*SolarMass,
+ i, Npart_g[i], Mass[i]*Msun,
CM[3*i+0] / RootBoxSize[0] + leftEdge[0],
CM[3*i+1] / RootBoxSize[1] + leftEdge[1],
CM[3*i+2] / RootBoxSize[2] + leftEdge[2]);

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/FOF_Initialize.C
--- a/src/enzo/FOF_Initialize.C
+++ b/src/enzo/FOF_Initialize.C
@@ -33,6 +33,7 @@
#include "Hierarchy.h"
#include "LevelHierarchy.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"

#include "FOF_allvars.h"
#include "FOF_nrutil.h"
@@ -82,13 +83,13 @@
D.Time = a / (1 + InitialRedshift);
}
else {
- D.BoxSize = LengthUnits / 3.086e21;
+ D.BoxSize = LengthUnits / kpc_cm;
D.Time = 1.0;
}

// Critical density in units of Msun / kpc^3
D.RhoCritical0 = 1.4775867e31 *
- ((3 * pow(100 * HubbleConstantNow / 3.086e19, 2)) / (8 * M_PI * GRAVITY));
+ ((3.0 * POW(100 * HubbleConstantNow / 3.086e19, 2)) / (8.0 * pi * GRAVITY));
//D.RhoCritical /= pow(D.Time, 3);

// Sometimes MassUnits is infinite (in cgs) when using single

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/FOF_ngbtree.C
--- a/src/enzo/FOF_ngbtree.C
+++ b/src/enzo/FOF_ngbtree.C
@@ -12,7 +12,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_ngbtree.h"

@@ -54,7 +54,7 @@
break;
}

- sr = th->len * pow((3.0/(4*M_PI)*SECFACTOR) * desngb /
+ sr = th->len * pow((3.0/(4*pi)*SECFACTOR) * desngb /
((float)(th->count)), 1.0/3);

} // ENDELSE

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/FOF_properties.C
--- a/src/enzo/FOF_properties.C
+++ b/src/enzo/FOF_properties.C
@@ -6,7 +6,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_proto.h"

@@ -134,7 +134,7 @@
// kpc^3, as is rho_crit.
e2 = D.OmegaLambda + D.Omega * pow(D.Time, -3.0);
// 1/a^3 factor converts comoving radius to proper.
- factor = 1e10 / (4*M_PI/3.0) * pow(D.Time, -3.0);
+ factor = 1e10 / (4*pi/3.0) * pow(D.Time, -3.0);
rho178 = 178.0 * D.RhoCritical0 * e2;
len4 = len/4;


diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/FSProb_ComputeOpacityLW.C
--- a/src/enzo/FSProb_ComputeOpacityLW.C
+++ b/src/enzo/FSProb_ComputeOpacityLW.C
@@ -21,12 +21,12 @@
************************************************************************/
#ifdef TRANSFER
#include "FSProb.h"
+#include "phys_constants.h"

int FSProb::ComputeOpacityLW(float *H2Density)
{

const float H2ISigma = 3.71e-18; // cm^-2
- const float mh = 1.673e-24;
double factor = H2ISigma * (DenUnits/mh);

int i, j, k, dim, size;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/FreeExpansionInitialize.C
--- a/src/enzo/FreeExpansionInitialize.C
+++ b/src/enzo/FreeExpansionInitialize.C
@@ -33,6 +33,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -128,7 +129,7 @@
&VelocityUnits, 0.0);

PressureUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
for (dim = 0; dim < MAX_DIMENSION; dim++)
FreeExpansionBField[dim] /= MagneticUnits;


diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/FreezeRateData.C
--- a/src/enzo/FreezeRateData.C
+++ b/src/enzo/FreezeRateData.C
@@ -30,6 +30,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "Hierarchy.h"
+#include "phys_constants.h"

/* function prototypes */
int GetUnits(float *DensityUnits, float *LengthUnits,
@@ -120,9 +121,7 @@

// put the temperature together
float mu = rho/num_density;
- float mp = 1.67262171e-24; // proton mass [g]
- float kb = 1.3806504e-16; // boltzmann constant [erg/K]
- float Temp = max((Gamma-1.0)*mu*mp*eint/kb, 1.0);
+ float Temp = max((Gamma-1.0)*mu*mh*eint/kboltz, 1.0);

// find temperature bin
float lamT = 3.15614e5/Temp;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -2394,6 +2394,41 @@

int FTStochasticForcing(int FieldDim); // WS

+
+ /* START Subgrid-scale modeling framework by P. Grete */
+
+ // Jacobians to be used in SGS model
+ float *JacVel[MAX_DIMENSION][MAX_DIMENSION];
+ float *JacB[MAX_DIMENSION][MAX_DIMENSION];
+
+ float *FilteredFields[7]; // filtered fields: rho, xyz-vel, Bxyz
+
+ // the scale-similarity model needs mixed filtered quantities
+ float *FltrhoUU[6];
+ float *FltBB[6];
+ float *FltUB[3];
+
+ int SGSUtil_ComputeJacobian(float *Jac[][MAX_DIMENSION],float* field1,float* field2,float* field3);
+ int SGSUtil_ComputeMixedFilteredQuantities();
+ int SGSUtil_FilterFields();
+
+ // the general functions that add the SGS terms to the dynamic eqns.
+ int SGS_AddEMFTerms(float **dU);
+ int SGS_AddMomentumTerms(float **dU);
+
+ // the different SGS models
+ void SGS_AddEMF_eddy_resistivity(float **EMF);
+ void SGS_AddEMF_nonlinear_compressive(float **EMF);
+ void SGS_AddMom_nonlinear_kinetic(float **Tau);
+ void SGS_AddMom_nonlinear_kinetic_scaled(float **Tau);
+ void SGS_AddMom_nonliner_magnetic(float **Tau);
+ void SGS_AddMom_eddy_viscosity_scaled(float **Tau);
+ void SGS_AddMom_scale_similarity_kinetic(float **Tau);
+ void SGS_AddMom_scale_similarity_magnetic(float **Tau);
+ void SGS_AddEMF_scale_similarity(float **EMF);
+
+ /* END Subgrid-scale modeling framework by P. Grete */
+
/* Comoving coordinate expansion terms. */

int ComovingExpansionTerms();

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_AddExternalAcceleration.C
--- a/src/enzo/Grid_AddExternalAcceleration.C
+++ b/src/enzo/Grid_AddExternalAcceleration.C
@@ -25,7 +25,8 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "SphericalInfall.h"
-
+#include "phys_constants.h"
+
/* function prototypes */

int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
@@ -39,7 +40,7 @@

if (ProblemType == 24) {

- float Pi = 3.14159, accel;
+ float accel;
FLOAT a = 1, dadt, rcubed, xpos, ypos = 0, zpos = 0;
int i, j, k, n = 0;

@@ -82,7 +83,7 @@
(i.e. 1/a^2 * a = 1/a). */

accel = GravitationalConstant*SphericalInfallFixedMass/
- (4.0*Pi*rcubed*a);
+ (4.0*pi*rcubed*a);

/* Apply force. */


diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_AddFeedbackSphere.C
--- a/src/enzo/Grid_AddFeedbackSphere.C
+++ b/src/enzo/Grid_AddFeedbackSphere.C
@@ -28,6 +28,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"

#ifdef CONFIG_BFLOAT_4
#define TOLERANCE 1e-06
@@ -51,8 +52,6 @@
{

const float WhalenMaxVelocity = 35; // km/s
- const double Msun = 1.989e33;
- const double c = 3.0e10;

int dim, i, j, k, index;
int sx, sy, sz;
@@ -117,7 +116,7 @@
if (MetalNum > 0 && SNColourNum > 0 && cstar->type == PopIII)
MetalNum = SNColourNum;

- float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
+ float BubbleVolume = (4.0 * pi / 3.0) * radius * radius * radius;

/***********************************************************************
SUPERNOVAE
@@ -144,7 +143,7 @@
// Correct if the volume with 27 cells is larger than the energy bubble volume
#ifdef UNUSED
float BoxVolume = 27 * CellWidth[0][0] * CellWidth[0][0] * CellWidth[0][0];
- float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
+ float BubbleVolume = (4.0 * pi / 3.0) * radius * radius * radius;
//printf("BoxVolume = %lg, BubbleVolume = %lg\n", BoxVolume, BubbleVolume);
if (BoxVolume > BubbleVolume) {
//printf("Reducing ejecta density by %g\n", BubbleVolume / BoxVolume);
@@ -174,7 +173,7 @@
// printf("grid::AFS: before: cstar->Mass = %lf\n", cstar->Mass);
if (cstar->FeedbackFlag != SUPERNOVA) {
float old_mass = (float)(cstar->Mass);
- cstar->Mass -= EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / Msun;
+ cstar->Mass -= EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / SolarMass;
float frac = old_mass / cstar->Mass;
cstar->vel[0] *= frac;
cstar->vel[1] *= frac;
@@ -466,7 +465,7 @@
metallicity_inside = 0.0, colour_inside = 0.0, rho_inside, rho_metal_inside, rho_colour_inside;
float m_cell_edge = 0.0, metal_cell_edge = 0.0, colour_cell_edge = 0.0,
metallicity_edge = 0.0, colour_edge = 0.0, rho_jet, rho_metal_jet, rho_colour_jet;
- float L_x, L_y, L_z, L_s, nx_L = 0.0, ny_L = 0.0, nz_L = 0.0, costheta = cos(3.1415926/3.9);
+ float L_x, L_y, L_z, L_s, nx_L = 0.0, ny_L = 0.0, nz_L = 0.0, costheta = cos(pi/3.9);
float EjectaMass, EjectaMetalMass, MBHJetsVelocity;

if (cstar->FeedbackFlag == MBH_JETS) {
@@ -492,7 +491,7 @@
/* Calculate mass that has accumulated thus far; since EjectaDensity is calculated
for isotropic MBH_THERMAL, we use it only to compute EjectaMass, not apply it directly. */

- cstar->NotEjectedMass += EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / Msun;
+ cstar->NotEjectedMass += EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / SolarMass;

// fprintf(stdout, "d1, d2, d3, i, j, k = %d %d %d / %d %d %d\n",
// GridDimension[0], GridDimension[1], GridDimension[2], i,j,k);
@@ -526,7 +525,7 @@

/* Find ejecta mass */

- EjectaMass = cstar->NotEjectedMass * Msun / DensityUnits / pow(LengthUnits,3.0);
+ EjectaMass = cstar->NotEjectedMass * SolarMass / DensityUnits / pow(LengthUnits,3.0);
EjectaMetalMass = EjectaMass * MBHFeedbackMetalYield;
OutputWhenJetsHaveNotEjected = FALSE;

@@ -561,8 +560,8 @@

//add random noise to theta and phi
srand(time(NULL));
- theta += MaximumNoiseAngle * M_PI / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
- phi += MaximumNoiseAngle * M_PI / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
+ theta += MaximumNoiseAngle * pi / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
+ phi += MaximumNoiseAngle * pi / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);

//get back to Cartesian coordinate; some tricks needed to preserve the signs of nx_L and ny_L
nx_L = sign(nx_L) * fabs(cos(theta))*sin(phi);
@@ -620,7 +619,7 @@
} // END jj-direction
} // END kk-direction

-// printf("EjectaM in Msun = %g, EjectaM = %g, EjectaMetalM = %g, m_cell_edge = %g, n_cell_edge = %d\n",
+// printf("EjectaM in SolarMass = %g, EjectaM = %g, EjectaMetalM = %g, m_cell_edge = %g, n_cell_edge = %d\n",
// cstar->NotEjectedMass, EjectaMass, EjectaMetalMass, m_cell_edge, n_cell_edge);

/* Calculate the jet density */
@@ -654,13 +653,13 @@
Note that EjectaThermalEnergy is never used; MBHFeedbackEnergyCoupling
should now be calculated considering gravitational redshift (Kim et al. 2010) */

- MBHJetsVelocity = c * sqrt( 2 * MBHFeedbackEnergyCoupling * MBHFeedbackRadiativeEfficiency
+ MBHJetsVelocity = clight * sqrt( 2 * MBHFeedbackEnergyCoupling * MBHFeedbackRadiativeEfficiency
/ MBHFeedbackMassEjectionFraction ) / VelocityUnits;

- if (MBHJetsVelocity * VelocityUnits > 0.99*c) {
+ if (MBHJetsVelocity * VelocityUnits > 0.99*clight) {
ENZO_VFAIL("grid::AddFS: MBHJetsVelocity is ultra-relativistic! (%g/ %g/ %g/ %g c)\n",
MBHFeedbackEnergyCoupling, MBHFeedbackRadiativeEfficiency,
- MBHFeedbackMassEjectionFraction, MBHJetsVelocity * VelocityUnits / c);
+ MBHFeedbackMassEjectionFraction, MBHJetsVelocity * VelocityUnits / clight);
}

/* Finally, add the jet feedback at the edges (outer part of the supercell) */
@@ -757,7 +756,7 @@
// old_mass, cstar->Mass, old_vel1, cstar->vel[1]);

fprintf(stdout, "grid::AddFS: jets injected (EjectaM = %g Ms, JetsVelocity = %g, grid v = %g, rho_jet = %g) along n_L = (%g, %g, %g)\n",
- EjectaMass * DensityUnits * pow(LengthUnits,3.0) / Msun,
+ EjectaMass * DensityUnits * pow(LengthUnits,3.0) / SolarMass,
MBHJetsVelocity, BaryonField[Vel3Num][n_cell_edge-1], rho_jet, nx_L, ny_L, nz_L);

} // END MBH_JETS
@@ -854,7 +853,7 @@

float kph, kheat;
sigma_HI *= (double) TimeUnits / ((double)LengthUnits * (double)LengthUnits)
- / (4.0 * M_PI);
+ / (4.0 * pi);
kph = (float) (Q_HI * sigma_HI);
kheat = kph * deltaE;

@@ -988,7 +987,7 @@
/* Now it's done, unmark. */

//cstar->FeedbackFlag = NO_FEEDBACK;
-
+

return SUCCESS;


diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_AddH2DissociationFromSources.C
--- a/src/enzo/Grid_AddH2DissociationFromSources.C
+++ b/src/enzo/Grid_AddH2DissociationFromSources.C
@@ -29,12 +29,12 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"
#ifdef FUTUREAP
#include "ActiveParticle.h"
#include "ActiveParticle_RadiationParticle.h"
#include "ActiveParticle_SmartStar.h"
#endif
-#include "phys_constants.h"

#define THRESHOLD_DENSITY_DB36 1e14
#define THRESHOLD_DENSITY_DB37 5e14
@@ -59,6 +59,7 @@
int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num, B1Num, B2Num, B3Num;
int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
DINum, DIINum, HDINum;
+
double l_char = 0.0, N_H2 = 0.0;
double H2ISigma = 3.71e-18;

@@ -130,7 +131,7 @@
// where the gas is fully molecular. If the solver still breaks,
// then this parameter can be safely increased to ~100 AU or turning
// on H2 self-shielding with RadiationShield = 2.
- float dilutionRadius = 4.848e-6 * pc / (double) LengthUnits; // 1 AU
+ float dilutionRadius = 4.848e-6 * pc_cm / (double) LengthUnits; // 1 AU
float dilRadius2 = dilutionRadius * dilutionRadius;
float LightTravelDist = TimeUnits * clight / LengthUnits;

@@ -196,9 +197,9 @@
int TemperatureField = 0;
/* Pre-compute some quantities to speed things up */
TemperatureField = this->GetTemperatureFieldNumberForH2Shield();
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
- kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * M_PI));
- kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+ kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * pi));
+ kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * pi));
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -291,7 +292,8 @@

double radius2, radius2_yz;

- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -365,9 +367,9 @@
int TemperatureField = 0;
/* Pre-compute some quantities to speed things up */
TemperatureField = this->GetTemperatureFieldNumberForH2Shield();
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
- kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * M_PI));
- kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+ kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * pi));
+ kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * pi));
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -455,7 +457,7 @@
{
float jeans_length = 0.0;

- jeans_length = 15*kboltz*T/(4.0*M_PI*GravConst*mh*dens*density_units);
+ jeans_length = 15*kboltz*T/(4.0*pi*GravConst*mh*dens*density_units);
return sqrt(jeans_length);
}


diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_AddH2DissociationFromTree.C
--- a/src/enzo/Grid_AddH2DissociationFromTree.C
+++ b/src/enzo/Grid_AddH2DissociationFromTree.C
@@ -23,6 +23,7 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"

#define MIN_OPENING_ANGLE 0.2 // 0.2 = arctan(11.3 deg)

@@ -38,8 +39,6 @@
int grid::AddH2DissociationFromTree(void)
{

- const double pc = 3.086e18, clight = 3e10;
-
int i, j, k, index, dim, ci;
FLOAT pos[MAX_DIMENSION];
FLOAT radius2;
@@ -77,7 +76,7 @@
H2ISigma *= (double)TimeUnits / ((double)LengthUnits * (double)LengthUnits);

// Dilution factor (prevent breaking in the rate solver near the star)
- float dilutionRadius = 10.0 * pc / (double) LengthUnits;
+ float dilutionRadius = 10.0 * pc_cm / (double) LengthUnits;
float dilRadius2 = dilutionRadius * dilutionRadius;

// Convert from #/s to RT units
@@ -87,7 +86,7 @@
/* Find sources in the tree that contribute to the cells */

SuperSourceEntry *Leaf;
- float factor = LConv_inv * H2ISigma / (4.0 * M_PI);
+ float factor = LConv_inv * H2ISigma / (4.0 * pi);
float angle;

Leaf = SourceClusteringTree;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_AddMagneticSupernovaeToList.C
--- a/src/enzo/Grid_AddMagneticSupernovaeToList.C
+++ b/src/enzo/Grid_AddMagneticSupernovaeToList.C
@@ -56,22 +56,22 @@
// below based on the resolution of the grid at the highest refinement level
if (UseMagneticSupernovaFeedback > 1){
sn_duration = 5.0 * dtFixed;
- MagneticSupernovaDuration = sn_duration * TimeUnits / SecondsPerYear;
+ MagneticSupernovaDuration = sn_duration * TimeUnits / yr_s;
sn_radius = 1.5 * this->CellWidth[0][0];
- MagneticSupernovaRadius = sn_radius * LengthUnits / pc;
+ MagneticSupernovaRadius = sn_radius * LengthUnits / pc_cm;
}
else {
// Converting time from years to seconds, then internal units
- sn_duration = MagneticSupernovaDuration * SecondsPerYear / TimeUnits;
+ sn_duration = MagneticSupernovaDuration * yr_s / TimeUnits;
if(sn_duration < 5.0 * dtFixed) {
printf("WARNING: Magnetic supernova feedback duration is less than the recommended minimum of 5 x dtFixed\n");
- printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / SecondsPerYear);
+ printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / yr_s);
}
// Converting radius from parsecs to cm, then internal units
- sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
+ sn_radius = MagneticSupernovaRadius * pc_cm / LengthUnits;
if(sn_radius < 1.5 * this->CellWidth[0][0]){
printf("WARNING: Magnetic supernova feedback radius is less than the recommended minimum of 1.5 x CellWidth\n");
- printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc);
+ printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc_cm);
}

}
@@ -93,7 +93,7 @@

random_u = (float)(mt_random()%32768)/32768.0; // random variable from 0 to 1
random_v = (float)(mt_random()%32768)/32768.0;
- random_phi = 2*M_PI*random_u; // 0 to 2pi
+ random_phi = 2*pi*random_u; // 0 to 2pi
random_theta = acos(2*random_v-1); // 0 to pi

// Setting up randomly oriented magnetic feedback of supernova

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_AddRadiationImpulse.C
--- a/src/enzo/Grid_AddRadiationImpulse.C
+++ b/src/enzo/Grid_AddRadiationImpulse.C
@@ -23,6 +23,7 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"

int FindField(int f, int farray[], int n);
int GetUnits(float *DensityUnits, float *LengthUnits,
@@ -33,8 +34,6 @@
FLOAT BirthTime, FLOAT* pos)
{

- const double pc = 3.086e19, clight = 3e10;
-
int i, j, k, index, dim, FieldNum;
float kdiss_r2;
FLOAT delx, dely, delz, DomainWidth[MAX_DIMENSION];
@@ -62,7 +61,7 @@
DomainWidth[dim] = DomainRightEdge[dim] - DomainLeftEdge[dim];

// Dilution factor (prevent breaking in the rate solver near the star)
- float dilutionRadius = 10.0 * pc / (double) LengthUnits;
+ float dilutionRadius = 10.0 * pc_cm / (double) LengthUnits;
float dilRadius2 = dilutionRadius * dilutionRadius;
float LightTravelDist = TimeUnits * clight / LengthUnits;

@@ -73,7 +72,7 @@
/* Loop over cells */

index = 0;
- kdiss_r2 = (float) (Luminosity * sigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (Luminosity * sigma / (4.0 * pi));
for (k = 0; k < GridDimension[2]; k++) {
delz = fabs(CellLeftEdge[2][k] + 0.5*CellWidth[2][k] - pos[2]);
delz = min(delz, DomainWidth[2]-delz);

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_CRShockTubesInitializeGrid.C
--- a/src/enzo/Grid_CRShockTubesInitializeGrid.C
+++ b/src/enzo/Grid_CRShockTubesInitializeGrid.C
@@ -10,6 +10,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -216,7 +217,7 @@
for (i = 0; i < GridDimension[0]; i++) {
x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
x = x - x0; // translate gaussian to center
- BaryonField[iCRD][i] = 1.0/sqrt(4.0*3.14159*CRkappa*t0)
+ BaryonField[iCRD][i] = 1.0/sqrt(4.0*pi*CRkappa*t0)
* PEXP( -x*x/(4.0*CRkappa*t0));
} // end i for
} // end if

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_ChangeParticleTypeBeforeSN.C
--- a/src/enzo/Grid_ChangeParticleTypeBeforeSN.C
+++ b/src/enzo/Grid_ChangeParticleTypeBeforeSN.C
@@ -23,6 +23,7 @@
#include "GridList.h"
#include "ExternalBoundary.h"
#include "Grid.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -40,7 +41,6 @@
if (Stars == NULL)
return SUCCESS;

- const float pc = 3.086e18, mh = 1.673e-24, Msun = 1.989e33;
const float PISNLowerMass = 140.0, PISNUpperMass = 260.0;
const float StartRefineAtTime = 0.99; // Percentage of stellar lifetime
const float EndRefineAtTime = 1.0;
@@ -98,7 +98,7 @@
if (factor < 1) {

Diameter = 2 * BufferZone * PopIIISupernovaRadius *
- (pc/LengthUnits);
+ (pc_cm/LengthUnits);

}

@@ -166,9 +166,6 @@

/* ASSUMES POP III PAIR-INSTABILITY SUPERNOVA FOR ENERGY */

- const float pc = 3.086e18, mh = 1.673e-24, Msun = 1.989e33;
- const float kb = 1.38e-16;
-
FLOAT StartTime, SoundCrossingTime;
float HeliumCoreMass, TransitionTime, ShockVelocity, STradius;
float BlastWaveRadius, SoundSpeed, SNTemperature;
@@ -183,14 +180,14 @@

// With a constant shock velocity, calculate time it takes for
// blastwave to travel from r=0 to r=PopIIISupernovaRadius
- ShockVelocity = 1.165 * sqrt(2.0 * SNenergy / (Mass * Msun));
- StartTime = PopIIISupernovaRadius * pc / ShockVelocity;
+ ShockVelocity = 1.165 * sqrt(2.0 * SNenergy / (Mass * SolarMass));
+ StartTime = PopIIISupernovaRadius * pc_cm / ShockVelocity;

// Because we inject thermal energy, the blastwave is delayed by a
// sound crossing time.
- SNTemperature = min(double(Mass*Msun) / double(mh) / kb, 1e8);
- SoundSpeed = sqrt(kb * SNTemperature / (0.6*mh));
- SoundCrossingTime = PopIIISupernovaRadius * pc / SoundSpeed;
+ SNTemperature = min(double(Mass*SolarMass) / double(mh) / kboltz, 1e8);
+ SoundSpeed = sqrt(kboltz * SNTemperature / (0.6*mh));
+ SoundCrossingTime = PopIIISupernovaRadius * pc_cm / SoundSpeed;

// units in cm
STradius = 3.62e19 * powf(Mass/100., 1./3.) * powf(n0, -1./3.);
@@ -212,7 +209,7 @@
// it slows by a factor of 2 near the transition to the ST phase)
if (Time < TransitionTime) {

- BlastWaveRadius = PopIIISupernovaRadius * pc + ShockVelocity * Time;
+ BlastWaveRadius = PopIIISupernovaRadius * pc_cm + ShockVelocity * Time;
printf("Free expansion\n");

} // ENDIF Free expansion phase

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_ClusterInitializeGrid.C
--- a/src/enzo/Grid_ClusterInitializeGrid.C
+++ b/src/enzo/Grid_ClusterInitializeGrid.C
@@ -21,6 +21,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -95,8 +96,8 @@

/* Set various units. */

- const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
- pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16, keV=1.1604e7;
+ const double keV=1.1604e7;
+
float DensityUnits, LengthUnits, TemperatureUnits = 1, TimeUnits,
VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6;
double MassUnits = 1;
@@ -110,7 +111,7 @@
BoxLength = ComovingBoxSize*ExpansionFactor/HubbleConstantNow; // in Mpc
} else {
CriticalDensity = 2.78e11*POW(0.74,2); // in Msolar/Mpc^3 for h=0.74
- BoxLength = LengthUnits / 3.086e24;
+ BoxLength = LengthUnits / Mpc_cm;
HubbleConstantNow = 1.0;
OmegaMatterNow = 1.0;
}
@@ -161,32 +162,32 @@
x1 = NFWRadius[i]/SphereCoreRadius[sphere];
NFWDensity[i] = SphereDensity[sphere]/(x1*(1.0+x1)*(1.0+x1)); // DM Density
if (SphereType[sphere]>=6 && SphereType[sphere] <= 8) { //aka 6, 7, 8: Perseus Cluster
- rkpc=NFWRadius[i]*LengthUnits/(1.0e-3*Mpc);
+ rkpc=NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm);
/* Initial Temperature */
if (rkpc > 300.0){
NFWTemp[i]=7.0594*1.3*POW((1.0+1.5*NFWRadius[i]/SphereRadius[sphere]),-1.6)*keV; // in K
} else{
- NFWTemp[i]=8.12e7*(1.0+POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc*71),3))/(2.3 + pow(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc*71),3));
+ NFWTemp[i]=8.12e7*(1.0+POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm*71),3))/(2.3 + pow(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm*71),3));
if (SphereType[sphere]==8)
NFWTemp[i]=8.12e7;
}
/* Set Gravity. NFW Dark Matter, BCG+BH */
Allg[i]=GravConst*PointSourceGravityConstant*SolarMass*
((log(1.0+x1)-x1/(1.0+x1)) /(log(1.0+1.0)-1.0/(1.0+1.0)))/POW(NFWRadius[i]*LengthUnits, 2.0) +
- ClusterSMBHBCG*POW((POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) +
+ ClusterSMBHBCG*POW((POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) +
GravConst*SolarMass* ClusterSMBHMass / POW(NFWRadius[i]*LengthUnits, 2) ;
}//end Perseus
else{ //Elliptical galaxies
if (SphereType[sphere]==1) { //NGC 4472
- NFWTemp[i]=1.16059e7*(1.17-(1.17-0.6)*exp(-NFWRadius[i]*LengthUnits/(2.0*5*1.0e-3*Mpc)));
+ NFWTemp[i]=1.16059e7*(1.17-(1.17-0.6)*exp(-NFWRadius[i]*LengthUnits/(2.0*5*1.0e-3*Mpc_cm)));
}
if (SphereType[sphere]==2) { //NGC 6482
- NFWTemp[i]=1.16059e7*(0.4+0.4*exp(-NFWRadius[i]*LengthUnits/(2.0*8.0*1.0e-3*Mpc)));
+ NFWTemp[i]=1.16059e7*(0.4+0.4*exp(-NFWRadius[i]*LengthUnits/(2.0*8.0*1.0e-3*Mpc_cm)));
}
Allg[i]=GravConst*PointSourceGravityConstant*SolarMass*
((log(1.0+x1)-x1/(1.0+x1)) /(log(1.0+1.0)-1.0/(1.0+1.0)))/POW(NFWRadius[i]*LengthUnits, 2.0) +
- GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(NFWRadius[i]*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc/1.8153, 2) + //ClusterSMBHBCG is M_* here
+ GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(NFWRadius[i]*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc_cm/1.8153, 2) + //ClusterSMBHBCG is M_* here
GravConst*SolarMass* ClusterSMBHMass / POW(NFWRadius[i]*LengthUnits, 2) ;
}

@@ -197,7 +198,7 @@
}
dpdr_old=dpdr;
if (i>0) {
- NFWPressure[i]= (NFWPressure[i-1] + (NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc*0.5*dpdr_old)/(1.0+(NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc*0.5*Allg[i]/(kboltz*NFWTemp[i]/(mu * mh)));
+ NFWPressure[i]= (NFWPressure[i-1] + (NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc_cm*0.5*dpdr_old)/(1.0+(NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc_cm*0.5*Allg[i]/(kboltz*NFWTemp[i]/(mu * mh)));
GasDensity[i]=NFWPressure[i]/(kboltz * NFWTemp[i]/(mu * mh));
dpdr = -Allg[i]* GasDensity[i];
}
@@ -275,8 +276,8 @@
temp1 = NFWTemp[m] + (NFWTemp[m-1] - NFWTemp[m])*
(r - NFWRadius[m])/(NFWRadius[m-1] - NFWRadius[m]); // in K
if (SphereType[sphere] == 6) {
- gas_dens1 =mh*(0.0192/(1.0+POW(r*LengthUnits/(18.0e-3*Mpc),3.0))+0.046/pow((1.0+pow(r*LengthUnits/(57.0e-3*Mpc), 2.0)), 1.8)+
- 0.0048/POW((1.0+pow(r*LengthUnits/(200.0e-3*Mpc), 2.0)), 1.1))/DensityUnits/0.88;
+ gas_dens1 =mh*(0.0192/(1.0+POW(r*LengthUnits/(18.0e-3*Mpc_cm),3.0))+0.046/pow((1.0+pow(r*LengthUnits/(57.0e-3*Mpc_cm), 2.0)), 1.8)+
+ 0.0048/POW((1.0+pow(r*LengthUnits/(200.0e-3*Mpc_cm), 2.0)), 1.1))/DensityUnits/0.88;
}
else{
gas_dens1 = (GasDensity[m] + (GasDensity[m-1] - GasDensity[m])*

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_ClusterSMBHEachGridGasMass.C
--- a/src/enzo/Grid_ClusterSMBHEachGridGasMass.C
+++ b/src/enzo/Grid_ClusterSMBHEachGridGasMass.C
@@ -65,7 +65,7 @@

int dim = 0;
float DiskRadius; //ClusterSMBHDiskRadius = 0.5; ClusterSMBHDiskRadiu now a parameter.
- DiskRadius = ClusterSMBHDiskRadius*kpc/LengthUnits; //from kpc to codeunits
+ DiskRadius = ClusterSMBHDiskRadius*kpc_cm/LengthUnits; //from kpc to codeunits
for (dim = 0; dim < GridRank; dim++) {
DiskCenter[dim] = PointSourceGravityPosition[dim];
DiskLeftCorner[dim] = PointSourceGravityPosition[dim]- DiskRadius;

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_ClusterSMBHFeedback.C
--- a/src/enzo/Grid_ClusterSMBHFeedback.C
+++ b/src/enzo/Grid_ClusterSMBHFeedback.C
@@ -85,7 +85,7 @@
float JetVelocity, FastJetVelocity; // Jet Velocity in km/s (should make parameter)-- gets value from parameter ClusterSMBHJetVelocity
JetScaleRadius = ClusterSMBHJetRadius/2.0; //JetScaleRadius is half the radius of the jet launch region in cellwidths
float DiskRadius; //ClusterSMBHDiskRadius = 0.5; //ClusterSMBHDiskRadiu now a parameter
- DiskRadius = ClusterSMBHDiskRadius*kpc/LengthUnits; //from kpc to codeunits
+ DiskRadius = ClusterSMBHDiskRadius*kpc_cm/LengthUnits; //from kpc to codeunits

for (dim = 0; dim < GridRank; dim++) {
JetCenter[dim] = PointSourceGravityPosition[dim];
@@ -100,8 +100,8 @@
DiskRightCorner[dim] = PointSourceGravityPosition[dim] + DiskRadius;
}

- JetLeftCorner[jet_dim] -= ClusterSMBHJetLaunchOffset*kpc/LengthUnits;
- JetRightCorner[jet_dim] += ClusterSMBHJetLaunchOffset*kpc/LengthUnits; //from kpc to codeunits
+ JetLeftCorner[jet_dim] -= ClusterSMBHJetLaunchOffset*kpc_cm/LengthUnits;
+ JetRightCorner[jet_dim] += ClusterSMBHJetLaunchOffset*kpc_cm/LengthUnits; //from kpc to codeunits


/* Compute indices of jet launch region. */
@@ -143,7 +143,7 @@
} // end: loop over dim

/* Compute mass and momentum to be put into cells in code units if Jet is on this grid. */
- JetMdot = (ClusterSMBHJetMdot*SolarMass/3.1557e7)/(MassUnits/TimeUnits); // from M_sun/yr to code units
+ JetMdot = (ClusterSMBHJetMdot*SolarMass/yr_s)/(MassUnits/TimeUnits); // from M_sun/yr to code units
if (JetOnGrid == true){
float JetNormalization = 0.0, density_normalization, radius, Tramp;
if (jet_dim == 2){
@@ -179,9 +179,9 @@

/* Convert to code units. */
density_normalization = (JetMdot/JetNormalization)*dtFixed/POW(CellWidth[0][0], 3);
- Tramp = ClusterSMBHTramp*1.0e6*3.1557e7/TimeUnits; // from Myr to code units
+ Tramp = ClusterSMBHTramp*1.0e6*yr_s/TimeUnits; // from Myr to code units

- JetVelocity = sqrt((ClusterSMBHJetEdot*1.0e44*ClusterSMBHKineticFraction*2)/(ClusterSMBHJetMdot*SolarMass/3.1557e7))/VelocityUnits;
+ JetVelocity = sqrt((ClusterSMBHJetEdot*1.0e44*ClusterSMBHKineticFraction*2)/(ClusterSMBHJetMdot*SolarMass/yr_s))/VelocityUnits;
JetVelocity *= min((Time-ClusterSMBHStartTime)/Tramp, 1.0); //linear ramp

/* Clip edge of jet launching disk so we don't set cell off the edge of the grid. */
@@ -212,7 +212,7 @@
JetVelocity_x = JetVelocity; // mutiplied by sincos later
JetVelocity_y = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001 || ClusterSMBHJetPrecessionPeriod < -0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_z = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
@@ -273,7 +273,7 @@
JetVelocity_z = JetVelocity; // mutiplied by sincos later
JetVelocity_y = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_x = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
@@ -335,7 +335,7 @@
JetVelocity_z = JetVelocity; // mutiplied by sincos later
JetVelocity_x = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_y = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));

diff -r ecced5aa9f43 -r 48fbcd4ec204 src/enzo/Grid_CollapseTestInitializeGrid.C
--- a/src/enzo/Grid_CollapseTestInitializeGrid.C
+++ b/src/enzo/Grid_CollapseTestInitializeGrid.C
@@ -24,6 +24,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"

#define NTHETA 1000
#define NR 1000
@@ -157,8 +158,6 @@

/* Set various units. */

- const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
- pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16, LightSpeed = 2.9979e10;
float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits,
VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6;

@@ -168,11 +167,11 @@
if (ComovingCoordinates) {
CosmologyComputeExpansionFactor(Time, &a, &dadt);
ExpansionFactor = a/(1.0+InitialRedshift);
- CriticalDensity = 2.78e11*pow(HubbleConstantNow, 2); // in Msolar/Mpc^3
+ CriticalDensity = 2.78e11*POW(HubbleConstantNow, 2); // in Msolar/Mpc^3
BoxLength = ComovingBoxSize*ExpansionFactor/HubbleConstantNow; // in Mpc
} else {
- CriticalDensity = 2.78e11*pow(0.74,2); // in Msolar/Mpc^3 for h=0.74
- BoxLength = LengthUnits / 3.086e24;
+ CriticalDensity = 2.78e11*POW(0.74,2); // in Msolar/Mpc^3 for h=0.74
+ BoxLength = LengthUnits / Mpc_cm;
HubbleConstantNow = 1.0;
OmegaMatterNow = 1.0;
}
@@ -193,20 +192,20 @@
NFWPressure[0] = 1.0 * kboltz * InitialTemperature / (mu * mh);
FILE *fptr = fopen("NFWProfile.out", "w");
for (i = 0; i < NFW_POINTS; i++) {
- NFWRadius[i] = SphereRadius[sphere]*pow(10, -3*(float(i)/NFW_POINTS));
+ NFWRadius[i] = SphereRadius[sphere]*POW(10, -3*(float(i)/NFW_POINTS));
x1 = NFWRadius[i]/SphereCoreRadius[sphere];
NFWDensity[i] = SphereDensity[sphere]/(x1*(1.0+x1)*(1.0+x1));
NFWMass[i] = 4.0*pi*SphereDensity[sphere]*
- (CriticalDensity/pow(ExpansionFactor, 3)) *
- pow(SphereCoreRadius[sphere]*BoxLength, 3) *
+ (CriticalDensity/POW(ExpansionFactor, 3)) *
+ POW(SphereCoreRadius[sphere]*BoxLength, 3) *
(log(1.0+x1) - x1/(x1+1.0)); // in Msolar
dpdr_old = dpdr;
dpdr = GravConst * NFWMass[i] * SolarMass *
NFWDensity[i] /
- pow(NFWRadius[i]*BoxLength*Mpc, 2);
+ POW(NFWRadius[i]*BoxLength*Mpc_cm, 2);
if (i > 0)
NFWPressure[i] = NFWPressure[i-1] -
- 0.5*(dpdr+dpdr_old)*(NFWRadius[i]-NFWRadius[i-1])*BoxLength*Mpc;
+ 0.5*(dpdr+dpdr_old)*(NFWRadius[i]-NFWRadius[i-1])*BoxLength*Mpc_cm;
NFWTemp[i] = NFWPressure[i]*mu*mh/(kboltz*NFWDensity[i]); // in K
NFWSigma[i] = sqrt(kboltz * NFWTemp[i] / (mu * mh)); // in cm/s
float mean_overdensity = 3.0*SphereDensity[sphere] / (x1*x1*x1) *
@@ -318,7 +317,7 @@
switch (SphereType[sphere]) {

case 1:
- SphereMass = (4*pi/3)*pow((SphereRadius[sphere]*LengthUnits), 3) *
+ SphereMass = (4*pi/3.0)*POW((SphereRadius[sphere]*LengthUnits), 3) *
(SphereDensity[sphere]*DensityUnits);
printf("mass = %"GSYM", lunit = %"GSYM", dunit = %"GSYM", rho = %"GSYM", r = %"GSYM"\n",
SphereMass, LengthUnits, DensityUnits, SphereDensity[sphere],
@@ -332,32 +331,32 @@
case 4:
// Integral of a Gaussian
term1 = 0.25 * sqrt(pi) * 0.3536 * // 2^(-3/2) = 0.3536
- pow((SphereCoreRadius[sphere]), 3) *
+ POW((SphereCoreRadius[sphere]), 3) *
ERF(SphereRadius[sphere]/SphereCoreRadius[sphere]);
term2 = 0.25 * SphereRadius[sphere] *
- pow(SphereCoreRadius[sphere], 2) *
- PEXP(-0.5 * pow((SphereRadius[sphere]/SphereCoreRadius[sphere]), 2));
+ POW(SphereCoreRadius[sphere], 2) *
+ PEXP(-0.5 * POW((SphereRadius[sphere]/SphereCoreRadius[sphere]), 2));
SphereMass = (4*pi*SphereDensity[sphere]*DensityUnits) *
- pow(LengthUnits, 3) * (term1 - term2);
+ POW(LengthUnits, 3) * (term1 - term2);
break;

case 5:
SphereCoreDens = (SphereDensity[sphere]*DensityUnits) *
- pow(SphereCoreRadius[sphere] / SphereRadius[sphere], -2);
- SphereCoreMass = (4*pi/3) *
- pow((SphereCoreRadius[sphere]*LengthUnits),3) * (SphereCoreDens);
+ POW(SphereCoreRadius[sphere] / SphereRadius[sphere], -2);
+ SphereCoreMass = (4*pi/3.0) *
+ POW((SphereCoreRadius[sphere]*LengthUnits),3) * (SphereCoreDens);
SphereMass = (4*pi*SphereDensity[sphere]*DensityUnits *
- pow(SphereRadius[sphere]*LengthUnits, 2) *
+ POW(SphereRadius[sphere]*LengthUnits, 2) *
((SphereRadius[sphere] - SphereCoreRadius[sphere]) *
LengthUnits)) + SphereCoreMass;
break;

case 6:
- SphereMass = pow(VelocitySound[sphere],3) /
- (sqrt(4*pi*pow((GravConst)/(4*pi),3) *
+ SphereMass = POW(VelocitySound[sphere],3) /
+ (sqrt(4*pi*POW((GravConst)/(4*pi),3) *
(SphereDensity[sphere]))) *
q(SphereCutOff[sphere]);
- SphereMass = SphereMass*DensityUnits*pow(LengthUnits,3);
+ SphereMass = SphereMass*DensityUnits*POW(LengthUnits,3);
break;

case 7:
@@ -376,7 +375,7 @@
break;

case 8:
- SphereMass = (4*pi/3)*pow((SphereRadius[sphere]*LengthUnits), 3) *
+ SphereMass = (4*pi/3.0)*POW((SphereRadius[sphere]*LengthUnits), 3) *
(SphereDensity[sphere]*DensityUnits);
ComputeRadialVelocity(SphereDensity[sphere], SphereMass,
SphereRadius[sphere], VelocityUnits,
@@ -402,15 +401,15 @@

BHMass = SphereMass / SolarMass; // in solar masses
SchwarzschildRadius = 2.0 * GravConst * SphereMass /
- (LightSpeed*LightSpeed);
- CavityRadius = 117.0 * SchwarzschildRadius * pow((BHMass/1e6), 0.08);
- InnerDensity = 4.31e-10 * pow((BHMass/1e8), -0.8) *
- pow((CavityRadius/SchwarzschildRadius/1e3), -0.6) / DensityUnits;
- InnerTemperature = 1.7e6 * pow((BHMass/1e6), -0.28);
- InnerScaleHeight = 0.46 * CavityRadius * pow(BHMass/1e6, -0.12) /
+ (clight*clight);
+ CavityRadius = 117.0 * SchwarzschildRadius * POW((BHMass/1e6), 0.08);
+ InnerDensity = 4.31e-10 * POW((BHMass/1e8), -0.8) *
+ POW((CavityRadius/SchwarzschildRadius/1e3), -0.6) / DensityUnits;
+ InnerTemperature = 1.7e6 * POW((BHMass/1e6), -0.28);
+ InnerScaleHeight = 0.46 * CavityRadius * POW(BHMass/1e6, -0.12) /
LengthUnits;
ThickenTransitionRadius = 1.9e3 * SchwarzschildRadius *
- pow((BHMass/1e6), 2.0/21) / LengthUnits;
+ POW((BHMass/1e6), 2.0/21) / LengthUnits;
CavityRadius /= LengthUnits;
printf("cgs: %"GSYM" %"GSYM" %"GSYM" %"GSYM" %"GSYM" %"GSYM"\n", SchwarzschildRadius,
CavityRadius*LengthUnits, InnerDensity*DensityUnits,
@@ -508,9 +507,9 @@
else if (xpos > 0 && ypos < 0)
theta = 2*pi + atan(ypos/xpos);
} else if (xpos == 0 && ypos > 0)
- theta = 3.14159 / 2.0;
+ theta = pi / 2.0;
else if (xpos == 0 && ypos < 0)
- theta = (3*3.14159) / 2.0;
+ theta = (3*pi) / 2.0;
else
theta = 0.0;

@@ -529,7 +528,7 @@
RotVelocity[2] = 2*pi*(xpos*sin(a)) /
SphereRotationalPeriod[sphere];
if (SphereType[sphere] == 7 || SphereType[sphere] == 9)
- GasRotVelocityCorrection = pow(rcyl/SphereRadius[sphere],-1.5);
+ GasRotVelocityCorrection = POW(rcyl/SphereRadius[sphere],-1.5);
for (dim = 0; dim < MAX_DIMENSION; dim++) {
Velocity[dim] = GasRotVelocityCorrection * RotVelocity[dim];
DMVelocity[dim] = DMRotVelocityCorrection * RotVelocity[dim];
@@ -566,10 +565,10 @@
if (SphereType[sphere] == 1)
dens1 = SphereDensity[sphere];

- /* 2) r^-2 power law */
+ /* 2) r^-2 POWer law */

if (SphereType[sphere] == 2)
- dens1 = SphereDensity[sphere]*pow(r/SphereRadius[sphere], -2);
+ dens1 = SphereDensity[sphere]*POW(r/SphereRadius[sphere], -2);

/* 3) NFW profile (use look-up table for temperature and
velocity dispersion)*/
@@ -591,17 +590,17 @@

if (SphereType[sphere] == 4) {
dens1 = SphereDensity[sphere]*
- PEXP(-0.5*pow(r/SphereCoreRadius[sphere], 2));
+ PEXP(-0.5*POW(r/SphereCoreRadius[sphere], 2));
}

- /* 5) r^-2 power law with core radius */
+ /* 5) r^-2 POWer law with core radius */

if (SphereType[sphere] == 5) {
if (r < SphereCoreRadius[sphere]) {
- dens1 = SphereDensity[sphere]*pow(SphereCoreRadius[sphere]/
+ dens1 = SphereDensity[sphere]*POW(SphereCoreRadius[sphere]/
SphereRadius[sphere], -2);
} else {
- dens1 = SphereDensity[sphere]*pow(r/SphereRadius[sphere], -2);
+ dens1 = SphereDensity[sphere]*POW(r/SphereRadius[sphere], -2);
}
}

@@ -638,14 +637,14 @@
/* 9) Circumbinary BH accretion disk (Lippai et al. 2008) */

if (SphereType[sphere] == 9) {
- MidplaneDensity = InnerDensity * pow(rcyl/CavityRadius, -0.6);
+ MidplaneDensity = InnerDensity * POW(rcyl/CavityRadius, -0.6);
MidplaneTemperature =
- InnerTemperature * pow(rcyl/CavityRadius, -0.9);
+ InnerTemperature * POW(rcyl/CavityRadius, -0.9);
if (rcyl < ThickenTransitionRadius)
ScaleHeight = InnerScaleHeight;
else
ScaleHeight = InnerScaleHeight *
- pow(rcyl/ThickenTransitionRadius, 1.05);
+ POW(rcyl/ThickenTransitionRadius, 1.05);
ScaleHeight = max(ScaleHeight, 2*CellWidth[0][0]);

// printf("r=%"FSYM", z=%"FSYM", h=%"GSYM", rho=%"GSYM", T=%"GSYM"\n",
@@ -708,7 +707,7 @@

if (dim == 0)
dens1 = SphereDensity[sphere]*PEXP(-drad/ScaleHeightR)/
- pow(cosh(zheight/max(ScaleHeightz, CellWidth[0][0])), 2);
+ POW(cosh(zheight/max(ScaleHeightz, CellWidth[0][0])), 2);

// if (dens1 < density)
// break;
@@ -723,11 +722,11 @@
float accel = PointSourceGravityConstant;
if (PointSourceGravity == 1)
accel = PointSourceGravityConstant/
- (pow(drad,3) + pow(PointSourceGravityCoreRadius, 3));
+ (POW(drad,3) + POW(PointSourceGravityCoreRadius, 3));
if (PointSourceGravity == 2) {
x1 = drad/PointSourceGravityCoreRadius;
accel = PointSourceGravityConstant*(log(1+x1)-x1/(1+x1))/
- pow(drad, 3);
+ POW(drad, 3);
}

float vel = sqrt(accel);
@@ -748,10 +747,10 @@

} // end: disk

- /* 11) stellar wind, with r^-2 power law density*/
+ /* 11) stellar wind, with r^-2 POWer law density*/
if (SphereType[sphere] == 11) {
radial_velocity = StellarWindSpeed/VelocityUnits;
- dens1 = StellarWindDensity*pow(r/StellarWindRadius, -2);
+ dens1 = StellarWindDensity*POW(r/StellarWindRadius, -2);
Velocity[0] += radial_velocity * xpos / r;
Velocity[1] += radial_velocity * ypos / r;
Velocity[2] += radial_velocity * zpos / r;
@@ -902,13 +901,13 @@

if (MultiSpecies > 1) {
BaryonField[HMNum][n] = CollapseTestInitialFractionHM *
- BaryonField[HIINum][n]* pow(temperature,float(0.88));
+ BaryonField[HIINum][n]* POW(temperature,float(0.88));
BaryonField[H2IINum][n] = CollapseTestInitialFractionH2II *
- 2.0*BaryonField[HIINum][n]* pow(temperature,float(1.8));
+ 2.0*BaryonField[HIINum][n]* POW(temperature,float(1.8));
if (ComovingCoordinates)
BaryonField[H2INum][n] = H2I_Fraction *
- BaryonField[0][n]*CoolData.HydrogenFractionByMass*pow(301.0,5.1)*
- pow(OmegaMatterNow, float(1.5))/
+ BaryonField[0][n]*CoolData.HydrogenFractionByMass*POW(301.0,5.1)*
+ POW(OmegaMatterNow, float(1.5))/
(OmegaMatterNow*BaryonMeanDensity)/
// CosmologySimulationOmegaBaryonNow/
HubbleConstantNow*2.0;
@@ -968,7 +967,7 @@

if (HydroMethod != Zeus_Hydro)
for (dim = 0; dim < GridRank; dim++)
- BaryonField[1][n] += 0.5*pow(BaryonField[ivel+dim][n], 2);
+ BaryonField[1][n] += 0.5*POW(BaryonField[ivel+dim][n], 2);

/* Set particles if being used (generate a number of particle
proportional to density). */
@@ -977,11 +976,11 @@
if (i >= GridStartIndex[0] && i <= GridEndIndex[0] &&
j >= GridStartIndex[1] && j <= GridEndIndex[1] &&
k >= GridStartIndex[2] && k <= GridEndIndex[2] ) {
- ParticleCount += density/pow(float(RefineBy), GridRank*level);
+ ParticleCount += density/POW(float(RefineBy), GridRank*level);
while (ParticleCount > 1) {
if (SetupLoopCount > 0) {
ParticleMass[npart] = ParticleMeanDensity*
- pow(float(RefineBy), GridRank*level);
+ POW(float(RefineBy), GridRank*level);
ParticleNumber[npart] = CollapseTestParticleCount++;
ParticleType[npart] = PARTICLE_TYPE_DARK_MATTER;

@@ -1065,8 +1064,8 @@
double BE(double r)
{
double factor;
- factor = 4.8089e-04*pow(r,5) - 1.0173e-02*pow(r,4) + 7.7899e-02*pow(r,3) -
- 2.3299e-01*pow(r,2) + 1.4721e-02*r + 1.0008e+00;
+ factor = 4.8089e-04*POW(r,5) - 1.0173e-02*POW(r,4) + 7.7899e-02*POW(r,3) -
+ 2.3299e-01*POW(r,2) + 1.4721e-02*r + 1.0008e+00;
return factor;
}

@@ -1075,8 +1074,8 @@
double q(double r)
{
double factor;
- factor = 0.0015970*pow(r,5) - 0.0229113*pow(r,4) + 0.0386709*pow(r,3) +
- 0.7350457*pow(r,2) - 0.5490283*r + 0.0872061;
+ factor = 0.0015970*POW(r,5) - 0.0229113*POW(r,4) + 0.0386709*POW(r,3) +
+ 0.7350457*POW(r,2) - 0.5490283*r + 0.0872061;
return factor;
}

@@ -1092,13 +1091,10 @@
double Maxwellian(double c_tilda, double vel_unit, double mu, double gamma)
{
// Set constants
- double mh = 1.67e-24;
- double kboltz = 1.38e-16;
- double pi = 3.14159;

// Compute temperature in cgs units
double c = c_tilda*vel_unit;
- double T = (pow(c,2)*mu*mh)/(kboltz*gamma);
+ double T = (POW(c,2)*mu*mh)/(kboltz*gamma);

// Compute random Maxwellian velocity
double mean = 0;
@@ -1114,8 +1110,8 @@

double ERF(double x)
{
- return (2.0 / sqrt(M_PI)) *
- (x - pow(x,3)/3.0 + pow(x,5)/10.0 - pow(x,7)/42.0 + pow(x,9)/216.0);
+ return (2.0 / sqrt(pi)) *
+ (x - POW(x,3)/3.0 + POW(x,5)/10.0 - POW(x,7)/42.0 + POW(x,9)/216.0);
}

int ComputeRadialVelocity(float density, double mass, float r_init,
@@ -1124,8 +1120,7 @@
double exterior_rho[], int Npts)
{

- const float theta0 = 0.5*M_PI, theta1 = 1.9*M_PI;
- const float Grav = 6.673e-8, Mpc = 3.086e24, yr = 3.1557e7, Msun = 1.989e33;
+ const float theta0 = 0.5*pi, theta1 = 1.9*pi;
const float delta_i = 0.5;
float dtheta, delta_t[NTHETA], Theta[NTHETA], Lambda[NTHETA], Chi;
float beta[NTHETA], little_d[NTHETA];
@@ -1138,44 +1133,44 @@
dtheta = (theta1 - theta0) / (NTHETA-1.0);
for (i = 0; i < NTHETA; i++) {
Theta[i] = theta0 + i*dtheta;
- Lambda[i] = pow(sin(Theta[i]/2), 2) *
- pow((Theta[i] - sin(Theta[i])) / M_PI, -8.0/9);
- delta_t[i] = 4.5 * pow((Theta[i] - sin(Theta[i])), 2) /
- pow((1.0 - cos(Theta[i])), 3) - 1.0;
+ Lambda[i] = POW(sin(Theta[i]/2), 2) *
+ POW((Theta[i] - sin(Theta[i])) / pi, -8.0/9);
+ delta_t[i] = 4.5 * POW((Theta[i] - sin(Theta[i])), 2) /
+ POW((1.0 - cos(Theta[i])), 3) - 1.0;
beta[i] = sin(0.5*Theta[i]) * sin(0.5*Theta[i]);
little_d[i] = 0.75 * (Theta[i] - sin(Theta[i]));
if (delta_t[i] < density)
iinit = i;
- if (Theta[i] < M_PI) // Turnaround
+ if (Theta[i] < pi) // Turnaround
iturn = i;
}

- t_init = 5.38e8 * yr * pow((1+InitialRedshift) / 10.0, -1.5);
- z_vir = pow((Theta[iinit] - sin(Theta[iinit])) / (2*M_PI), 2.0/3) *
+ t_init = 5.38e8 * yr_s * POW((1+InitialRedshift) / 10.0, -1.5);
+ z_vir = POW((Theta[iinit] - sin(Theta[iinit])) / (2*pi), 2.0/3.0) *
(1 + InitialRedshift) - 1.0;
- t_vir = 5.38e8 * yr * pow((1+z_vir) / 10.0, -1.5);
+ t_vir = 5.38e8 * yr_s * POW((1+z_vir) / 10.0, -1.5);
t_ta = 0.5 * t_vir;
- t_i = t_ta * pow(delta_i, 1.5) / (3*M_PI/4);
- rho_ci = 1.0 / (6 * M_PI * Grav * t_i*t_i);
- r_i = pow(mass / (4*M_PI/3 * rho_ci), 1.0/3);
+ t_i = t_ta * POW(delta_i, 1.5) / (3*pi/4);
+ rho_ci = 1.0 / (6 * pi * GravConst * t_i*t_i);
+ r_i = POW(mass / (4*pi/3.0 * rho_ci), 1.0/3.0);

- dm_rotvel_corr = 2 * pow(2*r_init/r_i, 2);
+ dm_rotvel_corr = 2 * POW(2*r_init/r_i, 2);
// dm_rotvel_corr = 1;

- r_ta_now = pow((3*M_PI/4), -8.0/9) * pow(delta_i, 1.0/3) * r_i
- * pow(t_init / t_i, 8.0/9);
+ r_ta_now = POW((3*pi/4), -8.0/9) * POW(delta_i, 1.0/3.0) * r_i
+ * POW(t_init / t_i, 8.0/9);

for (i = 0; i < NTHETA; i++) {
local_Vr[NTHETA-i-1] = (Lambda[i] * r_ta_now / t_init) *
sin(Theta[i]) * (Theta[i] - sin(Theta[i])) /
- pow(1.0 - cos(Theta[i]), 2);
+ POW(1.0 - cos(Theta[i]), 2);
Vl = local_Vr[NTHETA-i-1] / (r_ta_now / t_init);

// Subtract Hubble flow and convert to code units
- local_Vr[NTHETA-i-1] = (Vl - (2.0/3)*Lambda[i]) * r_ta_now / t_init;
+ local_Vr[NTHETA-i-1] = (Vl - (2.0/3.0)*Lambda[i]) * r_ta_now / t_init;
local_Vr[NTHETA-i-1] /= VelocityUnits;
Chi = 1.0 - 1.5 * Vl / Lambda[i];
- local_rho[NTHETA-i-1] = pow(little_d[i],2) / pow(beta[i],3) /
+ local_rho[NTHETA-i-1] = POW(little_d[i],2) / POW(beta[i],3) /
(1.0 + 3*Chi);
local_radius[NTHETA-i-1] = Lambda[i] * r_ta_now / LengthUnits;
if (local_radius[NTHETA-i-1] > r_init)
@@ -1211,7 +1206,7 @@
exterior_rho[i] = density;
}
// printf("%"ISYM": r = %"GSYM" pc, v_r = %"GSYM" km/s, rho = %"GSYM"\n",
-// i, radius_vr[i]*LengthUnits/3.086e18, Vr[i]*VelocityUnits/1e5,
+// i, radius_vr[i]*LengthUnits/pc, Vr[i]*VelocityUnits/1e5,
// exterior_rho[i]);
}


This diff is so big that we needed to truncate the remainder.

https://bitbucket.org/enzo/enzo-dev/commits/ff41cf66296c/
Changeset: ff41cf66296c
Branch: week-of-code
User: jwise77
Date: 2019-02-27 16:50:47+00:00
Summary: Adding option for bulk gas velocity in the ZeldovichPancake problem.
Added documentation for initial B-field, which was missing.
Affected #: 4 files

diff -r 48fbcd4ec204 -r ff41cf66296c doc/manual/source/parameters/index.rst
--- a/doc/manual/source/parameters/index.rst
+++ b/doc/manual/source/parameters/index.rst
@@ -4075,6 +4075,11 @@
x-axis)
``ZeldovichPancakeInitialTemperature`` (external)
Initial gas temperature. Units: degrees Kelvin. Default: 100
+``ZeldovichPancakeInitialGasVelocity`` (external)
+ Initial bulk gas velocity in the direction of the pancake
+ collapse. Units: km/s. Default: 0.0
+``ZeldovichPancakeInitialUniformBField`` (external)
+ Initial magnetic field. Units: Gauss. Default: 0.0 0.0 0.0
``ZeldovichPancakeOmegaBaryonNow`` (external)
Omega Baryon at redshift z=0; standard setting. Default: 1.0
``ZeldovichPancakeOmegaCDMNow`` (external)

diff -r 48fbcd4ec204 -r ff41cf66296c src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -1947,6 +1947,7 @@
float ZeldovichPancakeOmegaCDMNow,
float ZeldovichPancakeCollapseRedshift,
float ZeldovichPancakeInitialTemperature,
+ float ZeldovichPancakeInitialGasVelocity,
float ZeldovichPancakeInitialUniformBField[]);

/* 1D Pressureless Collapse: initialize grid. */

diff -r 48fbcd4ec204 -r ff41cf66296c src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
--- a/src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
+++ b/src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
@@ -40,13 +40,14 @@
float ZeldovichPancakeOmegaCDMNow,
float ZeldovichPancakeCollapseRedshift,
float ZeldovichPancakeInitialTemperature,
+ float ZeldovichPancakeInitialGasVelocity,
float ZeldovichPancakeInitialUniformBField[]
)
{
/* declarations */

float Amplitude, AmplitudeVel, kx, xLagrange, xEulerian, xEulerianOld;
- int dim, field, i, index;
+ int dim, field, i, j, k, index;

/* error check */

@@ -132,8 +133,10 @@
AmplitudeVel = -sqrt(2.0/3.0)*(1.0+ZeldovichPancakeCollapseRedshift) /
((1.0+InitialRedshift) * 2.0*pi);
kx = 2*pi/float(NumberOfZones);
-
+
/* set density, total energy and velocity in problem dimension */
+
+ float bulkv = ZeldovichPancakeInitialGasVelocity*km_cm / VelocityUnits;

for (i = 0; i < size; i++) {

@@ -174,16 +177,16 @@
if (DualEnergyFormalism)
BaryonField[iTE+1][i] = BaryonField[iTE][i];
if (HydroMethod != Zeus_Hydro)
- BaryonField[iTE][i] += 0.5*POW(AmplitudeVel*sin(kx*xEulerian),2);
+ BaryonField[iTE][i] += 0.5*POW(AmplitudeVel*sin(kx*xEulerian) + bulkv,2);

/* Set velocities (some may be set to zero later -- see below). */

if (HydroMethod == Zeus_Hydro) xEulerian -= 0.5; // only approximate
- BaryonField[vel][i] = AmplitudeVel * sin(kx*xEulerian);
+ BaryonField[vel][i] = AmplitudeVel * sin(kx*xEulerian) + bulkv;
if ( MaxVelocityIndex > 1)
- BaryonField[vel+1][i] = AmplitudeVel * sin(kx*xEulerian);
+ BaryonField[vel+1][i] = AmplitudeVel * sin(kx*xEulerian) + bulkv;
if ( MaxVelocityIndex > 2)
- BaryonField[vel+2][i] = AmplitudeVel * sin(kx*xEulerian);
+ BaryonField[vel+2][i] = AmplitudeVel * sin(kx*xEulerian) + bulkv;

if ( UseMHD ) {
BaryonField[iBx ][i] = ZeldovichPancakeInitialUniformBField[0];
@@ -236,6 +239,8 @@
ParticlePosition[dim][ipart] = pos0[dim] +
BaryonField[vel+dim][index] * InitialTimeInCodeUnits;
ParticleVelocity[dim][ipart] = BaryonField[vel+dim][index];
+ if (dim == ZeldovichPancakeDirection)
+ ParticleVelocity[dim][ipart] -= bulkv;
}
ipart++;
} // ENDFOR i

diff -r 48fbcd4ec204 -r ff41cf66296c src/enzo/ZeldovichPancakeInitialize.C
--- a/src/enzo/ZeldovichPancakeInitialize.C
+++ b/src/enzo/ZeldovichPancakeInitialize.C
@@ -79,6 +79,7 @@
float ZeldovichPancakeOmegaCDMNow = 0.0; // no dark matter
float ZeldovichPancakeCollapseRedshift = 1.0; // free parameter
float ZeldovichPancakeInitialTemperature = 100; // whatever
+ float ZeldovichPancakeInitialGasVelocity = 0.0; // km/s
float ZeldovichPancakeInitialUniformBField[MAX_DIMENSION]; // in Gauss

for (int dim = 0; dim < MAX_DIMENSION; dim++) {
@@ -106,6 +107,8 @@
&ZeldovichPancakeCollapseRedshift);
ret += sscanf(line, "ZeldovichPancakeInitialTemperature = %"FSYM,
&ZeldovichPancakeInitialTemperature);
+ ret += sscanf(line, "ZeldovichPancakeInitialGasVelocity = %"FSYM,
+ &ZeldovichPancakeInitialGasVelocity);
ret += sscanf(line, "ZeldovichPancakeInitialUniformBField = %"FSYM" %"FSYM" %"FSYM,
ZeldovichPancakeInitialUniformBField,
ZeldovichPancakeInitialUniformBField+1,
@@ -144,6 +147,7 @@
ZeldovichPancakeOmegaCDMNow,
ZeldovichPancakeCollapseRedshift,
ZeldovichPancakeInitialTemperature,
+ ZeldovichPancakeInitialGasVelocity,
ZeldovichPancakeInitialUniformBField
) == FAIL) {
ENZO_FAIL("Error in ZeldovichPancakeInitializeGrid.\n");


https://bitbucket.org/enzo/enzo-dev/commits/314b00183fe3/
Changeset: 314b00183fe3
Branch: week-of-code
User: jwise77
Date: 2019-02-27 17:00:47+00:00
Summary: Adding new ZeldovichPancake test problem with DM particles and a bulk
gas velocity.
Affected #: 4 files

diff -r ff41cf66296c -r 314b00183fe3 run/Cosmology/AMRZeldovichPancake_Streaming/AMRZeldovichPancake_Streaming.enzo
--- /dev/null
+++ b/run/Cosmology/AMRZeldovichPancake_Streaming/AMRZeldovichPancake_Streaming.enzo
@@ -0,0 +1,62 @@
+#
+# AMR PROBLEM DEFINITION FILE: Zeldovich Pancake (AMR version)
+#
+# define problem
+#
+ProblemType = 20 // Zeldovich pancake
+TopGridRank = 1
+TopGridDimensions = 16
+SelfGravity = 1 // gravity on
+TopGridGravityBoundary = 0 // Periodic BC for gravity
+LeftFaceBoundaryCondition = 3 // same for fluid
+RightFaceBoundaryCondition = 3
+
+
+StopCycle = 5
+#
+# problem parameters
+#
+ZeldovichPancakeCentralOffset = 0
+ZeldovichPancakeCollapseRedshift = 1
+ZeldovichPancakeInitialGasVelocity = 100
+ZeldovichPancakeOmegaBaryonNow = 0.15
+ZeldovichPancakeOmegaCDMNow = 0.85
+#
+# define cosmology parameters
+#
+ComovingCoordinates = 1 // Expansion ON
+CosmologyHubbleConstantNow = 0.5
+CosmologyComovingBoxSize = 64.0 // 64 Mpc/h
+CosmologyMaxExpansionRate = 0.01 //
+CosmologyInitialRedshift = 20 // start at z=20
+CosmologyFinalRedshift = 0
+CosmologyOmegaMatterNow = 1
+CosmologyOmegaLambdaNow = 0
+GravitationalConstant = 1 // this must be true for cosmology
+
+#
+# set I/O and stop/start parameters
+#
+dtDataDump = 100.0
+#
+# set hydro parameters
+#
+Gamma = 1.6667
+CourantSafetyNumber = 0.5
+PPMDiffusionParameter = 0 // diffusion off
+DualEnergyFormalism = 1 // use total & internal energy
+ConservativeInterpolation = 0
+FluxCorrection = 1
+InterpolationMethod = 1
+#
+# set grid refinement parameters
+#
+StaticHierarchy = 0 // dynamic hierarchy
+MaximumRefinementLevel = 2 // use up to 2 levels
+RefineBy = 4 // refinement factor
+CellFlaggingMethod = 2 3 // use mass & shock criteria for refinement
+MinimumOverDensityForRefinement = 1.5 // times the initial overdensity
+#
+# set some global parameters
+#
+MinimumEfficiency = 0.4 // better value for 1d than 0.2

diff -r ff41cf66296c -r 314b00183fe3 run/Cosmology/AMRZeldovichPancake_Streaming/AMRZeldovichPancake_Streaming.enzotest
--- /dev/null
+++ b/run/Cosmology/AMRZeldovichPancake_Streaming/AMRZeldovichPancake_Streaming.enzotest
@@ -0,0 +1,14 @@
+name = 'AMRZeldovichPancake_Streaming'
+answer_testing_script = None
+nprocs = 1
+runtime = 'short'
+hydro = True
+gravity = True
+AMR = True
+cosmology = True
+dimensionality = 1
+author = 'John Wise'
+max_time_minutes = 1
+fullsuite = True
+pushsuite = True
+quicksuite = True

diff -r ff41cf66296c -r 314b00183fe3 run/Cosmology/AMRZeldovichPancake_Streaming/notes.txt
--- /dev/null
+++ b/run/Cosmology/AMRZeldovichPancake_Streaming/notes.txt
@@ -0,0 +1,23 @@
+AMRZeldovichPancake
+----------
+
+(John Wise, July 2010)
+
+This test simulates a collapsing sinusoidal cosmological pertubation
+in one-dimension, commonly known as a Zel'dovich pancake. This
+problem tests both the hydrodynamics and gravity solvers and the
+implementation of cosmological expansion. The system will form a
+caustic in the center of the domain with a density and temperature
+peak. There should be a small dip in temperature at the center of the
+broad peak. In flat cosmology, there exists an analytical solution in
+the linear phase of collapse and is given in Zel'dovich (1970).
+
+This test runs to completion and creates 2 outputs -- the initial
+(z=20) and final (z=0) states. There are two levels of refinement by
+factors of 4. The finest resolution element is the same as the
+non-AMR version of this test. There are some small differences between
+Enzo v1.5 and v2.0 at the parent-child grid boundaries.
+
+(John Wise, February 2019)
+
+Modified to add DM particles and a bulk gas velocity

diff -r ff41cf66296c -r 314b00183fe3 run/Cosmology/AMRZeldovichPancake_Streaming/scripts.py
--- /dev/null
+++ b/run/Cosmology/AMRZeldovichPancake_Streaming/scripts.py
@@ -0,0 +1,21 @@
+from matplotlib import use; use('Agg')
+import matplotlib.pyplot as plt
+from yt.mods import *
+
+pf = load('DD0001/data0001')
+ray = pf.h.ortho_ray(0,(0,0))
+
+plt.subplot(311)
+plt.semilogy(ray['x'], ray['Density'], c='k', ls='-', marker='.', ms=5)
+plt.ylabel(r'Density (cm$^{-3}$)')
+
+plt.subplot(312)
+plt.semilogy(ray['x'], ray['Temperature'], c='k', ls='-', marker='.', ms=5)
+plt.ylabel('Temperature (K)')
+
+plt.subplot(313)
+plt.plot(ray['x'], 1e-5*ray['x-velocity'], c='k', ls='-', marker='.', ms=5)
+plt.ylabel('Velocity (km/s)')
+plt.xlabel('x')
+
+plt.savefig('AMRZeldovichPancake.png')


https://bitbucket.org/enzo/enzo-dev/commits/e63161949122/
Changeset: e63161949122
Branch: week-of-code
User: jwise77
Date: 2019-02-28 12:32:12+00:00
Summary: Cell width not calculated correctly for Zeldovich particles in >1 dimensions
Affected #: 1 file

diff -r 314b00183fe3 -r e63161949122 src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
--- a/src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
+++ b/src/enzo/Grid_ZeldovichPancakeInitializeGrid.C
@@ -221,8 +221,8 @@

int ipart = 0;
FLOAT pos0[MAX_DIMENSION], dr[MAX_DIMENSION];
- for (dim = 0 ; dim < GridRank; dim++)
- dr[dim] = (GridRank > dim) ? CellWidth[dim][0] :
+ for (dim = 0 ; dim < MAX_DIMENSION; dim++)
+ dr[dim] = (dim < GridRank) ? CellWidth[dim][0] :
DomainRightEdge[dim] - DomainLeftEdge[dim];

for (k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {


https://bitbucket.org/enzo/enzo-dev/commits/f96d4827b734/
Changeset: f96d4827b734
Branch: week-of-code
User: jwise77
Date: 2019-03-28 14:55:29+00:00
Summary: Adding parameter in new ZP test problem parameter file output
Affected #: 2 files

diff -r e63161949122 -r f96d4827b734 src/enzo/ZeldovichPancakeInitialize.C
--- a/src/enzo/ZeldovichPancakeInitialize.C
+++ b/src/enzo/ZeldovichPancakeInitialize.C
@@ -195,8 +195,10 @@
ZeldovichPancakeOmegaCDMNow);
fprintf(Outfptr, "ZeldovichPancakeCollapseRedshift = %"FSYM"\n",
ZeldovichPancakeCollapseRedshift);
- fprintf(Outfptr, "ZeldovichPancakeInitialTemperature = %"FSYM"\n\n",
+ fprintf(Outfptr, "ZeldovichPancakeInitialTemperature = %"FSYM"\n",
ZeldovichPancakeInitialTemperature);
+ fprintf(Outfptr, "ZeldovichPancakeInitialGasVelocity = %"FSYM"\n",
+ ZeldovichPancakeInitialGasVelocity);

for (int dim = 0; dim < MAX_DIMENSION; dim++)
ZeldovichPancakeInitialUniformBField[dim] *= MagneticUnits;

diff -r e63161949122 -r f96d4827b734 src/enzo/libconfig/libconfig.c
--- a/src/enzo/libconfig/libconfig.c
+++ b/src/enzo/libconfig/libconfig.c
@@ -33,9 +33,11 @@

#include <locale.h>

-#define HAVE_XLOCALE_H
+//#define HAVE_XLOCALE_H
#ifdef HAVE_XLOCALE_H
#include <xlocale.h>
+#else
+#include <locale.h>
#endif

#include <stdlib.h>


https://bitbucket.org/enzo/enzo-dev/commits/98b7c039fbd3/
Changeset: 98b7c039fbd3
Branch: week-of-code
User: jwise77
Date: 2019-03-29 13:45:46+00:00
Summary: Merging to tip
Affected #: 18 files

diff -r f96d4827b734 -r 98b7c039fbd3 .hgtags
--- a/.hgtags
+++ b/.hgtags
@@ -23,3 +23,4 @@
1bd78f067e9168391ac76d6234d5af1e2bd9273b gold-standard-v7
cc6d4f93e6ee6df2c5fdd8107f135a12410a826e gold-standard-v8
d4cee1981c1fb6e9f4d812ac0ee99bdcba7be7f5 gold-standard-v9
+69a7950b8040128440aa68daa79802aa8e21a5b1 gold-standard-v10

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/Grid_PrepareGreensFunction.C
--- a/src/enzo/Grid_PrepareGreensFunction.C
+++ b/src/enzo/Grid_PrepareGreensFunction.C
@@ -53,13 +53,13 @@

/* Set the constant to be used. */

- float GravitationalConstant;
+ float GravConst_factor;
if (GridRank == 3)
- GravitationalConstant = -GravConst/(4.0*pi);
+ GravConst_factor = -GravitationalConstant/(4.0*pi);
if (GridRank == 2)
- GravitationalConstant = -GravConst*0.5/pi;
+ GravConst_factor = -GravitationalConstant*0.5/pi;
if (GridRank == 1)
- GravitationalConstant = -GravConst*0.5;
+ GravConst_factor = -GravitationalConstant*0.5;

/* Set Greens' function. */

@@ -78,12 +78,12 @@
r = max(r, GravitatingMassFieldCellSize);
r *= GravitatingMassFieldCellSize;
if (GridRank == 3)
- PotentialField[n] = GravitationalConstant/r;
+ PotentialField[n] = GravConst_factor/r;
if (GridRank == 2)
- PotentialField[n] = GravitationalConstant*log(r);
+ PotentialField[n] = GravConst_factor*log(r);
if (GridRank == 1)

- PotentialField[n] = GravitationalConstant*r;
+ PotentialField[n] = GravConst_factor*r;

}
}

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/PrepareIsolatedGreensFunction.C
--- a/src/enzo/PrepareIsolatedGreensFunction.C
+++ b/src/enzo/PrepareIsolatedGreensFunction.C
@@ -91,13 +91,13 @@

/* Set the constant to be used. */

- float GravitationalConstant;
+ float GravConst_factor;
if (GridRank == 3)
- GravitationalConstant = -GravConst*RealCellVolume/(4.0*pi);
+ GravConst_factor = -GravitationalConstant*RealCellVolume/(4.0*pi);
if (GridRank == 2)
- GravitationalConstant = GravConst*RealCellVolume/(2.0*pi);
+ GravConst_factor = GravitationalConstant*RealCellVolume/(2.0*pi);
if (GridRank == 1)
- GravitationalConstant = GravConst*RealCellVolume/2.0;
+ GravConst_factor = GravitationalConstant*RealCellVolume/2.0;

/* Set Greens' function. */

@@ -122,11 +122,11 @@
r = max(r, 0.38*RealCellWidth[0]);
// r *= GravitatingMassFieldCellSize;
if (GridRank == 3)
- GreensFunction->Data[n] = GravitationalConstant/r;
+ GreensFunction->Data[n] = GravConst_factor/r;
if (GridRank == 2)
- GreensFunction->Data[n] = GravitationalConstant*log(r);
+ GreensFunction->Data[n] = GravConst_factor*log(r);
if (GridRank == 1)
- GreensFunction->Data[n] = GravitationalConstant*r;
+ GreensFunction->Data[n] = GravConst_factor*r;
}
}
}

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/pop3_maker.F
--- a/src/enzo/pop3_maker.F
+++ b/src/enzo/pop3_maker.F
@@ -149,7 +149,8 @@
c 4) t_cool < t_free-fall (if T < 1.1e4 skip this check)
c
dtot = ( d(i,j,k) + dm(i,j,k) )*d1
- tdyn = sqrt(3._RKIND*pi_val/32._RKIND/GravConst*dtot)/t1
+ tdyn = sqrt(3._RKIND*pi_val/32._RKIND/
+ & GravConst*dtot)/t1
c if (tdyn .lt. cooltime(i,j,k) .and.
c & temp(i,j,k) .gt. 1.1e4)
c & write(40,*) '1',tdyn,cooltime(i,j,k),temp(i,j,k)

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/projplane.F
--- a/src/enzo/projplane.F
+++ b/src/enzo/projplane.F
@@ -345,7 +345,7 @@
if (ifield .eq. 3) plane(i,j) = plane(i,j) +
& sample1**1.5_RKIND *grid2(i0,n,k0)**2 * weight1
if (ifield .eq. 4) plane(i,j) = plane(i,j) +
- & sample1 * grid2(n,j0,k0) * weight1
+ & sample1 * grid2(i0,n,k0) * weight1
c
enddo
enddo
@@ -381,7 +381,7 @@
if (ifield .eq. 3) plane(i,j) = plane(i,j) +
& sample1**1.5_RKIND *grid2(i0,j0,n)**2 * weight1
if (ifield .eq. 4) plane(i,j) = plane(i,j) +
- & sample1 * grid2(n,j0,k0) * weight1
+ & sample1 * grid2(i0,j0,n) * weight1
c
enddo
enddo

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/solve_rate.F
--- a/src/enzo/solve_rate.F
+++ b/src/enzo/solve_rate.F
@@ -127,7 +127,8 @@
tbase1 = utim
xbase1 = uxyz/(aye*uaye) ! uxyz is [x]*a = [x]*[a]*a' '
dbase1 = urho*(aye*uaye)**3 ! urho is [dens]/a^3 = [dens]/([a]*a')^3 '
- coolunit = (uaye**5 * xbase1**2 * mass_h**2) / (tbase1**3 * dbase1)
+ coolunit = (uaye**5 * xbase1**2 * mass_h**2) /
+ & (tbase1**3 * dbase1)
c
c Set log values of start and end of lookup tables
c

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/solve_rate_cool.F
--- a/src/enzo/solve_rate_cool.F
+++ b/src/enzo/solve_rate_cool.F
@@ -270,7 +270,8 @@
tbase1 = utim
xbase1 = uxyz/(aye*uaye) ! uxyz is [x]*a = [x]*[a]*a' '
dbase1 = urho*(aye*uaye)**3 ! urho is [dens]/a^3 = [dens]/([a]*a')^3 '
- coolunit = (uaye**5 * xbase1**2 * mass_h**2) / (tbase1**3 * dbase1)
+ coolunit = (uaye**5 * xbase1**2 * mass_h**2) /
+ & (tbase1**3 * dbase1)
uvel = uxyz / utim
c chunit = (7.17775e-12_RKIND)/(2._RKIND*uvel*uvel*mass_h) ! 4.5 eV per H2 formed
chunit = (ev2erg)/(2._RKIND*uvel*uvel*mass_h) ! 1 eV per H2 formed

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_feedback_pn_snia.F
--- a/src/enzo/star_feedback_pn_snia.F
+++ b/src/enzo/star_feedback_pn_snia.F
@@ -180,7 +180,8 @@
c
c EdotSNII = eSNII/0.150081 * MdotSNII*(c_light/v1)**2
if (imetalSNIa .eq. 1) then
- EdotSNIa = eSNIa/1.234e-3_RKIND * MdotSNIa*(c_light/v1)**2
+ EdotSNIa = eSNIa/1.234e-3_RKIND * MdotSNIa*
+ & (c_light/v1)**2
else
EdotSNIa = 0._RKIND
endif

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker1.F
--- a/src/enzo/star_maker1.F
+++ b/src/enzo/star_maker1.F
@@ -147,7 +147,8 @@
c 4) t_cool < t_free-fall (if T < 1.1e4 skip this check)
c
dtot = ( d(i,j,k) + dm(i,j,k) )*d1
- tdyn = sqrt(3._RKIND*pi_val/32._RKIND/GravConst/dtot)/t1
+ tdyn = sqrt(3._RKIND*pi_val/32._RKIND/
+ & GravConst/dtot) / t1

if (tdyn .lt. cooltime(i,j,k) .and.
& temp(i,j,k) .gt. 1.1e4_RKIND) goto 10
@@ -160,8 +161,10 @@
c
bmass = d(i,j,k)*dble(d1)*dble(x1*dx)**3 / SolarMass
isosndsp2 = sndspdC * temp(i,j,k)
- jeanmass = pi_val/(6._RKIND*sqrt(d(i,j,k)*dble(d1))) *
- & dble(pi_val * isosndsp2 / GravConst)**1.5_RKIND / SolarMass
+ jeanmass = pi_val/(6._RKIND*
+ & sqrt(d(i,j,k)*dble(d1)))*
+ & dble(pi_val * isosndsp2 /
+ & GravConst)**1.5_RKIND / SolarMass

if (bmass .lt. jeanmass) goto 10
c

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker10.F
--- a/src/enzo/star_maker10.F
+++ b/src/enzo/star_maker10.F
@@ -146,7 +146,8 @@
c 4) t_cool < t_free-fall (if T < 1.1e4 skip this check)
c
dtot = ( d(i,j,k) + dm(i,j,k) )*d1
- tdyn = sqrt(3.0*pi_val/32.0/GravConst/dtot)/t1
+ tdyn = sqrt(3.0*pi_val/32.0/
+ & GravConst/dtot)/t1

if (tdyn .lt. cooltime(i,j,k) .and.
& temp(i,j,k) .gt. 1.1e4) goto 10
@@ -159,8 +160,10 @@
c
bmass = d(i,j,k)*dble(d1)*dble(x1*dx)**3 / SolarMass
isosndsp2 = sndspdC * temp(i,j,k)
- jeanmass = pi_val/(6.0*sqrt(d(i,j,k)*dble(d1))) *
- & dble(pi_val * isosndsp2 / GravConst)**1.5 / SolarMass
+ jeanmass = pi_val/(6.0*
+ & sqrt(d(i,j,k)*dble(d1))) *
+ & dble(pi_val * isosndsp2 /
+ & GravConst)**1.5 / SolarMass

c
c THIS IS COMMENTED OUT - NO JEANS MASS CRITERION IN THIS ALGORITHM!!!
@@ -457,7 +460,7 @@
c Calculate how much of the star formation in this
c timestep would have gone into supernova energy.
c
- energy = sn_param * mform * (c_light/v1)**2 / distcells
+ energy = sn_param * mform * (c_light/v1)**2/distcells
c
c (aug 7/00 addition: share ejected energy between gas
c being ejected and that remaining in star particle)

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker2.F
--- a/src/enzo/star_maker2.F
+++ b/src/enzo/star_maker2.F
@@ -153,7 +153,8 @@
c 4) t_cool < t_free-fall (if T < 1.1e4 skip this check)
c
dtot = ( d(i,j,k) + dm(i,j,k) )*d1
- tdyn = sqrt(3._RKIND*pi_val/32._RKIND/GravConst/dtot)/t1
+ tdyn = sqrt(3._RKIND*pi_val/32._RKIND/
+ & GravConst/dtot)/t1

if (tdyn .lt. cooltime(i,j,k) .and.
& temp(i,j,k) .gt. 1.1e4_RKIND) goto 10
@@ -167,9 +168,10 @@
bmass = d(i,j,k)*dble(d1)*dble(x1*dx)**3 / SolarMass
isosndsp2 = sndspdC * temp(i,j,k)
if (usejeans .eq. 1) then
- jeanmass = pi_val/(6._RKIND*sqrt(d(i,j,k)*dble(d1)))*
- & dble(pi_val * isosndsp2 / GravConst)**1.5_RKIND /
- & SolarMass
+ jeanmass = pi_val/(6._RKIND*
+ & sqrt(d(i,j,k)*dble(d1)))*
+ & dble(pi_val * isosndsp2 /
+ & GravConst)**1.5_RKIND / SolarMass

if (bmass .lt. jeanmass) goto 10
endif
@@ -462,7 +464,7 @@
c Calculate how much of the star formation in this
c timestep would have gone into supernova energy.
c
- energy = sn_param * mform * (c_light/v1)**2 / distcells
+ energy = sn_param * mform * (c_light/v1)**2/distcells
c
c (aug 7/00 addition: share ejected energy between gas
c being ejected and that remaining in star particle)

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker3.F
--- a/src/enzo/star_maker3.F
+++ b/src/enzo/star_maker3.F
@@ -154,7 +154,8 @@
c 4) t_cool < t_free-fall (if T < 1.1e4 skip this check)
c
dtot = ( d(i,j,k) + dm(i,j,k) )*d1
- tdyn = sqrt(3._RKIND*pi_val/32._RKIND/GravConst/dtot)/t1
+ tdyn = sqrt(3._RKIND*pi_val/32._RKIND/
+ & GravConst/dtot)/t1

if (tdyn .lt. cooltime(i,j,k) .and.
& temp(i,j,k) .gt. 1.1e4_RKIND) goto 10
@@ -167,8 +168,10 @@
c
bmass = d(i,j,k)*dble(d1)*dble(x1*dx)**3 / SolarMass
isosndsp2 = sndspdC * temp(i,j,k)
- jeanmass = pi_val/(6._RKIND*sqrt(d(i,j,k)*dble(d1))) *
- & dble(pi_val * isosndsp2 / GravConst)**1.5_RKIND / SolarMass
+ jeanmass = pi_val/(6._RKIND*
+ & sqrt(d(i,j,k)*dble(d1))) *
+ & dble(pi_val * isosndsp2 /
+ & GravConst)**1.5_RKIND / SolarMass

c
c THIS IS COMMENTED OUT - NO JEANS MASS CRITERION IN THIS ALGORITHM!!!
@@ -474,7 +477,7 @@
c Calculate how much of the star formation in this
c timestep would have gone into supernova energy.
c
- energy = sn_param * mform * (c_light/v1)**2 / distcells
+ energy = sn_param * mform * (c_light/v1)**2/distcells
c
c (aug 7/00 addition: share ejected energy between gas
c being ejected and that remaining in star particle)

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker3mom.F
--- a/src/enzo/star_maker3mom.F
+++ b/src/enzo/star_maker3mom.F
@@ -139,7 +139,8 @@
c 4) t_cool < t_free-fall (if T < 1.1e4 skip this check)
c
dtot = ( d(i,j,k) + dm(i,j,k) )*dunits
- tdyn = sqrt(3._RKIND*pi_val/32._RKIND/GravConst/dtot)/t1
+ tdyn = sqrt(3._RKIND*pi_val/32._RKIND/
+ & GravConst/dtot)/t1

if (tdyn .lt. cooltime(i,j,k) .and.
& temp(i,j,k) .gt. 1.1e4_RKIND) goto 10
@@ -152,8 +153,10 @@
c
bmass = d(i,j,k)*dble(dunits)*dble(x1*dx)**3 / SolarMass
isosndsp2 = sndspdC * temp(i,j,k)
- jeanmass = pi_val/(6._RKIND*sqrt(d(i,j,k)*dble(dunits))) *
- & dble(pi_val * isosndsp2 / GravConst)**1.5_RKIND / SolarMass
+ jeanmass = pi_val/(6._RKIND*
+ & sqrt(d(i,j,k)*dble(dunits))) *
+ & dble(pi_val * isosndsp2 /
+ & GravConst)**1.5_RKIND / SolarMass

c
c THIS IS COMMENTED OUT - NO JEANS MASS CRITERION IN THIS ALGORITHM!!!

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker5.F
--- a/src/enzo/star_maker5.F
+++ b/src/enzo/star_maker5.F
@@ -163,8 +163,8 @@
c coolrate is negative, gas is cooling (which is the opposite of how
c it's defined in springel & Hernquist 2003, MNRAS, 339, 289, eqn. 16
y = -1._RKIND*tstar*coolrate(i,j,k)/(d(i,j,k)*dble(d1))/
- & (beta*usn-(1._RKIND-beta)*1.5_RKIND*kboltz*temp(i,j,k)/
- & 0.6_RKIND/mass_h) !Eq(16)
+ & (beta*usn-(1._RKIND-beta)*1.5_RKIND*kboltz*
+ & temp(i,j,k)/ 0.6_RKIND/mass_h) !Eq(16)

c Calculate the fraction of mass in cold clouds
if (y .le. 0._RKIND) then

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker7.F
--- a/src/enzo/star_maker7.F
+++ b/src/enzo/star_maker7.F
@@ -191,7 +191,8 @@
c 4) t_cool < t_free-fall (if T < 1.1e4 skip this check)
c
dtot = ( d(i,j,k) + dm(i,j,k) )*d1
- tdyn = sqrt(3._RKIND*pi_val/32._RKIND/GravConst/dtot)/t1
+ tdyn = sqrt(3._RKIND*pi_val/32._RKIND/
+ & GravConst/dtot)/t1

if (tdyn .lt. cooltime(i,j,k) .and.
& temp(i,j,k) .gt. 1.1e4_RKIND) goto 10
@@ -207,8 +208,10 @@
c
bmass = d(i,j,k)*dble(d1)*dble(x1*dx)**3 / SolarMass
isosndsp2 = sndspdC * temp(i,j,k)
- jeanmass = pi_val/(6._RKIND*sqrt(d(i,j,k)*dble(d1))) *
- & dble(pi_val * isosndsp2 / GravConst)**1.5_RKIND / SolarMass
+ jeanmass = pi_val/(6._RKIND*
+ & sqrt(d(i,j,k)*dble(d1))) *
+ & dble(pi_val * isosndsp2 /
+ & GravConst)**1.5_RKIND / SolarMass
c
c if (bmass .lt. jeanmass) goto 10
c

diff -r f96d4827b734 -r 98b7c039fbd3 src/enzo/star_maker_ssn.F
--- a/src/enzo/star_maker_ssn.F
+++ b/src/enzo/star_maker_ssn.F
@@ -139,11 +139,14 @@
c

dtot = ( d(i,j,k) + dm(i,j,k) )*d1
- tdyn = sqrt(3._RKIND*pi_val/32._RKIND/GravConst/dtot)/t1
+ tdyn = sqrt(3._RKIND*pi_val/32._RKIND/
+ & GravConst/dtot)/t1
bmass = d(i,j,k)*dble(d1)*dble(x1*dx)**3 / SolarMass
isosndsp2 = sndspdC * temp(i,j,k)
- jeanmass = pi_val/(6._RKIND*sqrt(d(i,j,k)*dble(d1))) *
- & dble(pi_val * isosndsp2 / GravConst)**1.5_RKIND / SolarMass
+ jeanmass = pi_val/(6._RKIND*
+ & sqrt(d(i,j,k)*dble(d1))) *
+ & dble(pi_val * isosndsp2 /
+ & GravConst)**1.5_RKIND / SolarMass

if (jeanmass > bmass) then
goto 10
@@ -423,13 +426,15 @@

num = d(i,j,k) * d1 / mu / mass_h
stromgren_radius = (
- & (3._RKIND*s49_tot(n)*1d49)/(4._RKIND*pi_val*alpha*num**2)
- & )**(1._RKIND/3._RKIND)
- stromgren_volume =
- & (4._RKIND/3._RKIND)*pi_val*stromgren_radius**3
+ & (3._RKIND*s49_tot(n)*1d49)/
+ & (4._RKIND*pi_val*alpha*num**2)
+ & )**(1._RKIND/3._RKIND)
+ stromgren_volume = (4._RKIND/3._RKIND)*
+ & pi_val*stromgren_radius**3
cell_volume = dx*x1*dx*x1*dx*x1

- ionized = kboltz*1e4 / mu / mass_h / x1**2 * t1**2
+ ionized = kboltz*1e4 / mu / mass_h /
+ & x1**2 * t1**2
if (stromgren_volume .le. cell_volume) then
ionized = ionized * stromgren_volume/cell_volume
endif

diff -r f96d4827b734 -r 98b7c039fbd3 src/performance_tools/performance_tools.py
--- a/src/performance_tools/performance_tools.py
+++ b/src/performance_tools/performance_tools.py
@@ -40,6 +40,10 @@
"""
Checks to see if an object is listlike (but not a string)
"""
+ try:
+ basestring
+ except NameError:
+ basestring = str
return not isinstance(obj, basestring)

def preserve_extrema(extrema, xdata, ydata):
@@ -124,19 +128,19 @@
>>> x_smooth = smooth(x_noisy)
"""
if x.ndim != 1:
- raise ValueError, "smooth only accepts 1 dimension arrays."
+ raise ValueError("smooth only accepts 1 dimension arrays.")

if x.size < window_len:
- raise ValueError, "Input vector needs to be bigger than window size."
+ raise ValueError("Input vector needs to be bigger than window size.")

if window_len<3:
return x

if window_len%2 == 0:
- raise ValueError, "smooth requires an odd-valued window_len."
+ raise ValueError("smooth requires an odd-valued window_len.")

if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
- raise ValueError, "Window is one of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
+ raise ValueError("Window is one of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'")

s = np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]
if window == 'flat': # moving average
@@ -273,7 +277,8 @@
line_value = np.array([cycle] + line_list[1:],dtype='float64')
line_value = np.nan_to_num(line_value) # error checking

- data[line_key][i] = line_value
+ # new numpy requires this to be cast as a tuple
+ data[line_key][i] = to_tuple(line_value)

### Make sure all cycles are set for all keys, even those that
### didn't output every cycle
@@ -647,7 +652,7 @@

### Filter out field_label's for which the data is all zeros
### (e.g. Group_WriteAllData if no cycles had outputs)
- field_label = filter(lambda x: sum(data[x]["Min Time"] + data[x]["Max Time"]) > 0.0, field_label)
+ field_label = [x for x in field_label if sum(data[x]["Min Time"] + data[x]["Max Time"]) > 0.0]

num_fields = len(field_label)
field_label.sort()
@@ -866,7 +871,7 @@
### there are including any that were defined in the original
### field_label argument.
if repeated_field:
- key_list = data.keys()
+ key_list = list(data.keys())
if repeated_field == "All":
field_label = key_list

@@ -887,7 +892,7 @@

### Filter out field_label's for which the data is all zeros
### (e.g. Group_WriteAllData if no cycles had outputs)
- field_label = filter(lambda x: sum(data[x]["Min Time"] + data[x]["Max Time"]) > 0.0, field_label)
+ field_label = [x for x in field_label if sum(data[x]["Min Time"] + data[x]["Max Time"]) > 0.0]

### Loop through the y datasets to figure out the extrema
for i in range(len(field_label)):
@@ -987,6 +992,14 @@
pl.savefig(filename)
pl.clf()

+def to_tuple(a):
+ """
+ Convert an array to a tuple
+ """
+ try:
+ return tuple(to_tuple(i) for i in a)
+ except TypeError:
+ return a

### *** DEFAULT COMMAND LINE BEHAVIOR ***


diff -r f96d4827b734 -r 98b7c039fbd3 test.sh
--- a/test.sh
+++ b/test.sh
@@ -33,7 +33,7 @@
export MATPLOTLIBRC=$HOME

export ENZOTEST_DIR=$HOME/enzo_test
-export GOLD_STANDARD_TAG="gold-standard-v9"
+export GOLD_STANDARD_TAG="gold-standard-v10"

# Build the gold standard version.
cd $BITBUCKET_CLONE_DIR

Repository URL: https://bitbucket.org/enzo/enzo-dev/

--

This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.