diff --git a/configure.ac b/configure.ac
index 960b098a6f5917caf428c84684415ee17a44f81a..6558696ba25ed72f1789f5149b0ae022427c4610 100644
--- a/configure.ac
+++ b/configure.ac
@@ -38,47 +38,6 @@ AX_CHECK_ENABLE_DEBUG
AC_PROG_CC
AM_PROG_CC_C_O
-
-# Master subgrid options
-# If you add a restriction (e.g. no cooling, chemistry or hydro)
-# you will need to check for overwrite after reading the additional options.
-# As an example for this, see the call to AC_ARG_WITH for cooling.
-AC_ARG_WITH([subgrid],
- [AS_HELP_STRING([--with-subgrid=<subgrid>],
- [Master switch for subgrid methods. Inexperienced user should start from here @<:@none, GEAR, EAGLE default: none@:>@]
- )],
- [with_subgrid="$withval"],
- [with_subgrid=none]
-)
-
-# Default values
-with_subgrid_cooling=none
-with_subgrid_chemistry=none
-with_subgrid_hydro=none
-
-case "$with_subgrid" in
- yes)
- AC_MSG_ERROR([Invalid option. A subgrid model must be chosen.])
- ;;
- none)
- ;;
- GEAR)
- with_subgrid_cooling=grackle
- with_subgrid_chemistry=GEAR
- with_subgrid_hydro=gadget2
- ;;
- EAGLE)
- with_subgrid_cooling=EAGLE
- with_subgrid_chemistry=EAGLE
- with_subgrid_hydro=gadget2
- ;;
- *)
- AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
- ;;
-esac
-
-
-
# If debug is selected then we also define SWIFT_DEVELOP_MODE to control
# any developer code options.
if test "x$ax_enable_debug" != "xno"; then
@@ -594,42 +553,6 @@ AH_VERBATIM([__STDC_FORMAT_MACROS],
#define __STDC_FORMAT_MACROS 1
#endif])
-# Check for grackle.
-have_grackle="no"
-AC_ARG_WITH([grackle],
- [AS_HELP_STRING([--with-grackle=PATH],
- [root directory where grackle is installed @<:@yes/no@:>@]
- )],
- [with_grackle="$withval"],
- [with_grackle="no"]
-)
-if test "x$with_grackle" != "xno"; then
- AC_PROG_FC
- AC_FC_LIBRARY_LDFLAGS
- if test "x$with_grackle" != "xyes" -a "x$with_grackle" != "x"; then
- GRACKLE_LIBS="-L$with_grackle/lib -lgrackle"
- GRACKLE_INCS="-I$with_grackle/include"
- else
- GRACKLE_LIBS="-lgrackle"
- GRACKLE_INCS=""
- fi
-
- have_grackle="yes"
-
- AC_CHECK_LIB(
- [grackle],
- [initialize_chemistry_data],
- [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.])
- AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
- ],
- [AC_MSG_ERROR(Cannot find grackle library!)],
- [$GRACKLE_LIBS $GRACKLE_INCS $FCLIBS]
- )
-fi
-AC_SUBST([GRACKLE_LIBS])
-AC_SUBST([GRACKLE_INCS])
-AM_CONDITIONAL([HAVEGRACKLE],[test -n "$GRACKLE_LIBS"])
-
# Check for FFTW. We test for this in the standard directories by default,
# and only disable if using --with-fftw=no or --without-fftw. When a value
# is given GSL must be found.
@@ -989,6 +912,46 @@ fi
# Various package configuration options.
+# Master subgrid options
+# If you add a restriction (e.g. no cooling, chemistry or hydro)
+# you will need to check for overwrite after reading the additional options.
+# As an example for this, see the call to AC_ARG_WITH for cooling.
+AC_ARG_WITH([subgrid],
+ [AS_HELP_STRING([--with-subgrid=<subgrid>],
+ [Master switch for subgrid methods. Inexperienced user should start from here @<:@none, GEAR, EAGLE default: none@:>@]
+ )],
+ [with_subgrid="$withval"],
+ [with_subgrid=none]
+)
+
+# Default values
+with_subgrid_cooling=none
+with_subgrid_chemistry=none
+with_subgrid_hydro=none
+
+case "$with_subgrid" in
+ yes)
+ AC_MSG_ERROR([Invalid option. A subgrid model must be chosen.])
+ ;;
+ none)
+ ;;
+ GEAR)
+ with_subgrid_cooling=grackle
+ with_subgrid_chemistry=GEAR
+ with_subgrid_hydro=gadget2
+ ;;
+ EAGLE)
+ with_subgrid_cooling=EAGLE
+ with_subgrid_chemistry=EAGLE
+ with_subgrid_hydro=gadget2
+ ;;
+ *)
+ AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
+ ;;
+esac
+
+
+
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
@@ -1193,6 +1156,42 @@ case "$with_riemann" in
;;
esac
+# Check for grackle.
+have_grackle="no"
+AC_ARG_WITH([grackle],
+ [AS_HELP_STRING([--with-grackle=PATH],
+ [root directory where grackle is installed @<:@yes/no@:>@]
+ )],
+ [with_grackle="$withval"],
+ [with_grackle="no"]
+)
+if test "x$with_grackle" != "xno"; then
+ AC_PROG_FC
+ AC_FC_LIBRARY_LDFLAGS
+ if test "x$with_grackle" != "xyes" -a "x$with_grackle" != "x"; then
+ GRACKLE_LIBS="-L$with_grackle/lib -lgrackle"
+ GRACKLE_INCS="-I$with_grackle/include"
+ else
+ GRACKLE_LIBS="-lgrackle"
+ GRACKLE_INCS=""
+ fi
+
+ have_grackle="yes"
+
+ AC_CHECK_LIB(
+ [grackle],
+ [initialize_chemistry_data],
+ [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.])
+ AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
+ ],
+ [AC_MSG_ERROR(Cannot find grackle library!)],
+ [$GRACKLE_LIBS $GRACKLE_INCS $FCLIBS]
+ )
+fi
+AC_SUBST([GRACKLE_LIBS])
+AC_SUBST([GRACKLE_INCS])
+AM_CONDITIONAL([HAVEGRACKLE],[test -n "$GRACKLE_LIBS"])
+
# Cooling function
AC_ARG_WITH([cooling],
[AS_HELP_STRING([--with-cooling=<function>],
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 2b177a8a0d6abff8b2d7e83cb23b2950e3da2efe..cba52250ccc37f50ed130c70d8a5039d8c786474 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -1991,6 +1991,7 @@ INCLUDE_FILE_PATTERNS =
PREDEFINED = "__attribute__(x)= "
PREDEFINED += HAVE_HDF5
+PREDEFINED += HAVE_FFTW
PREDEFINED += WITH_MPI
PREDEFINED += WITH_VECTORIZATION
PREDEFINED += __GNUC__
diff --git a/doc/RTD/source/GettingStarted/parameter_file.rst b/doc/RTD/source/GettingStarted/parameter_file.rst
index 32a4f1220d0dc1c65074cc14b53d2c4edd666a67..cecf2fa085b23946b8bcbebd324372b01f4940c2 100644
--- a/doc/RTD/source/GettingStarted/parameter_file.rst
+++ b/doc/RTD/source/GettingStarted/parameter_file.rst
@@ -1,6 +1,8 @@
.. Parameter File
Loic Hausammann, 1 june 2018
+.. _Parameter_File_label:
+
Parameter File
==============
diff --git a/doc/RTD/source/InitialConditions/index.rst b/doc/RTD/source/InitialConditions/index.rst
index e9684ac4ffde5886f7a110de2cd7fb0fbb572a5e..eba438c722fbf4ffd78984aa55d6bfa5efcd71ad 100644
--- a/doc/RTD/source/InitialConditions/index.rst
+++ b/doc/RTD/source/InitialConditions/index.rst
@@ -4,45 +4,63 @@
Initial Conditions
==================
-To run anything more than examples from our suite, you will need to be able to
-produce your own initial conditions for SWIFT. We use the same initial conditions
-as the popular GADGET-2 code, which uses the HDF5 file format.
+To run anything more than examples from our suite, you will need to be able to
+produce your own initial conditions for SWIFT. We use the same initial
+conditions format as the popular `GADGET-2
+<https://wwwmpa.mpa-garching.mpg.de/~volker/gadget/>`_ code, which uses HDF5 for
+its type 3 format. Note that we do not support the GADGET-2 types 1 and 2
+formats.
+
+The original GADGET-2 file format only contains 2 types of particles: gas
+particles and 5 sorts of collisionless particles that allow users to run with 5
+separate particle masses and softenings. In SWIFT, we expand on this by using
+two of these types for stars and black holes.
+
+GADGET-2 can have initial conditions split over many files. This allow multiple
+ones to be read in parallel and is the only way the code can handle more than
+2^31 particles. This limitation is not in place in SWIFT. A single file can
+contain any number of particles (well... up to 2^64...) and the file is read in
+parallel by HDF5 when running on more than one compute node.
As the original documentation for the GADGET-2 initial conditions format is
-quite sparse, we lay out here all of the necessary components. If you are generating
-your initial conditions from python, we recommend you use the h5py package. We
-provide a writing wrapper for this for our initial conditions in
+quite sparse, we lay out here all of the necessary components. If you are
+generating your initial conditions from python, we recommend you use the h5py
+package. We provide a writing wrapper for this for our initial conditions in
``examples/KeplerianRing/write_gadget.py``.
-You can find out more about the HDF5 format on their webpages, here:
-https://support.hdfgroup.org/HDF5/doc/H5.intro.html
+You can find out more about the HDF5 format on their `webpages
+<https://support.hdfgroup.org/HDF5/doc/H5.intro.html>`_.
Structure of the File
---------------------
-There are several groups that contain 'auxilliary' information, such as ``Header``.
-Particle data is placed in groups that signify particle type.
-
-+---------------------+------------------------+
-| Group Name | Physical Particle Type |
-+=====================+========================+
-| ``PartType0`` | Gas |
-+---------------------+------------------------+
-| ``PartType1`` | Dark Matter |
-+---------------------+------------------------+
-| ``PartType2`` | Ignored |
-+---------------------+------------------------+
-| ``PartType3`` | Ignored |
-+---------------------+------------------------+
-| ``PartType4`` | Stars |
-+---------------------+------------------------+
-| ``PartType5`` | Black Holes |
-+---------------------+------------------------+
-
-Currently, not all of these particle types are included in SWIFT. Note that the
-only particles that have hydrodynamical forces calculated between them are those
-in ``PartType0``.
+There are several groups that contain 'auxilliary' information, such as
+``Header``. Particle data is placed in separate groups depending of the type of
+the particles. Some types are currently ignored by SWIFT but are kept in the
+file format for compatibility reasons.
+
++---------------------+------------------------+----------------------------+
+| HDF5 Group Name | Physical Particle Type | In code ``enum part_type`` |
++=====================+========================+============================+
+| ``/PartType0/`` | Gas | ``swift_type_gas`` |
++---------------------+------------------------+----------------------------+
+| ``/PartType1/`` | Dark Matter | ``swift_type_dark_matter`` |
++---------------------+------------------------+----------------------------+
+| ``/PartType2/`` | Ignored | |
++---------------------+------------------------+----------------------------+
+| ``/PartType3/`` | Ignored | |
++---------------------+------------------------+----------------------------+
+| ``/PartType4/`` | Stars | ``swift_type_star`` |
++---------------------+------------------------+----------------------------+
+| ``/PartType5/`` | Black Holes | ``swift_type_black_hole`` |
++---------------------+------------------------+----------------------------+
+
+The last column in the table gives the ``enum`` value from ``part_type.h``
+corresponding to a given entry in the files.
+
+Note that the only particles that have hydrodynamical forces calculated between
+them are those in ``PartType0``.
Necessary Components
@@ -55,27 +73,34 @@ script.
Header
~~~~~~
-In ``Header``, the following attributes are required:
-
-+ ``BoxSize``, a floating point number or N-dimensional (usually 3) array
- that describes the size of the box.
+In the ``/Header/`` group, the following attributes are required:
+
++ ``Dimension``, an integer indicating the dimensionality of the ICs (1,2 or 3).
+ Note that this parameter is an addition to the GADGET-2 format and will be
+ ignored by GADGET. SWIFT will use this value to verify that the dimensionality
+ of the code matches the ICs. If this parameter is not provided, it defaults
+ to 3.
++ ``BoxSize``, a floating point number or N-dimensional (usually 3) array that
+ describes the size of the box. If only one number is provided (as per the
+ GADGET-2 standard) then the box is assumed have the same size along all the
+ axis. In cosmological runs, this is the comoving box-size expressed in the
+ units specified in the ``/Units`` group (see below). Note that, unlike GADGET,
+ we express all quantities in "h-free" units. So that, for instance, we express
+ the box side-length in ``Mpc`` and not ``Mpc/h``.
++ ``NumPart_Total``, a length 6 array of integers that tells the code how many
+ particles of each type are in the initial conditions file. Unlike traditional
+ GADGET-2 files, these can be >2^31.
++ ``NumPart_Total_HighWord``, a historical length-6 array that tells the code
+ the number of high word particles in the initial conditions there are. If you
+ are unsure, just set this to ``[0, 0, 0, 0, 0, 0]``. This does have to be
+ present but can be a set of 0s unless you have more than 2^31 particles and
+ want to be fully compliant with GADGET-2. Note that, as SWIFT supports
+ ``NumPart_Total`` to be >2^31, the use of ``NumPart_Total_HighWord`` is only
+ here for compatibility reasons.
+ ``Flag_Entropy_ICs``, a historical value that tells the code if you have
included entropy or internal energy values in your intial conditions files.
- Acceptable values are 0 or 1.
-+ ``NumPart_Total``, a length 6 array of integers that tells the code how many
- particles are of each type are in the initial conditions file.
-+ ``NumPart_Total_HighWord``, a historical length-6 array that tells the code
- the number of high word particles in the initial conditions there are. If
- you are unsure, just set this to ``[0, 0, 0, 0, 0, 0]``. This does have to be
- present, but unlike GADGET-2, this can be a set of 0s unless you have more than
- 2^31 particles.
-+ ``NumFilesPerSnapshot``, again a historical integer value that tells the code
- how many files there are per snapshot. You will probably want to set this to 1
- and simply have a single HDF5 file for your initial conditions; SWIFT can
- leverage parallel-HDF5 to read from this single file in parallel.
-+ ``NumPart_ThisFile``, a length 6 array of integers describing the number of
- particles in this file. If you have followed the above advice, this will be
- exactly the same as the ``NumPart_Total`` array.
+ Acceptable values are 0 or 1. We recommend using internal energies over
+ entropy in the ICs and hence have this flag set to 0.
You may want to include the following for backwards-compatibility with many
GADGET-2 based analysis programs:
@@ -83,63 +108,102 @@ GADGET-2 based analysis programs:
+ ``MassTable``, an array of length 6 which gives the masses of each particle
type. SWIFT ignores this and uses the individual particle masses, but some
programs will crash if it is not included.
-+ ``Time``, the internal code time of the start (set this to 0).
-
++ ``NumPart_ThisFile``, a length 6 array of integers describing the number of
+ particles in this file. If you have followed the above advice, this will be
+ exactly the same as the ``NumPart_Total`` array. As SWIFT only uses ICs
+ contained in a single file, this is not necessary for SWIFT-only ICs.
++ ``NumFilesPerSnapshot``, again a historical integer value that tells the code
+ how many files there are per snapshot. You will probably want to set this to 1.
++ ``Time``, time of the start of the simulation in internal units or expressed
+ as a scale-factor for cosmological runs. SWIFT ignores this and reads it from
+ the parameter file.
+
RuntimePars
~~~~~~~~~~~
-In ``RuntimePars``, the following attributes are required:
+In the ``/RuntimePars/``, the following attributes are required:
+ ``PeriodicBoundaryConditionsOn``, a flag to tell the code whether or not you
have periodic boundaries switched on. Again, this is historical; it should be
set to 1 (default) if you have the code running in periodic mode, or 0 otherwise.
-Units
-~~~~~
-
-In ``Units``, you will need to specify what units your initial conditions are
-in. If these are not present, the code assumes that you are using the same
-units for your initial conditions as are in your parameterfile, but it is best
-to include them to be on the safe side. You will need:
-
-+ ``Unit current in cgs (U_I)``
-+ ``Unit length in cgs (U_L)``
-+ ``Unit mass in cgs (U_M)``
-+ ``Unit temperature in cgs (U_T)``
-+ ``Unit time in cgs (U_t)``
-
-These are all floating point numbers.
-
-
Particle Data
~~~~~~~~~~~~~
Now for the interesting part! You can include particle data groups for each
-individual particle type (e.g. ``PartType0``) that have the following _datasets_:
+individual particle type (e.g. ``/PartType0/``) that have the following *datasets*:
+ ``Coordinates``, an array of shape (N, 3) where N is the number of particles
- of that type, that are the cartesian co-ordinates of the particles. Co-ordinates
- must be positive, but will be wrapped on reading to be within the periodic box.
-+ ``Velocities``, an array of shape (N, 3) that is the cartesian velocities
- of the particles.
+ of that type, that are the cartesian co-ordinates of the
+ particles. Co-ordinates must be within the box so, in the case of a cube
+ within [0, L)^3 where L is the side-length of the simulation volume. In the
+ case of cosmological simulations, these are the co-moving positions.
++ ``Velocities``, an array of shape (N, 3) that is the cartesian velocities of
+ the particles. When running cosmological simulations, these are the peculiar
+ velocities. Note that this is different from GADGET which uses peculiar
+ velocities divided by ``sqrt(a)`` (see below for a fix).
+ ``ParticleIDs``, an array of length N that are unique identifying numbers for
each particle. Note that these have to be unique to a particle, and cannot be
- the same even between particle types. Please ensure that your IDs are positive
- integer numbers.
+ the same even between particle types. The **IDs must be >1**. 0 or negative
+ IDs will be rejected by the code.
+ ``Masses``, an array of length N that gives the masses of the particles.
-For ``PartType0`` (i.e. particles that interact through hydrodynamics), you will
+For ``PartType0`` (i.e. particles that interact through hydro-dynamics), you will
need the following auxilliary items:
-+ ``InternalEnergy``, an array of length N that gives the internal energies of
- the particles. For PressureEntropy, you can specify ``Entropy`` instead.
-+ ``SmoothingLength``, the smoothing lenghts of the particles. These will be
- tidied up a bit, but it is best if you provide accurate numbers.
++ ``SmoothingLength``, the smoothing lengths of the particles. These will be
+ tidied up a bit, but it is best if you provide accurate numbers. In
+ cosmological runs, these are the co-moving smoothing lengths.
++ ``InternalEnergy``, an array of length N that gives the internal energies per
+ unit mass of the particles. If the hydro-scheme used in the code is based on
+ another thermodynamical quantity (entropy or total energy, etc.), the
+ conversion will happen inside the code. In cosmological runs, this is the
+ **physical** internal energy per unit mass. This has the dimension of velocity
+ squared.
+
+
+Note that for cosmological runs, all quantities have to be expressed in "h-free"
+dimensions. This means ``Mpc`` and not ``Mpc/h`` for instance. If the ICs have
+been generated for GADGET (where h-full values are expected), the parameter
+``InitialConditions:cleanup_h_factors`` can be set to ``1`` in the
+:ref:`Parameter_File_label` to make SWIFT convert the quantities read in to
+h-free quantities. Switching this parameter on will also affect the box size
+read from the ``/Header/`` group (see above).
+
+Similarly, GADGET cosmological ICs have traditionally used velocities expressed
+as peculiar velocities divided by ``sqrt(a)``. This can be undone by swicthing
+on the parameter ``InitialConditions:cleanup_velocity_factors`` in the
+:ref:`Parameter_File_label`.
+
+
+Optional Components
+-------------------
+
+In the ``/Units/`` HDF5 group, you cans specify what units your initial conditions are
+in. If this group is not present, the code assumes that you are using the same
+units for your initial conditions as in your :ref:`Parameter_File_label`
+(i.e. as the internal units system used by the code), but it is best to include
+them to be on the safe side. You will need:
+
++ ``Unit length in cgs (U_L)``
++ ``Unit mass in cgs (U_M)``
++ ``Unit time in cgs (U_t)``
++ ``Unit current in cgs (U_I)``
++ ``Unit temperature in cgs (U_T)``
+
+These are all floating point numbers. Note that we specify the time units and
+not the velocity units.
+
+If the units specified in the initial conditions are different from the internal
+units (specified in the parameter file), SWIFT will perform a conversion of all
+the quantities when reading in the ICs. This includes a conversion of the box
+size read from the ``/Header/`` group.
+
Summary
-~~~~~~~
+-------
You should have an HDF5 file with the following structure:
@@ -147,11 +211,9 @@ You should have an HDF5 file with the following structure:
Header/
BoxSize=[x, y, z]
- Flag_Entropy_ICs=1
- NumPart_Total=[0, 1, 2, 3, 4, 5]
+ Flag_Entropy_ICs=0
+ NumPart_Total=[0, 1, 0, 0, 4, 5]
NumPart_Total_HighWord=[0, 0, 0, 0, 0, 0]
- NumFilesPerSnapshot=1
- NumPart_ThisFile=[0, 1, 2, 3, 4, 5]
RuntimePars/
PeriodicBoundariesOn=1
Units/
diff --git a/examples/AgoraDisk/agora_disk.yml b/examples/AgoraDisk/agora_disk.yml
new file mode 100644
index 0000000000000000000000000000000000000000..598f1db2858c3e420c55ae31c9ca10d37c2f48b7
--- /dev/null
+++ b/examples/AgoraDisk/agora_disk.yml
@@ -0,0 +1,69 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.085e21 # kpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+Scheduler:
+ max_top_level_cells: 8
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.5 # The end time of the simulation (in internal units).
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: agora_disk # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
+ compression: 4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.08 # Comoving softening length (in internal units).
+ max_physical_softening: 0.08 # Physical softening length (in internal units).
+ mesh_side_length: 32 # Number of cells along each axis for the periodic gravity mesh.
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 10 # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./agora_disk.hdf5 # The file to read
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
+ shift: [674.1175, 674.1175, 674.1175] # (Optional) A shift to apply to all particles read from the ICs (in internal units).
+
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+ lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
+ minimum_temperature: 1.0e2 # Minimal temperature (Kelvin)
+ mean_molecular_weight: 0.59 # Mean molecular weight
+ hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
+ cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
+
+# Cooling with Grackle 2.0
+GrackleCooling:
+ CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ WithUVbackground: 1 # Enable or not the UV background
+ Redshift: 0 # Redshift to use (-1 means time based redshift)
+ WithMetalCooling: 1 # Enable or not the metal cooling
+ ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
+ ProvideSpecificHeatingRates: 0 # User provide specific heating rates
+ SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
+ OutputMode: 1 # Write in output corresponding primordial chemistry mode
+ MaxSteps: 1000
+ ConvergenceLimit: 1e-2
diff --git a/examples/AgoraDisk/changeType.py b/examples/AgoraDisk/changeType.py
new file mode 100644
index 0000000000000000000000000000000000000000..9a0eba9def3a07b01f18f727de220d8f7193858f
--- /dev/null
+++ b/examples/AgoraDisk/changeType.py
@@ -0,0 +1,92 @@
+#!/usr/bin/env python3
+
+from h5py import File
+from sys import argv
+import numpy as np
+
+"""
+Change particle types in order to match the implemented types
+"""
+
+# number of particle type
+N_type = 6
+
+debug = 0
+
+
+def getOption():
+ if len(argv) != 2:
+ raise IOError("You need to provide a filename")
+
+ # get filename and read it
+ filename = argv[-1]
+
+ return filename
+
+
+def groupName(part_type):
+ return "PartType%i" % part_type
+
+
+def changeType(f, old, new):
+ # check if directory exists
+ old_group = groupName(old)
+ if old_group not in f:
+ raise IOError("Cannot find group '%s'" % old)
+ old = f[old_group]
+
+ new_group = groupName(new)
+ if new_group not in f:
+ f.create_group(new_group)
+ new = f[new_group]
+
+ for name in old:
+ if debug:
+ print("Moving '%s' from '%s' to '%s'"
+ % (name, old_group, new_group))
+
+ tmp = old[name][:]
+ del old[name]
+ if name in new:
+ new_tmp = new[name][:]
+ if debug:
+ print("Found previous data:", tmp.shape, new_tmp.shape)
+ tmp = np.append(tmp, new_tmp, axis=0)
+ del new[name]
+
+ if debug:
+ print("With new shape:", tmp.shape)
+ new.create_dataset(name, tmp.shape)
+ new[name][:] = tmp
+
+ del f[old_group]
+
+
+def countPart(f):
+ npart = []
+ for i in range(N_type):
+ name = groupName(i)
+ if name in f:
+ grp = f[groupName(i)]
+ N = grp["Masses"].shape[0]
+ else:
+ N = 0
+ npart.append(N)
+
+ f["Header"].attrs["NumPart_ThisFile"] = npart
+ f["Header"].attrs["NumPart_Total"] = npart
+ f["Header"].attrs["NumPart_Total_HighWord"] = [0]*N_type
+
+
+if __name__ == "__main__":
+ filename = getOption()
+
+ f = File(filename)
+
+ changeType(f, 2, 1)
+ changeType(f, 3, 1)
+ changeType(f, 4, 1)
+
+ countPart(f)
+
+ f.close()
diff --git a/examples/AgoraDisk/cleanupSwift.py b/examples/AgoraDisk/cleanupSwift.py
new file mode 100644
index 0000000000000000000000000000000000000000..e3c364081fc82986968727695e1142690781bb67
--- /dev/null
+++ b/examples/AgoraDisk/cleanupSwift.py
@@ -0,0 +1,28 @@
+#!/usr/bin/env python3
+
+# ./translate_particles.py filename output_name
+from h5py import File
+import sys
+from shutil import copyfile
+
+NPartType = 1
+filename = sys.argv[-2]
+out = sys.argv[-1]
+
+copyfile(filename, out)
+
+f = File(out)
+
+name = "PartType0/ElementAbundance"
+if name in f:
+ del f[name]
+
+for i in range(NPartType):
+ name = "PartType%i" % i
+ if name not in f:
+ continue
+
+ grp = f[name + "/SmoothingLength"]
+ grp[:] *= 1.823
+
+f.close()
diff --git a/examples/AgoraDisk/getIC.sh b/examples/AgoraDisk/getIC.sh
new file mode 100644
index 0000000000000000000000000000000000000000..620a751bedaf6c646119247270fad6dd3f740fde
--- /dev/null
+++ b/examples/AgoraDisk/getIC.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+if [ "$#" -ne 1 ]; then
+ echo "You need to provide the resolution (e.g. ./getIC.sh low)."
+ echo "The possible options are low, med and high."
+ exit
+fi
+
+wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/$1
diff --git a/examples/AgoraDisk/getSolution.sh b/examples/AgoraDisk/getSolution.sh
new file mode 100644
index 0000000000000000000000000000000000000000..41fbab800098bfaa381ca2e83cab0210c8dcf924
--- /dev/null
+++ b/examples/AgoraDisk/getSolution.sh
@@ -0,0 +1,36 @@
+#!/bin/bash
+
+OPTIND=1
+
+with_cooling=0
+
+function show_help {
+ echo "Valid options are:"
+ echo "\t -h \t Show this help"
+ echo "\t -C \t Download solution with cooling"
+}
+
+while getopts "h?C" opt; do
+ case "$opt" in
+ h|\?)
+ show_help
+ exit
+ ;;
+ C)
+ with_cooling=1
+ ;;
+ esac
+done
+
+# cleanup work space
+rm snapshot_0000.hdf5 snapshot_0500.hdf5
+
+if [ $with_cooling -eq 1 ]; then
+ wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0000.hdf5
+ wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0500.hdf5
+else
+ wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0000.hdf5
+ wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0500.hdf5
+fi
+
+
diff --git a/examples/AgoraDisk/plotSolution.py b/examples/AgoraDisk/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c3de1728eb6f22bee72361db1e7f0c77613eb25
--- /dev/null
+++ b/examples/AgoraDisk/plotSolution.py
@@ -0,0 +1,2634 @@
+#!/usr/bin/env python2
+#######################################################################
+#
+# UNIFIED ANALYSIS SCRIPT FOR DISK SIMULATION FOR THE AGORA PROJECT
+#
+# FOR SCRIPT HISTORY SEE VERSION CONTROL CHANGELOG
+#
+# Note: This script requires yt-3.2 or yt-dev. Older versions may
+# yield incorrect results.
+#
+#######################################################################
+# This script is a copy of the AGORA project (https://bitbucket.org/mornkr/agora-analysis-script/)
+# modified in order to take into account SWIFT
+import matplotlib
+matplotlib.use('Agg')
+import sys
+import math
+import copy
+import csv
+import numpy as np
+import matplotlib.colorbar as cb
+import matplotlib.lines as ln
+import yt.utilities.physical_constants as constants
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from yt.mods import *
+from yt.units.yt_array import YTQuantity
+from mpl_toolkits.axes_grid1 import AxesGrid
+from matplotlib.offsetbox import AnchoredText
+from matplotlib.ticker import MaxNLocator
+from yt.fields.particle_fields import add_volume_weighted_smoothed_field
+from yt.fields.particle_fields import add_nearest_neighbor_field
+from yt.analysis_modules.star_analysis.api import StarFormationRate
+from yt.analysis_modules.halo_analysis.api import *
+from yt.data_objects.particle_filters import add_particle_filter
+from scipy.stats import kde
+from subprocess import call
+#mylog.setLevel(1)
+
+import sys
+
+if "-C" in sys.argv:
+ with_cooling = 1
+else:
+ with_cooling = 0
+
+draw_density_map = 1#1 # 0/1 = OFF/ON
+draw_temperature_map = 1#1 # 0/1 = OFF/ON
+draw_cellsize_map = 2#2 # 0/1/2/3 = OFF/ON/ON now with resolution map where particle_size is defined as [M/rho]^(1/3) for SPH/ON with both cell_size and resolution maps
+draw_elevation_map = 1#1 # 0/1 = OFF/ON
+draw_metal_map = 0#0 # 0/1/2 = OFF/ON/ON with total metal mass in the simulation annotated (this will be inaccurate for SPH)
+draw_zvel_map = 0#0 # 0/1 = OFF/ON
+draw_star_map = 0#0 # 0/1 = OFF/ON
+draw_star_clump_stats = 0#0 # 0/1/2/3 = OFF/ON/ON with additional star_map with annotated clumps/ON with addtional star_map + extra clump mass function with GIZMO-ps2
+draw_SFR_map = 0###0 # 0/1/2 = OFF/ON/ON with local K-S plot using patches
+draw_PDF = 0###0 # 0/1/2/3 = OFF/ON/ON with constant pressure/entropy lines/ON with additional annotations such as 1D profile from a specific code (CHANGA)
+draw_pos_vel_PDF = 0##0 # 0/1/2/3/4 = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
+draw_star_pos_vel_PDF = 0##0 # 0/1/2/3/4 = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
+draw_rad_height_PDF = 0##0 # 0/1/2/3 = OFF/ON/ON with 1D profile/ON with analytic ftn subtracted
+draw_metal_PDF = 0##0 # 0/1 = OFF/ON
+draw_density_DF = 0#0 # 0/1/2 = OFF/ON/ON with difference plot between 1st and 2nd datasets (when 2, dataset_num should be set to 2)
+draw_radius_DF = 0#0 # 0/1 = OFF/ON
+draw_star_radius_DF = 0#0 # 0/1/2 = OFF/ON/ON with SFR profile and K-S plot (when 2, this automatically turns on draw_radius_DF)
+draw_height_DF = 0#0 # 0/1 = OFF/ON
+draw_SFR = 0#0 # 0/1/2 = OFF/ON/ON with extra line with GIZMO-ps2
+draw_cut_through = 0#0 # 0/1 = OFF/ON
+add_nametag = 1#1 # 0/1 = OFF/ON
+add_mean_fractional_dispersion = 0#0 # 0/1 = OFF/ON
+
+dataset_num = 2 # 1/2 = 1st dataset(Grackle+noSF)/2nd dataset(Grackle+SF+ThermalFbck) for AGORA Paper 4
+yt_version_pre_3_2_3 = 0 # 0/1 = NO/YES to "Is the yt version being used pre yt-dev-3.2.3?"
+times = [0, 500] # in Myr
+figure_width = 30 # in kpc
+n_ref = 64 # 8; for SPH codes
+over_refine_factor = 1 # 2; for SPH codes
+aperture_size_SFR_map = 750 # in pc = Used if draw_SFR_map = 1, 750 matches Bigiel et al. (2008)
+young_star_cutoff_SFR_map = 20 # in Myr = Used if draw_SFR_map = 1
+young_star_cutoff_star_radius_DF = 20 # in Myr = Used if draw_star_radius_DF = 1
+mean_dispersion_radius_range = [2, 10] # in kpc = Used if add_mean_fractional_dispersion = 1, for draw_pos_vel_PDF, draw_star_pos_vel_PDF, draw_rad_height_PDF, draw_radius_DF, draw_star_radius_DF
+mean_dispersion_height_range = [0, 0.6] # in kpc = Used if add_mean_fractional_dispersion = 1, for draw_height_DF
+mean_dispersion_time_range = [50, 500] # in Myr = Used if add_mean_fractional_dispersion = 1, for draw_SFR
+mean_dispersion_density_range = [1e-25, 1e-22]# in g/cm3 = Used if add_mean_fractional_dispersion = 1, for draw_density_DF
+disk_normal_vector = [0., 0., 1.]
+
+gadget_default_unit_base = {'UnitLength_in_cm' : 3.08568e+21,
+ 'UnitMass_in_g' : 1.989e+43,
+ 'UnitVelocity_in_cm_per_s' : 100000}
+color_names = ['red', 'magenta', 'orange', 'gold', 'green', 'cyan', 'blue', 'blueviolet', 'black']
+linestyle_names = ['-']
+marker_names = ['s', 'o', 'p', 'v', '^', '<', '>', 'h', '*']
+
+# file_location = ['/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository-for-use/Grackle+noSF/',
+# '/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository-for-use/Grackle+SF+ThermalFbck/']
+# file_location = ['/global/project/projectdirs/agora/AGORA-DISK-repository-for-use/Grackle+noSF/',
+# '/global/project/projectdirs/agora/AGORA-DISK-repository-for-use/Grackle+SF+ThermalFbck/']
+# codes = ['ART-I', 'ART-II', 'ENZO', 'RAMSES', 'CHANGA', 'GASOLINE', 'GADGET-3', 'GEAR', 'GIZMO']
+# filenames = [[[file_location[0]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[0]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d'],
+# [file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000000.art', file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000098.art'],
+# [file_location[0]+'ENZO/DD0000/DD0000', file_location[0]+'ENZO/DD0100/DD0100'],
+# [file_location[0]+'RAMSES/output_00001/info_00001.txt', file_location[0]+'RAMSES/output_00068/info_00068.txt'],
+# [file_location[0]+'CHANGA/disklow/disklow.000000', file_location[0]+'CHANGA/disklow/disklow.000500'],
+# [file_location[0]+'GASOLINE/LOW_dataset1.00001', file_location[0]+'GASOLINE/LOW_dataset1.00335'],
+# [file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_000.hdf5', file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_010.hdf5'],
+# [file_location[0]+'GEAR/snapshot_0000', file_location[0]+'GEAR/snapshot_0500'],
+# [file_location[0]+'GIZMO/snapshot_temp_000', file_location[0]+'GIZMO/snapshot_temp_100']],
+# [[file_location[1]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[1]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d'],
+# [file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000000.art', file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000311.art'],
+# [file_location[1]+'ENZO/DD0000/DD0000', file_location[1]+'ENZO/DD0050/DD0050'],
+# [file_location[1]+'RAMSES/output_00001/info_00001.txt', file_location[1]+'RAMSES/output_00216/info_00216.txt'],
+# [file_location[1]+'CHANGA/disklow/disklow.000000', file_location[1]+'CHANGA/disklow/disklow.000500'],
+# [file_location[1]+'GASOLINE/LOW_dataset1.00001', file_location[1]+'GASOLINE/LOW_dataset2.00335'],
+# [file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_000.hdf5', file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_010.hdf5'],
+# [file_location[1]+'GEAR/snapshot_0000', file_location[1]+'GEAR/snapshot_0500'],
+# [file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
+codes = ['SWIFT', 'GEAR']
+filenames = [[["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"],
+ ["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]],
+ [["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"],
+ ["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]]]
+
+# codes = ["SWIFT"]
+# filenames = [[["./agora_disk_0000.hdf5", "./agora_disk_0050.hdf5"]],
+# [["./agora_disk_0000.hdf5", "./agora_disk_0050.hdf5"]]]
+# codes = ['ART-I']
+# filenames = [[[file_location[0]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[0]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d']],
+# [[file_location[1]+'ART-I/IC/AGORA_Galaxy_LOW.d', file_location[1]+'ART-I/t0.5Gyr/10MpcBox_csf512_02350.d']]]
+# codes = ['ART-II']
+# filenames = [[[file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000000.art', file_location[0]+'ART-II/noSF_def_2p/OUT/AGORA_LOW_000098.art']],
+# [[file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000000.art', file_location[1]+'ART-II/SF_FBth_def_2p/OUT/AGORA_LOW_000311.art']]]
+# codes = ['ENZO']
+# filenames = [[[file_location[0]+'ENZO/DD0000/DD0000', file_location[0]+'ENZO/DD0100/DD0100']],
+# [[file_location[1]+'ENZO/DD0000/DD0000', file_location[1]+'ENZO/DD0050/DD0050']]]
+# codes = ['RAMSES']
+# filenames = [[[file_location[0]+'RAMSES/output_00001/info_00001.txt', file_location[0]+'RAMSES/output_00068/info_00068.txt']],
+# [[file_location[1]+'RAMSES/output_00001/info_00001.txt', file_location[1]+'RAMSES/output_00216/info_00216.txt']]]
+# codes = ['CHANGA']
+# filenames = [[[file_location[0]+'CHANGA/disklow/disklow.000000', file_location[0]+'CHANGA/disklow/disklow.000500']],
+# [[file_location[1]+'CHANGA/disklow/disklow.000000', file_location[1]+'CHANGA/disklow/disklow.000500']]]
+# codes = ['GASOLINE']
+# filenames = [[[file_location[0]+'GASOLINE/LOW_dataset1.00001', file_location[0]+'GASOLINE/LOW_dataset1.00335']],
+# [[file_location[1]+'GASOLINE/LOW_dataset1.00001', file_location[1]+'GASOLINE/LOW_dataset2.00335']]]
+# codes = ['GADGET-3']
+# filenames = [[[file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_000.hdf5', file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_010.hdf5']],
+# [[file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_000.hdf5', file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_010.hdf5']]]
+# codes = ['GEAR']
+# filenames = [[['snapshot_0000', 'snapshot_0500']],
+# [['snapshot_0000', 'snapshot_0500']]]
+# codes = ['GIZMO']
+# filenames = [[[file_location[0]+'GIZMO/snapshot_temp_000', file_location[0]+'GIZMO/snapshot_temp_100']],
+# [[file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
+
+def load_dataset(dataset_num, time, code, codes, filenames_entry):
+ if codes[code] == 'ART-I': # ART frontend doesn't find the accompanying files, so we specify them; see http://yt-project.org/docs/dev/examining/loading_data.html#art-data
+ dirnames = filenames_entry[code][time][:filenames_entry[code][time].rfind('/')+1]
+ if time == 0:
+ timestamp = ''
+ else:
+ timestamp = filenames_entry[code][time][filenames_entry[code][time].rfind('_'):filenames_entry[code][time].rfind('.')]
+ pf = load(filenames_entry[code][time], file_particle_header=dirnames+'PMcrd'+timestamp+'.DAT', file_particle_data=dirnames+'PMcrs0'+timestamp+'.DAT', file_particle_stars=dirnames+'stars'+timestamp+'.dat')
+ elif codes[code] == 'CHANGA' or codes[code] == 'GASOLINE': # For TIPSY frontend, always make sure to place your parameter file in the same directory as your datasets
+ pf = load(filenames_entry[code][time], n_ref=n_ref, over_refine_factor=over_refine_factor)
+ elif codes[code] == 'GADGET-3': # For GADGET-3 2nd dataset, there somehow exist particles very far from the center; so we use [-2000, 2000] for a bounding_box
+ pf = load(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-2000.0, 2000.0], [-2000.0, 2000.0], [-2000.0, 2000.0]], n_ref=n_ref, over_refine_factor=over_refine_factor)
+ elif codes[code] == "SWIFT":
+ pf = load(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-2000.0, 2000.0], [-2000.0, 2000.0], [-2000.0, 2000.0]], n_ref=n_ref, over_refine_factor=over_refine_factor)
+ elif codes[code] == 'GEAR':
+ pf = load(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-2000.0, 2000.0], [-2000.0, 2000.0], [-2000.0, 2000.0]], n_ref=n_ref, over_refine_factor=over_refine_factor)
+ # from yt.frontends.gadget.definitions import gadget_header_specs
+ # from yt.frontends.gadget.definitions import gadget_ptype_specs
+ # from yt.frontends.gadget.definitions import gadget_field_specs
+ # if with_cooling:
+ # gadget_header_specs["chemistry"] = (('h1', 4, 'c'),('h2', 4, 'c'),('empty', 256, 'c'),)
+ # header_spec = "default+chemistry"
+ # else:
+ # header_spec = "default"
+ # gear_ptype_specs = ("Gas", "Stars", "Halo", "Disk", "Bulge", "Bndry")
+ # if dataset_num == 1:
+ # pf = GadgetDataset(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-1000.0, 1000.0], [-1000.0, 1000.0], [-1000.0, 1000.0]], header_spec=header_spec,
+ # ptype_spec=gear_ptype_specs, n_ref=n_ref, over_refine_factor=over_refine_factor)
+ # elif dataset_num == 2:
+ # # For GEAR 2nd dataset (Grackle+SF+ThermalFbck), "Metals" is acutally 10-species field; check how metallicity field is added below.
+ # if with_cooling:
+ # agora_gear = ( "Coordinates",
+ # "Velocities",
+ # "ParticleIDs",
+ # "Mass",
+ # ("InternalEnergy", "Gas"),
+ # ("Density", "Gas"),
+ # ("SmoothingLength", "Gas"),
+ # ("Metals", "Gas"),
+ # ("StellarFormationTime", "Stars"),
+ # ("StellarInitMass", "Stars"),
+ # ("StellarIDs", "Stars"),
+ # ("StellarDensity", "Stars"),
+ # ("StellarSmoothingLength", "Stars"),
+ # ("StellarMetals", "Stars"),
+ # ("Opt1", "Stars"),
+ # ("Opt2", "Stars"),
+ # )
+ # else:
+ # agora_gear = ( "Coordinates",
+ # "Velocities",
+ # "ParticleIDs",
+ # "Mass",
+ # ("InternalEnergy", "Gas"),
+ # ("Density", "Gas"),
+ # ("SmoothingLength", "Gas"),
+ # )
+ # gadget_field_specs["agora_gear"] = agora_gear
+ # pf = GadgetDataset(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-1000.0, 1000.0], [-1000.0, 1000.0], [-1000.0, 1000.0]], header_spec=header_spec,
+ # ptype_spec=gear_ptype_specs, field_spec="agora_gear", n_ref=n_ref, over_refine_factor=over_refine_factor)
+ elif codes[code] == 'GIZMO':
+ from yt.frontends.gadget.definitions import gadget_field_specs
+ if dataset_num == 1:
+ agora_gizmo = ( "Coordinates",
+ "Velocities",
+ "ParticleIDs",
+ "Mass",
+ ("Temperature", "Gas"),
+ ("Density", "Gas"),
+ ("Electron_Number_Density", "Gas"),
+ ("HI_NumberDensity", "Gas"),
+ ("SmoothingLength", "Gas"),
+ )
+ elif dataset_num == 2:
+ agora_gizmo = ( "Coordinates",
+ "Velocities",
+ "ParticleIDs",
+ "Mass",
+ ("Temperature", "Gas"),
+ ("Density", "Gas"),
+ ("ElectronAbundance", "Gas"),
+ ("NeutralHydrogenAbundance", "Gas"),
+ ("SmoothingLength", "Gas"),
+ ("StarFormationRate", "Gas"),
+ ("StellarFormationTime", "Stars"),
+ ("Metallicity", ("Gas", "Stars")),
+ ("DelayTime", "Stars"),
+ ("StellarInitMass", "Stars"),
+ )
+ gadget_field_specs["agora_gizmo"] = agora_gizmo
+ pf = load(filenames_entry[code][time], unit_base = gadget_default_unit_base, bounding_box=[[-1000.0, 1000.0], [-1000.0, 1000.0], [-1000.0,1000.0]], field_spec="agora_gizmo",
+ n_ref=n_ref, over_refine_factor=over_refine_factor)
+ else:
+ pf = load(filenames_entry[code][time])
+ return pf
+
+fig_density_map = []
+fig_temperature_map = []
+fig_cellsize_map = []
+fig_cellsize_map_2 = []
+fig_elevation_map = []
+fig_metal_map = []
+fig_zvel_map = []
+fig_star_map = []
+fig_star_map_2 = []
+fig_star_map_3 = []
+fig_degr_sfr_map = []
+fig_degr_density_map = []
+fig_PDF = []
+fig_pos_vel_PDF = []
+fig_star_pos_vel_PDF = []
+fig_rad_height_PDF = []
+fig_metal_PDF = []
+grid_density_map = []
+grid_temperature_map = []
+grid_cellsize_map = []
+grid_cellsize_map_2 = []
+grid_elevation_map = []
+grid_metal_map = []
+grid_zvel_map = []
+grid_star_map = []
+grid_star_map_2 = []
+grid_star_map_3 = []
+grid_degr_sfr_map = []
+grid_degr_density_map = []
+grid_PDF = []
+grid_pos_vel_PDF = []
+grid_star_pos_vel_PDF = []
+grid_rad_height_PDF = []
+grid_metal_PDF = []
+star_clump_masses_hop = []
+star_clump_masses_fof = []
+star_clump_masses_hop_ref = []
+star_clump_masses_fof_ref = []
+pos_vel_xs = []
+pos_vel_profiles = []
+pos_disp_xs = []
+pos_disp_profiles = []
+pos_disp_vert_xs = []
+pos_disp_vert_profiles = []
+star_pos_vel_xs = []
+star_pos_vel_profiles = []
+star_pos_disp_xs = []
+star_pos_disp_profiles = []
+star_pos_disp_vert_xs = []
+star_pos_disp_vert_profiles = []
+rad_height_xs = []
+rad_height_profiles = []
+density_DF_xs = []
+density_DF_profiles = []
+density_DF_1st_xs = []
+density_DF_1st_profiles = []
+radius_DF_xs = []
+radius_DF_profiles = []
+surface_density = []
+star_radius_DF_xs = []
+star_radius_DF_profiles = []
+star_surface_density = []
+sfr_radius_DF_xs = []
+sfr_radius_DF_profiles = []
+sfr_surface_density = []
+height_DF_xs = []
+height_DF_profiles = []
+height_surface_density = []
+sfr_ts = []
+sfr_cum_masses = []
+sfr_sfrs = []
+surf_dens_SFR_map = []
+sfr_surf_dens_SFR_map = []
+cut_through_zs = []
+cut_through_zvalues = []
+cut_through_xs = []
+cut_through_xvalues = []
+
+for time in range(len(times)):
+ if draw_density_map == 1:
+ fig_density_map += [plt.figure(figsize=(100,20))]
+ grid_density_map += [AxesGrid(fig_density_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ if draw_temperature_map == 1:
+ fig_temperature_map += [plt.figure(figsize=(100,20))]
+ grid_temperature_map += [AxesGrid(fig_temperature_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ if draw_cellsize_map >= 1:
+ fig_cellsize_map += [plt.figure(figsize=(100,20))]
+ grid_cellsize_map += [AxesGrid(fig_cellsize_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ fig_cellsize_map_2 += [plt.figure(figsize=(100,20))]
+ grid_cellsize_map_2 += [AxesGrid(fig_cellsize_map_2[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ if draw_elevation_map == 1:
+ fig_elevation_map += [plt.figure(figsize=(100,20))]
+ grid_elevation_map += [AxesGrid(fig_elevation_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (1, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "6%", cbar_pad = 0.02)]
+ if draw_metal_map >= 1:
+ fig_metal_map += [plt.figure(figsize=(100,20))]
+ grid_metal_map += [AxesGrid(fig_metal_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ if draw_zvel_map == 1:
+ fig_zvel_map += [plt.figure(figsize=(100,20))]
+ grid_zvel_map += [AxesGrid(fig_zvel_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ if draw_star_map == 1:
+ fig_star_map += [plt.figure(figsize=(100,20))]
+ grid_star_map += [AxesGrid(fig_star_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ if draw_star_clump_stats >= 1:
+ star_clump_masses_hop.append([])
+ star_clump_masses_fof.append([])
+ fig_star_map_2 += [plt.figure(figsize=(100,20))]
+ grid_star_map_2 += [AxesGrid(fig_star_map_2[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ fig_star_map_3 += [plt.figure(figsize=(100,20))]
+ grid_star_map_3 += [AxesGrid(fig_star_map_3[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.02)]
+ if draw_SFR_map >= 1:
+ fig_degr_density_map += [plt.figure(figsize=(100,20))]
+ grid_degr_density_map+= [AxesGrid(fig_degr_density_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (1, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "6%", cbar_pad = 0.02)]
+ fig_degr_sfr_map += [plt.figure(figsize=(100,20))]
+ grid_degr_sfr_map += [AxesGrid(fig_degr_sfr_map[time], (0.01,0.01,0.99,0.99), nrows_ncols = (1, len(codes)), axes_pad = 0.02, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "6%", cbar_pad = 0.02)]
+ surf_dens_SFR_map.append([])
+ sfr_surf_dens_SFR_map.append([])
+ if draw_SFR_map >= 2 or draw_star_radius_DF >= 2:
+ def import_text(filename, separator):
+ for line in csv.reader(open(filename), delimiter=separator, skipinitialspace=True):
+ if line:
+ yield line
+ if draw_PDF >= 1:
+ fig_PDF += [plt.figure(figsize=(50, 80))]
+ grid_PDF += [AxesGrid(fig_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
+ if draw_pos_vel_PDF >= 1:
+ fig_pos_vel_PDF += [plt.figure(figsize=(50, 80))]
+ grid_pos_vel_PDF += [AxesGrid(fig_pos_vel_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
+ pos_vel_xs.append([])
+ pos_vel_profiles.append([])
+ pos_disp_xs.append([])
+ pos_disp_profiles.append([])
+ pos_disp_vert_xs.append([])
+ pos_disp_vert_profiles.append([])
+ if draw_star_pos_vel_PDF >= 1:
+ fig_star_pos_vel_PDF += [plt.figure(figsize=(50, 80))]
+ grid_star_pos_vel_PDF+= [AxesGrid(fig_star_pos_vel_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
+ star_pos_vel_xs.append([])
+ star_pos_vel_profiles.append([])
+ star_pos_disp_xs.append([])
+ star_pos_disp_profiles.append([])
+ star_pos_disp_vert_xs.append([])
+ star_pos_disp_vert_profiles.append([])
+ if draw_rad_height_PDF >= 1:
+ fig_rad_height_PDF += [plt.figure(figsize=(50, 80))]
+ grid_rad_height_PDF += [AxesGrid(fig_rad_height_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
+ rad_height_xs.append([])
+ rad_height_profiles.append([])
+ if draw_metal_PDF == 1:
+ fig_metal_PDF += [plt.figure(figsize=(50, 80))]
+ grid_metal_PDF += [AxesGrid(fig_metal_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
+ label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
+ if draw_density_DF >= 1:
+ density_DF_xs.append([])
+ density_DF_profiles.append([])
+ density_DF_1st_xs.append([])
+ density_DF_1st_profiles.append([])
+ if draw_density_DF == 2 and dataset_num == 1:
+ print("This won't work; resetting draw_density_DF to 1...")
+ draw_density_DF == 1
+ if draw_star_radius_DF >= 1:
+ star_radius_DF_xs.append([])
+ star_radius_DF_profiles.append([])
+ star_surface_density.append([])
+ sfr_radius_DF_xs.append([])
+ sfr_radius_DF_profiles.append([])
+ sfr_surface_density.append([])
+ if draw_star_radius_DF == 2:
+ draw_radius_DF = 1
+ if draw_radius_DF == 1:
+ radius_DF_xs.append([])
+ radius_DF_profiles.append([])
+ surface_density.append([])
+ if draw_height_DF == 1:
+ height_DF_xs.append([])
+ height_DF_profiles.append([])
+ height_surface_density.append([])
+ if draw_SFR >= 1:
+ sfr_ts.append([])
+ sfr_cum_masses.append([])
+ sfr_sfrs.append([])
+ if draw_cut_through == 1:
+ cut_through_zs.append([])
+ cut_through_zvalues.append([])
+ cut_through_xs.append([])
+ cut_through_xvalues.append([])
+
+for time in range(len(times)):
+
+ for code in range(len(codes)):
+
+ if (dataset_num != 1 and dataset_num != 2) or (filenames[dataset_num-1][code] == []):
+ continue
+
+ ####################################
+ # PRE-ANALYSIS STEPS #
+ ####################################
+
+ # LOAD DATASETS
+ pf = load_dataset(dataset_num, time, code, codes, filenames[dataset_num-1])
+
+ # PARTICLE FILED NAMES FOR SPH CODES, AND STELLAR PARTICLE FILTERS
+ PartType_Gas_to_use = "Gas"
+ PartType_Star_to_use = "Stars"
+ PartType_StarBeforeFiltered_to_use = "Stars"
+ MassType_to_use = "Mass"
+ MetallicityType_to_use = "Metallicity"
+ FormationTimeType_to_use = "StellarFormationTime" # for GADGET/GEAR/GIZMO, this field has to be added in frontends/sph/fields.py, in which only "FormationTime" can be recognized
+ SmoothingLengthType_to_use = "SmoothingLength"
+
+ if codes[code] == 'CHANGA' or codes[code] == 'GASOLINE':
+ MetallicityType_to_use = "Metals"
+ PartType_Star_to_use = "NewStars"
+ PartType_StarBeforeFiltered_to_use = "Stars"
+ FormationTimeType_to_use = "FormationTime"
+ SmoothingLengthType_to_use = "smoothing_length"
+ def NewStars(pfilter, data): # see http://yt-project.org/docs/dev/analyzing/filtering.html#filtering-particle-fields
+ return (data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0)
+ add_particle_filter(PartType_Star_to_use, function=NewStars, filtered_type=PartType_StarBeforeFiltered_to_use, requires=[FormationTimeType_to_use])
+ pf.add_particle_filter(PartType_Star_to_use)
+ pf.periodicity = (True, True, True) # this is needed especially when bPeriodic = 0 in GASOLINE, to avoid RuntimeError in geometry/selection_routines.pyx:855
+ elif codes[code] == 'ART-I':
+ PartType_StarBeforeFiltered_to_use = "stars"
+ FormationTimeType_to_use = "particle_creation_time"
+ def Stars(pfilter, data):
+ return (data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0)
+# return ((data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0) & (data[(pfilter.filtered_type, "particle_index")] >= 212500)) # without the fix metioned above, the above line won't work because all "stars"="specie1" have the same wrong particle_creation_time of 6.851 Gyr; so you will have to use this quick and dirty patch
+ add_particle_filter(PartType_Star_to_use, function=Stars, filtered_type=PartType_StarBeforeFiltered_to_use, requires=[FormationTimeType_to_use, "particle_index"])
+ pf.add_particle_filter(PartType_Star_to_use)
+ elif codes[code] == 'ART-II':
+ PartType_StarBeforeFiltered_to_use = "STAR"
+ if time != 0: # BIRTH_TIME field exists in ART-II but as a dimensionless quantity for some reason in frontends/artio/fields.py; so we create StellarFormationTime field
+ def _FormationTime(field, data):
+ return pf.arr(data["STAR", "BIRTH_TIME"].d, 'code_time')
+ pf.add_field(("STAR", FormationTimeType_to_use), function=_FormationTime, particle_type=True, take_log=False, units="code_time")
+ def Stars(pfilter, data):
+ return (data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0)
+ add_particle_filter(PartType_Star_to_use, function=Stars, filtered_type=PartType_StarBeforeFiltered_to_use, requires=[FormationTimeType_to_use])
+ pf.add_particle_filter(PartType_Star_to_use)
+ elif codes[code] == 'ENZO':
+ PartType_StarBeforeFiltered_to_use = "all"
+ FormationTimeType_to_use = "creation_time"
+ def Stars(pfilter, data):
+ return ((data[(pfilter.filtered_type, "particle_type")] == 2) & (data[(pfilter.filtered_type, FormationTimeType_to_use)] > 0))
+ add_particle_filter(PartType_Star_to_use, function=Stars, filtered_type=PartType_StarBeforeFiltered_to_use, requires=["particle_type", FormationTimeType_to_use])
+ pf.add_particle_filter(PartType_Star_to_use)
+ elif codes[code] == "GADGET-3":
+ PartType_Gas_to_use = "PartType0"
+ PartType_Star_to_use = "PartType4"
+ PartType_StarBeforeFiltered_to_use = "PartType4"
+ MassType_to_use = "Masses"
+ elif codes[code] == "SWIFT":
+ PartType_Gas_to_use = "PartType0"
+ PartType_Star_to_use = "PartType2"
+ PartType_StarBeforeFiltered_to_use = "PartType2"
+ MassType_to_use = "Masses"
+ elif codes[code] == "GEAR":
+ PartType_Gas_to_use = "PartType0"
+ PartType_Star_to_use = "PartType1"
+ PartType_StarBeforeFiltered_to_use = "PartType1"
+ MassType_to_use = "Masses"
+ elif codes[code] == 'RAMSES':
+ PartType_StarBeforeFiltered_to_use = "all"
+ pf.current_time = pf.arr(pf.parameters['time'], 'code_time') # reset pf.current_time because it is incorrectly set up in frontends/ramses/data_structure.py, and I don't wish to mess with units there
+ FormationTimeType_to_use = "particle_age" # particle_age field actually means particle creation time, at least for this particular dataset, so the new field below is not needed
+ # if time != 0: # Only particle_age field exists in RAMSES (only for new stars + IC stars), so we create StellarFormationTime field
+ # def _FormationTime(field, data):
+ # return pf.current_time - data["all", "particle_age"].in_units("s")
+ # pf.add_field(("all", FormationTimeType_to_use), function=_FormationTime, particle_type=True, take_log=False, units="code_time")
+ def Stars(pfilter, data):
+ return (data[(pfilter.filtered_type, "particle_age")] > 0)
+ add_particle_filter(PartType_Star_to_use, function=Stars, filtered_type=PartType_StarBeforeFiltered_to_use, requires=["particle_age"])
+ pf.add_particle_filter(PartType_Star_to_use)
+
+ # AXIS SWAP FOR PLOT COLLECTION
+ pf.coordinates.x_axis[1] = 0
+ pf.coordinates.y_axis[1] = 2
+ pf.coordinates.x_axis['y'] = 0
+ pf.coordinates.y_axis['y'] = 2
+
+ # ADDITIONAL FIELDS I: GLOBALLY USED FIELDS
+ def _density_squared(field, data):
+ return data[("gas", "density")]**2
+ pf.add_field(("gas", "density_squared"), function=_density_squared, units="g**2/cm**6")
+
+ # ADDITIONAL FIELDS II: FOR CELL SIZE AND RESOLUTION MAPS
+ def _CellSizepc(field,data):
+ return (data[("index", "cell_volume")].in_units('pc**3'))**(1/3.)
+ pf.add_field(("index", "cell_size"), function=_CellSizepc, units='pc', display_name="Resolution $\Delta$ x", take_log=True )
+ def _Inv2CellVolumeCode(field,data):
+ return data[("index", "cell_volume")]**-2
+ pf.add_field(("index", "cell_volume_inv2"), function=_Inv2CellVolumeCode, units='code_length**(-6)', display_name="Inv2CellVolumeCode", take_log=True)
+ if draw_cellsize_map >= 2:
+ if codes[code] == 'CHANGA' or codes[code] == 'GEAR' or codes[code] == 'GADGET-3' or codes[code] == 'GASOLINE' or codes[code] == 'GIZMO' or codes[code] == "SWIFT":
+ def _ParticleSizepc(field, data):
+ return (data[(PartType_Gas_to_use, MassType_to_use)]/data[(PartType_Gas_to_use, "Density")])**(1./3.)
+ pf.add_field((PartType_Gas_to_use, "particle_size"), function=_ParticleSizepc, units="pc", display_name="Resolution $\Delta$ x", particle_type=True, take_log=True)
+ def _Inv2ParticleVolumepc(field, data):
+ return (data[(PartType_Gas_to_use, MassType_to_use)]/data[(PartType_Gas_to_use, "Density")])**(-2.)
+ pf.add_field((PartType_Gas_to_use, "particle_volume_inv2"), function=_Inv2ParticleVolumepc, units="pc**(-6)", display_name="Inv2ParticleVolumepc", particle_type=True, take_log=True)
+ # Also creating smoothed field following an example in yt-project.org/docs/dev/cookbook/calculating_information.html; use hardcoded num_neighbors as in frontends/gadget/fields.py
+ fn = add_volume_weighted_smoothed_field(PartType_Gas_to_use, "Coordinates", MassType_to_use, SmoothingLengthType_to_use, "Density", "particle_size", pf.field_info, nneighbors=64)
+ fn = add_volume_weighted_smoothed_field(PartType_Gas_to_use, "Coordinates", MassType_to_use, SmoothingLengthType_to_use, "Density", "particle_volume_inv2", pf.field_info, nneighbors=64)
+
+ # Alias doesn't work -- e.g. pf.field_info.alias(("gas", "particle_size"), fn[0]) -- check alias below; so I simply add ("gas", "particle_size")
+ def _ParticleSizepc_2(field, data):
+ return data["deposit", PartType_Gas_to_use+"_smoothed_"+"particle_size"]
+ pf.add_field(("gas", "particle_size"), function=_ParticleSizepc_2, units="pc", force_override=True, display_name="Resolution $\Delta$ x", particle_type=False, take_log=True)
+ def _Inv2ParticleVolumepc_2(field, data):
+ return data["deposit", PartType_Gas_to_use+"_smoothed_"+"particle_volume_inv2"]
+ pf.add_field(("gas", "particle_volume_inv2"), function=_Inv2ParticleVolumepc_2, units="pc**(-6)", force_override=True, display_name="Inv2ParticleVolumepc", particle_type=False, take_log=True)
+
+ # ADDITIONAL FIELDS III: TEMPERATURE
+ if codes[code] == 'GEAR' or codes[code] == 'GADGET-3' or codes[code] == 'RAMSES' or codes[code] == "SWIFT":
+ # From grackle/src/python/utilities/convenience.py: Calculate a tabulated approximation to mean molecular weight (valid for data that used Grackle 2.0 or below)
+ def calc_mu_table_local(temperature):
+ tt = np.array([1.0e+01, 1.0e+02, 1.0e+03, 1.0e+04, 1.3e+04, 2.1e+04, 3.4e+04, 6.3e+04, 1.0e+05, 1.0e+09])
+ mt = np.array([1.18701555, 1.15484424, 1.09603514, 0.9981496, 0.96346395, 0.65175895, 0.6142901, 0.6056833, 0.5897776, 0.58822635])
+ logttt= np.log(temperature)
+ logmu = np.interp(logttt,np.log(tt),np.log(mt)) # linear interpolation in log-log space
+ return np.exp(logmu)
+ temperature_values = []
+ mu_values = []
+ T_over_mu_values = []
+ current_temperature = 1e1
+ final_temperature = 1e7
+ dlogT = 0.1
+ while current_temperature < final_temperature:
+ temperature_values.append(current_temperature)
+ current_mu = calc_mu_table_local(current_temperature)
+ mu_values.append(current_mu)
+ T_over_mu_values.append(current_temperature/current_mu)
+ current_temperature = np.exp(np.log(current_temperature)+dlogT)
+ def convert_T_over_mu_to_T(T_over_mu):
+ logT_over_mu = np.log(T_over_mu)
+ logT = np.interp(logT_over_mu, np.log(T_over_mu_values), np.log(temperature_values)) # linear interpolation in log-log space
+ return np.exp(logT)
+ if codes[code] == 'GEAR' or codes[code] == 'GADGET-3' or codes[code] == "SWIFT":
+ def _Temperature_3(field, data):
+ gamma = 5.0/3.0
+ T_over_mu = (data[PartType_Gas_to_use, "InternalEnergy"] * (gamma-1) * constants.mass_hydrogen_cgs / constants.boltzmann_constant_cgs).in_units('K').d # T/mu
+ return YTArray(convert_T_over_mu_to_T(T_over_mu), 'K') # now T
+ pf.add_field((PartType_Gas_to_use, "Temperature"), function=_Temperature_3, particle_type=True, force_override=True, units="K")
+ elif codes[code] == 'RAMSES':
+ # The pressure field includes the artificial pressure support term, so one needs to be careful (compare with the exsiting frontends/ramses/fields.py)
+ def _temperature_3(field, data):
+ T_J = 1800.0 # in K
+ n_J = 8.0 # in H/cc
+ gamma_0 = 2.0
+ x_H = 0.76
+ mH = 1.66e-24 # from pymses/utils/constants/__init__.py (vs. in yt, mass_hydrogen_cgs = 1.007947*amu_cgs = 1.007947*1.660538921e-24 = 1.6737352e-24)
+ kB = 1.3806504e-16 # from pymses/utils/constants/__init__.py (vs. in yt, boltzmann_constant_cgs = 1.3806488e-16)
+ if time != 0:
+ T_over_mu = data["gas", "pressure"].d/data["gas", "density"].d * mH / kB - T_J * (data["gas", "density"].d * x_H / mH / n_J)**(gamma_0 - 1.0) # T/mu = T2 in Ramses
+ else:
+ T_over_mu = data["gas", "pressure"].d/data["gas", "density"].d * mH / kB # IC: no pressure support
+ return YTArray(convert_T_over_mu_to_T(T_over_mu), 'K') # now T
+ pf.add_field(("gas", "temperature"), function=_temperature_3, force_override=True, units="K")
+
+ # ADDITIONAL FIELDS IV: METALLICITY (IN MASS FRACTION, NOT IN ZSUN)
+ if draw_metal_map >= 1 or draw_metal_PDF == 1:
+ if codes[code] == 'ART-I': # metallicity field in ART-I has a different meaning (see frontends/art/fields.py), and metallicity field in ART-II is missing
+ def _metallicity_2(field, data):
+ return (data["gas", "metal_ii_density"] + data["gas", "metal_ia_density"]) / data["gas", "density"] # metal_ia_density needs to be added to account for initial 0.5 Zsun, even though we don't have SNe Ia
+ pf.add_field(("gas", "metallicity"), function=_metallicity_2, force_override=True, display_name="Metallicity", take_log=True, units="")
+ elif codes[code] == 'ART-II':
+ def _metallicity_2(field, data):
+ return data["gas", "metal_ii_density"] / data["gas", "density"]
+ pf.add_field(("gas", "metallicity"), function=_metallicity_2, force_override=True, display_name="Metallicity", take_log=True, units="")
+ elif codes[code] == 'ENZO': # metallicity field in ENZO is in Zsun, so we create a new field
+ def _metallicity_2(field, data):
+ return data["gas", "metal_density"] / data["gas", "density"]
+ pf.add_field(("gas", "metallicity"), function=_metallicity_2, force_override=True, display_name="Metallicity", take_log=True, units="")
+ elif codes[code] == 'GEAR': # "Metals" in GEAR is 10-species field ([:,9] is the total metal fraction), so requires a change in _vector_fields in frontends/gadget/io.py: added ("Metals", 10)
+ def _metallicity_2(field, data):
+ if len(data[PartType_Gas_to_use, "Metals"].shape) == 1:
+ return data[PartType_Gas_to_use, "Metals"]
+ else:
+ return data[PartType_Gas_to_use, "Metals"][:,9].in_units("") # in_units("") turned out to be crucial!; otherwise code_metallicity will be used and it will mess things up
+ # We are creating ("Gas", "Metallicity") here, different from ("Gas", "metallicity") which is auto-generated by yt but doesn't work properly
+ pf.add_field((PartType_Gas_to_use, MetallicityType_to_use), function=_metallicity_2, display_name="Metallicity", particle_type=True, take_log=True, units="")
+ # Also creating smoothed field following an example in yt-project.org/docs/dev/cookbook/calculating_information.html; use hardcoded num_neighbors as in frontends/gadget/fields.py
+ fn = add_volume_weighted_smoothed_field(PartType_Gas_to_use, "Coordinates", MassType_to_use, "SmoothingLength", "Density", MetallicityType_to_use, pf.field_info, nneighbors=64)
+ # Alias doesn't work -- e.g. pf.field_info.alias(("gas", "metallicity"), fn[0]) -- probably because pf=GadgetDataset()?, not load()?; so I add and replace existing ("gas", "metallicity")
+ def _metallicity_3(field, data):
+ return data["deposit", PartType_Gas_to_use+"_smoothed_"+MetallicityType_to_use]
+ pf.add_field(("gas", "metallicity"), function=_metallicity_3, force_override=True, display_name="Metallicity", particle_type=False, take_log=True, units="")
+ if draw_metal_map >= 2:
+ def _metal_mass(field, data):
+ return data["gas", "cell_mass"] * data["gas", "metallicity"]
+ pf.add_field(("gas", "metal_mass"), function=_metal_mass, force_override=True, display_name="Metal Mass", take_log=True, units="Msun")
+
+ # ADDITIONAL FIELDS V: FAKE PARTICLE FIELDS, CYLINDRICAL COORDINATES, etc.
+ def rho_agora_disk(r, z):
+ r_d = YTArray(3.432, 'kpc')
+ z_d = 0.1*r_d
+ M_d = YTArray(4.297e10, 'Msun')
+ f_gas = 0.2
+ r_0 = (f_gas*M_d / (4.0*np.pi*r_d**2*z_d)).in_units('g/cm**3')
+ return r_0*numpy.exp(-r/r_d)*numpy.exp(-z/z_d)
+
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ def _cylindrical_z_abs(field, data):
+ return numpy.abs(data[("index", "cylindrical_z")])
+ pf.add_field(("index", "cylindrical_z_abs"), function=_cylindrical_z_abs, take_log=False, particle_type=False, units="cm")
+ def _density_minus_analytic(field, data):
+ return data[("gas", "density")] - rho_agora_disk(data[("index", "cylindrical_r")], data[("index", "cylindrical_z_abs")])
+ pf.add_field(("gas", "density_minus_analytic"), function=_density_minus_analytic, take_log=False, particle_type=False, display_name="density residual abs", units="g/cm**3")
+ else:
+ def _particle_position_cylindrical_z_abs(field, data):
+ return numpy.abs(data[(PartType_Gas_to_use, "particle_position_cylindrical_z")])
+ pf.add_field((PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), function=_particle_position_cylindrical_z_abs, take_log=False, particle_type=True, units="cm")
+ # particle_type=False doesn't make sense, but is critical for PhasePlot/ProfilePlot to work for particle fields
+ # requires a change in data_objects/data_container.py: remove raise YTFieldTypeNotFound(ftype)
+ def _Density_2(field, data):
+ return data[(PartType_Gas_to_use, "Density")].in_units('g/cm**3')
+ pf.add_field((PartType_Gas_to_use, "Density_2"), function=_Density_2, take_log=True, particle_type=False, display_name="Density", units="g/cm**3")
+ def _Temperature_2(field, data):
+ return data[(PartType_Gas_to_use, "Temperature")].in_units('K')
+ pf.add_field((PartType_Gas_to_use, "Temperature_2"), function=_Temperature_2, take_log=True, particle_type=False, display_name="Temperature", units="K")
+ def _Mass_2(field, data):
+ return data[(PartType_Gas_to_use, MassType_to_use)].in_units('Msun')
+ pf.add_field((PartType_Gas_to_use, "Mass_2"), function=_Mass_2, take_log=True, particle_type=False, display_name="Mass", units="Msun")
+ if draw_metal_map >= 1 or draw_metal_PDF == 1:
+ def _Metallicity_2(field, data):
+ return data[(PartType_Gas_to_use, MetallicityType_to_use)]
+ pf.add_field((PartType_Gas_to_use, "Metallicity_2"), function=_Metallicity_2, take_log=True, particle_type=False, display_name="Metallicity", units="")
+ def _Density_2_minus_analytic(field, data):
+ return data[(PartType_Gas_to_use, "density")] - rho_agora_disk(data[(PartType_Gas_to_use, "particle_position_cylindrical_radius")], data[(PartType_Gas_to_use, "particle_position_cylindrical_z_abs")])
+ pf.add_field((PartType_Gas_to_use, "Density_2_minus_analytic"), function=_Density_2_minus_analytic, take_log=False, particle_type=False, display_name="Density Residual", units="g/cm**3")
+
+ ####################################
+ # MAIN ANALYSIS ROUTINES #
+ ####################################
+
+ # FIND CENTER AND PROJ_REGION
+ v, cen = pf.h.find_max(("gas", "density")) # find the center to keep the galaxy at the center of all the images; here we assume that the gas disk is no bigger than 30 kpc in radius
+ sp = pf.sphere(cen, (30.0, "kpc"))
+ cen2 = sp.quantities.center_of_mass(use_gas=True, use_particles=False).in_units("kpc")
+ sp2 = pf.sphere(cen2, (1.0, "kpc"))
+ cen3 = sp2.quantities.max_location(("gas", "density"))
+ center = pf.arr([cen3[1].d, cen3[2].d, cen3[3].d], 'code_length') # naive usage such as YTArray([cen3[1], cen3[2], cen3[3]]) doesn't work somehow for ART-II data
+ if yt_version_pre_3_2_3 == 1:
+ center = pf.arr([cen3[2].d, cen3[3].d, cen3[4].d], 'code_length') # for yt-3.2.3 or before
+ if codes[code] == "GASOLINE" and dataset_num == 2 and time == 1:
+ #center = pf.arr([2.7912903206399826e+20, 1.5205303149849894e+21, 1.5398968883245956e+21], 'cm') # a temporary hack for GASOLINE (center of the most massive stellar clump)
+ center = pf.arr([ 0.09045956, 0.49277032, 0.49904659], 'code_length')
+ # if codes[code] == "GADGET-3" and dataset_num == 2 and time == 1:
+ # #center = pf.arr([6.184520935812296e+21, 4.972678132728082e+21, 6.559067311284074e+21], 'cm') # a temporary hack for GADGET-3 (center of the most massive stellar clump)
+ # center = pf.arr([ 2.00426673674, 1.61153523084, 2.12564895041], 'code_length')
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ proj_region = pf.box(center - YTArray([figure_width, figure_width, figure_width], 'kpc'),
+ center + YTArray([figure_width, figure_width, figure_width], 'kpc')) # projected images made using a (2*figure_width)^3 box for AMR codes
+ else:
+ proj_region = pf.all_data()
+
+ # DENSITY MAPS
+ if draw_density_map == 1:
+ my_cmap = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap.set_bad(my_cmap(0)) # cmap range [0, 1)
+ my_cmap.set_under(my_cmap(0))
+ for ax in range(1, 3):
+ p = ProjectionPlot(pf, ax, ("gas", "density"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = None, fontsize=9)
+ p.set_zlim(("gas", "density"), 1e-5, 1e-1)
+ p.set_cmap(("gas", "density"), my_cmap)
+ plot = p.plots[("gas", "density")]
+
+ plot.figure = fig_density_map[time]
+ plot.axes = grid_density_map[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot.cax = grid_density_map[time].cbar_axes[0]
+ p._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_density_map[time][code].axes.add_artist(at)
+
+ # TEMPERATURE MAPS
+ if draw_temperature_map == 1:
+ my_cmap = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap.set_bad(my_cmap(0.6)) # (log(1e4) - log(1e1)) / (log(1e6) - log(1e1)) = 0.6
+ my_cmap.set_under(my_cmap(0))
+ for ax in range(1, 3):
+ p2 = ProjectionPlot(pf, ax, ("gas", "temperature"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = ("gas", "density_squared"), fontsize=9)
+ p2.set_zlim(("gas", "temperature"), 1e1, 1e6)
+ p2.set_cmap(("gas", "temperature"), my_cmap)
+ plot2 = p2.plots[("gas", "temperature")]
+
+ plot2.figure = fig_temperature_map[time]
+ plot2.axes = grid_temperature_map[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot2.cax = grid_temperature_map[time].cbar_axes[0]
+ p2._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_temperature_map[time][code].axes.add_artist(at)
+
+ # CELL-SIZE MAPS
+ if draw_cellsize_map >= 1:
+ my_cmap = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap.set_bad(my_cmap(0.99999)) # 1 doesn't work!
+ my_cmap.set_under(my_cmap(0))
+ if draw_cellsize_map == 1 or draw_cellsize_map == 3:
+ for ax in range(1, 3):
+ p25 = ProjectionPlot(pf, ax, ("index", "cell_size"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = ("index", "cell_volume_inv2"), fontsize=9)
+ p25.set_zlim(("index", "cell_size"), 10, 1e3)
+ p25.set_cmap(("index", "cell_size"), my_cmap)
+ plot25 = p25.plots[("index", "cell_size")]
+
+ plot25.figure = fig_cellsize_map[time]
+ plot25.axes = grid_cellsize_map[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot25.cax = grid_cellsize_map[time].cbar_axes[0]
+ p25._setup_plots()
+
+ # Create another map with a different resolution definition if requested
+ if draw_cellsize_map == 2 or draw_cellsize_map == 3:
+ for ax in range(1, 3):
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p251 = ProjectionPlot(pf, ax, ("index", "cell_size"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), \
+ weight_field = ("index", "cell_volume_inv2"), fontsize=9)
+ p251.set_zlim(("index", "cell_size"), 10, 1e3)
+ p251.set_cmap(("index", "cell_size"), my_cmap)
+ plot251 = p251.plots[("index", "cell_size")]
+ else:
+ p251 = ProjectionPlot(pf, ax, ("gas", "particle_size"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), \
+ weight_field = ("gas", "particle_volume_inv2"), fontsize=9)
+ p251.set_zlim(("gas", "particle_size"), 10, 1e3)
+ p251.set_cmap(("gas", "particle_size"), my_cmap)
+ plot251 = p251.plots[("gas", "particle_size")]
+
+ plot251.figure = fig_cellsize_map_2[time]
+ plot251.axes = grid_cellsize_map_2[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot251.cax = grid_cellsize_map_2[time].cbar_axes[0]
+ p251._setup_plots()
+
+ if add_nametag == 1:
+ if draw_cellsize_map == 1 or draw_cellsize_map == 3:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_cellsize_map[time][code].axes.add_artist(at)
+ if draw_cellsize_map == 2 or draw_cellsize_map == 3:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_cellsize_map_2[time][code].axes.add_artist(at)
+
+ # ELEVATION MAPS
+ if draw_elevation_map == 1:
+ def _CellzElevationpc(field,data):
+ return ((data[("index", "z")] - center[2]).in_units('pc'))
+ pf.add_field(("index", "z_elevation"), function=_CellzElevationpc, units='kpc', display_name="$z$ Elevation", take_log=False)
+ my_cmap = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap.set_bad(my_cmap(0.5))
+ my_cmap.set_under(my_cmap(0))
+ for ax in range(2, 3):
+ p255 = ProjectionPlot(pf, ax, ("index", "z_elevation"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = ("gas", "density"), fontsize=9)
+ p255.set_zlim(("index", "z_elevation"), -1, 1)
+ p255.set_cmap(("index", "z_elevation"), my_cmap)
+ plot255 = p255.plots[("index", "z_elevation")]
+
+ plot255.figure = fig_elevation_map[time]
+ plot255.axes = grid_elevation_map[time][(ax-2)*len(codes)+code].axes
+ if code == 0: plot255.cax = grid_elevation_map[time].cbar_axes[0]
+ p255._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_elevation_map[time][code].axes.add_artist(at)
+
+ # Z-VELOCITY MAPS
+ if draw_zvel_map == 1:
+ my_cmap = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap.set_bad(my_cmap(0.5))
+ my_cmap.set_under(my_cmap(0))
+ for ax in range(1, 3):
+ p26 = ProjectionPlot(pf, ax, ("gas", "velocity_z"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = ("gas", "density_squared"), fontsize=9)
+ p26.set_zlim(("gas", "velocity_z"), -1e7, 1e7)
+ p26.set_cmap(("gas", "velocity_z"), my_cmap)
+ plot26 = p26.plots[("gas", "velocity_z")]
+
+ plot26.figure = fig_zvel_map[time]
+ plot26.axes = grid_zvel_map[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot26.cax = grid_zvel_map[time].cbar_axes[0]
+ p26._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_zvel_map[time][code].axes.add_artist(at)
+
+ # METAL MAPS
+ if draw_metal_map >= 1:
+ my_cmap = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap.set_bad(my_cmap(0.347)) # (0.02041 - 0.01) / (0.04 - 0.01) = 0.347
+ my_cmap.set_under(my_cmap(0))
+ for ax in range(1, 3):
+ pf.field_info[("gas", "metallicity")].take_log = False
+ p26 = ProjectionPlot(pf, ax, ("gas", "metallicity"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = ("gas", "density_squared"), fontsize=9)
+ p26.set_zlim(("gas", "metallicity"), 0.01, 0.04)
+ p26.set_cmap(("gas", "metallicity"), my_cmap)
+ plot26 = p26.plots[("gas", "metallicity")]
+
+ plot26.figure = fig_metal_map[time]
+ plot26.axes = grid_metal_map[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot26.cax = grid_metal_map[time].cbar_axes[0]
+ p26._setup_plots()
+
+ if add_nametag == 1:
+ if draw_metal_map == 2:
+ sp = pf.sphere(center, (1, "Mpc"))
+ total_metal_mass = sp.quantities.total_quantity([("gas", "metal_mass")])
+ at = AnchoredText("%s - %e" % (codes[code], total_metal_mass), loc=2, prop=dict(size=6), frameon=True)
+ else:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_metal_map[time][code].axes.add_artist(at)
+
+ # STELLAR MAPS
+ if draw_star_map == 1 and time != 0:
+ for ax in range(1, 3):
+ p27 = ParticleProjectionPlot(pf, ax, (PartType_Star_to_use, "particle_mass"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = None, fontsize=9)
+ p27.set_unit((PartType_Star_to_use, "particle_mass"), 'Msun')
+ p27.set_zlim((PartType_Star_to_use, "particle_mass"), 1e4, 1e7)
+ p27.set_cmap((PartType_Star_to_use, "particle_mass"), 'algae')
+ p27.set_buff_size(400) # default is 800
+ p27.set_colorbar_label((PartType_Star_to_use, "particle_mass"), "Stellar Mass Per Pixel ($\mathrm{M}_{\odot}$)")
+ plot27 = p27.plots[(PartType_Star_to_use, "particle_mass")]
+ # p27 = ParticleProjectionPlot(pf, ax, (PartType_Star_to_use, "particle_velocity_cylindrical_theta"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = (PartType_Star_to_use, "particle_mass"), fontsize=9)
+ # p27.set_unit((PartType_Star_to_use, "particle_velocity_cylindrical_theta"), 'km/s')
+ # p27.set_buff_size(400)
+ # p27.set_colorbar_label((PartType_Star_to_use, "particle_velocity_cylindrical_theta"), "Rotational Velocity (km/s)")
+ # plot27 = p27.plots[(PartType_Star_to_use, "particle_velocity_cylindrical_theta")]
+
+ plot27.figure = fig_star_map[time]
+ plot27.axes = grid_star_map[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot27.cax = grid_star_map[time].cbar_axes[0]
+ p27._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_star_map[time][code].axes.add_artist(at)
+
+ # STELLAR CLUMP STATISTICS AND ANNOTATED STELLAR MAPS
+ if draw_star_clump_stats >= 1 and time != 0:
+ # For make yt's HaloCatalog to work with non-cosmological dataset, a fix needed to be applied to analysis_modules/halo_finding/halo_objects.py: self.period = ds.arr([3.254, 3.254, 3.254], 'Mpc')
+ pf.hubble_constant = 0.71; pf.omega_lambda = 0.73; pf.omega_matter = 0.27; pf.omega_curvature = 0.0; pf.current_redshift = 0.0 # another trick to make HaloCatalog work especially for ART-I dataset
+ if os.path.exists("./halo_catalogs/hop_%s_%05d/hop_%s_%05d.0.h5" % (codes[code], times[time], codes[code], times[time])) == False:
+ hc = HaloCatalog(data_ds=pf, finder_method='hop', output_dir="./halo_catalogs/hop_%s_%05d" % (codes[code], times[time]), \
+ finder_kwargs={'threshold': 2e8, 'dm_only': False, 'ptype': PartType_Star_to_use})
+# finder_kwargs={'threshold': 2e5, 'dm_only': False, 'ptype': "all"})
+ hc.add_filter('quantity_value', 'particle_mass', '>', 2.6e6, 'Msun') # exclude halos with less than 30 particles
+ hc.add_filter('quantity_value', 'particle_mass', '<', 2.6e8, 'Msun') # exclude the most massive halo (threshold 1e8.4 is hand-picked, so one needs to be careful!)
+ hc.create()
+ halo_ds = load("./halo_catalogs/hop_%s_%05d/hop_%s_%05d.0.h5" % (codes[code], times[time], codes[code], times[time]))
+ hc = HaloCatalog(halos_ds=halo_ds, output_dir="./halo_catalogs/hop_%s_%05d" % (codes[code], times[time]))
+ hc.load()
+
+ halo_ad = hc.halos_ds.all_data()
+ star_clump_masses_hop[time].append(np.log10(halo_ad['particle_mass'][:].in_units("Msun")))
+
+ if os.path.exists("./halo_catalogs/fof_%s_%05d/fof_%s_%05d.0.h5" % (codes[code], times[time], codes[code], times[time])) == False:
+ hc2 = HaloCatalog(data_ds=pf, finder_method='fof', output_dir="./halo_catalogs/fof_%s_%05d" % (codes[code], times[time]), \
+ finder_kwargs={'link': 0.0025, 'dm_only': False, 'ptype': PartType_Star_to_use})
+# finder_kwargs={'link': 0.02, 'dm_only': False, 'ptype': "all"})
+ hc2.add_filter('quantity_value', 'particle_mass', '>', 2.6e6, 'Msun') # exclude halos with less than 30 particles
+ hc2.add_filter('quantity_value', 'particle_mass', '<', 2.6e8, 'Msun') # exclude the most massive halo (threshold 1e8.4 is hand-picked, so one needs to be careful!)
+ hc2.create()
+ halo_ds2 = load("./halo_catalogs/fof_%s_%05d/fof_%s_%05d.0.h5" % (codes[code], times[time], codes[code], times[time]))
+ hc2 = HaloCatalog(halos_ds=halo_ds2, output_dir="./halo_catalogs/fof_%s_%05d" % (codes[code], times[time]))
+ hc2.load()
+
+ halo_ad2 = hc2.halos_ds.all_data()
+ star_clump_masses_fof[time].append(np.log10(halo_ad2['particle_mass'][:].in_units("Msun")))
+ # most_massive = halo_ad['particle_mass'].max().in_units('Msun')
+ # cen4 = [halo_ad['particle_position_x'][halo_ad['particle_mass'] == most_massive][0].in_units('code_length').d,
+ # halo_ad['particle_position_y'][halo_ad['particle_mass'] == most_massive][0].in_units('code_length').d,
+ # halo_ad['particle_position_z'][halo_ad['particle_mass'] == most_massive][0].in_units('code_length').d]
+
+
+ # Add additional star_map with annotated clumps if requested
+ if draw_star_clump_stats >= 2:
+ for ax in range(1, 3):
+ p271 = ParticleProjectionPlot(pf, ax, (PartType_Star_to_use, "particle_mass"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = None, fontsize=9)
+ p271.set_unit((PartType_Star_to_use, "particle_mass"), 'Msun')
+ p271.set_zlim((PartType_Star_to_use, "particle_mass"), 1e4, 1e7)
+ p271.set_cmap((PartType_Star_to_use, "particle_mass"), 'algae')
+ p271.set_buff_size(400) # default is 800
+ p271.set_colorbar_label((PartType_Star_to_use, "particle_mass"), "Stellar Mass Per Pixel ($\mathrm{M}_{\odot}$)")
+ plot271 = p271.plots[(PartType_Star_to_use, "particle_mass")]
+ if ax == 2:
+ p271.annotate_halos(hc, factor=1.0, circle_args={'linewidth':0.7, 'alpha':0.8, 'facecolor':'none', 'edgecolor':'k'})
+
+ plot271.figure = fig_star_map_2[time]
+ plot271.axes = grid_star_map_2[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot271.cax = grid_star_map_2[time].cbar_axes[0]
+ p271._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_star_map_2[time][code].axes.add_artist(at)
+
+ for ax in range(1, 3):
+ p272 = ParticleProjectionPlot(pf, ax, (PartType_Star_to_use, "particle_mass"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = None, fontsize=9)
+ p272.set_unit((PartType_Star_to_use, "particle_mass"), 'Msun')
+ p272.set_zlim((PartType_Star_to_use, "particle_mass"), 1e4, 1e7)
+ p272.set_cmap((PartType_Star_to_use, "particle_mass"), 'algae')
+ p272.set_buff_size(400) # default is 800
+ p272.set_colorbar_label((PartType_Star_to_use, "particle_mass"), "Stellar Mass Per Pixel ($\mathrm{M}_{\odot}$)")
+ plot272 = p272.plots[(PartType_Star_to_use, "particle_mass")]
+ # p272 = ParticleProjectionPlot(pf, ax, ("all", "particle_mass"), center = center, data_source=proj_region, width = (figure_width, 'kpc'), weight_field = None, fontsize=9)
+ # p272.set_unit(("all", "particle_mass"), 'Msun')
+ # p272.set_zlim(("all", "particle_mass"), 1e4, 1e9)
+ # p272.set_cmap(("all", "particle_mass"), 'algae')
+ # p272.set_buff_size(400) # default is 800
+ # p272.set_colorbar_label(("all", "particle_mass"), "Mass Per Pixel ($\mathrm{M}_{\odot}$)")
+ # plot272 = p272.plots[("all", "particle_mass")]
+ if ax == 2:
+ p272.annotate_halos(hc2, factor=1.0, circle_args={'linewidth':0.7, 'alpha':0.8, 'facecolor':'none', 'edgecolor':'k'})#, annotate_field='particle_mass')
+
+ plot272.figure = fig_star_map_3[time]
+ plot272.axes = grid_star_map_3[time][(ax-1)*len(codes)+code].axes
+ if code == 0: plot272.cax = grid_star_map_3[time].cbar_axes[0]
+ p272._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_star_map_3[time][code].axes.add_artist(at)
+
+ # SFR MAPS
+ if draw_SFR_map >= 1 and time != 0:
+ sp = pf.sphere(center, (1.0*figure_width, "kpc"))
+ xy_bins = int(float(figure_width) / (float(aperture_size_SFR_map/1000.)))
+ my_cmap = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap.set_bad(my_cmap(0))
+ my_cmap.set_under(my_cmap(0))
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ def _position_x_from_center(field, data):
+ return (data[("index", "x")] - center[0]).in_units('kpc')
+ pf.add_field(("index", "position_x_from_center"), function=_position_x_from_center, take_log=False, particle_type=False, units="kpc")
+ def _position_y_from_center(field, data):
+ return (data[("index", "y")] - center[1]).in_units('kpc')
+ pf.add_field(("index", "position_y_from_center"), function=_position_y_from_center, take_log=False, particle_type=False, units="kpc")
+ def _cell_mass_surface_density_SFR_map(field, data):
+ trans = np.zeros(data[("gas", "cell_mass")].shape)
+ ind = np.where(data[("gas", "cell_mass")] > 0)
+ trans[ind] = data[("gas", "cell_mass")][ind].in_units('Msun').d/(float(aperture_size_SFR_map))**2
+ return data.ds.arr(trans, "Msun/pc**2").in_base(data.ds.unit_system.name)
+ pf.add_field(("gas", "cell_mass_surface_density_SFR_map"), function=_cell_mass_surface_density_SFR_map, \
+ display_name="$\mathrm{\Sigma}_\mathrm{gas}$", units='Msun/pc**2', particle_type=False, take_log=True)
+
+ p28 = PhasePlot(sp, ("index", "position_x_from_center"), ("index", "position_y_from_center"), \
+ ("gas", "cell_mass_surface_density_SFR_map"), weight_field=None, fontsize=9, x_bins=xy_bins, y_bins=xy_bins)
+ p28.set_zlim(("gas", "cell_mass_surface_density_SFR_map"), 1e0, 1e3)
+ p28.set_cmap(("gas", "cell_mass_surface_density_SFR_map"), my_cmap)
+ plot28 = p28.plots[("gas", "cell_mass_surface_density_SFR_map")]
+ else:
+ def _particle_position_x_from_center(field, data):
+ return (data[(PartType_Gas_to_use, "particle_position_x")] - center[0]).in_units('kpc')
+ pf.add_field((PartType_Gas_to_use, "particle_position_x_from_center"), function=_particle_position_x_from_center, take_log=False, particle_type=True, units="kpc")
+ def _particle_position_y_from_center(field, data):
+ return (data[(PartType_Gas_to_use, "particle_position_y")] - center[0]).in_units('kpc')
+ pf.add_field((PartType_Gas_to_use, "particle_position_y_from_center"), function=_particle_position_y_from_center, take_log=False, particle_type=True, units="kpc")
+ def _particle_mass_surface_density_SFR_map(field, data):
+ trans = np.zeros(data[(PartType_Gas_to_use, "particle_mass")].shape)
+ ind = np.where(data[(PartType_Gas_to_use, "particle_mass")] > 0)
+ trans[ind] = data[(PartType_Gas_to_use, "particle_mass")][ind].in_units('Msun').d/(float(aperture_size_SFR_map))**2
+ return data.ds.arr(trans, "Msun/pc**2").in_base(data.ds.unit_system.name)
+ pf.add_field((PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"), function=_particle_mass_surface_density_SFR_map, \
+ display_name="$\mathrm{\Sigma}_\mathrm{gas}$", units='Msun/pc**2', particle_type=True, take_log=True)
+
+ p28 = ParticlePhasePlot(sp, (PartType_Gas_to_use, "particle_position_x_from_center"), (PartType_Gas_to_use, "particle_position_y_from_center"), \
+ (PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"), weight_field=None, deposition='cic', fontsize=9, x_bins=xy_bins, y_bins=xy_bins)
+ p28.set_zlim((PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"), 1e0, 1e3)
+ p28.set_cmap((PartType_Gas_to_use, "particle_mass_surface_density_SFR_map"), my_cmap)
+ plot28 = p28.plots[(PartType_Gas_to_use, "particle_mass_surface_density_SFR_map")]
+
+ p28.set_xlabel("x (kpc)")
+ p28.set_ylabel("y (kpc)")
+ p28.set_xlim(-15, 15)
+ p28.set_ylim(-15, 15)
+ plot28.figure = fig_degr_density_map[time]
+ plot28.axes = grid_degr_density_map[time][code].axes
+ if code == 0: plot28.cax = grid_degr_density_map[time].cbar_axes[0]
+ p28._setup_plots()
+
+ def _particle_position_x_from_center_2(field, data):
+ return (data[(PartType_StarBeforeFiltered_to_use, "particle_position_x")] - center[0]).in_units('kpc')
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_position_x_from_center"), function=_particle_position_x_from_center_2, take_log=False, particle_type=True, units="kpc")
+ pf.add_particle_filter(PartType_Star_to_use) # This is needed for a filtered particle type PartType_Star_to_use to work, because we have just created new particle fields.
+ def _particle_position_y_from_center_2(field, data):
+ return (data[(PartType_StarBeforeFiltered_to_use, "particle_position_y")] - center[0]).in_units('kpc')
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_position_y_from_center"), function=_particle_position_y_from_center_2, take_log=False, particle_type=True, units="kpc")
+ pf.add_particle_filter(PartType_Star_to_use)
+ def _particle_mass_young_stars_SFR_surface_density_SFR_map(field, data):
+ trans = np.zeros(data[(PartType_StarBeforeFiltered_to_use, "particle_mass")].shape)
+ ind = np.where(data[(PartType_StarBeforeFiltered_to_use, FormationTimeType_to_use)].in_units('Myr') > (pf.current_time.in_units('Myr').d - young_star_cutoff_SFR_map)) # mass for young stars only
+ trans[ind] = data[(PartType_StarBeforeFiltered_to_use, "particle_mass")][ind].in_units('Msun').d/(young_star_cutoff_SFR_map*1e6)/(float(aperture_size_SFR_map)/1000.)**2
+ return data.ds.arr(trans, "Msun/yr/kpc**2").in_base(data.ds.unit_system.name)
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_mass_young_stars_SFR_surface_density_SFR_map"), function=_particle_mass_young_stars_SFR_surface_density_SFR_map, \
+ display_name="$\mathrm{\Sigma}_\mathrm{SFR}$", units='Msun/yr/kpc**2', particle_type=True, take_log=True)
+ pf.add_particle_filter(PartType_Star_to_use)
+
+ p281 = ParticlePhasePlot(sp, (PartType_Star_to_use, "particle_position_x_from_center"), (PartType_Star_to_use, "particle_position_y_from_center"), \
+ (PartType_Star_to_use, "particle_mass_young_stars_SFR_surface_density_SFR_map"), deposition='cic', weight_field=None, fontsize=9, x_bins=xy_bins, y_bins=xy_bins)
+ p281.set_zlim((PartType_Star_to_use, "particle_mass_young_stars_SFR_surface_density_SFR_map"), 3e-4, 3e-1)
+ p281.set_cmap((PartType_Star_to_use, "particle_mass_young_stars_SFR_surface_density_SFR_map"), my_cmap)
+ plot281 = p281.plots[(PartType_Star_to_use, "particle_mass_young_stars_SFR_surface_density_SFR_map")]
+
+ p281.set_xlabel("x (kpc)")
+ p281.set_ylabel("y (kpc)")
+ p281.set_xlim(-15, 15)
+ p281.set_ylim(-15, 15)
+ plot281.figure = fig_degr_sfr_map[time]
+ plot281.axes = grid_degr_sfr_map[time][code].axes
+ if code == 0: plot281.cax = grid_degr_sfr_map[time].cbar_axes[0]
+ p281._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_degr_density_map[time][code].axes.add_artist(at)
+ at2 = AnchoredText("%s" % codes[code], loc=2, prop=dict(size=6), frameon=True)
+ grid_degr_sfr_map[time][code].axes.add_artist(at2)
+
+ # Record degraded maps for the later use in a local K-S plot
+ if draw_SFR_map == 2 and time != 0:
+ fn = p28.profile.save_as_dataset() # see http://yt-project.org/docs/dev/reference/api/generated/yt.data_objects.profiles.ProfileND.save_as_dataset.html
+ phaseplot_ds = load(fn)
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ phaseplot_surf_dens = phaseplot_ds.data[('data', "cell_mass_surface_density_SFR_map")]
+ else:
+ phaseplot_surf_dens = phaseplot_ds.data[('data', "particle_mass_surface_density_SFR_map")]
+ fn_1 = p281.profile.save_as_dataset()
+ phaseplot_ds_1 = load(fn_1)
+ phaseplot_sfr_surf_dens = phaseplot_ds_1.data[('data', "particle_mass_young_stars_SFR_surface_density_SFR_map")]
+ xy_shape = phaseplot_surf_dens.shape[0]
+ surf_dens_SFR_map[time].append(phaseplot_surf_dens.reshape(1, xy_shape**2)[0])
+ sfr_surf_dens_SFR_map[time].append(phaseplot_sfr_surf_dens.reshape(1, xy_shape**2)[0])
+ call(["rm", "%s_%s.h5" % (str(pf), 'ParticleProfile')]) # Remove unwanted .h5 files saved
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ call(["rm", "%s_%s.h5" % (str(pf), 'Profile2D')])
+
+ # DENSITY-TEMPERATURE PDF
+ if draw_PDF >= 1:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ my_cmap2 = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap2.set_under(my_cmap2(0))
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p3 = PhasePlot(sp, ("gas", "density"), ("gas", "temperature"), ("gas", "cell_mass"), weight_field=None, fontsize=12, x_bins=300, y_bins=300)
+ p3.set_unit("cell_mass", 'Msun')
+ p3.set_zlim(("gas", "cell_mass"), 1e3, 1e8)
+ p3.set_cmap(("gas", "cell_mass"), my_cmap2)
+ p3.set_colorbar_label(("gas", "cell_mass"), "Mass ($\mathrm{M}_{\odot}$)")
+ plot3 = p3.plots[("gas", "cell_mass")]
+ else:
+ # Because ParticlePhasePlot doesn't yet work for a log-log PDF for some reason, I will do the following trick.
+ p3 = PhasePlot(sp, (PartType_Gas_to_use, "Density_2"), (PartType_Gas_to_use, "Temperature_2"), (PartType_Gas_to_use, "Mass_2"), weight_field=None, fontsize=12, x_bins=300, y_bins=300)
+ p3.set_zlim((PartType_Gas_to_use, "Mass_2"), 1e3, 1e8)
+ p3.set_cmap((PartType_Gas_to_use, "Mass_2"), my_cmap2)
+ plot3 = p3.plots[(PartType_Gas_to_use, "Mass_2")]
+
+ p3.set_xlim(1e-29, 1e-21)
+ p3.set_ylim(10, 1e7)
+
+ plot3.figure = fig_PDF[time]
+ plot3.axes = grid_PDF[time][code].axes
+ if code == 0: plot3.cax = grid_PDF[time].cbar_axes[0]
+ p3._setup_plots()
+
+ # Add constant pressure and constant entropy lines to guide the eyes
+ if draw_PDF >= 2:
+ dummy_density = np.logspace(-29, -21, num=4)
+ for constant_i in range(1, 100):
+ constant = 1e-100*100**constant_i
+ line = ln.Line2D(dummy_density, constant/dummy_density, linestyle=":", linewidth=0.6, color='k', alpha=0.7) # constant pressure P ~ n*T
+ grid_PDF[time][code].axes.add_line(line)
+ line2 = ln.Line2D(dummy_density, constant*dummy_density**(2./3.), linestyle="-.", linewidth=0.6, color='k', alpha=0.7) # constant entropy S ~ n**(-2/3)*T
+ grid_PDF[time][code].axes.add_line(line2)
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=3, prop=dict(size=10), frameon=True)
+ grid_PDF[time][code].axes.add_artist(at)
+
+ # POSITION-VELOCITY PDF FOR GAS
+ if draw_pos_vel_PDF >= 1:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ sp.set_field_parameter("normal", disk_normal_vector)
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ sp_dense = sp.cut_region(["obj['gas', 'density'].in_units('g/cm**3') > 1.e-25"]) # For AMR codes, consider only the cells that are dense enough
+ pf.field_info[("index", "cylindrical_r")].take_log = False
+ pf.field_info[("gas", "cylindrical_tangential_velocity")].take_log = False
+ p4 = PhasePlot(sp_dense, ("index", "cylindrical_r"), ("gas", "cylindrical_tangential_velocity"), ("gas", "cell_mass"), weight_field=None, fontsize=12, x_bins=300, y_bins=300)
+ p4.set_unit("cylindrical_r", 'kpc')
+ p4.set_unit("cylindrical_tangential_velocity", 'km/s')
+ p4.set_unit("cell_mass", 'Msun')
+ p4.set_zlim(("gas", "cell_mass"), 1e3, 1e8)
+ p4.set_cmap(("gas", "cell_mass"), 'algae')
+ p4.set_colorbar_label(("gas", "cell_mass"), "Mass ($\mathrm{M}_{\odot}$)")
+ plot4 = p4.plots[("gas", "cell_mass")]
+ else:
+ pf.field_info[(PartType_Gas_to_use, "particle_position_cylindrical_radius")].take_log = False
+ pf.field_info[(PartType_Gas_to_use, "particle_velocity_cylindrical_theta")].take_log = False
+ p4 = ParticlePhasePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_cylindrical_theta"), \
+ (PartType_Gas_to_use, MassType_to_use), deposition='ngp', weight_field=None, fontsize=12, x_bins=300, y_bins=300)
+ p4.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p4.set_unit("particle_velocity_cylindrical_theta", 'km/s')
+ p4.set_unit(MassType_to_use, 'Msun')
+ p4.set_zlim((PartType_Gas_to_use, MassType_to_use), 1e3, 1e8)
+ p4.set_cmap((PartType_Gas_to_use, MassType_to_use), 'algae')
+ p4.set_colorbar_label((PartType_Gas_to_use, MassType_to_use), "Mass ($\mathrm{M}_{\odot}$)")
+ plot4 = p4.plots[(PartType_Gas_to_use, MassType_to_use)]
+
+ p4.set_xlabel("Cylindrical Radius (kpc)")
+ p4.set_ylabel("Rotational Velocity (km/s)")
+ p4.set_xlim(0, 14)
+ p4.set_ylim(-100, 350)
+
+ plot4.figure = fig_pos_vel_PDF[time]
+ plot4.axes = grid_pos_vel_PDF[time][code].axes
+ if code == 0: plot4.cax = grid_pos_vel_PDF[time].cbar_axes[0]
+ p4._setup_plots()
+
+ # Add 1D profile line if requested
+ if draw_pos_vel_PDF >= 2 and time != 0:
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p5 = ProfilePlot(sp_dense, ("index", "cylindrical_r"), ("gas", "cylindrical_tangential_velocity"), \
+ weight_field=("gas", "cell_mass"), n_bins=50, x_log=False)
+ p5.set_log("cylindrical_tangential_velocity", False)
+ p5.set_log("cylindrical_r", False)
+ p5.set_unit("cylindrical_r", 'kpc')
+ p5.set_xlim(1e-3, 14)
+ p5.set_ylim("cylindrical_tangential_velocity", -100, 350)
+ line = ln.Line2D(p5.profiles[0].x.in_units('kpc'), p5.profiles[0]["cylindrical_tangential_velocity"].in_units('km/s'), linestyle="-", linewidth=2, color='k', alpha=0.7)
+ pos_vel_xs[time].append(p5.profiles[0].x.in_units('kpc').d)
+ pos_vel_profiles[time].append(p5.profiles[0]["cylindrical_tangential_velocity"].in_units('km/s').d)
+ else:
+ p5 = ProfilePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_cylindrical_theta"), \
+ weight_field=(PartType_Gas_to_use, MassType_to_use), n_bins=50, x_log=False)
+ p5.set_log("particle_velocity_cylindrical_theta", False)
+ p5.set_log("particle_position_cylindrical_radius", False)
+ p5.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p5.set_xlim(1e-3, 14)
+ p5.set_ylim("particle_velocity_cylindrical_theta", -100, 350)
+ line = ln.Line2D(p5.profiles[0].x.in_units('kpc'), p5.profiles[0]["particle_velocity_cylindrical_theta"].in_units('km/s'), linestyle="-", linewidth=2, color='k', alpha=0.7)
+ pos_vel_xs[time].append(p5.profiles[0].x.in_units('kpc').d)
+ pos_vel_profiles[time].append(p5.profiles[0]["particle_velocity_cylindrical_theta"].in_units('km/s').d)
+ grid_pos_vel_PDF[time][code].axes.add_line(line)
+
+ # Add velocity dispersion profile if requested
+ if draw_pos_vel_PDF >= 3 and time != 0:
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ # To calculate velocity dispersion -- residual velocity components other than the rotational velocity found above -- we will need a mass-weighted average of [v_i - v_rot_local]**2
+ # Here we assumed that the disk is on x-y plane; i.e. the variable disk_normal_vector above needs to be [0.0, 0.0, 1.0]
+ def _local_rotational_velocity_x(field, data):
+ trans = np.zeros(data[("gas", "velocity_x")].shape)
+ dr = 0.5*(pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0])
+ for radius, v_rot_local in zip(pos_vel_xs[time][code], pos_vel_profiles[time][code]):
+ ind = np.where((data[("index", "cylindrical_r")].in_units("kpc") >= (radius - dr)) & (data[("index", "cylindrical_r")].in_units("kpc") < (radius + dr)))
+ trans[ind] = -np.sin(data["index", 'cylindrical_theta'][ind]) * v_rot_local * 1e5 # in cm/s
+ return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
+ pf.add_field(("gas", "local_rotational_velocity_x"), function=_local_rotational_velocity_x, take_log=False, particle_type=False, units="cm/s")
+ def _local_rotational_velocity_y(field, data):
+ trans = np.zeros(data[("gas", "velocity_y")].shape)
+ dr = 0.5*(pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0])
+ for radius, v_rot_local in zip(pos_vel_xs[time][code], pos_vel_profiles[time][code]):
+ ind = np.where((data[("index", "cylindrical_r")].in_units("kpc") >= (radius - dr)) & (data[("index", "cylindrical_r")].in_units("kpc") < (radius + dr)))
+ trans[ind] = np.cos(data["index", 'cylindrical_theta'][ind]) * v_rot_local * 1e5
+ return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
+ pf.add_field(("gas", "local_rotational_velocity_y"), function=_local_rotational_velocity_y, take_log=False, particle_type=False, units="cm/s")
+ def _velocity_minus_local_rotational_velocity_squared(field, data):
+ return (data[("gas", "velocity_x")] - data[("gas", "local_rotational_velocity_x")])**2 + \
+ (data[("gas", "velocity_y")] - data[("gas", "local_rotational_velocity_y")])**2 + \
+ (data[("gas", "velocity_z")])**2
+ pf.add_field(("gas", "velocity_minus_local_rotational_velocity_squared"), function=_velocity_minus_local_rotational_velocity_squared,
+ take_log=False, particle_type=False, units="cm**2/s**2")
+# slc = SlicePlot(pf, 'z', ("gas", "local_rotational_velocity_y"), center = center, width = (figure_width, 'kpc'))
+# slc = SlicePlot(pf, 'z', ("gas", "velocity_minus_local_rotational_velocity_squared"), center = center, width = (figure_width, 'kpc')) # this should give zeros everywhere for IC at t=0
+# slc.save()
+
+ p55 = ProfilePlot(sp_dense, ("index", "cylindrical_r"), ("gas", "velocity_minus_local_rotational_velocity_squared"), \
+ weight_field=("gas", "cell_mass"), n_bins=50, x_log=False)
+ p55.set_log("cylindrical_r", False)
+ p55.set_unit("cylindrical_r", 'kpc')
+ p55.set_xlim(1e-3, 14)
+ pos_disp_xs[time].append(p55.profiles[0].x.in_units('kpc').d)
+ pos_disp_profiles[time].append(np.sqrt(p55.profiles[0]["velocity_minus_local_rotational_velocity_squared"]).in_units('km/s').d)
+ else:
+ def _particle_local_rotational_velocity_x(field, data):
+ trans = np.zeros(data[(PartType_Gas_to_use, "particle_velocity_x")].shape)
+ dr = 0.5*(pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0])
+ for radius, v_rot_local in zip(pos_vel_xs[time][code], pos_vel_profiles[time][code]):
+ ind = np.where((data[(PartType_Gas_to_use, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
+ (data[(PartType_Gas_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+ trans[ind] = -np.sin(data[(PartType_Gas_to_use, "particle_position_cylindrical_theta")][ind]) * v_rot_local * 1e5
+ return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
+ pf.add_field((PartType_Gas_to_use, "particle_local_rotational_velocity_x"), function=_particle_local_rotational_velocity_x, take_log=False, particle_type=True, units="cm/s")
+ def _particle_local_rotational_velocity_y(field, data):
+ trans = np.zeros(data[(PartType_Gas_to_use, "particle_velocity_y")].shape)
+ dr = 0.5*(pos_vel_xs[time][code][1] - pos_vel_xs[time][code][0])
+ for radius, v_rot_local in zip(pos_vel_xs[time][code], pos_vel_profiles[time][code]):
+ ind = np.where((data[(PartType_Gas_to_use, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
+ (data[(PartType_Gas_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+ trans[ind] = np.cos(data[(PartType_Gas_to_use, "particle_position_cylindrical_theta")][ind]) * v_rot_local * 1e5
+ return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
+ pf.add_field((PartType_Gas_to_use, "particle_local_rotational_velocity_y"), function=_particle_local_rotational_velocity_y, take_log=False, particle_type=True, units="cm/s")
+ def _particle_velocity_minus_local_rotational_velocity_squared(field, data):
+ return (data[(PartType_Gas_to_use, "particle_velocity_x")] - data[(PartType_Gas_to_use, "particle_local_rotational_velocity_x")])**2 + \
+ (data[(PartType_Gas_to_use, "particle_velocity_y")] - data[(PartType_Gas_to_use, "particle_local_rotational_velocity_y")])**2 + \
+ (data[(PartType_Gas_to_use, "particle_velocity_z")])**2
+ pf.add_field((PartType_Gas_to_use, "particle_velocity_minus_local_rotational_velocity_squared"), function=_particle_velocity_minus_local_rotational_velocity_squared,
+ take_log=False, particle_type=True, units="cm**2/s**2")
+
+ p55 = ProfilePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_minus_local_rotational_velocity_squared"), \
+ weight_field=(PartType_Gas_to_use, MassType_to_use), n_bins=50, x_log=False)
+ p55.set_log("particle_position_cylindrical_radius", False)
+ p55.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p55.set_xlim(1e-3, 14)
+ pos_disp_xs[time].append(p55.profiles[0].x.in_units('kpc').d)
+ pos_disp_profiles[time].append(np.sqrt(p55.profiles[0]["particle_velocity_minus_local_rotational_velocity_squared"]).in_units('km/s').d)
+
+ # Add vertical velocity dispersion profile if requested
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ def _velocity_z_squared(field, data):
+ return (data[("gas", "velocity_z")])**2
+ pf.add_field(("gas", "velocity_z_squared"), function=_velocity_z_squared, take_log=False, particle_type=False, units="cm**2/s**2")
+
+ p551 = ProfilePlot(sp_dense, ("index", "cylindrical_r"), ("gas", "velocity_z_squared"), weight_field=("gas", "cell_mass"), n_bins=50, x_log=False)
+ p551.set_log("cylindrical_r", False)
+ p551.set_unit("cylindrical_r", 'kpc')
+ p551.set_xlim(1e-3, 14)
+ pos_disp_vert_xs[time].append(p551.profiles[0].x.in_units('kpc').d)
+ pos_disp_vert_profiles[time].append(np.sqrt(p551.profiles[0]["velocity_z_squared"]).in_units('km/s').d)
+ else:
+ def _particle_velocity_z_squared(field, data):
+ return (data[(PartType_Gas_to_use, "particle_velocity_z")])**2
+ pf.add_field((PartType_Gas_to_use, "particle_velocity_z_squared"), function=_particle_velocity_z_squared, take_log=False, particle_type=True, units="cm**2/s**2")
+
+ p551 = ProfilePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_velocity_z_squared"), \
+ weight_field=(PartType_Gas_to_use, MassType_to_use), n_bins=50, x_log=False)
+ p551.set_log("particle_position_cylindrical_radius", False)
+ p551.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p551.set_xlim(1e-3, 14)
+ pos_disp_vert_xs[time].append(p551.profiles[0].x.in_units('kpc').d)
+ pos_disp_vert_profiles[time].append(np.sqrt(p551.profiles[0]["particle_velocity_z_squared"]).in_units('km/s').d)
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=4, prop=dict(size=10), frameon=True)
+ grid_pos_vel_PDF[time][code].axes.add_artist(at)
+
+ # POSITION-VELOCITY PDF FOR NEW STARS
+ if draw_star_pos_vel_PDF >= 1 and time != 0:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ sp.set_field_parameter("normal", disk_normal_vector)
+ pf.field_info[(PartType_Star_to_use, "particle_position_cylindrical_radius")].take_log = False
+ pf.field_info[(PartType_Star_to_use, "particle_velocity_cylindrical_theta")].take_log = False
+ pf.field_info[(PartType_Star_to_use, "particle_mass")].output_units = 'code_mass' # this turned out to be crucial!; otherwise wrong output_unit 'g' is assumed in ParticlePhasePlot->create_profile in visualizaiton/particle_plots.py for ART-I/ENZO/RAMSES
+ p41 = ParticlePhasePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_cylindrical_theta"), \
+ (PartType_Star_to_use, "particle_mass"), deposition='ngp', weight_field=None, fontsize=12, x_bins=300, y_bins=300)
+ p41.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p41.set_unit("particle_velocity_cylindrical_theta", 'km/s')
+ p41.set_unit("particle_mass", 'Msun') # requires a change in set_unit in visualization/profile_plotter.py: remove self.plots[field].zmin, self.plots[field].zmax = (None, None)
+# p41.set_unit((PartType_Star_to_use, "particle_mass"), 'Msun') # Neither this nor above works without such change
+ p41.set_zlim((PartType_Star_to_use, "particle_mass"), 1e3, 1e7)
+ p41.set_cmap((PartType_Star_to_use, "particle_mass"), 'algae')
+
+ p41.set_colorbar_label((PartType_Star_to_use, "particle_mass"), "Newly Formed Stellar Mass ($\mathrm{M}_{\odot}$)")
+ plot41 = p41.plots[(PartType_Star_to_use, "particle_mass")]
+
+ p41.set_xlabel("Cylindrical Radius (kpc)")
+ p41.set_ylabel("Rotational Velocity (km/s)")
+ p41.set_xlim(0, 14)
+ p41.set_ylim(-100, 350)
+
+ plot41.figure = fig_star_pos_vel_PDF[time]
+ plot41.axes = grid_star_pos_vel_PDF[time][code].axes
+ if code == 0: plot41.cax = grid_star_pos_vel_PDF[time].cbar_axes[0]
+ p41._setup_plots()
+
+ # Add 1D profile line if requested
+ if draw_star_pos_vel_PDF >= 2 and time != 0:
+ p51 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_cylindrical_theta"), \
+ weight_field=(PartType_Star_to_use, "particle_mass"), n_bins=50, x_log=False)
+ p51.set_log((PartType_Star_to_use, "particle_velocity_cylindrical_theta"), False)
+ p51.set_log((PartType_Star_to_use, "particle_position_cylindrical_radius"), False)
+ p51.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p51.set_xlim(1e-3, 14)
+ p51.set_ylim("particle_velocity_cylindrical_theta", -100, 350)
+ line = ln.Line2D(p51.profiles[0].x.in_units('kpc'), p51.profiles[0]["particle_velocity_cylindrical_theta"].in_units('km/s'), linestyle="-", linewidth=2, color='k', alpha=0.7)
+ star_pos_vel_xs[time].append(p51.profiles[0].x.in_units('kpc').d)
+ star_pos_vel_profiles[time].append(p51.profiles[0]["particle_velocity_cylindrical_theta"].in_units('km/s').d)
+ grid_star_pos_vel_PDF[time][code].axes.add_line(line)
+
+ # Add velocity dispersion profile if requested
+ if draw_star_pos_vel_PDF >= 3 and time != 0:
+ def _particle_local_rotational_velocity_x(field, data):
+ trans = np.zeros(data[(PartType_StarBeforeFiltered_to_use, "particle_velocity_x")].shape)
+ dr = 0.5*(star_pos_vel_xs[time][code][1] - star_pos_vel_xs[time][code][0])
+ for radius, v_rot_local in zip(star_pos_vel_xs[time][code], star_pos_vel_profiles[time][code]):
+ ind = np.where((data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
+ (data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+ trans[ind] = -np.sin(data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_theta")][ind]) * v_rot_local * 1e5
+ return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_local_rotational_velocity_x"), function=_particle_local_rotational_velocity_x, take_log=False, particle_type=True, units="cm/s")
+ def _particle_local_rotational_velocity_y(field, data):
+ trans = np.zeros(data[(PartType_StarBeforeFiltered_to_use, "particle_velocity_y")].shape)
+ dr = 0.5*(star_pos_vel_xs[time][code][1] - star_pos_vel_xs[time][code][0])
+ for radius, v_rot_local in zip(star_pos_vel_xs[time][code], star_pos_vel_profiles[time][code]):
+ ind = np.where((data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_radius")].in_units("kpc") >= (radius - dr)) & \
+ (data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_radius")].in_units("kpc") < (radius + dr)))
+ trans[ind] = np.cos(data[(PartType_StarBeforeFiltered_to_use, "particle_position_cylindrical_theta")][ind]) * v_rot_local * 1e5
+ return data.ds.arr(trans, "cm/s").in_base(data.ds.unit_system.name)
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_local_rotational_velocity_y"), function=_particle_local_rotational_velocity_y, take_log=False, particle_type=True, units="cm/s")
+ def _particle_velocity_minus_local_rotational_velocity_squared(field, data):
+ return (data[(PartType_StarBeforeFiltered_to_use, "particle_velocity_x")] - data[(PartType_StarBeforeFiltered_to_use, "particle_local_rotational_velocity_x")])**2 + \
+ (data[(PartType_StarBeforeFiltered_to_use, "particle_velocity_y")] - data[(PartType_StarBeforeFiltered_to_use, "particle_local_rotational_velocity_y")])**2 + \
+ (data[(PartType_StarBeforeFiltered_to_use, "particle_velocity_z")])**2
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_velocity_minus_local_rotational_velocity_squared"), function=_particle_velocity_minus_local_rotational_velocity_squared,
+ take_log=False, particle_type=True, units="cm**2/s**2")
+ pf.add_particle_filter(PartType_Star_to_use) # This is needed for a filtered particle type PartType_Star_to_use to work, because we have just created new particle fields.
+
+ p515 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_minus_local_rotational_velocity_squared"), \
+ weight_field=(PartType_Star_to_use, "particle_mass"), n_bins=50, x_log=False)
+ p515.set_log((PartType_Star_to_use, "particle_position_cylindrical_radius"), False)
+ p515.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p515.set_xlim(1e-3, 14)
+ star_pos_disp_xs[time].append(p515.profiles[0].x.in_units('kpc').d)
+ star_pos_disp_profiles[time].append(np.sqrt(p515.profiles[0]["particle_velocity_minus_local_rotational_velocity_squared"]).in_units('km/s').d)
+
+ # Add vertical velocity dispersion profile if requested
+ if draw_star_pos_vel_PDF >= 4 and time != 0:
+ def _particle_velocity_z_squared(field, data):
+ return (data[(PartType_StarBeforeFiltered_to_use, "particle_velocity_z")])**2
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_velocity_z_squared"), function=_particle_velocity_z_squared,
+ take_log=False, particle_type=True, units="cm**2/s**2")
+ pf.add_particle_filter(PartType_Star_to_use) # This is needed for a filtered particle type PartType_Star_to_use to work, because we have just created new particle fields.
+
+ p5151 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_velocity_z_squared"), \
+ weight_field=(PartType_Star_to_use, "particle_mass"), n_bins=50, x_log=False)
+ p5151.set_log((PartType_Star_to_use, "particle_position_cylindrical_radius"), False)
+ p5151.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p5151.set_xlim(1e-3, 14)
+ star_pos_disp_vert_xs[time].append(p5151.profiles[0].x.in_units('kpc').d)
+ star_pos_disp_vert_profiles[time].append(np.sqrt(p5151.profiles[0]["particle_velocity_z_squared"]).in_units('km/s').d)
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=4, prop=dict(size=10), frameon=True)
+ grid_star_pos_vel_PDF[time][code].axes.add_artist(at)
+
+ # RADIUS-HEIGHT PDF
+ if draw_rad_height_PDF >= 1:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ if draw_rad_height_PDF == 1 or draw_rad_height_PDF == 2:
+ p55 = PhasePlot(sp, ("index", "cylindrical_r"), ("index", "cylindrical_z_abs"), ("gas", "density"), weight_field=("gas", "cell_mass"), fontsize=12, x_bins=200, y_bins=200)
+ p55.set_zlim(("gas", "density"), 1e-26, 1e-21)
+ p55.set_cmap(("gas", "density"), 'algae')
+ p55.set_log("cylindrical_r", False)
+ p55.set_log("cylindrical_z_abs", False)
+ p55.set_unit("cylindrical_r", 'kpc')
+ p55.set_unit("cylindrical_z_abs", 'kpc')
+ plot55 = p55.plots[("gas", "density")]
+ elif draw_rad_height_PDF == 3:
+ p55 = PhasePlot(sp, ("index", "cylindrical_r"), ("index", "cylindrical_z_abs"), ("gas", "density_minus_analytic"), weight_field=("gas", "cell_mass"), fontsize=12, x_bins=200, y_bins=200)
+ p55.set_zlim(("gas", "density_minus_analytic"), -1e-24, 1e-24)
+ p55.set_cmap(("gas", "density_minus_analytic"), 'algae')
+ p55.set_log("cylindrical_r", False)
+ p55.set_log("cylindrical_z_abs", False)
+ p55.set_unit("cylindrical_r", 'kpc')
+ p55.set_unit("cylindrical_z_abs", 'kpc')
+ plot55 = p55.plots[("gas", "density_minus_analytic")]
+ else:
+ # Because ParticlePhasePlot doesn't work for these newly created fields for some reason, I will do the following trick.
+ if draw_rad_height_PDF == 1 or draw_rad_height_PDF == 2:
+ p55 = PhasePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), (PartType_Gas_to_use, "Density_2"), weight_field=(PartType_Gas_to_use, "Mass_2"), fontsize=12, x_bins=200, y_bins=200)
+ p55.set_zlim((PartType_Gas_to_use, "Density_2"), 1e-26, 1e-21)
+ p55.set_cmap((PartType_Gas_to_use, "Density_2"), 'algae')
+ p55.set_log("particle_position_cylindrical_radius", False)
+ p55.set_log("particle_position_cylindrical_z_abs", False)
+ p55.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p55.set_unit("particle_position_cylindrical_z_abs", 'kpc')
+ plot55 = p55.plots[(PartType_Gas_to_use, "Density_2")]
+ elif draw_rad_height_PDF == 3:
+ p55 = PhasePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), (PartType_Gas_to_use, "Density_2_minus_analytic"), weight_field=(PartType_Gas_to_use, "Mass_2"), fontsize=12, x_bins=200, y_bins=200)
+ p55.set_zlim((PartType_Gas_to_use, "Density_2_minus_analytic"), -1e-24, 1e-24)
+ p55.set_cmap((PartType_Gas_to_use, "Density_2_minus_analytic"), 'algae')
+ p55.set_log("particle_position_cylindrical_radius", False)
+ p55.set_log("particle_position_cylindrical_z_abs", False)
+ p55.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p55.set_unit("particle_position_cylindrical_z_abs", 'kpc')
+ plot55 = p55.plots[(PartType_Gas_to_use, "Density_2_minus_analytic")]
+
+ p55.set_xlabel("Cylindrical Radius (kpc)")
+ p55.set_ylabel("Vertical Height (kpc)")
+ p55.set_xlim(0, 14)
+ p55.set_ylim(0, 1.4)
+
+ plot55.figure = fig_rad_height_PDF[time]
+ plot55.axes = grid_rad_height_PDF[time][code].axes
+ if code == 0: plot55.cax = grid_rad_height_PDF[time].cbar_axes[0]
+ p55._setup_plots()
+
+ # Add 1D profile line if requested
+ if draw_rad_height_PDF == 2:
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p56 = ProfilePlot(sp, ("index", "cylindrical_r"), ("index", "cylindrical_z_abs"), \
+ weight_field=("gas", "cell_mass"), n_bins=50, x_log=False)
+ p56.set_log("cylindrical_r", False)
+ p56.set_log("cylindrical_z_abs", False)
+ p56.set_unit("cylindrical_r", 'kpc')
+ p56.set_unit("cylindrical_z_abs", 'kpc')
+ p56.set_xlim(0, 14)
+ p56.set_ylim("cylindrical_z_abs", 0, 1.4)
+ line = ln.Line2D(p56.profiles[0].x.in_units('kpc'), p56.profiles[0]["cylindrical_z_abs"].in_units('kpc'), linestyle="-", linewidth=2, color='k', alpha=0.7)
+ rad_height_xs[time].append(p56.profiles[0].x.in_units('kpc').d)
+ rad_height_profiles[time].append(p56.profiles[0]["cylindrical_z_abs"].in_units('kpc').d)
+ else:
+ p56 = ProfilePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), \
+ weight_field=(PartType_Gas_to_use, MassType_to_use), n_bins=50, x_log=False)
+ p56.set_log("particle_position_cylindrical_radius", False)
+ p56.set_log("particle_position_cylindrical_z_abs", False)
+ p56.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p56.set_unit("particle_position_cylindrical_z_abs", 'kpc')
+ p56.set_xlim(0, 14)
+ p56.set_ylim("particle_position_cylindrical_z_abs", 0, 1.4)
+ line = ln.Line2D(p56.profiles[0].x.in_units('kpc'), p56.profiles[0]["particle_position_cylindrical_z_abs"].in_units('kpc'), linestyle="-", linewidth=2, color='k', alpha=0.7)
+ rad_height_xs[time].append(p56.profiles[0].x.in_units('kpc').d)
+ rad_height_profiles[time].append(p56.profiles[0]["particle_position_cylindrical_z_abs"].in_units('kpc').d)
+ grid_rad_height_PDF[time][code].axes.add_line(line)
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=1, prop=dict(size=10), frameon=True)
+ grid_rad_height_PDF[time][code].axes.add_artist(at)
+
+ # DENSITY-TEMPERATURE-METALLICITY PDF
+ if draw_metal_PDF == 1:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ my_cmap2 = copy.copy(matplotlib.cm.get_cmap('algae'))
+ my_cmap2.set_under('w')
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ pf.field_info[("gas", "metallicity")].take_log = False
+ p3 = PhasePlot(sp, ("gas", "density"), ("gas", "temperature"), ("gas", "metallicity"), weight_field=("gas", "cell_mass"), fontsize=12, x_bins=300, y_bins=300)
+ p3.set_zlim(("gas", "metallicity"), 0.01, 0.04)
+ p3.set_cmap(("gas", "metallicity"), my_cmap2)
+ p3.set_colorbar_label(("gas", "metallicity"), "Metallicity (Mass-weighted average of mass fraction)")
+ plot3 = p3.plots[("gas", "metallicity")]
+ else:
+ # Because ParticlePhasePlot doesn't yet work for a log-log PDF for some reason, I will do the following trick.
+ pf.field_info[(PartType_Gas_to_use, "Metallicity_2")].take_log = False
+ p3 = PhasePlot(sp, (PartType_Gas_to_use, "Density_2"), (PartType_Gas_to_use, "Temperature_2"), (PartType_Gas_to_use, "Metallicity_2"), weight_field=(PartType_Gas_to_use, "Mass_2"), fontsize=12, x_bins=300, y_bins=300)
+ p3.set_zlim((PartType_Gas_to_use, "Metallicity_2"), 0.01, 0.04)
+ p3.set_cmap((PartType_Gas_to_use, "Metallicity_2"), my_cmap2)
+ plot3 = p3.plots[(PartType_Gas_to_use, "Metallicity_2")]
+
+ p3.set_xlim(1e-29, 1e-21)
+ p3.set_ylim(10, 1e7)
+
+ plot3.figure = fig_metal_PDF[time]
+ plot3.axes = grid_metal_PDF[time][code].axes
+ if code == 0: plot3.cax = grid_metal_PDF[time].cbar_axes[0]
+ p3._setup_plots()
+
+ if add_nametag == 1:
+ at = AnchoredText("%s" % codes[code], loc=3, prop=dict(size=10), frameon=True)
+ grid_metal_PDF[time][code].axes.add_artist(at)
+
+ # DENSITY DF (DISTRIBUTION FUNCTION)
+ if draw_density_DF >= 1:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p6 = ProfilePlot(sp, ("gas", "density"), ("gas", "cell_mass"), weight_field=None, n_bins=50, x_log=True, accumulation=False)
+# p6 = ProfilePlot(sp, ("gas", "density"), ("gas", "cell_mass"), weight_field=None, n_bins=50, x_log=True, accumulation=True)
+ p6.set_log("cell_mass", True)
+ p6.set_xlim(1e-29, 1e-21)
+ density_DF_xs[time].append(p6.profiles[0].x.in_units('g/cm**3').d)
+ density_DF_profiles[time].append(p6.profiles[0]["cell_mass"].in_units('Msun').d)
+ else:
+ # Because ParticleProfilePlot doesn't exist, I will do the following trick.
+ p6 = ProfilePlot(sp, (PartType_Gas_to_use, "Density_2"), (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=True, accumulation=False)
+# p6 = ProfilePlot(sp, (PartType_Gas_to_use, "Density_2"), (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=True, accumulation=True)
+ p6.set_log("Mass_2", True)
+ p6.set_xlim(1e-29, 1e-21)
+ density_DF_xs[time].append(p6.profiles[0].x.in_units('g/cm**3').d)
+ density_DF_profiles[time].append(p6.profiles[0]["Mass_2"].in_units('Msun').d)
+
+ # Add difference plot between 1st and 2nd datasets, if requested
+ if draw_density_DF == 2 and time != 0:
+ if dataset_num == 2:
+ pf_1st = load_dataset(1, time, code, codes, filenames[0]) # load 1st datasets
+ v, cen = pf_1st.h.find_max(("gas", "density"))
+ sp = pf_1st.sphere(cen, (30.0, "kpc"))
+ cen2 = sp.quantities.center_of_mass(use_gas=True, use_particles=False).in_units("kpc")
+ sp2 = pf_1st.sphere(cen2, (1.0, "kpc"))
+ cen3 = sp2.quantities.max_location(("gas", "density"))
+ center_1st = pf_1st.arr([cen3[1].d, cen3[2].d, cen3[3].d], 'code_length')
+ if yt_version_pre_3_2_3 == 1:
+ center_1st = pf_1st.arr([cen3[2].d, cen3[3].d, cen3[4].d], 'code_length') # for yt-3.2.3 or before
+ sp_1st = pf_1st.sphere(center_1st, (0.5*figure_width, "kpc"))
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p61 = ProfilePlot(sp_1st, ("gas", "density"), ("gas", "cell_mass"), weight_field=None, n_bins=50, x_log=True, accumulation=False)
+ p61.set_log("cell_mass", True)
+ p61.set_xlim(1e-29, 1e-21)
+ density_DF_1st_xs[time].append(p61.profiles[0].x.in_units('g/cm**3').d)
+ density_DF_1st_profiles[time].append(p61.profiles[0]["cell_mass"].in_units('Msun').d)
+ else:
+ def _Density_2(field, data):
+ return data[(PartType_Gas_to_use, "Density")].in_units('g/cm**3')
+ pf_1st.add_field((PartType_Gas_to_use, "Density_2"), function=_Density_2, take_log=True, particle_type=False, display_name="Density", units="g/cm**3")
+ def _Mass_2(field, data):
+ return data[(PartType_Gas_to_use, MassType_to_use)].in_units('Msun')
+ pf_1st.add_field((PartType_Gas_to_use, "Mass_2"), function=_Mass_2, take_log=True, particle_type=False, display_name="Mass", units="Msun")
+ p61 = ProfilePlot(sp_1st, (PartType_Gas_to_use, "Density_2"), (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=True, accumulation=False)
+ p61.set_log("Mass_2", True)
+ p61.set_xlim(1e-29, 1e-21)
+ density_DF_1st_xs[time].append(p61.profiles[0].x.in_units('g/cm**3').d)
+ density_DF_1st_profiles[time].append(p61.profiles[0]["Mass_2"].in_units('Msun').d)
+ else:
+ print("This won't work; consider setting dataset_num to 2...")
+ continue
+
+ # CYLINDRICAL RADIUS DF + RADIALLY-BINNED GAS SURFACE DENSITY
+ if draw_radius_DF == 1:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ sp.set_field_parameter("normal", disk_normal_vector)
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p7 = ProfilePlot(sp, ("index", "cylindrical_r"), ("gas", "cell_mass"), weight_field=None, n_bins=50, x_log=False, accumulation=False)
+ p7.set_log("cell_mass", True)
+ p7.set_log("cylindrical_r", False)
+ p7.set_unit("cylindrical_r", 'kpc')
+ p7.set_xlim(1e-3, 15)
+ radius_DF_xs[time].append(p7.profiles[0].x.in_units('kpc').d)
+ radius_DF_profiles[time].append(p7.profiles[0]["cell_mass"].in_units('Msun').d)
+ else:
+ p7 = ProfilePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_radius"), (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=50, x_log=False, accumulation=False)
+ p7.set_log("Mass_2", True)
+ p7.set_log("particle_position_cylindrical_radius", False)
+ p7.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p7.set_xlim(1e-3, 15)
+ radius_DF_xs[time].append(p7.profiles[0].x.in_units('kpc').d)
+ radius_DF_profiles[time].append(p7.profiles[0]["Mass_2"].in_units('Msun').d)
+
+ # CYLINDRICAL RADIUS DF + RADIALLY-BINNED SURFACE DENSITY FOR NEW STARS
+ if draw_star_radius_DF >= 1 and time != 0:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ sp.set_field_parameter("normal", disk_normal_vector)
+ pf.field_info[(PartType_Star_to_use, "particle_mass")].take_log = True
+ pf.field_info[(PartType_Star_to_use, "particle_mass")].output_units = 'code_mass' # this turned out to be crucial!; check output_units above
+ pf.field_info[(PartType_Star_to_use, "particle_position_cylindrical_radius")].take_log = False
+ p71 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_mass"), weight_field=None, n_bins=50, x_log=False, accumulation=False)
+ p71.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p71.set_xlim(1e-3, 15)
+ star_radius_DF_xs[time].append(p71.profiles[0].x.in_units('kpc').d)
+ star_radius_DF_profiles[time].append(p71.profiles[0]["particle_mass"].in_units('Msun').d)
+
+ # Add RADIALLY-BINNED SFR SURFACE DENSITY PROFILE if requested (SFR estimated using stars younger than 20 Myrs old)
+ if draw_star_radius_DF == 2 and time != 0:
+ def _particle_mass_young_stars_star_radius_DF(field, data):
+ trans = np.zeros(data[(PartType_StarBeforeFiltered_to_use, "particle_mass")].shape)
+ ind = np.where(data[(PartType_StarBeforeFiltered_to_use, FormationTimeType_to_use)].in_units('Myr') > (pf.current_time.in_units('Myr').d - young_star_cutoff_star_radius_DF))
+ trans[ind] = data[(PartType_StarBeforeFiltered_to_use, "particle_mass")][ind].in_units('code_mass')
+ return data.ds.arr(trans, "code_mass").in_base(data.ds.unit_system.name)
+ pf.add_field((PartType_StarBeforeFiltered_to_use, "particle_mass_young_stars_star_radius_DF"), function=_particle_mass_young_stars_star_radius_DF, \
+ display_name="Young Stellar Mass", units='code_mass', particle_type=True, take_log=True)
+ pf.add_particle_filter(PartType_Star_to_use)
+
+ pf.field_info[(PartType_Star_to_use, "particle_mass_young_stars_star_radius_DF")].take_log = True
+ pf.field_info[(PartType_Star_to_use, "particle_mass_young_stars_star_radius_DF")].output_units = 'code_mass' # this turned out to be crucial!; check output_units above
+ pf.field_info[(PartType_Star_to_use, "particle_position_cylindrical_radius")].take_log = False
+ p72 = ProfilePlot(sp, (PartType_Star_to_use, "particle_position_cylindrical_radius"), (PartType_Star_to_use, "particle_mass_young_stars_star_radius_DF"), \
+ weight_field=None, n_bins=50, x_log=False, accumulation=False)
+ p72.set_unit("particle_position_cylindrical_radius", 'kpc')
+ p72.set_xlim(1e-3, 15)
+ sfr_radius_DF_xs[time].append(p72.profiles[0].x.in_units('kpc').d)
+ sfr_radius_DF_profiles[time].append(p72.profiles[0]["particle_mass_young_stars_star_radius_DF"].in_units('Msun').d/young_star_cutoff_star_radius_DF/1e6) # in Msun/yr
+
+ # VERTICAL HEIGHT DF + VERTICALLY-BINNED GAS SURFACE DENSITY
+ if draw_height_DF == 1:
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ sp.set_field_parameter("normal", disk_normal_vector)
+ if codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES":
+ p8 = ProfilePlot(sp, ("index", "cylindrical_z_abs"), ("gas", "cell_mass"), weight_field=None, n_bins=10, x_log=False, accumulation=False)
+ p8.set_log("cell_mass", True)
+ p8.set_log("cylindrical_z_abs", False)
+ p8.set_unit("cylindrical_z_abs", 'kpc')
+ p8.set_xlim(1e-3, 1.4)
+ height_DF_xs[time].append(p8.profiles[0].x.in_units('kpc').d)
+ height_DF_profiles[time].append(p8.profiles[0]["cell_mass"].in_units('Msun').d)
+ else:
+ p8 = ProfilePlot(sp, (PartType_Gas_to_use, "particle_position_cylindrical_z_abs"), (PartType_Gas_to_use, "Mass_2"), weight_field=None, n_bins=10, x_log=False, accumulation=False)
+ p8.set_log("Mass_2", True)
+ p8.set_log("particle_position_cylindrical_z_abs", False)
+ p8.set_unit("particle_position_cylindrical_z_abs", 'kpc')
+ p8.set_xlim(1e-3, 1.4)
+ height_DF_xs[time].append(p8.profiles[0].x.in_units('kpc').d)
+ height_DF_profiles[time].append(p8.profiles[0]["Mass_2"].in_units('Msun').d)
+
+ # STAR FORMATION RATE + CUMULATIVE STELLAR MASS GROWTH IN TIME
+ if draw_SFR >= 1 and time != 0:
+ from yt.units.dimensions import length # Below are tricks to make StarFormationRate() work, particularly with "volume" argument, as it currently works only with comoving datasets
+ pf.unit_registry.add('pccm', pf.unit_registry.lut['pc'][0], length, "\\rm{pc}/(1+z)")
+ pf.hubble_constant = 0.71; pf.omega_lambda = 0.73; pf.omega_matter = 0.27; pf.omega_curvature = 0.0
+
+ sp = pf.sphere(center, (0.5*figure_width, "kpc"))
+ draw_SFR_mass = sp[(PartType_Star_to_use, "particle_mass")].in_units('Msun')
+ draw_SFR_ct = sp[(PartType_Star_to_use, FormationTimeType_to_use)].in_units('Myr')
+ sfr = StarFormationRate(pf, star_mass = draw_SFR_mass, star_creation_time = draw_SFR_ct,
+ volume = sp.volume(), bins = 25) # 25 bins hardcoded; see: http://yt-project.org/docs/dev/analyzing/analysis_modules/star_analysis.html
+
+ sfr_ts[time].append(sfr.time.in_units('Myr')) # in Myr
+ sfr_cum_masses[time].append(sfr.Msol_cumulative) # in Msun
+ sfr_sfrs[time].append(sfr.Msol_yr) # in Msun/yr
+
+ # DENSITY ALONG THE ORTHO-RAY OBJECT CUTTING THROUGH THE CENTER
+ if draw_cut_through == 1:
+ ray = pf.ortho_ray(2, (center[0].in_units('code_length'), center[1].in_units('code_length'))) # see: http://yt-project.org/doc/visualizing/manual_plotting.html#line-plots
+ srt = np.argsort(ray['z'])
+ cut_through_zs[time].append(np.array(ray['z'][srt].in_units('kpc').d - center[2].in_units('kpc').d))
+ cut_through_zvalues[time].append(np.array(ray[("gas", "density")][srt].in_units('g/cm**3').d))
+
+ ray = pf.ortho_ray(0, (center[1].in_units('code_length'), center[2].in_units('code_length')))
+ srt = np.argsort(ray['x'])
+ cut_through_xs[time].append(np.array(ray['x'][srt].in_units('kpc').d - center[2].in_units('kpc').d))
+ cut_through_xvalues[time].append(np.array(ray[("gas", "density")][srt].in_units('g/cm**3').d))
+
+ ####################################
+ # POST-ANALYSIS STEPS #
+ ####################################
+
+ # SAVE FIGURES
+ if draw_density_map == 1:
+ fig_density_map[time].savefig("Sigma_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_temperature_map == 1:
+ fig_temperature_map[time].savefig("Temp_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_cellsize_map == 1 or draw_cellsize_map == 3:
+ fig_cellsize_map[time].savefig("Cell_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_cellsize_map == 2 or draw_cellsize_map == 3:
+ fig_cellsize_map_2[time].savefig("Resolution_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_elevation_map == 1:
+ fig_elevation_map[time].savefig("Elevation_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_metal_map >= 1:
+ fig_metal_map[time].savefig("Metal_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_zvel_map == 1:
+ fig_zvel_map[time].savefig("z-vel_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_star_map == 1 and time != 0:
+ fig_star_map[time].savefig("Star_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_star_clump_stats >= 1 and time != 0:
+ if draw_star_clump_stats >= 2:
+ fig_star_map_2[time].savefig("Star_with_clumps_HOP_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ fig_star_map_3[time].savefig("Star_with_clumps_FOF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_star_clump_stats == 3:
+ pf_clump_ref = load_dataset(2, 1, 0, ['GIZMO'],
+ [['dummy', '/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository/Grackle+SF+ThermalFbck/GIZMO/AGORA_disk_second_ps1.5/snapshot_temp_100_last']])
+ PartType_Star_to_use = "Stars"
+ if os.path.exists("./halo_catalogs/hop_%s_%05d_high_floor/hop_%s_%05d_high_floor.0.h5" % ('GIZMO', 500, 'GIZMO', 500)) == False:
+ hc_clump_ref = HaloCatalog(data_ds=pf_clump_ref, finder_method='hop', output_dir="./halo_catalogs/hop_%s_%05d_high_floor" % ('GIZMO', 500), \
+ finder_kwargs={'threshold': 2e8, 'dm_only': False, 'ptype': PartType_Star_to_use})
+ hc_clump_ref.add_filter('quantity_value', 'particle_mass', '>', 2.6e6, 'Msun')
+ hc_clump_ref.add_filter('quantity_value', 'particle_mass', '<', 2.6e8, 'Msun')
+ hc_clump_ref.create()
+ halo_ds_clump_ref = load("./halo_catalogs/hop_%s_%05d_high_floor/hop_%s_%05d_high_floor.0.h5" % ('GIZMO', 500, 'GIZMO', 500))
+ hc_clump_ref = HaloCatalog(halos_ds=halo_ds_clump_ref, output_dir="./halo_catalogs/hop_%s_%05d_high_floor" % ('GIZMO', 500))
+ hc_clump_ref.load()
+
+ halo_ad_clump_ref = hc_clump_ref.halos_ds.all_data()
+ star_clump_masses_hop_ref.append(np.log10(halo_ad_clump_ref['particle_mass'][:].in_units("Msun")))
+
+ if os.path.exists("./halo_catalogs/fof_%s_%05d_high_floor/fof_%s_%05d_high_floor.0.h5" % ('GIZMO', 500, 'GIZMO', 500)) == False:
+ hc2_clump_ref = HaloCatalog(data_ds=pf_clump_ref, finder_method='fof', output_dir="./halo_catalogs/fof_%s_%05d_high_floor" % ('GIZMO', 500), \
+ finder_kwargs={'link': 0.0025, 'dm_only': False, 'ptype': PartType_Star_to_use})
+ hc2_clump_ref.add_filter('quantity_value', 'particle_mass', '>', 2.6e6, 'Msun')
+ hc2_clump_ref.add_filter('quantity_value', 'particle_mass', '<', 2.6e8, 'Msun')
+ hc2_clump_ref.create()
+ halo_ds2_clump_ref = load("./halo_catalogs/fof_%s_%05d_high_floor/fof_%s_%05d_high_floor.0.h5" % ('GIZMO', 500, 'GIZMO', 500))
+ hc2_clump_ref = HaloCatalog(halos_ds=halo_ds2_clump_ref, output_dir="./halo_catalogs/fof_%s_%05d_high_floor" % ('GIZMO', 500))
+ hc2_clump_ref.load()
+
+ halo_ad2_clump_ref = hc2_clump_ref.halos_ds.all_data()
+ star_clump_masses_fof_ref.append(np.log10(halo_ad2_clump_ref['particle_mass'][:].in_units("Msun")))
+ plt.clf()
+ fig = plt.figure(figsize=(8, 4))
+ gridspec.GridSpec(1, 2)
+ for star_clump_stats_i in range(1,3,1):
+ codes_plotted = []
+# plt.subplot(1,2,star_clump_stats_i, aspect=0.25)
+ plt.subplot2grid((1, 2), (0, star_clump_stats_i-1))
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ if (star_clump_stats_i == 1 and (codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES")) or \
+ (star_clump_stats_i == 2 and (codes[code] == "CHANGA" or codes[code] == "GASOLINE" or codes[code] == "GADGET-3" or codes[code] == "GEAR" or codes[code] == "GIZMO" or codes[code] == "SWIFT")):
+ hist = np.histogram(star_clump_masses_hop[time][code], bins=10, range=(6., 8.5))
+ dbin = 0.5*(hist[1][1] - hist[1][0])
+ lines = plt.plot(hist[1][:-1]+dbin, hist[0], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ codes_plotted.append(codes[code])
+ plt.xlim(6, 8.5)
+ plt.ylim(-0.1, 15)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$")
+ if star_clump_stats_i == 1:
+ plt.ylabel("$\mathrm{Stellar\ Clump\ Counts, \ \ N_{clump}(M)}$")
+ plt.legend(codes_plotted, loc=1, frameon=True, ncol=1, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='xx-small')
+ plt.savefig("star_clump_stats_HOP_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ # Reiterate for cumulative plots
+ fig = plt.figure(figsize=(8, 4))
+ gridspec.GridSpec(1, 2)
+ for star_clump_stats_i in range(1,3,1):
+ codes_plotted = []
+# plt.subplot(1,2,star_clump_stats_i, aspect=0.15)
+ plt.subplot2grid((1, 2), (0, star_clump_stats_i-1))
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ if (star_clump_stats_i == 1 and (codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES")) or \
+ (star_clump_stats_i == 2 and (codes[code] == "CHANGA" or codes[code] == "GASOLINE" or codes[code] == "GADGET-3" or codes[code] == "GEAR" or codes[code] == "GIZMO" or codes[code] == "SWIFT")):
+ hist = np.histogram(star_clump_masses_hop[time][code], bins=10, range=(6., 8.5))
+ dbin = 0.5*(hist[1][1] - hist[1][0])
+ lines = plt.plot(hist[1][:-1]+dbin, np.cumsum(hist[0][::-1])[::-1], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], \
+ marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ codes_plotted.append(codes[code])
+ if draw_star_clump_stats == 3 and star_clump_stats_i == 2:
+ hist_clump_ref = np.histogram(star_clump_masses_hop_ref, bins=10, range=(6., 8.5))
+ dbin = 0.5*(hist[1][1] - hist[1][0])
+ lines = plt.plot(hist_clump_ref[1][:-1]+dbin, np.cumsum(hist_clump_ref[0][::-1])[::-1], color='k', linestyle='--', marker='*', markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ codes_plotted.append('GIZMO-ps2')
+ plt.xlim(6, 8.5)
+# plt.ylim(-0.1, 30)
+ plt.ylim(0.9, 50)
+ plt.semilogy()
+ plt.grid(True)
+ plt.xlabel("$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$")
+ if star_clump_stats_i == 1:
+ plt.ylabel("$\mathrm{Cumulative\ Clump\ Counts, \ \ N_{clump}(>M)}}$")
+ plt.legend(codes_plotted, loc=1, frameon=True, ncol=1, fancybox=True, labelspacing=0.15, borderpad=0.3)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='xx-small')
+ plt.savefig("star_clump_stats_HOP_cumulative_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+
+ fig = plt.figure(figsize=(8, 4))
+ gridspec.GridSpec(1, 2)
+ for star_clump_stats_i in range(1,3,1):
+ codes_plotted = []
+# plt.subplot(1,2,star_clump_stats_i, aspect=0.25)
+ plt.subplot2grid((1, 2), (0, star_clump_stats_i-1))
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ if (star_clump_stats_i == 1 and (codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES")) or \
+ (star_clump_stats_i == 2 and (codes[code] == "CHANGA" or codes[code] == "GASOLINE" or codes[code] == "GADGET-3" or codes[code] == "GEAR" or codes[code] == "GIZMO" or codes[code] == "SWIFT")):
+ hist = np.histogram(star_clump_masses_fof[time][code], bins=10, range=(6., 8.5))
+ dbin = 0.5*(hist[1][1] - hist[1][0])
+ lines = plt.plot(hist[1][:-1]+dbin, hist[0], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ codes_plotted.append(codes[code])
+ plt.xlim(6, 8.5)
+ plt.ylim(-0.1, 15)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$")
+ if star_clump_stats_i == 1:
+ plt.ylabel("$\mathrm{Stellar\ Clump\ Counts, \ \ N_{clump}(M)}$")
+ plt.legend(codes_plotted, loc=1, frameon=True, ncol=1, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='xx-small')
+ plt.savefig("star_clump_stats_FOF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ # Reiterate for cumulative plots
+ fig = plt.figure(figsize=(8, 4))
+ gridspec.GridSpec(1, 2)
+ for star_clump_stats_i in range(1,3,1):
+ codes_plotted = []
+# plt.subplot(1,2,star_clump_stats_i, aspect=0.15)
+ plt.subplot2grid((1, 2), (0, star_clump_stats_i-1))
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ if (star_clump_stats_i == 1 and (codes[code] == "ART-I" or codes[code] == "ART-II" or codes[code] == "ENZO" or codes[code] == "RAMSES")) or \
+ (star_clump_stats_i == 2 and (codes[code] == "CHANGA" or codes[code] == "GASOLINE" or codes[code] == "GADGET-3" or codes[code] == "GEAR" or codes[code] == "GIZMO" or codes[code] == "SWIFT")):
+ hist = np.histogram(star_clump_masses_fof[time][code], bins=10, range=(6., 8.5))
+ dbin = 0.5*(hist[1][1] - hist[1][0])
+ lines = plt.plot(hist[1][:-1]+dbin, np.cumsum(hist[0][::-1])[::-1], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], \
+ marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ codes_plotted.append(codes[code])
+ if draw_star_clump_stats == 3 and star_clump_stats_i == 2:
+ hist_clump_ref = np.histogram(star_clump_masses_fof_ref, bins=10, range=(6., 8.5))
+ dbin = 0.5*(hist[1][1] - hist[1][0])
+ lines = plt.plot(hist_clump_ref[1][:-1]+dbin, np.cumsum(hist_clump_ref[0][::-1])[::-1], color='k', linestyle='--', marker='*', markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ codes_plotted.append('GIZMO-ps2')
+ plt.xlim(6, 8.5)
+# plt.ylim(-0.1, 30)
+ plt.ylim(0.9, 50)
+ plt.semilogy()
+ plt.grid(True)
+ plt.xlabel("$\mathrm{log[Newly\ Formed\ Stellar\ Clump\ Mass\ (M_{\odot})]}$")
+ if star_clump_stats_i == 1:
+ plt.ylabel("$\mathrm{Cumulative\ Clump\ Counts, \ \ N_{clump}(>M)}}$")
+ plt.legend(codes_plotted, loc=1, frameon=True, ncol=1, fancybox=True, labelspacing=0.15, borderpad=0.3)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='xx-small')
+ plt.savefig("star_clump_stats_FOF_cumulative_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_SFR_map >= 1 and time != 0:
+ fig_degr_density_map[time].savefig("degraded_Sigma_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ fig_degr_sfr_map[time].savefig("degraded_SFR_map_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_SFR_map == 2 and time != 0:
+ # If requested, draw the local K-S plot using the degraded maps created above
+ plt.clf()
+# plt.subplot(111)
+ fig = plt.figure(figsize=(8, 7))
+ Bigiel_surface_density = []
+ Bigiel_sfr_surface_density = []
+ for code in range(len(codes)):
+ # Remove bins where SFR surface density is zero
+ KS_x = np.array(surf_dens_SFR_map[time][code])
+ KS_x[np.where(np.array(sfr_surf_dens_SFR_map[time][code]) < 1e-10)] = 0
+ KS_x = np.log10(KS_x)
+ KS_y = np.log10(np.array(sfr_surf_dens_SFR_map[time][code]))
+ KS_x = KS_x[~np.isinf(KS_x)]
+ KS_y = KS_y[~np.isinf(KS_y)]
+ if codes[code] == "CHANGA":
+ plt.scatter(KS_x, KS_y, color='k', edgecolor='k', s=20, linewidth=0.7, marker=marker_names[code], alpha=0.1)
+ # Draw a 80% contour rather than scattering all the datapoints; see http://stackoverflow.com/questions/19390320/scatterplot-contours-in-matplotlib
+ Gaussian_density_estimation_nbins = 20
+ kernel = kde.gaussian_kde(np.vstack([KS_x, KS_y]))
+ xi, yi = np.mgrid[KS_x.min():KS_x.max():Gaussian_density_estimation_nbins*1j, KS_y.min():KS_y.max():Gaussian_density_estimation_nbins*1j]
+ zi = np.reshape(kernel(np.vstack([xi.flatten(), yi.flatten()])), xi.shape)
+ contours = plt.contour(xi, yi, zi, np.array([0.2]), colors=color_names[code], linewidths=1.2, alpha=1.0) # 80% percentile contour
+# contours = plt.contour(xi, yi, zi, np.array([0.3173]), colors=color_names[code], linewidths=1.2, alpha=1.0) # 68.27% percentile contour (1-sigma)
+ contours.collections[0].set_label(codes[code]) # setting names for legend
+ plt.clabel(contours, fmt=codes[code], inline=True, fontsize=7) # setting labels
+ plt.xlim(0, 3)
+ plt.ylim(-4, 1)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{log[Gas\ Surface\ Density\ (M_{\odot}/pc^2)]}$")
+ plt.ylabel("$\mathrm{log[Star\ Formation\ Rate\ Surface\ Density\ (M_{\odot}/yr/kpc^2)]}$")
+ plt.legend(loc=2, frameon=True, ncol=2, fancybox=True) # note difference from others since we want the legends to list contours, not scattered points
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ t = np.arange(-2, 5, 0.01)
+ for line in import_text("./Bigieletal2008_Fig8_Contour.txt", " "):
+ Bigiel_surface_density.append(float(line[0]) + np.log10(1.36)) # for the factor 1.36, see Section 2.3.1 of Bigiel et al. 2008
+ Bigiel_sfr_surface_density.append(float(line[1]))
+ plt.fill(Bigiel_surface_density, Bigiel_sfr_surface_density, fill=True, color='b', alpha = 0.1, hatch='\\') # contour by Bigiel et al. 2008 (Fig 8), or Feldmann et al. 2012 (Fig 1)
+ plt.plot(t, 1.37*t - 3.78, 'k--', linewidth = 2, alpha = 0.7) # observational fit by Kennicutt et al. 2007
+# plt.axhline(y = np.log10(8.593e4/young_star_cutoff_SFR_map/1e6/(float(aperture_size_SFR_map)/1000.)**2),
+# color='k', linestyle ='-.', linewidth=1, alpha=0.7) # SFR surface density cutoff due to limited star particle mass resolution
+ plt.savefig("K-S_local_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_PDF >= 1:
+ if draw_PDF == 3 and time != 0:
+ # If requested, find a 1D profile line from a specific snapshot as a reference (in this case ENZO or CHANGA noSF run at 500 Myr), then we add it to all panels
+ # pf_profile_ref = load_dataset(1, 1, 0, ['ENZO'], [['dummy', file_location[0]+'ENZO/DD0100/DD0100']])
+ pf_profile_ref = load_dataset(1, 1, 0, ['CHANGA'], [['dummy', file_location[0]+'CHANGA/disklow/disklow.000500']])
+ def _Density_2(field, data):
+ return data[(PartType_Gas_to_use, "Density")].in_units('g/cm**3')
+ pf_profile_ref.add_field((PartType_Gas_to_use, "Density_2"), function=_Density_2, take_log=True, particle_type=False, display_name="Density", units="g/cm**3")
+ def _Temperature_2(field, data):
+ return data[(PartType_Gas_to_use, "Temperature")].in_units('K')
+ pf_profile_ref.add_field((PartType_Gas_to_use, "Temperature_2"), function=_Temperature_2, take_log=True, particle_type=False, display_name="Temperature", units="K")
+ def _Mass_2(field, data):
+ return data[(PartType_Gas_to_use, MassType_to_use)].in_units('Msun')
+ pf_profile_ref.add_field((PartType_Gas_to_use, "Mass_2"), function=_Mass_2, take_log=True, particle_type=False, display_name="Mass", units="Msun")
+
+ v, cen = pf_profile_ref.h.find_max(("gas", "density"))
+ sp = pf_profile_ref.sphere(cen, (30.0, "kpc"))
+ cen2 = sp.quantities.center_of_mass(use_gas=True, use_particles=False).in_units("kpc")
+ sp2 = pf_profile_ref.sphere(cen2, (1.0, "kpc"))
+ cen3 = sp2.quantities.max_location(("gas", "density"))
+ center_profile_ref = pf_profile_ref.arr([cen3[1].d, cen3[2].d, cen3[3].d], 'code_length')
+ if yt_version_pre_3_2_3 == 1:
+ center_profile_ref = pf_profile_ref.arr([cen3[2].d, cen3[3].d, cen3[4].d], 'code_length')
+ sp_profile_ref = pf_profile_ref.sphere(center_profile_ref, (0.5*figure_width, "kpc"))
+
+ # p35 = ProfilePlot(sp_profile_ref, ("gas", "density"), ("gas", "temperature"), weight_field=("gas", "cell_mass"), n_bins=30)
+ # p35.set_xlim(1e-29, 1e-21)
+ # p35.set_ylim("temperature", 10, 1e7)
+ p35 = ProfilePlot(sp_profile_ref, (PartType_Gas_to_use, "Density_2"), (PartType_Gas_to_use, "Temperature_2"), weight_field=(PartType_Gas_to_use, "Mass_2"), n_bins=30)
+ p35.set_xlim(1e-26, 1e-21)
+ p35.set_ylim("Temperature_2", 10, 1e7)
+ for code in range(len(codes)):
+# line_profile_ref = ln.Line2D(p35.profiles[0].x.in_units('g/cm**3'), p35.profiles[0]["temperature"].in_units('K'), linestyle="--", linewidth=2, color='k', alpha=0.7)
+ line_profile_ref = ln.Line2D(p35.profiles[0].x.in_units('g/cm**3'), p35.profiles[0]["Temperature_2"].in_units('K'), linestyle="--", linewidth=2, color='k', alpha=0.7)
+ grid_PDF[time][code].axes.add_line(line_profile_ref)
+ fig_PDF[time].savefig("PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_pos_vel_PDF >= 1:
+ fig_pos_vel_PDF[time].savefig("pos_vel_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_pos_vel_PDF >= 2 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=0.04)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(pos_vel_xs[time][code], pos_vel_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0, 270)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Rotational\ Velocity\ (km/s)}$", fontsize='large')
+ plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(pos_vel_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(pos_vel_xs[time][code], (pos_vel_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-0.5, 0.5)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{v} - \mathrm{\bar v})/\mathrm{\bar v}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(pos_vel_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < pos_vel_xs[time][0]) & (pos_vel_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+# plt.text(12, 0.3, "mean dispersion = %.3f or %.3f dex " % (mean_fractional_dispersion, mean_fractional_dispersion_in_dex), ha="center", va="center", size=8, bbox=dict(boxstyle="round", fc="w", ec="0.5", alpha=0.9))
+ plt.savefig("pos_vel_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_pos_vel_PDF >= 3 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=0.064)
+ fig = plt.figure(figsize=(8, 10))
+ gridspec.GridSpec(5, 1)
+ plt.subplot2grid((5, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(pos_disp_xs[time][code], pos_disp_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0, 170)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Velocity\ Dispersion\ (km/s)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((5, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(pos_disp_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(pos_disp_xs[time][code], (pos_disp_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(pos_disp_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < pos_disp_xs[time][0]) & (pos_disp_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("pos_disp_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.subplot2grid((5, 1), (4, 0), rowspan=1)
+ for code in range(len(codes)):
+ lines = plt.plot(pos_disp_xs[time][code], pos_disp_vert_profiles[time][code]/pos_disp_profiles[time][code], color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Ratio,}\/\mathrm{\sigma}_{\perp}/\mathrm{\sigma}$", fontsize='large')
+ plt.savefig("pos_disp_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_pos_vel_PDF == 4 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=0.064)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(pos_disp_vert_xs[time][code], pos_disp_vert_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0, 170)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Vertical\ Velocity\ Dispersion\ (km/s)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(pos_disp_vert_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(pos_disp_vert_xs[time][code], (pos_disp_vert_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(pos_disp_vert_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < pos_disp_vert_xs[time][0]) & (pos_disp_vert_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("pos_disp_vert_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_star_pos_vel_PDF >= 1 and time != 0:
+ fig_star_pos_vel_PDF[time].savefig("star_pos_vel_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_star_pos_vel_PDF >= 2 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=0.04)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(star_pos_vel_xs[time][code], star_pos_vel_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0, 270)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Rotational\ Velocity\ (km/s)}$", fontsize='large')
+ plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(star_pos_vel_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(star_pos_vel_xs[time][code], (star_pos_vel_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-0.5, 0.5)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{v} - \mathrm{\bar v})/\mathrm{\bar v}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(star_pos_vel_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < star_pos_vel_xs[time][0]) & (star_pos_vel_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("star_pos_vel_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_star_pos_vel_PDF >= 3 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=0.064)
+ fig = plt.figure(figsize=(8, 10))
+ gridspec.GridSpec(5, 1)
+ plt.subplot2grid((5, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(star_pos_disp_xs[time][code], star_pos_disp_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0, 170)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Velocity\ Dispersion\ (km/s)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((5, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(star_pos_disp_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(star_pos_disp_xs[time][code], (star_pos_disp_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(star_pos_disp_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < star_pos_disp_xs[time][0]) & (star_pos_disp_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("star_pos_disp_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.subplot2grid((5, 1), (4, 0), rowspan=1)
+ for code in range(len(codes)):
+ lines = plt.plot(star_pos_disp_xs[time][code], star_pos_disp_vert_profiles[time][code]/star_pos_disp_profiles[time][code], color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Ratio,}\/\mathrm{\sigma}_{\perp}/\mathrm{\sigma}$", fontsize='large')
+ plt.savefig("star_pos_disp_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_star_pos_vel_PDF == 4 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=0.064)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(star_pos_disp_vert_xs[time][code], star_pos_disp_vert_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0, 170)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Vertical\ Velocity\ Dispersion\ (km/s)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(star_pos_disp_vert_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(star_pos_disp_vert_xs[time][code], (star_pos_disp_vert_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\sigma} - \mathrm{\bar \sigma})/\mathrm{\bar \sigma}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(star_pos_disp_vert_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < star_pos_disp_vert_xs[time][0]) & (star_pos_disp_vert_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("star_pos_disp_vert_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_rad_height_PDF >= 1:
+ fig_rad_height_PDF[time].savefig("rad_height_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_rad_height_PDF == 2 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=18)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(rad_height_xs[time][code], rad_height_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(0, 0.45)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Average\ Vertical\ Height\ (kpc)}$", fontsize='large')
+ plt.legend(codes, loc=2, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(rad_height_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(rad_height_xs[time][code], (rad_height_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{h} - \mathrm{\bar h})/\mathrm{\bar h}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(rad_height_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < rad_height_xs[time][0]) & (rad_height_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("rad_height_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_metal_PDF == 1:
+ fig_metal_PDF[time].savefig("metal_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ if draw_density_DF >= 1:
+ plt.clf()
+# plt.subplot(111, aspect=1)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(density_DF_xs[time][code], density_DF_profiles[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogx()
+ plt.semilogy()
+ plt.xlim(1e-28, 1e-21)
+ plt.ylim(1e4, 1e9) # accumulation=False
+ #plt.ylim(1e5, 2e10) # accumulation=True
+ if time != 0 and dataset_num == 2:
+ plt.axvline(x = 1.67e-24*10, color='k', linestyle ='--', linewidth=2, alpha=0.7) # SF threshold density
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Mass,}\/\mathrm{d}M\mathrm{/dlog}\/\mathrm{\rho}\/\mathrm{(M_{\odot})}$", fontsize='large')
+ plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(density_DF_profiles[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(density_DF_xs[time][code], (density_DF_profiles[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogx()
+ plt.xlim(1e-28, 1e-21)
+ plt.ylim(-1.0, 3.0)
+ plt.grid(True)
+ plt.locator_params(axis='y',nbins=4)
+ plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{m} - \mathrm{\bar m})/\mathrm{\bar m}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(density_DF_profiles[time]), axis=0)/ave_profiles)[(mean_dispersion_density_range[0] < density_DF_xs[time][0]) & (density_DF_xs[time][0] < mean_dispersion_density_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %.1e < rho < %.1e = %.3f (%.3f dex) " % (mean_dispersion_density_range[0], mean_dispersion_density_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.27, 0.84), size=10, xycoords='axes fraction')
+ plt.savefig("density_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_density_DF == 2 and time != 0 and dataset_num == 2:
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(16, 1)
+ plt.subplot2grid((16, 1), (0, 0), rowspan=9)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(density_DF_xs[time][code], (density_DF_profiles[time][code] - density_DF_1st_profiles[time][code])/1e8, color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xscale('log')
+ plt.xlim(1e-28, 1e-21)
+ plt.ylim(-7, 3)
+ plt.axvline(x = 1.67e-24*10, color='k', linestyle ='--', linewidth=2, alpha=0.7) # SF threshold density
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Mass\ Change,}\/\mathrm{\Delta}(\mathrm{d}M\mathrm{/dlog}\/\mathrm{\rho})\/\mathrm{(10^8 M_{\odot})}$", fontsize='large')
+ plt.legend(codes, loc=3, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((16, 1), (9, 0), rowspan=7)
+ for code in range(len(codes)):
+ lines = plt.plot(density_DF_xs[time][code], (density_DF_profiles[time][code] - density_DF_1st_profiles[time][code])/1e8, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize='large')
+ plt.ylabel(r"$\mathtt{symlog}[\/\mathrm{\Delta}(\mathrm{d}M\mathrm{/dlog}\/\mathrm{\rho})\/\mathrm{(10^8 M_{\odot})}\/]$", fontsize='large')
+ plt.xscale('log')
+ plt.yscale('symlog', linthreshy=0.01)
+ plt.xlim(1e-28, 1e-21)
+ plt.ylim(-10, 10)
+ plt.axvline(x = 1.67e-24*10, color='k', linestyle ='--', linewidth=2, alpha=0.7) # SF threshold density
+ plt.grid(True)
+ plt.savefig("density_DF_change_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_radius_DF == 1:
+ plt.clf()
+# plt.subplot(111, aspect=1)
+ fig = plt.figure(figsize=(8, 6))
+ for code in range(len(codes)):
+ lines = plt.plot(radius_DF_xs[time][code], np.add.accumulate(radius_DF_profiles[time][code]), color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(0, 14)
+ plt.ylim(1e7, 2e10)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Mass\ (M_{\odot})}$", fontsize='large')
+ plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ plt.savefig("radius_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+
+ plt.clf()
+# plt.subplot(111, aspect=1)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ temp = []
+ dr = 0.5*(radius_DF_xs[time][code][1] - radius_DF_xs[time][code][0]) # Here we assume that ProfilePlot was made with linearly binned radius_DF_xs
+ for radius in range(len(radius_DF_profiles[time][code])):
+ surface_area = np.pi*(((radius_DF_xs[time][code][radius]+dr)*1e3)**2 - ((radius_DF_xs[time][code][radius]-dr)*1e3)**2)
+ temp.append(radius_DF_profiles[time][code][radius] / surface_area)
+ surface_density[time].append(temp)
+ lines = plt.plot(radius_DF_xs[time][code], surface_density[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(0, 14)
+ plt.ylim(1e-1, 2e3)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Surface\ Density\ (M_{\odot}/pc^2)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(surface_density[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(radius_DF_xs[time][code], (surface_density[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(surface_density[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < radius_DF_xs[time][0]) & (radius_DF_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("gas_surface_density_radial_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_star_radius_DF >= 1 and time != 0:
+ plt.clf()
+# plt.subplot(111, aspect=1)
+ fig = plt.figure(figsize=(8, 6))
+ for code in range(len(codes)):
+ lines = plt.plot(star_radius_DF_xs[time][code], np.add.accumulate(star_radius_DF_profiles[time][code]), color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(0, 14)
+ plt.ylim(1e7, 2e10)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Newly\ Formed\ Stellar\ Mass\ (M_{\odot})}$", fontsize='large')
+ plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ plt.savefig("star_radius_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+
+ plt.clf()
+# plt.subplot(111, aspect=1)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ temp = []
+ dr = 0.5*(star_radius_DF_xs[time][code][1] - star_radius_DF_xs[time][code][0])
+ for radius in range(len(star_radius_DF_profiles[time][code])):
+ surface_area = np.pi*(((star_radius_DF_xs[time][code][radius]+dr)*1e3)**2 - ((star_radius_DF_xs[time][code][radius]-dr)*1e3)**2)
+ temp.append(star_radius_DF_profiles[time][code][radius] / surface_area)
+ star_surface_density[time].append(temp)
+ lines = plt.plot(star_radius_DF_xs[time][code], star_surface_density[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(0, 14)
+ plt.ylim(1e-1, 2e3)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Newly\ Formed\ Stellar\ Surface\ Density\ (M_{\odot}/pc^2)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(star_surface_density[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(star_radius_DF_xs[time][code], (star_surface_density[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 3.0)
+ plt.grid(True)
+ plt.locator_params(axis='y',nbins=4)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(star_surface_density[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < star_radius_DF_xs[time][0]) & (star_radius_DF_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.02, 0.84), size=10, xycoords='axes fraction')
+ plt.savefig("star_surface_density_radial_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+
+ if draw_star_radius_DF == 2 and time != 0:
+# plt.subplot(111, aspect=1)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ temp = []
+ dr = 0.5*(sfr_radius_DF_xs[time][code][1] - sfr_radius_DF_xs[time][code][0])
+ for radius in range(len(sfr_radius_DF_profiles[time][code])):
+ surface_area = np.pi*((sfr_radius_DF_xs[time][code][radius]+dr)**2 - (sfr_radius_DF_xs[time][code][radius]-dr)**2)
+ temp.append(sfr_radius_DF_profiles[time][code][radius] / surface_area)
+ sfr_surface_density[time].append(temp)
+ lines = plt.plot(sfr_radius_DF_xs[time][code], sfr_surface_density[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(0, 14)
+ plt.ylim(1e-4, 1e1)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Star\ Formation\ Rate\ Surface\ Density\ (M_{\odot}/yr/kpc^2)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(sfr_surface_density[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(sfr_radius_DF_xs[time][code], (sfr_surface_density[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 14)
+ plt.ylim(-1.0, 3.0)
+ plt.grid(True)
+ plt.locator_params(axis='y',nbins=4)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$", fontsize='large')
+ # if add_mean_fractional_dispersion == 1:
+ # mean_fractional_dispersion = (np.std(np.array(sfr_surface_density[time]), axis=0)/ave_profiles)[(mean_dispersion_radius_range[0] < sfr_radius_DF_xs[time][0]) & (sfr_radius_DF_xs[time][0] < mean_dispersion_radius_range[1])].mean()
+ # mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ # plt.annotate("Mean fractional dispersion for %d < r < %d kpc = %.3f (%.3f dex) " % (mean_dispersion_radius_range[0], mean_dispersion_radius_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.02, 0.84), size=10, xycoords='axes fraction')
+ plt.savefig("sfr_surface_density_radial_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+
+ # Draw K-S plot; below assumes that surface_density (or radius_DF_profiles) and sfr_surface_density (or sfr_radius_DF_profiles) have the identical size (n_bins in ProfilePlot)
+# plt.subplot(111)
+ fig = plt.figure(figsize=(8, 7))
+ KS_fit_t1 = []
+ KS_fit_t2 = []
+ Bigiel_surface_density = []
+ Bigiel_sfr_surface_density = []
+ for code in range(len(codes)):
+ # Remove bins where SFR surface density is zero
+ KS_x = np.array(surface_density[time][code])
+ KS_x[np.where(np.array(sfr_surface_density[time][code]) < 1e-10)] = 0
+ KS_x = np.log10(KS_x)
+ KS_y = np.log10(np.array(sfr_surface_density[time][code]))
+ KS_x = KS_x[~np.isinf(KS_x)]
+ KS_y = KS_y[~np.isinf(KS_y)]
+ t1, t2 = np.polyfit(KS_x, KS_y, 1)
+ KS_fit_t1.append(t1)
+ KS_fit_t2.append(t2)
+ plt.scatter(KS_x, KS_y, color=color_names[code], edgecolor=color_names[code], s=30, linewidth=0.7, marker=marker_names[code], alpha=0.8)
+ plt.xlim(0, 3)
+ plt.ylim(-4, 1)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{log[Gas\ Surface\ Density\ (M_{\odot}/pc^2)]}$")
+ plt.ylabel("$\mathrm{log[Star\ Formation\ Rate\ Surface\ Density\ (M_{\odot}/yr/kpc^2)]}$")
+ plt.legend(codes, loc=2, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ t = np.arange(-2, 5, 0.01)
+ for line in import_text("./Bigieletal2008_Fig8_Contour.txt", " "):
+ Bigiel_surface_density.append(float(line[0]) + np.log10(1.36)) # for the factor 1.36, see Section 2.3.1 of Bigiel et al. 2008
+ Bigiel_sfr_surface_density.append(float(line[1]))
+ plt.fill(Bigiel_surface_density, Bigiel_sfr_surface_density, fill=True, color='b', alpha = 0.1, hatch='\\') # contour by Bigiel et al. 2008 (Fig 8), or Feldmann et al. 2012 (Fig 1)
+# plt.plot(Bigiel_surface_density, Bigiel_sfr_surface_density, 'b-', linewidth = 0.7, alpha = 0.4)
+ plt.plot(t, 1.37*t - 3.78, 'k--', linewidth = 2, alpha = 0.7) # observational fit by Kennicutt et al. 2007
+ plt.savefig("K-S_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ for code in range(len(codes)):
+ plt.plot(t, np.polyval([KS_fit_t1[code], KS_fit_t2[code]], t), color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], alpha = 0.7) # linear fits
+ plt.savefig("K-S_with_fits_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_height_DF == 1:
+ plt.clf()
+# plt.subplot(111, aspect=0.5)
+ fig = plt.figure(figsize=(8, 6))
+ for code in range(len(codes)):
+ lines = plt.plot(height_DF_xs[time][code], np.add.accumulate(height_DF_profiles[time][code]), color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(0, 1.4)
+ plt.ylim(1e9, 2e10)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Mass\ (M_{\odot})}$", fontsize='large')
+ plt.legend(codes, loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ plt.savefig("height_DF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+
+ plt.clf()
+# plt.subplot(111, aspect=0.8)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ temp = []
+ dh = height_DF_xs[time][code][1] - height_DF_xs[time][code][0]
+ for height in range(len(height_DF_profiles[time][code])):
+ surface_area = 2 * dh*1e3 * figure_width*1e3 # surface_area = 2*d(height)*figure_width in pc^2
+ temp.append(height_DF_profiles[time][code][height] / surface_area)
+ height_surface_density[time].append(temp)
+ lines = plt.plot(height_DF_xs[time][code], height_surface_density[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(0, 1.4)
+ plt.ylim(1e-1, 2e3)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Surface\ Density\ (M_{\odot}/pc^2)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(height_surface_density[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(height_DF_xs[time][code], (height_surface_density[time][code] - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, 1.4)
+ plt.ylim(-1.0, 3.0)
+ plt.grid(True)
+ plt.locator_params(axis='y',nbins=4)
+ plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\Sigma} - \mathrm{\bar \Sigma})/\mathrm{\bar \Sigma}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(height_surface_density[time]), axis=0)/ave_profiles)[(mean_dispersion_height_range[0] < height_DF_xs[time][0]) & (height_DF_xs[time][0] < mean_dispersion_height_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %.1f < z < %.1f kpc = %.3f (%.3f dex) " % (mean_dispersion_height_range[0], mean_dispersion_height_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.02, 0.84), size=10, xycoords='axes fraction')
+ plt.savefig("gas_surface_density_vertical_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_SFR >= 1 and time != 0:
+ plt.clf()
+# plt.subplot(111)#, aspect=1e-7)
+ fig = plt.figure(figsize=(8, 6))
+ for code in range(len(codes)):
+ lines = plt.plot(sfr_ts[time][code], sfr_cum_masses[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, times[time])
+ plt.ylim(0, 2.5e9)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Time\ (Myr)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Stellar\ Mass\ (M_{\odot})}$", fontsize='large')
+ plt.legend(codes, loc=2, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ plt.savefig("Stellar_mass_evolution_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+# plt.subplot(111, aspect=55)
+ fig = plt.figure(figsize=(8, 8))
+ gridspec.GridSpec(4, 1)
+ plt.subplot2grid((4, 1), (0, 0), rowspan=3)
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ for code in range(len(codes)):
+ lines = plt.plot(sfr_ts[time][code], sfr_sfrs[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, times[time])
+ plt.ylim(0, 8)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Time\ (Myr)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Star\ Formation\ Rate\ (M_{\odot}/yr)}$", fontsize='large')
+ if draw_SFR == 2:
+ pf_sfr_ref = load_dataset(2, 1, 0, ['GIZMO'],
+ [['dummy', '/lustre/ki/pfs/mornkr/080112_CHaRGe/pfs-hyades/AGORA-DISK-repository/Grackle+SF+ThermalFbck/GIZMO/AGORA_disk_second_ps1.5/snapshot_temp_100_last']])
+ pf_sfr_ref.unit_registry.add('pccm', pf_sfr_ref.unit_registry.lut['pc'][0], length, "\\rm{pc}/(1+z)")
+ pf_sfr_ref.hubble_constant = 0.71; pf_sfr_ref.omega_lambda = 0.73; pf_sfr_ref.omega_matter = 0.27; pf_sfr_ref.omega_curvature = 0.0
+ sp_sfr_ref = pf_sfr_ref.sphere(center, (0.5*figure_width, "kpc"))
+ draw_SFR_mass_ref = sp_sfr_ref[(PartType_Star_to_use, "particle_mass")].in_units('Msun')
+ draw_SFR_ct_ref = sp_sfr_ref[(PartType_Star_to_use, FormationTimeType_to_use)].in_units('Myr')
+ sfr_ref = StarFormationRate(pf_sfr_ref, star_mass = draw_SFR_mass_ref, star_creation_time = draw_SFR_ct_ref,
+ volume = sp_sfr_ref.volume(), bins = 25)
+ lines = plt.plot(sfr_ref.time.in_units('Myr'), sfr_ref.Msol_yr, color='k', linestyle='--', marker='*', markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.legend(codes+['GIZMO-ps2'], loc=2, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+ ave_profiles = np.mean(np.array(sfr_sfrs[time]), axis=0)
+ for code in range(len(codes)):
+ lines = plt.plot(sfr_ts[time][code], (sfr_sfrs[time][code].d - ave_profiles)/ave_profiles, color=color_names[code], linestyle='none', marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.xlim(0, times[time])
+ plt.ylim(-1.0, 1.0)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Time\ (Myr)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Residual,}\/(\mathrm{\dot{\rho}_\star} - \mathrm{\bar{\dot{\rho}_\star}})/\mathrm{\bar{\dot{\rho}_\star}}$", fontsize='large')
+ if add_mean_fractional_dispersion == 1:
+ mean_fractional_dispersion = (np.std(np.array(sfr_sfrs[time]), axis=0)/ave_profiles)[(mean_dispersion_time_range[0] < sfr_ts[time][0]) & (sfr_ts[time][0] < mean_dispersion_time_range[1])].mean()
+ mean_fractional_dispersion_in_dex = np.log10(1.+mean_fractional_dispersion)
+ plt.annotate("Mean fractional dispersion for %d < t < %d Myr = %.3f (%.3f dex) " % (mean_dispersion_time_range[0], mean_dispersion_time_range[1], mean_fractional_dispersion, mean_fractional_dispersion_in_dex), xy=(0.35, 0.08), size=10, xycoords='axes fraction')
+ plt.savefig("SFR_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+ if draw_cut_through == 1:
+ plt.clf()
+# plt.subplot(111, aspect=0.5)
+ fig = plt.figure(figsize=(8, 6))
+ for code in range(len(codes)):
+ lines = plt.plot(cut_through_zs[time][code], cut_through_zvalues[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(-1.4, 1.4)
+ plt.ylim(1e-26, 1e-22)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Vertical\ Height\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Density\ (g/cm^3)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ z = np.arange(-1.4, 1.4, 0.05)
+ plt.plot(z, rho_agora_disk(0, np.abs(z)), linestyle="--", linewidth=2, color='k', alpha=0.7)
+ plt.savefig("cut_through_z_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
+# plt.subplot(111, aspect=0.5)
+ fig = plt.figure(figsize=(8, 6))
+ for code in range(len(codes)):
+ lines = plt.plot(cut_through_xs[time][code], cut_through_xvalues[time][code], color=color_names[code], linestyle=linestyle_names[np.mod(code, len(linestyle_names))], marker=marker_names[code], markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogy()
+ plt.xlim(-14, 14)
+ plt.ylim(1e-26, 1e-22)
+ plt.grid(True)
+ plt.xlabel("$\mathrm{Cylindrical\ Radius\ (kpc)}$", fontsize='large')
+ plt.ylabel("$\mathrm{Density\ (g/cm^3)}$", fontsize='large')
+ plt.legend(codes, loc=1, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ x = np.arange(-14, 14, 0.05)
+ plt.plot(x, rho_agora_disk(np.abs(x), 0), linestyle="--", linewidth=2, color='k', alpha=0.7)
+ plt.savefig("cut_through_x_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
+ plt.clf()
diff --git a/examples/AgoraDisk/run.sh b/examples/AgoraDisk/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..d7e284db52c2e6750fd713b3607a7f423bac7769
--- /dev/null
+++ b/examples/AgoraDisk/run.sh
@@ -0,0 +1,52 @@
+#!/bin/bash
+
+# This example is based on the AGORA disk article (DOI: 10.3847/1538-4357/833/2/202)
+
+# currently only the low resolution is available
+sim=low
+
+# enable cooling or not
+cooling=0
+
+# make run.sh fail if a subcommand fails
+set -e
+
+# define flags
+flag=
+if [ $cooling -eq 1 ]
+then
+ flag=-C
+fi
+
+# Generate the initial conditions if they are not present.
+if [ ! -e $sim.hdf5 ]
+then
+ echo "Fetching initial glass file for the Sedov blast example..."
+ ./getIC.sh $sim.hdf5
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+ echo "Fetching the Cloudy tables required by Grackle..."
+ ../getCoolingTable.sh
+fi
+
+# copy the initial conditions
+cp $sim.hdf5 agora_disk.hdf5
+# Update the particle types
+python3 changeType.py agora_disk.hdf5
+
+# Run SWIFT
+#../swift $flag -s -G -t 4 agora_disk.yml 2>&1 | tee output.log
+
+
+echo "Changing smoothing length to be Gadget compatible"
+python3 cleanupSwift.py agora_disk_0000.hdf5 agora_disk_IC.hdf5
+python3 cleanupSwift.py agora_disk_0050.hdf5 agora_disk_500Myr.hdf5
+
+echo "Fetching GEAR solution..."
+./getSolution.sh $flag
+
+echo "Plotting..."
+python3 plotSolution.py $flag
diff --git a/examples/CoolingBox/run.sh b/examples/CoolingBox/run.sh
index 61106d1f4819145ad4fb3a017918ca0e074a87c2..30b2177a6e8bb95a20146397f8b6a5021161b27f 100755
--- a/examples/CoolingBox/run.sh
+++ b/examples/CoolingBox/run.sh
@@ -17,7 +17,7 @@ fi
if [ ! -e CloudyData_UVB=HM2012.h5 ]
then
echo "Fetching the Cloudy tables required by Grackle..."
- ./getCoolingTable.sh
+ ../getCoolingTable.sh
fi
# Run SWIFT
diff --git a/examples/KeplerianRing/README.md b/examples/KeplerianRing/README.md
index d486cbfe412c7ff9ca121c87bbc548d0ece078dd..1c361f275d60ef1ca46d696e2e9507bb749e531c 100644
--- a/examples/KeplerianRing/README.md
+++ b/examples/KeplerianRing/README.md
@@ -2,7 +2,8 @@ Keplerian Ring Test
===================
This test uses the '2D Keplerian Ring' to assess the accuracy of SPH
-calculations. For more information on the test, please see Cullen & Dehnen 2009
+calculations. It is known to fail with all schemes in SWIFT. For more
+information on the test, please see Cullen & Dehnen 2009
(http://arxiv.org/abs/1006.1524) Section 4.3.
It is key that for this test in particular that the initial conditions are
diff --git a/examples/KeplerianRing/run.sh b/examples/KeplerianRing/run.sh
index ee8f7b247f38b74a204f9caf2bd543b72c4f94fa..0195846a8839a27083594c20569b1fd4d49f4c16 100755
--- a/examples/KeplerianRing/run.sh
+++ b/examples/KeplerianRing/run.sh
@@ -10,3 +10,9 @@ fi
rm -rf keplerian_ring_*.hdf5
../swift -g -s -t 1 -v 1 keplerian_ring.yml 2>&1 | tee output.log
+
+echo
+echo
+echo "This test is expected to fail, please read the README.md file"
+echo
+echo
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 2c5376a678b9d816e9575509b094e0a4168a6e42..b9c601e039d3f7566e56b0220935e63e42994b28 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -56,19 +56,23 @@ swift_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_
swift_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(MPI_LIBS) $(EXTRA_LIBS)
# Scripts to generate ICs
-EXTRA_DIST = BigCosmoVolume/makeIC.py \
- BigPerturbedBox/makeIC_fcc.py \
- CosmoVolume/cosmoVolume.yml CosmoVolume/getIC.sh CosmoVolume/run.sh \
- CoolingBox/coolingBox.yml CoolingBox/energy_plot.py CoolingBox/makeIC.py CoolingBox/run.sh \
+EXTRA_DIST = CoolingBox/coolingBox.yml CoolingBox/energy_plot.py CoolingBox/makeIC.py CoolingBox/run.sh \
EAGLE_6/eagle_6.yml EAGLE_6/getIC.sh EAGLE_6/README EAGLE_6/run.sh \
EAGLE_12/eagle_12.yml EAGLE_12/getIC.sh EAGLE_12/README EAGLE_12/run.sh \
EAGLE_25/eagle_25.yml EAGLE_25/getIC.sh EAGLE_25/README EAGLE_25/run.sh \
EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \
EAGLE_100/eagle_100.yml EAGLE_100/getIC.sh EAGLE_100/README EAGLE_100/run.sh \
+ EAGLE_DMO_12/eagle_12.yml EAGLE_DMO_12/getIC.sh EAGLE_DMO_12/README EAGLE_DMO_12/run.sh \
+ EAGLE_DMO_25/eagle_25.yml EAGLE_DMO_25/getIC.sh EAGLE_DMO_25/README EAGLE_DMO_25/run.sh \
+ EAGLE_DMO_50/eagle_50.yml EAGLE_DMO_50/getIC.sh EAGLE_DMO_50/README EAGLE_DMO_50/run.sh \
+ EAGLE_DMO_100/eagle_100.yml EAGLE_DMO_100/getIC.sh EAGLE_DMO_100/README EAGLE_DMO_100/run.sh \
+ EvrardCollapse_3D/evrard.yml EvrardCollapse_3D/makeIC.py EvrardCollapse_3D/plotSolution.py EvrardCollapse_3D/run.sh EvrardCollapse_3D/getReference.sh \
ExternalPointMass/externalPointMass.yml ExternalPointMass/makeIC.py ExternalPointMass/run.sh ExternalPointMass/energy_plot.py \
GreshoVortex_2D/getGlass.sh GreshoVortex_2D/gresho.yml GreshoVortex_2D/makeIC.py GreshoVortex_2D/plotSolution.py GreshoVortex_2D/run.sh \
+ GreshoVortex_3D/getGlass.sh GreshoVortex_3D/gresho.yml GreshoVortex_3D/makeIC.py GreshoVortex_3D/plotSolution.py GreshoVortex_3D/run.sh \
HydrostaticHalo/README HydrostaticHalo/hydrostatic.yml HydrostaticHalo/makeIC.py HydrostaticHalo/run.sh \
HydrostaticHalo/density_profile.py HydrostaticHalo/velocity_profile.py HydrostaticHalo/internal_energy_profile.py HydrostaticHalo/test_energy_conservation.py \
+ InteractingBlastWaves_1D/run.sh InteractingBlastWaves_1D/makeIC.py InteractingBlastWaves_1D/plotSolution.py InteractingBlastWaves_1D/interactingBlastWaves.yml InteractingBlastWaves_1D/getReference.sh \
IsothermalPotential/README IsothermalPotential/run.sh IsothermalPotential/energy_plot.py IsothermalPotential/isothermal.yml IsothermalPotential/makeIC.py \
KelvinHelmholtz_2D/kelvinHelmholtz.yml KelvinHelmholtz_2D/makeIC.py KelvinHelmholtz_2D/plotSolution.py KelvinHelmholtz_2D/run.sh \
MultiTypes/makeIC.py MultiTypes/multiTypes.yml MultiTypes/run.sh \
@@ -83,13 +87,15 @@ EXTRA_DIST = BigCosmoVolume/makeIC.py \
SineWavePotential_1D/makeIC.py SineWavePotential_1D/plotSolution.py SineWavePotential_1D/run.sh SineWavePotential_1D/sineWavePotential.yml \
SineWavePotential_2D/makeIC.py SineWavePotential_2D/plotSolution.py SineWavePotential_2D/run.sh SineWavePotential_2D/sineWavePotential.yml \
SineWavePotential_3D/makeIC.py SineWavePotential_3D/plotSolution.py SineWavePotential_3D/run.sh SineWavePotential_3D/sineWavePotential.yml \
+ SmallCosmoVolume/README SmallCosmoVolume/getIC.sh SmallCosmoVolume/run.sh SmallCosmoVolume/small_cosmo_volume.yml \
SodShock_1D/makeIC.py SodShock_1D/plotSolution.py SodShock_1D/run.sh SodShock_1D/sodShock.yml \
SodShock_2D/getGlass.sh SodShock_2D/makeIC.py SodShock_2D/plotSolution.py SodShock_2D/run.sh SodShock_2D/sodShock.yml \
SodShock_3D/getGlass.sh SodShock_3D/makeIC.py SodShock_3D/plotSolution.py SodShock_3D/run.sh SodShock_3D/sodShock.yml \
SquareTest_2D/makeIC.py SquareTest_2D/plotSolution.py SquareTest_2D/run.sh SquareTest_2D/square.yml \
UniformBox_2D/makeIC.py UniformBox_2D/run.sh UniformBox_2D/uniformPlane.yml \
UniformBox_3D/makeICbig.py UniformBox_3D/makeIC.py UniformBox_3D/run.sh UniformBox_3D/uniformBox.yml \
- UniformDMBox/makeIC.py
+ UniformDMBox/makeIC.py \
+ ZeldovichPancake_3D/makeIC.py ZeldovichPancake_3D/zeldovichPancake.yml ZeldovichPancake_3D/run.sh ZeldovichPancake_3D/plotSolution.py
# Default parameter file
EXTRA_DIST += parameter_example.yml
@@ -106,3 +112,5 @@ EXTRA_DIST += analyse_threadpool_tasks.py \
EXTRA_DIST += plot_scaling_results.py \
plot_scaling_results_breakdown.py
+# Script for gravity accuracy
+EXTRA_DIST += plot_gravity_checks.py
diff --git a/examples/ZeldovichPancake_3D/makeIC.py b/examples/ZeldovichPancake_3D/makeIC.py
index 170be68a59552a8d346ac5ad346b5dc0f2d43399..15fb8bdef95f6e830e78f4d1b2c419051a6f00af 100644
--- a/examples/ZeldovichPancake_3D/makeIC.py
+++ b/examples/ZeldovichPancake_3D/makeIC.py
@@ -40,7 +40,6 @@ k_in_J_K = 1.38064852e-23
# Parameters
rho_0 = 1.8788e-26 # h^2 kg m^-3
H_0 = 1. / Mpc_in_m * 10**5 # s^-1
-#lambda_i = 1.975e24 # h^-1 m (= 64 h^-1 Mpc)
lambda_i = 64. / H_0 * 10**5 # h^-1 m (= 64 h^-1 Mpc)
x_min = -0.5 * lambda_i
x_max = 0.5 * lambda_i
diff --git a/examples/CoolingBox/getCoolingTable.sh b/examples/getCoolingTable.sh
old mode 100755
new mode 100644
similarity index 100%
rename from examples/CoolingBox/getCoolingTable.sh
rename to examples/getCoolingTable.sh
diff --git a/examples/plot_gravity_checks.py b/examples/plot_gravity_checks.py
old mode 100644
new mode 100755
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 025779f17496e7bf30fdf12353c4381c7d6292ce..d70d58c6ba508ba4282ac9dd32565478afb40692 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_SHADOWSWIFT_HYDRO_H
+#define SWIFT_SHADOWSWIFT_HYDRO_H
#include <float.h>
#include "adiabatic_index.h"
@@ -41,15 +43,23 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float CFL_condition = hydro_properties->CFL_condition;
- float R = get_radius_dimension_sphere(p->cell.volume);
-
- if (p->timestepvars.vmax == 0.) {
- /* vmax can be zero in vacuum. We force the time step to become the maximal
- time step in this case */
- return FLT_MAX;
- } else {
- return CFL_condition * R / fabsf(p->timestepvars.vmax);
+ float vrel[3];
+ vrel[0] = p->primitives.v[0] - xp->v_full[0];
+ vrel[1] = p->primitives.v[1] - xp->v_full[1];
+ vrel[2] = p->primitives.v[2] - xp->v_full[2];
+ float vmax =
+ sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
+ sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ vmax = max(vmax, p->timestepvars.vmax);
+
+ const float psize =
+ cosmo->a *
+ powf(p->cell.volume / hydro_dimension_unit_sphere, hydro_dimension_inv);
+ float dt = FLT_MAX;
+ if (vmax > 0.) {
+ dt = psize / vmax;
}
+ return CFL_condition * dt;
}
/**
@@ -102,7 +112,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[2] = mass * p->primitives.v[2];
#ifdef EOS_ISOTHERMAL_GAS
- p->conserved.energy = mass * const_isothermal_internal_energy;
+ p->conserved.energy = mass * gas_internal_energy_from_entropy(0.f, 0.f);
#else
p->conserved.energy *= mass;
#endif
@@ -117,12 +127,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->v[0] = 0.;
p->v[1] = 0.;
p->v[2] = 0.;
+#else
+ p->v[0] = p->primitives.v[0];
+ p->v[1] = p->primitives.v[1];
+ p->v[2] = p->primitives.v[2];
#endif
/* set the initial velocity of the cells */
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
+
+ /* ignore accelerations present in the initial condition */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
}
/**
@@ -137,7 +156,8 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part* p, const struct hydro_space* hs) {
- p->density.wcount = 0.0f;
+ /* make sure we don't enter the no neighbour case in runner.c */
+ p->density.wcount = 1.0f;
p->density.wcount_dh = 0.0f;
voronoi_cell_init(&p->cell, p->x, hs->anchor, hs->side);
@@ -180,12 +200,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
if (hnew < p->h) {
/* Iteration succesful: we accept, but manually set h to a smaller value
for the next time step */
- p->density.wcount = 1.0f;
+ const float hinvdim = pow_dimension(1.0f / p->h);
+ p->density.wcount = hinvdim;
p->h = 1.1f * hnew;
} else {
/* Iteration not succesful: we force h to become 1.1*hnew */
p->density.wcount = 0.0f;
- p->density.wcount_dh = 1.0f / (1.1f * hnew - p->h);
+ p->density.wcount_dh = p->h / (1.1f * hnew - p->h);
return;
}
volume = p->cell.volume;
@@ -286,8 +307,48 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.v_full[0] = xp->v_full[0];
p->force.v_full[1] = xp->v_full[1];
p->force.v_full[2] = xp->v_full[2];
+
+ p->conserved.flux.mass = 0.0f;
+ p->conserved.flux.momentum[0] = 0.0f;
+ p->conserved.flux.momentum[1] = 0.0f;
+ p->conserved.flux.momentum[2] = 0.0f;
+ p->conserved.flux.energy = 0.0f;
+}
+
+/**
+ * @brief Prepare a particle for the gradient calculation.
+ *
+ * This function is called after the density loop and before the gradient loop.
+ *
+ * We use it to set the physical timestep for the particle and to copy the
+ * actual velocities, which we need to boost our interfaces during the flux
+ * calculation. We also initialize the variables used for the time step
+ * calculation.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
+ struct part* restrict p, struct xpart* restrict xp,
+ const struct cosmology* cosmo) {
+
+ /* Initialize time step criterion variables */
+ p->timestepvars.vmax = 0.;
}
+/**
+ * @brief Resets the variables that are required for a gradient calculation.
+ *
+ * This function is called after hydro_prepare_gradient.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_gradient(
+ struct part* restrict p) {}
+
/**
* @brief Finishes the gradient calculation.
*
@@ -385,53 +446,36 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
const struct hydro_props* hydro_props) {
- float vcell[3];
-
/* Update the conserved variables. We do this here and not in the kick,
since we need the updated variables below. */
- p->conserved.mass += p->conserved.flux.mass;
- p->conserved.momentum[0] += p->conserved.flux.momentum[0];
- p->conserved.momentum[1] += p->conserved.flux.momentum[1];
- p->conserved.momentum[2] += p->conserved.flux.momentum[2];
- p->conserved.energy += p->conserved.flux.energy;
+ p->conserved.mass += p->conserved.flux.mass * dt;
+ p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
+ p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt;
+ p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt;
#ifdef EOS_ISOTHERMAL_GAS
/* reset the thermal energy */
- p->conserved.energy = p->conserved.mass * const_isothermal_internal_energy;
-
-#ifdef SHADOWFAX_TOTAL_ENERGY
- p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] +
- p->conserved.momentum[1] * p->primitives.v[1] +
- p->conserved.momentum[2] * p->primitives.v[2]);
+ p->conserved.energy =
+ p->conserved.mass * gas_internal_energy_from_entropy(0.f, 0.f);
+#else
+ p->conserved.energy += p->conserved.flux.energy * dt;
#endif
-#endif
+#if defined(SHADOWFAX_FIX_CELLS)
+ p->v[0] = 0.0f;
+ p->v[1] = 0.0f;
+ p->v[2] = 0.0f;
+#else
+ if (p->conserved.mass > 0.0f && p->primitives.rho > 0.0f) {
- /* reset fluxes */
- /* we can only do this here, since we need to keep the fluxes for inactive
- particles */
- p->conserved.flux.mass = 0.0f;
- p->conserved.flux.momentum[0] = 0.0f;
- p->conserved.flux.momentum[1] = 0.0f;
- p->conserved.flux.momentum[2] = 0.0f;
- p->conserved.flux.energy = 0.0f;
+ const float inverse_mass = 1.f / p->conserved.mass;
- if (p->conserved.mass > 0.) {
- /* We want the cell velocity to be as close as possible to the fluid
- velocity */
- vcell[0] = p->conserved.momentum[0] / p->conserved.mass;
- vcell[1] = p->conserved.momentum[1] / p->conserved.mass;
- vcell[2] = p->conserved.momentum[2] / p->conserved.mass;
- } else {
- vcell[0] = 0.;
- vcell[1] = 0.;
- vcell[2] = 0.;
- }
+ /* Normal case: set particle velocity to fluid velocity. */
+ p->v[0] = p->conserved.momentum[0] * inverse_mass;
+ p->v[1] = p->conserved.momentum[1] * inverse_mass;
+ p->v[2] = p->conserved.momentum[2] * inverse_mass;
#ifdef SHADOWFAX_STEER_CELL_MOTION
- /* To prevent stupid things like cell crossovers or generators that move
- outside their cell, we steer the motion of the cell somewhat */
- if (p->primitives.rho) {
float centroid[3], d[3];
float volume, csnd, R, vfac, fac, dnrm;
voronoi_get_centroid(&p->cell, centroid);
@@ -449,32 +493,29 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
} else {
vfac = csnd / dnrm;
}
- } else {
- vfac = 0.0f;
+ p->v[0] += vfac * d[0];
+ p->v[1] += vfac * d[1];
+ p->v[2] += vfac * d[2];
}
- vcell[0] += vfac * d[0];
- vcell[1] += vfac * d[1];
- vcell[2] += vfac * d[2];
- }
#endif
-#if defined(SHADOWFAX_FIX_CELLS)
- xp->v_full[0] = 0.;
- xp->v_full[1] = 0.;
- xp->v_full[2] = 0.;
+ } else {
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+ }
+#endif
- p->v[0] = 0.;
- p->v[1] = 0.;
- p->v[2] = 0.;
-#else
- xp->v_full[0] = vcell[0];
- xp->v_full[1] = vcell[1];
- xp->v_full[2] = vcell[2];
+ /* Now make sure all velocity variables are up to date. */
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
- p->v[0] = xp->v_full[0];
- p->v[1] = xp->v_full[1];
- p->v[2] = xp->v_full[2];
-#endif
+ if (p->gpart) {
+ p->gpart->v_full[0] = p->v[0];
+ p->gpart->v_full[1] = p->v[1];
+ p->gpart->v_full[2] = p->v[2];
+ }
}
/**
@@ -799,3 +840,5 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_density(
return cosmo->a3_inv * p->primitives.rho;
}
+
+#endif /* SWIFT_SHADOWSWIFT_HYDRO_H */
diff --git a/src/hydro/Shadowswift/hydro_gradients.h b/src/hydro/Shadowswift/hydro_gradients.h
index 285d889a1a6e10662a06979f69290aabd4206059..4e7a9911d8d4fc586fe7a56687dd4c4ae9ec8de2 100644
--- a/src/hydro/Shadowswift/hydro_gradients.h
+++ b/src/hydro/Shadowswift/hydro_gradients.h
@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
struct part* pi, struct part* pj, float hi, float hj, const float* dx,
- float r, float* xij_i, float* Wi, float* Wj, float mindt) {
+ float r, float* xij_i, float* Wi, float* Wj) {
float dWi[5], dWj[5];
float xij_j[3];
@@ -132,59 +132,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
- /* time */
- dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
- Wi[2] * pi->primitives.gradients.rho[1] +
- Wi[3] * pi->primitives.gradients.rho[2] +
- Wi[0] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
- dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
- Wi[2] * pi->primitives.gradients.v[0][1] +
- Wi[3] * pi->primitives.gradients.v[0][2] +
- pi->primitives.gradients.P[0] / Wi[0]);
- dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
- Wi[2] * pi->primitives.gradients.v[1][1] +
- Wi[3] * pi->primitives.gradients.v[1][2] +
- pi->primitives.gradients.P[1] / Wi[0]);
- dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
- Wi[2] * pi->primitives.gradients.v[2][1] +
- Wi[3] * pi->primitives.gradients.v[2][2] +
- pi->primitives.gradients.P[2] / Wi[0]);
- dWi[4] -=
- 0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
- Wi[2] * pi->primitives.gradients.P[1] +
- Wi[3] * pi->primitives.gradients.P[2] +
- hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
-
- dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
- Wj[2] * pj->primitives.gradients.rho[1] +
- Wj[3] * pj->primitives.gradients.rho[2] +
- Wj[0] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
- dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
- Wj[2] * pj->primitives.gradients.v[0][1] +
- Wj[3] * pj->primitives.gradients.v[0][2] +
- pj->primitives.gradients.P[0] / Wj[0]);
- dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
- Wj[2] * pj->primitives.gradients.v[1][1] +
- Wj[3] * pj->primitives.gradients.v[1][2] +
- pj->primitives.gradients.P[1] / Wj[0]);
- dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
- Wj[2] * pj->primitives.gradients.v[2][1] +
- Wj[3] * pj->primitives.gradients.v[2][2] +
- pj->primitives.gradients.P[2] / Wj[0]);
- dWj[4] -=
- 0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
- Wj[2] * pj->primitives.gradients.P[1] +
- Wj[3] * pj->primitives.gradients.P[2] +
- hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
-
Wi[0] += dWi[0];
Wi[1] += dWi[1];
Wi[2] += dWi[2];
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 9ac1debf3184c25603412867c41c62a1131345f3..eda8e3759d9e08dac8073ebed9fb36dd0c5b99f6 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -143,7 +143,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float vmax, dvdotdx;
float vi[3], vj[3], vij[3];
float Wi[5], Wj[5];
- float dti, dtj, mindt;
float n_unit[3];
A = voronoi_get_face(&pi->cell, pj->id, xij_i);
@@ -168,9 +167,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[3] = pj->primitives.v[2];
Wj[4] = pj->primitives.P;
- dti = pi->force.dt;
- dtj = pj->force.dt;
-
/* calculate the maximal signal velocity */
vmax = 0.0f;
if (Wi[0] > 0.) {
@@ -192,10 +188,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pj->timestepvars.vmax = fmaxf(pj->timestepvars.vmax, vmax);
}
- /* The flux will be exchanged using the smallest time step of the two
- * particles */
- mindt = fminf(dti, dtj);
-
/* compute the normal vector of the interface */
for (k = 0; k < 3; ++k) {
n_unit[k] = -dx[k] / r;
@@ -219,13 +211,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[2] -= vij[1];
Wj[3] -= vij[2];
- hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj, mindt);
+ hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
/* we don't need to rotate, we can use the unit vector in the Riemann problem
* itself (see GIZMO) */
if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) {
- printf("mindt: %g\n", mindt);
printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0],
pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P);
#ifdef USE_GRADIENTS
@@ -266,20 +257,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Update conserved variables */
/* eqn. (16) */
- pi->conserved.flux.mass -= mindt * A * totflux[0];
- pi->conserved.flux.momentum[0] -= mindt * A * totflux[1];
- pi->conserved.flux.momentum[1] -= mindt * A * totflux[2];
- pi->conserved.flux.momentum[2] -= mindt * A * totflux[3];
- pi->conserved.flux.energy -= mindt * A * totflux[4];
+ pi->conserved.flux.mass -= A * totflux[0];
+ pi->conserved.flux.momentum[0] -= A * totflux[1];
+ pi->conserved.flux.momentum[1] -= A * totflux[2];
+ pi->conserved.flux.momentum[2] -= A * totflux[3];
+ pi->conserved.flux.energy -= A * totflux[4];
#ifndef SHADOWFAX_TOTAL_ENERGY
float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
pi->primitives.v[1] * pi->primitives.v[1] +
pi->primitives.v[2] * pi->primitives.v[2]);
- pi->conserved.flux.energy += mindt * A * totflux[1] * pi->primitives.v[0];
- pi->conserved.flux.energy += mindt * A * totflux[2] * pi->primitives.v[1];
- pi->conserved.flux.energy += mindt * A * totflux[3] * pi->primitives.v[2];
- pi->conserved.flux.energy -= mindt * A * totflux[0] * ekin;
+ pi->conserved.flux.energy += A * totflux[1] * pi->primitives.v[0];
+ pi->conserved.flux.energy += A * totflux[2] * pi->primitives.v[1];
+ pi->conserved.flux.energy += A * totflux[3] * pi->primitives.v[2];
+ pi->conserved.flux.energy -= A * totflux[0] * ekin;
#endif
/* here is how it works:
@@ -295,20 +286,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/
if (mode == 1 || pj->force.active == 0) {
- pj->conserved.flux.mass += mindt * A * totflux[0];
- pj->conserved.flux.momentum[0] += mindt * A * totflux[1];
- pj->conserved.flux.momentum[1] += mindt * A * totflux[2];
- pj->conserved.flux.momentum[2] += mindt * A * totflux[3];
- pj->conserved.flux.energy += mindt * A * totflux[4];
+ pj->conserved.flux.mass += A * totflux[0];
+ pj->conserved.flux.momentum[0] += A * totflux[1];
+ pj->conserved.flux.momentum[1] += A * totflux[2];
+ pj->conserved.flux.momentum[2] += A * totflux[3];
+ pj->conserved.flux.energy += A * totflux[4];
#ifndef SHADOWFAX_TOTAL_ENERGY
ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
pj->primitives.v[1] * pj->primitives.v[1] +
pj->primitives.v[2] * pj->primitives.v[2]);
- pj->conserved.flux.energy -= mindt * A * totflux[1] * pj->primitives.v[0];
- pj->conserved.flux.energy -= mindt * A * totflux[2] * pj->primitives.v[1];
- pj->conserved.flux.energy -= mindt * A * totflux[3] * pj->primitives.v[2];
- pj->conserved.flux.energy += mindt * A * totflux[0] * ekin;
+ pj->conserved.flux.energy -= A * totflux[1] * pj->primitives.v[0];
+ pj->conserved.flux.energy -= A * totflux[2] * pj->primitives.v[1];
+ pj->conserved.flux.energy -= A * totflux[3] * pj->primitives.v[2];
+ pj->conserved.flux.energy += A * totflux[0] * ekin;
#endif
}
}
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index c5448f77353e1859c1f8853394bbefbe26d0a3a9..f79fd832248fba8fbc55bd9fcec57e645be93159 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -72,6 +72,7 @@ void hydro_props_init(struct hydro_props *p,
/* change the meaning of target_neighbours and delta_neighbours */
p->target_neighbours = 1.0f;
p->delta_neighbours = 0.0f;
+ p->eta_neighbours = 1.0f;
#endif
/* Maximal smoothing length */
@@ -122,7 +123,7 @@ void hydro_props_init(struct hydro_props *p,
/* Compute the minimal energy (Note the temp. read is in internal units) */
double u_min = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
- u_min *= p->initial_temperature;
+ u_min *= p->minimal_temperature;
u_min *= hydro_one_over_gamma_minus_one;
/* Correct for hydrogen mass fraction */
diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c
index be24117fd7947e467c0ceca36a68e339952f7815..49e66763f7f9ef47ca79d0024fa65b9ce7803fc1 100644
--- a/src/mesh_gravity.c
+++ b/src/mesh_gravity.c
@@ -38,6 +38,8 @@
#include "runner.h"
#include "space.h"
+#ifdef HAVE_FFTW
+
/**
* @brief Returns 1D index of a 3D NxNxN array using row-major style.
*
@@ -264,6 +266,8 @@ void mesh_to_gparts_CIC(struct gpart* gp, const double* pot, int N, double fac,
#endif
}
+#endif
+
/**
* @brief Compute the potential, including periodic correction on the mesh.
*
@@ -532,7 +536,7 @@ void pm_mesh_struct_restore(struct pm_mesh* mesh, FILE* stream) {
restart_read_blocks((void*)mesh, sizeof(struct pm_mesh), 1, stream, NULL,
"gravity props");
-
+#ifdef HAVE_FFTW
const int N = mesh->N;
/* Allocate the memory for the combined density and potential array */
@@ -540,7 +544,6 @@ void pm_mesh_struct_restore(struct pm_mesh* mesh, FILE* stream) {
if (mesh->potential == NULL)
error("Error allocating memory for the long-range gravity mesh.");
-#ifdef HAVE_FFTW
#else
error("No FFTW library found. Cannot compute periodic long-range forces.");
#endif
diff --git a/src/partition.c b/src/partition.c
index f59187c7aea2abb8521d932770e25d544341f798..f39bcb1e2baf7916df1e369a6681435560fe9d41 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -1072,7 +1072,7 @@ static void repart_edge_parmetis(int vweights, int eweights, int timebins,
struct task *t = &tasks[j];
/* Skip un-interesting tasks. */
- if (t->cost == 0) continue;
+ if (t->cost == 0.f) continue;
/* Get the task weight based on costs. */
double w = (double)t->cost;
diff --git a/src/queue.c b/src/queue.c
index b22313d71d733ae1a6b3419d492356f9d914f61f..3a9919163a4fb9d146c055e58859a175858c17eb 100644
--- a/src/queue.c
+++ b/src/queue.c
@@ -259,7 +259,7 @@ struct task *queue_gettask(struct queue *q, const struct task *prev,
int k = ind;
if (k < qcount) {
qtid[k] = qtid[qcount];
- int w = qtasks[qtid[k]].weight;
+ const float w = qtasks[qtid[k]].weight;
while (k > 0 && w > qtasks[qtid[(k - 1) / 2]].weight) {
int temp = q->tid[k];
q->tid[k] = q->tid[(k - 1) / 2];
diff --git a/src/scheduler.c b/src/scheduler.c
index 844939ba1643379ff0ab5be4355212f80098d37d..38313251b89f9f7dd6a6e8dd5d93a7eba3894f5a 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1185,12 +1185,14 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
/* Run through the tasks backwards and set their weights. */
for (int k = nr_tasks - 1; k >= 0; k--) {
struct task *t = &tasks[tid[k]];
- t->weight = 0;
+ t->weight = 0.f;
for (int j = 0; j < t->nr_unlock_tasks; j++)
if (t->unlock_tasks[j]->weight > t->weight)
t->weight = t->unlock_tasks[j]->weight;
- float cost = 0;
+ float cost = 0.f;
+#if defined(WITH_MPI) && defined(HAVE_PARMETIS)
int partcost = 1;
+#endif
const float count_i = (t->ci != NULL) ? t->ci->count : 0.f;
const float count_j = (t->cj != NULL) ? t->cj->count : 0.f;
@@ -1277,14 +1279,18 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
cost = wscale * count_i + wscale * gcount_i;
break;
case task_type_send:
+#if defined(WITH_MPI) && defined(HAVE_PARMETIS)
partcost = 0;
+#endif
if (count_i < 1e5)
cost = 10.f * (wscale * count_i) * count_i;
else
cost = 2e9;
break;
case task_type_recv:
+#if defined(WITH_MPI) && defined(HAVE_PARMETIS)
partcost = 0;
+#endif
if (count_i < 1e5)
cost = 5.f * (wscale * count_i) * count_i;
else
@@ -1295,12 +1301,10 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
break;
}
-/* Note if the cost is > 1e9 we cap it as we don't care. That's large
- * compared to other possible costs. */
-#if defined(WITH_MPI) && defined(HAVE_METIS)
- if (partcost) t->cost = (cost < 1e9) ? cost : 1e9;
+#if defined(WITH_MPI) && defined(HAVE_PARMETIS)
+ if (partcost) t->cost = cost;
#endif
- t->weight += (cost < 1e9) ? cost : 1e9;
+ t->weight += cost;
}
if (verbose)
diff --git a/src/single_io.c b/src/single_io.c
index be33bba4d7c95454892b02836f073c085e07b2ea..a0f02878b52c89beca94d15c09ef7d456ce0a4eb 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -347,8 +347,8 @@ void read_ic_single(const char* fileName,
/* GADGET has only cubic boxes (in cosmological mode) */
double boxSize[3] = {0.0, -1.0, -1.0};
/* GADGET has 6 particle types. We only keep the type 0 & 1 for now...*/
- int numParticles[swift_type_count] = {0};
- int numParticles_highWord[swift_type_count] = {0};
+ long long numParticles[swift_type_count] = {0};
+ long long numParticles_highWord[swift_type_count] = {0};
size_t N[swift_type_count] = {0};
int dimension = 3; /* Assume 3D if nothing is specified */
size_t Ndm = 0;
@@ -388,13 +388,12 @@ void read_ic_single(const char* fileName,
io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
*flag_entropy = flag_entropy_temp[0];
io_read_attribute(h_grp, "BoxSize", DOUBLE, boxSize);
- io_read_attribute(h_grp, "NumPart_Total", UINT, numParticles);
- io_read_attribute(h_grp, "NumPart_Total_HighWord", UINT,
+ io_read_attribute(h_grp, "NumPart_Total", LONGLONG, numParticles);
+ io_read_attribute(h_grp, "NumPart_Total_HighWord", LONGLONG,
numParticles_highWord);
for (int ptype = 0; ptype < swift_type_count; ++ptype)
- N[ptype] = ((long long)numParticles[ptype]) +
- ((long long)numParticles_highWord[ptype] << 32);
+ N[ptype] = (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
/* Get the box size if not cubic */
dim[0] = boxSize[0];
diff --git a/src/task.h b/src/task.h
index a345d4a29522a709afa5e3abd2a163f7550fb7ed..eb352f272e2f3599b57d8c8b4559fbb6236442e4 100644
--- a/src/task.h
+++ b/src/task.h
@@ -136,11 +136,11 @@ struct task {
int rank;
/*! Weight of the task */
- int weight;
+ float weight;
#if defined(WITH_MPI) && defined(HAVE_PARMETIS)
/*! Individual cost estimate for this task. */
- int cost;
+ float cost;
#endif
/*! Number of tasks unlocked by this one */
diff --git a/src/units.c b/src/units.c
index 6152d7f42f9a409bea9d057862c04d0ae9fb6a75..04e74bc4d7040ed1bde73184b125eec5d8a7fe97 100644
--- a/src/units.c
+++ b/src/units.c
@@ -50,6 +50,20 @@ void units_init_cgs(struct unit_system* us) {
us->UnitTemperature_in_cgs = 1.;
}
+/**
+ * @brief Initialises the unit_system structure with SI system
+ *
+ * @param us The unit_system to initialize
+ */
+void units_init_si(struct unit_system* us) {
+
+ us->UnitMass_in_cgs = 1000.;
+ us->UnitLength_in_cgs = 100.;
+ us->UnitTime_in_cgs = 1.;
+ us->UnitCurrent_in_cgs = 1.;
+ us->UnitTemperature_in_cgs = 1.;
+}
+
/**
* @brief Initialise the unit_system with values for the base units.
*
diff --git a/src/units.h b/src/units.h
index da2c209815b07d1d5597a598ee4a61f3132e39db..08b738c5303db8b40dfbe51799d67da8df3936ce 100644
--- a/src/units.h
+++ b/src/units.h
@@ -96,6 +96,7 @@ enum unit_conversion_factor {
};
void units_init_cgs(struct unit_system*);
+void units_init_si(struct unit_system*);
void units_init(struct unit_system* us, double U_M_in_cgs, double U_L_in_cgs,
double U_t_in_cgs, double U_C_in_cgs, double U_T_in_cgs);
void units_init_from_params(struct unit_system*, struct swift_params*,
diff --git a/tests/test27cells.c b/tests/test27cells.c
index ada1b782cfff3866bf26937391007947e9c9a175..f62c169486250ba940c22f21ba1556cca4060c1a 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -217,8 +217,20 @@ void clean_up(struct cell *ci) {
* @brief Initializes all particles field to be ready for a density calculation
*/
void zero_particle_fields(struct cell *c) {
+#ifdef SHADOWFAX_SPH
+ struct hydro_space hs;
+ hs.anchor[0] = 0.;
+ hs.anchor[1] = 0.;
+ hs.anchor[2] = 0.;
+ hs.side[0] = 1.;
+ hs.side[1] = 1.;
+ hs.side[2] = 1.;
+ struct hydro_space *hspointer = &hs;
+#else
+ struct hydro_space *hspointer = NULL;
+#endif
for (int pid = 0; pid < c->count; pid++) {
- hydro_init_part(&c->parts[pid], NULL);
+ hydro_init_part(&c->parts[pid], hspointer);
}
}
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
index 526a7661eaf67a63ca5ce3497c5ff366ffcbed8e..0ca3d0ca9aa2adfc5d6237a668ee93c7c4daba63 100644
--- a/tests/testPotentialPair.c
+++ b/tests/testPotentialPair.c
@@ -95,11 +95,6 @@ int main(int argc, char *argv[]) {
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
-/* Choke on FPEs */
-#ifdef HAVE_FE_ENABLE_EXCEPT
- feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
-#endif
-
/* Initialise a few things to get us going */
/* Non-truncated forces first */