Merge branch 'planetary_docs' into 'master'

Planetary documentation

See merge request !733

.gitignore

@@ -27,6 +27,7 @@ examples/swift_mpi
 examples/*/*.xmf
 examples/*/*.h5
 examples/*/*.png
+examples/*/*.pdf
 examples/*/*.mp4
 examples/*/*.txt
 examples/*/dependency_graph.csv

doc/RTD/source/EquationOfState/index.rst

-.. Equation of State
+.. Equations of State
    Loic Hausammann, 6th April 2018
+   Jacob Kegerreis, 3rd February 2019
 
 .. _equation_of_state:
 
-Equation of State
-=================
+Equations of State
+==================
 
-Currently (if the documentation was well updated), we have two different
-equation of states implemented: ideal gas and isothermal.  They describe the
-relations between our main thermodynamical variables: the internal energy
-(\\(u\\)), the density (\\(\\rho\\)), the entropy (\\(A\\)) and the pressure
-(\\(P\\)).
+Currently (if the documentation was well updated), we have two different gas
+equations of state (EoS) implemented: ideal and isothermal; as well as a variety  
+of EoS for "planetary" materials. 
+The EoS describe the relations between our main thermodynamical variables: 
+the internal energy (\\(u\\)), the density (\\(\\rho\\)), the entropy (\\(A\\)) 
+and the pressure (\\(P\\)).
 
-Equations
----------
+Gas EoS
+-------
 
 In the following section, the variables not yet defined are: \\(\\gamma\\) for
 the adiabatic index and \\( c_s \\) for the speed of sound.
@@ -37,12 +39,55 @@ the adiabatic index and \\( c_s \\) for the speed of sound.
    "\\( c_s\\)", "", "\\(\\sqrt{ u \\gamma \\left( \\gamma - 1 \\right) } \\)", ""
 
 
+
+Planetary EoS
+-------------
+Configuring SWIFT with the ``--with-equation-of-state=planetary`` and 
+``--with-hydro=planetary`` options enables the use of multiple EoS.
+Every SPH particle then requires and carries the additional ``MaterialID`` flag 
+from the initial conditions file. This flag indicates the particle's material 
+and which EoS it should use. 
+
+So far, we have implemented several Tillotson, SESAME, and Hubbard \& MacFarlane 
+(1980) materials, with more on their way.
+The material's ID is set by a base type ID (multiplied by 100), plus a minor 
+type:
+
++ Tillotson (Melosh, 2007): ``1``
+    + Iron: ``100``
+    + Granite: ``101``
+    + Water: ``102``
++ Hubbard \& MacFarlane (1980): ``2``
+    + Hydrogen-helium atmosphere: ``200``
+    + Ice H20-CH4-NH3 mix: ``201``
+    + Rock SiO2-MgO-FeS-FeO mix: ``202``
++ SESAME (and similar): ``3``
+    + Iron (2140): ``300``
+    + Basalt (7530): ``301``
+    + Water (7154): ``302``
+    + Senft \& Stewart (2008) water (in a SESAME-style table): ``303``
+
+Unlike the EoS for an ideal or isothermal gas, these more complicated materials 
+do not always include transformations between the internal energy, 
+temperature, and entropy. At the moment, we have only implemented 
+\\(P(\\rho, u)\\) and \\(c_s(\\rho, u)\\). 
+This is sufficient for the simple :ref:`planetary_sph` hydrodynamics scheme, 
+but makes these materials currently incompatible with other entropy-based 
+schemes.
+
+The Tillotson sound speed was derived using 
+\\(c_s^2 = \\left. \\dfrac{\\partial P}{\\partial \\rho} \\right|_S \\)
+as described in Kegerreis et al. (2019).
+The table files for the HM80 and SESAME-style EoS can be downloaded using 
+the ``examples/EoSTables/get_eos_tables.sh`` script.
+
+
 How to Implement a New Equation of State
 ----------------------------------------
 
 See :ref:`new_option` for a full list of required changes.
 
-You will need to provide a ``equation_of_state.h`` file containing: the
+You will need to provide an ``equation_of_state.h`` file containing: the
 definition of ``eos_parameters``, IO functions and transformations between the
 different variables: \\(u(\\rho, A)\\), \\(u(\\rho, P)\\), \\(P(\\rho,A)\\),
 \\(P(\\rho, u)\\), \\(A(\\rho, P)\\), \\(A(\\rho, u)\\), \\(c_s(\\rho, A)\\),

doc/RTD/source/HydroSchemes/index.rst

@@ -15,6 +15,7 @@ schemes available in SWIFT, as well as how to implement your own.
 
    traditional_sph
    minimal_sph
+   planetary
    hopkins_sph
    gizmo
    adding_your_own

doc/RTD/source/HydroSchemes/planetary.rst

+.. Planetary SPH
+    Jacob Kegerreis, 3rd February 2019
+
+.. _planetary_sph:
+
+Planetary (Density-Energy, Multi-Material) SPH
+==============================================
+
+.. toctree::
+   :maxdepth: 2
+   :hidden:
+   :caption: Contents:
+
+This scheme is the same as the Minimal SPH scheme but also allows multiple 
+materials, meaning that different SPH particles can be assigned different 
+:ref:`equation_of_state` (EoS).
+
+To use the planetary scheme and the corresponding planetary EoS, use 
+
+.. code-block:: bash
+
+    ./configure --with-hydro=planetary --with-equation-of-state=planetary
+
+Every SPH particle then requires and carries the additional ``MaterialID`` flag 
+from the initial conditions file. This flag indicates the particle's material 
+and which EoS it should use. 
\ No newline at end of file

doc/RTD/source/conf.py

@@ -18,7 +18,7 @@
 
 # -- Project information -----------------------------------------------------
 
-project = 'SWIFT: SPH WIth Fine-grained inter-dependent Tasking'
+project = 'SWIFT: SPH With Inter-dependent Fine-grained Tasking'
 copyright = '2018, SWIFT Collaboration'
 author = 'SWIFT Team'
 

doc/RTD/source/index.rst

@@ -3,7 +3,7 @@
    You can adapt this file completely to your liking, but it should at least
    contain the root `toctree` directive.
 
-Welcome to SWIFT: SPH WIth Fine-grained inter-dependent Tasking's documentation!
+Welcome to SWIFT: SPH With Inter-dependent Fine-grained Tasking's documentation!
 ================================================================================
 
 Want to get started using SWIFT? Check out the on-boarding guide available

examples/EoSTables/get_eos_tables.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/HM80_HHe.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/HM80_ice.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/HM80_rock.txt
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_basalt_7530.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_iron_2140.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_water_7154.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SS08_water.txt

examples/GiantImpactUranus/README.md

+An example planetary simulation of a giant impact onto the young Uranus with 
+~10^6 SPH particles, as described in Kegerreis et al. (2018), ApJ, 861, 52.

examples/GiantImpactUranus/configuration.yml

+with-hydro:             planetary 
+with-equation-of-state: planetary

examples/GiantImpactUranus/get_init_cond.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/uranus_1e6.hdf5

examples/GiantImpactUranus/output_list.txt

+# Time
+4000
+9000
+14000
+20000
+30000
+40000
\ No newline at end of file

examples/GiantImpactUranus/plot_solution.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2019 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
+# Plot the snapshots from the example giant impact on Uranus, showing the 
+# particles in a thin slice near z=0, coloured by their material, similarish
+# (but not identical) to Fig. 2 in Kegerreis et al. (2018).
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+# import swiftsimio as sw
+import h5py
+
+font_size   = 20
+params      = { 
+    'axes.labelsize'    : font_size,
+    'font.size'         : font_size,
+    'xtick.labelsize'   : font_size,
+    'ytick.labelsize'   : font_size,
+    'text.usetex'       : True,
+    'font.family'       : 'serif',
+    }
+matplotlib.rcParams.update(params)
+
+# Snapshot output times
+output_list = [4000, 9000, 14000, 20000, 30000, 40000]
+
+# Material IDs ( = type_id * type_factor + unit_id )
+type_factor = 100
+type_HM80   = 2
+id_body     = 10000
+# Name and ID
+Di_mat_id   = {
+    'HM80_HHe'      : type_HM80 * type_factor,      # Hydrogen-helium atmosphere
+    'HM80_ice'      : type_HM80 * type_factor + 1,  # H20-CH4-NH3 ice mix
+    'HM80_ice_2'    : type_HM80 * type_factor + 1 + id_body,
+    'HM80_rock'     : type_HM80 * type_factor + 2,  # SiO2-MgO-FeS-FeO rock mix
+    'HM80_rock_2'   : type_HM80 * type_factor + 2 + id_body,
+    }
+# ID and colour
+Di_id_colour    = {
+    Di_mat_id['HM80_HHe']       : '#33DDFF',
+    Di_mat_id['HM80_ice']       : 'lightsteelblue',
+    Di_mat_id['HM80_ice_2']     : '#A080D0',
+    Di_mat_id['HM80_rock']      : 'slategrey',
+    Di_mat_id['HM80_rock_2']    : '#706050',
+    }
+   
+def get_snapshot_slice(snapshot):
+    """ Load and select the particles to plot. """
+    # Load particle data
+    # data    = load("uranus_1e6_%06d.hdf5" % snapshot)
+    # id      = data.gas.particle_ids
+    # pos     = data.gas.coordinates
+    # mat_id  = data.gas.material
+    with h5py.File("uranus_1e6_%06d.hdf5" % snapshot, 'r') as f:
+        id      = f['PartType0/ParticleIDs'].value
+        pos     = (f['PartType0/Coordinates'].value
+                   - 0.5 * f['Header'].attrs['BoxSize'])
+        mat_id  = f['PartType0/MaterialID'].value
+                   
+    # Edit the material ID of particles in the impactor
+    num_in_target   = 869104
+    sel_id          = np.where(num_in_target < id)[0]
+    mat_id[sel_id]  += id_body
+    
+    # Select particles in a thin slice around z=0
+    z_min   = -0.1
+    z_max   = 0.1
+    sel_z   = np.where((z_min < pos[:, 2]) & (pos[:, 2] < z_max))[0]
+    pos     = pos[sel_z]
+    mat_id  = mat_id[sel_z]
+    
+    return pos, mat_id
+
+def plot_snapshot_slice(pos, mat_id):
+    """ Plot the particles, coloured by their material. """
+    colour  = np.empty(len(pos), dtype=object)
+    for id, c in Di_id_colour.items():
+        sel_c           = np.where(mat_id == id)[0]
+        colour[sel_c]   = c
+
+    ax.scatter(pos[:, 0], pos[:, 1], c=colour, edgecolors='none', marker='.', 
+               s=10, alpha=0.5, zorder=0)
+
+# Set up the figure
+fig     = plt.figure(figsize=(12, 8))
+gs      = matplotlib.gridspec.GridSpec(2, 3)
+axes    = [plt.subplot(gs[i_y, i_x]) for i_y in range(2) for i_x in range(3)]
+
+# Plot each snapshot
+for i_ax, ax in enumerate(axes):
+    plt.sca(ax)
+    ax.set_rasterization_zorder(1)
+    
+    # Load and select the particles to plot
+    pos, mat_id = get_snapshot_slice(output_list[i_ax])
+    
+    # Plot the particles, coloured by their material
+    plot_snapshot_slice(pos, mat_id)
+    
+    # Axes etc.
+    ax.set_aspect('equal')
+    ax.set_facecolor('k')
+    
+    ax.set_xlim(-13, 13)
+    ax.set_ylim(-13, 13)
+
+    if i_ax in [0, 3]:
+        ax.set_ylabel(r"y Postion $(R_\oplus)$")
+    else: 
+        ax.set_yticklabels([])
+    if 2 < i_ax:
+        ax.set_xlabel(r"x Postion $(R_\oplus)$")
+    else: 
+        ax.set_xticklabels([])
+        
+    # Corner time labels 
+    x   = ax.get_xlim()[0] + 0.04 * (ax.get_xlim()[1] - ax.get_xlim()[0])
+    y   = ax.get_ylim()[0] + 0.89 * (ax.get_ylim()[1] - ax.get_ylim()[0])
+    ax.text(x, y, "%.1f h" % (output_list[i_ax] / 60**2), color='w')
+
+plt.subplots_adjust(wspace=0, hspace=0)
+plt.tight_layout()
+
+# Save
+plt.savefig("uranus_1e6.pdf", dpi=200)
\ No newline at end of file

examples/GiantImpactUranus/run.sh

+#!/bin/bash
+
+# Get the initial conditions if they are not present.
+if [ ! -e uranus_1e6.hdf5 ]
+then
+    echo "Fetching initial conditions file for the Uranus impact example..."
+    ./get_init_cond.sh
+fi
+
+# Get the EoS tables if they are not present.
+cd ../EoSTables
+if [ ! -e HM80_HHe.txt ] || [ ! -e HM80_ice.txt ] || [ ! -e HM80_rock.txt ] 
+then
+    echo "Fetching equations of state tables for the Uranus impact example..."
+    ./get_eos_tables.sh
+fi
+cd ../GiantImpactUranus
+
+# Run SWIFT
+../swift -s -G -t 8 uranus_1e6.yml 2>&1 | tee output.log
+
+# Plot the solution
+python3 plot_solution.py

examples/GiantImpactUranus/system.yml

+input:
+  - uranus_1e6.yml
+  - output_list.txt
+
+output:
+  - uranus_1e6.pdf
+
+swift_parameters:   -s -G
+swift_threads:      8

examples/GiantImpactUranus/uranus_1e6.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+    UnitMass_in_cgs:        5.9724e27   # Grams
+    UnitLength_in_cgs:      6.371e8     # Centimeters
+    UnitVelocity_in_cgs:    6.371e8     # Centimeters per second
+    UnitCurrent_in_cgs:     1           # Amperes
+    UnitTemp_in_cgs:        1           # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:      
+    file_name:  uranus_1e6.hdf5         # The initial conditions file to read
+    periodic:   0                       # Are we running with periodic ICs?
+
+# Parameters governing the time integration
+TimeIntegration:
+    time_begin:     0                   # The starting time of the simulation (in internal units).
+    time_end:       40000               # The end time of the simulation (in internal units).
+    dt_min:         0.0001              # The minimal time-step size of the simulation (in internal units).
+    dt_max:         100                 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+    basename:           uranus_1e6      # Common part of the name of output files
+    time_first:         0               # Time of the first output (in internal units)
+    delta_time:         1000            # Time difference between consecutive outputs (in internal units)
+    int_time_label_on:  1               # Enable to label the snapshots using the time rounded to an integer (in internal units)
+    output_list_on:     1               # Enable the output list
+    output_list:        output_list.txt # File containing the output times (see documentation in "Parameter File" section)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+    time_first: 0                       # Time of the first output (in internal units)
+    delta_time: 1000                    # Time between statistics output
+
+# Parameters controlling restarts
+Restarts:
+    enable: 0                           # Whether to enable dumping restarts at fixed intervals.
+
+# 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).
+    delta_neighbours:   0.1             # The tolerance for the targetted number of neighbours.
+    CFL_condition:      0.2             # Courant-Friedrich-Levy condition for time integration.
+    h_max:              0.5             # Maximal allowed smoothing length (in internal units).
+    viscosity_alpha:    1.5             # Override for the initial value of the artificial viscosity.
+
+# 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.003       # Comoving softening length (in internal units).
+    max_physical_softening: 0.003       # Physical softening length (in internal units).
+
+# Parameters for the task scheduling
+Scheduler:
+    max_top_level_cells:    64          # Maximal number of top-level cells in any dimension. The nu
+
+# Parameters related to the equation of state
+EoS:
+    planetary_use_HM80:             1   # Whether to initialise the Hubbard & MacFarlane (1980) EOS
+    planetary_HM80_HHe_table_file:  ../EoSTables/HM80_HHe.txt
+    planetary_HM80_ice_table_file:  ../EoSTables/HM80_ice.txt
+    planetary_HM80_rock_table_file: ../EoSTables/HM80_rock.txt

examples/parameter_example.yml

@@ -184,13 +184,13 @@ EoS:
   planetary_use_ANEOS:  0   # (Optional) Whether to prepare the ANEOS EOS
   planetary_use_SESAME: 0   # (Optional) Whether to prepare the SESAME EOS
                             # (Optional) Table file paths
-  planetary_HM80_HHe_table_file:        ./equation_of_state/planetary_HM80_HHe.txt
-  planetary_HM80_ice_table_file:        ./equation_of_state/planetary_HM80_ice.txt
-  planetary_HM80_rock_table_file:       ./equation_of_state/planetary_HM80_rock.txt
-  planetary_SESAME_iron_table_file:     ./equation_of_state/planetary_SESAME_iron_2140.txt
-  planetary_SESAME_basalt_table_file:   ./equation_of_state/planetary_SESAME_basalt_7530.txt
-  planetary_SESAME_water_table_file:    ./equation_of_state/planetary_SESAME_water_7154.txt
-  planetary_SS08_water_table_file:      ./equation_of_state/planetary_SS08_water.txt
+  planetary_HM80_HHe_table_file:        ./EoSTables/planetary_HM80_HHe.txt
+  planetary_HM80_ice_table_file:        ./EoSTables/planetary_HM80_ice.txt
+  planetary_HM80_rock_table_file:       ./EoSTables/planetary_HM80_rock.txt
+  planetary_SESAME_iron_table_file:     ./EoSTables/planetary_SESAME_iron_2140.txt
+  planetary_SESAME_basalt_table_file:   ./EoSTables/planetary_SESAME_basalt_7530.txt
+  planetary_SESAME_water_table_file:    ./EoSTables/planetary_SESAME_water_7154.txt
+  planetary_SS08_water_table_file:      ./EoSTables/planetary_SS08_water.txt
 
 # Parameters related to external potentials --------------------------------------------