SmoothedMetallicity example done and works

examples/SmoothedMetallicity/getGlass.sh

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

examples/SmoothedMetallicity/makeIC.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@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/>.
+# 
+##############################################################################
+
+import h5py
+import numpy as np
+
+# Generates a swift IC file for the Sedov blast test in a periodic cubic box
+
+# Parameters
+gamma = 5./3.      # Gas adiabatic index
+rho0 = 1.          # Background density
+P0 = 1.e-6         # Background pressure
+low_metal = -6     # Low iron fraction
+high_metal = -5    # high iron fraction
+sigma_metal = 0.1  # relative standard deviation for the metallicities
+Nelem = 9          # Gear: 9, EAGLE: 9
+fileName = "smoothed_metallicity.hdf5"
+
+# ---------------------------------------------------
+glass = h5py.File("glassCube_32.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:, :]
+h = glass["/PartType0/SmoothingLength"][:]
+
+numPart = np.size(h)
+vol = 1.
+
+# Generate extra arrays
+v = np.zeros((numPart, 3))
+ids = np.linspace(1, numPart, numPart)
+m = np.zeros(numPart)
+u = np.zeros(numPart)
+r = np.zeros(numPart)
+Z = np.zeros((numPart, Nelem))
+
+r = np.sqrt((pos[:, 0] - 0.5)**2 + (pos[:, 1] - 0.5)**2 + (pos[:, 2] - 0.5)**2)
+m[:] = rho0 * vol / numPart
+u[:] = P0 / (rho0 * (gamma - 1))
+
+# set metallicities
+select = pos[:, 0] < 0.5
+nber = sum(select)
+sigma = abs(sigma_metal*low_metal)
+Z[select, :] = low_metal + np.random.normal(loc=0., scale=sigma,
+                                            size=(nber, Nelem))
+nber = numPart - nber
+sigma = abs(sigma_metal*high_metal)
+Z[np.logical_not(select), :] = high_metal + np.random.normal(
+    loc=0., scale=sigma, size=(nber, Nelem))
+
+# --------------------------------------------------
+
+# File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 1.]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+# Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+# Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset('Coordinates', data=pos, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset('ElementAbundance', data=Z, dtype='f')
+
+
+file.close()

examples/SmoothedMetallicity/plotSolution.py

+#!/usr/bin/env python
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    Matthieu Schaller (matthieu.schaller@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/>.
+# 
+##############################################################################
+
+# Computes the analytical solution of the 3D Smoothed Metallicity example.
+
+import h5py
+import sys
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+# Parameters
+high_metal = -5    # High metal abundance
+low_metal = -6     # low metal abundance
+sigma_metal = 0.1  # relative standard deviation for Z
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+Nelem = 4           # Number of element
+
+
+# Plot parameters
+params = {
+    'axes.labelsize': 10,
+    'axes.titlesize': 10,
+    'font.size': 12,
+    'legend.fontsize': 12,
+    'xtick.labelsize': 10,
+    'ytick.labelsize': 10,
+    'text.usetex': True,
+    'figure.figsize': (9.90, 6.45),
+    'figure.subplot.left': 0.045,
+    'figure.subplot.right': 0.99,
+    'figure.subplot.bottom': 0.05,
+    'figure.subplot.top': 0.99,
+    'figure.subplot.wspace': 0.15,
+    'figure.subplot.hspace': 0.12,
+    'lines.markersize': 6,
+    'lines.linewidth': 3.,
+    'text.latex.unicode': True
+}
+
+plt.rcParams.update(params)
+plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("smoothed_metallicity_%04d.hdf5" % snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+chemistry = sim["/SubgridScheme"].attrs["Chemistry Model"]
+git = sim["Code"].attrs["Git Revision"]
+
+pos = sim["/PartType0/Coordinates"][:, :]
+d = pos[:, 0] - boxSize / 2
+metal = sim["/PartType0/SmoothedElementAbundance"][:, :]
+h = sim["/PartType0/SmoothingLength"][:]
+h = np.mean(h)
+
+N = 1000
+d_a = np.linspace(-boxSize / 2., boxSize / 2., N)
+
+# Now, work our the solution....
+
+
+def calc_a(d, high_metal, low_metal, std_dev, h):
+    """
+    Compute analytical solution:
+                        ___ High Metallicity
+    Low Metallicity ___/
+
+    where the linear part length is given by the smoothing length.
+
+    Parameters
+    ----------
+
+    d: np.array
+        Position to compute
+    high_metal: float
+        Value on the high metallicity side
+
+    low_metal: float
+        Value on the low metallicity side
+
+    std_dev: float
+        Standard deviation of the noise added
+
+    h: float
+        Mean smoothing length
+    """
+
+    # solution
+    a = np.zeros([len(d), Nelem])
+    # function at +- 1 sigma
+    sigma = np.zeros([len(d), Nelem, 2])
+
+    for i in range(Nelem):
+        # compute low metallicity side
+        s = d < -h
+        a[s, i] = low_metal
+        # compute high metallicity side
+        s = d > h
+        a[s, i] = high_metal
+
+        # compute non constant parts
+        m = (high_metal - low_metal) / (2.0 * h)
+        c = (high_metal + low_metal) / 2.0
+        # compute left linear part
+        s = d < - boxSize / 2.0 + h
+        a[s, i] = - m * (d[s] + boxSize / 2.0) + c
+        # compute middle linear part
+        s = np.logical_and(d >= -h, d <= h)
+        a[s, i] = m * d[s] + c
+        # compute right linear part
+        s = d > boxSize / 2.0 - h
+        a[s, i] = - m * (d[s] - boxSize / 2.0) + c
+
+    sigma[:, :, 0] = a * (1 + std_dev)
+    sigma[:, :, 1] = a * (1 - std_dev)
+    return a, sigma
+
+
+# compute solution
+sol, sig = calc_a(d_a, high_metal, low_metal, sigma_metal, h)
+
+# Plot the interesting quantities
+plt.figure()
+
+# Metallicity --------------------------------
+e = 0
+plt.subplot(221)
+plt.plot(d, metal[:, e], '.', color='r', ms=0.5, alpha=0.2)
+plt.plot(d_a, sol[:, e], '--', color='b', alpha=0.8, lw=1.2)
+plt.fill_between(d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b",
+                 interpolate=True, alpha=0.5)
+plt.xlabel("${\\rm{Distance}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Metallicity}}~Z$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(low_metal-1, high_metal+1)
+
+# Metallicity --------------------------------
+e = 1
+plt.subplot(222)
+plt.plot(d, metal[:, e], '.', color='r', ms=0.5, alpha=0.2)
+plt.plot(d_a, sol[:, e], '--', color='b', alpha=0.8, lw=1.2)
+plt.fill_between(d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b",
+                 interpolate=True, alpha=0.5)
+plt.xlabel("${\\rm{Distance}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Metallicity}}~Z$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(low_metal-1, high_metal+1)
+
+# Metallicity --------------------------------
+e = 2
+plt.subplot(223)
+plt.plot(d, metal[:, e], '.', color='r', ms=0.5, alpha=0.2)
+plt.plot(d_a, sol[:, e], '--', color='b', alpha=0.8, lw=1.2)
+plt.fill_between(d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b",
+                 interpolate=True, alpha=0.5)
+plt.xlabel("${\\rm{Distance}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Metallicity}}~Z$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(low_metal-1, high_metal+1)
+
+# Information -------------------------------------
+plt.subplot(224, frameon=False)
+
+plt.text(-0.49, 0.9, "Smoothed Metallicity in 3D at $t=%.2f$" % time,
+         fontsize=10)
+plt.plot([-0.49, 0.1], [0.82, 0.82], 'k-', lw=1)
+plt.text(-0.49, 0.7, "$\\textsc{Swift}$ %s" % git, fontsize=10)
+plt.text(-0.49, 0.6, scheme, fontsize=10)
+plt.text(-0.49, 0.5, kernel, fontsize=10)
+plt.text(-0.49, 0.4, chemistry + "'s Chemistry", fontsize=10)
+plt.text(-0.49, 0.3, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+         fontsize=10)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.savefig("SmoothedMetallicity.png", dpi=200)

examples/SmoothedMetallicity/run.sh

+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+    echo "Fetching initial glass file for the SmoothedMetallicity example..."
+    ./getGlass.sh
+fi
+if [ ! -e smoothed_metallicity.hdf5 ]
+then
+    echo "Generating initial conditions for the SmoothedMetallicity example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 4 smoothed_metallicity.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 4

examples/SmoothedMetallicity/smoothed_metallicity.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   2e-4  # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # 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:            smoothed_metallicity # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          5e-5  # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-5 # Time between statistics output
+
+# 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.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./smoothed_metallicity.hdf5          # The file to read
+  h_scaling:           3.33
+

src/cell.c

@@ -49,6 +49,7 @@
 /* Local headers. */
 #include "active.h"
 #include "atomic.h"
+#include "chemistry.h"
 #include "drift.h"
 #include "engine.h"
 #include "error.h"
@@ -2423,6 +2424,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
       /* Get ready for a density calculation */
       if (part_is_active(p, e)) {
         hydro_init_part(p, &e->s->hs);
+	chemistry_init_part(p, e->chemistry);
       }
     }
 

src/chemistry/gear/chemistry.h

@@ -52,20 +52,6 @@ chemistry_get_element_name(enum chemistry_element elem) {
   return chemistry_element_names[elem];
 }
 
-/**
- * @brief Sets the chemistry properties of the (x-)particles to a valid start
- * state.
- *
- * Nothing to do here.
- *
- * @param p Pointer to the particle data.
- * @param xp Pointer to the extended particle data.
- * @param data The global chemistry information.
- */
-__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
-    const struct part* restrict p, struct xpart* restrict xp,
-    const struct chemistry_data* data) {}
-
 /**
  * @brief Initialises the chemistry properties.
  *
@@ -87,7 +73,7 @@ static INLINE void chemistry_init_backend(
  */
 static INLINE void chemistry_print_backend(const struct chemistry_data* data) {
 
-  message("Chemistry function is 'gear'.");
+  message("Chemistry function is 'Gear'.");
 }
 
 
@@ -105,8 +91,8 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
 
   struct chemistry_part_data *cpd = &p->chemistry_data;
 
-  for(size_t i=0; i < chemistry_element_count; i++) {
-    cpd->smoothed_metal_mass_fraction[i] = 0;
+  for(int i=0; i < chemistry_element_count; i++) {
+    cpd->smoothed_metal_mass_fraction[i] = 0.f;
   }
 }
 
@@ -121,24 +107,41 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
  *
  * @param p The particle to act upon
  */
-__attribute__((always_inline)) INLINE static void chemistry_end(
+__attribute__((always_inline)) INLINE static void chemistry_end_density(
     struct part *restrict p, const struct chemistry_data *cd) {
 
   /* Some smoothing length multiples. */
   const float h = p->h;
   const float h_inv = 1.0f / h;                       /* 1/h */
-  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
+  const float factor = pow_dimension(h_inv) / p->rho; /* 1 / h^d * rho */
   const float m = p->mass;
 
   struct chemistry_part_data *cpd = &p->chemistry_data;
 
-  for(size_t i=0; i < chemistry_element_count; i++) {
+  for(int i=0; i < chemistry_element_count; i++) {
     /* Final operation on the density (add self-contribution). */
     cpd->smoothed_metal_mass_fraction[i] += m * cpd->metal_mass_fraction[i]* kernel_root;
 
     /* Finish the calculation by inserting the missing h-factors */
-    cpd->smoothed_metal_mass_fraction[i] *= h_inv_dim;
+    cpd->smoothed_metal_mass_fraction[i] *= factor;
   }
 }
 
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ * @param data The global chemistry information.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
+    struct part* restrict p, struct xpart* restrict xp,
+    const struct chemistry_data* data) {
+  chemistry_init_part(p, data);
+}
+
 #endif /* SWIFT_CHEMISTRY_GEAR_H */

src/chemistry/gear/chemistry_iact.h

@@ -60,9 +60,9 @@ const struct chemistry_data *chem_data) {
   kernel_deval(uj, &wj, &wj_dx);
 
   /* Compute contribution to the smooth metallicity */
-  for(size_t i=0; i < chemistry_element_count; i++) {
-    chi->smoothed_metal_mass_fraction[i] += mj * chj->smoothed_metal_mass_fraction[i] * wi;
-    chj->smoothed_metal_mass_fraction[i] += mi * chi->smoothed_metal_mass_fraction[i] * wj;
+  for(int i=0; i < chemistry_element_count; i++) {
+    chi->smoothed_metal_mass_fraction[i] += mj * chj->metal_mass_fraction[i] * wi;
+    chj->smoothed_metal_mass_fraction[i] += mi * chi->metal_mass_fraction[i] * wj;
   }
 }
 
@@ -89,9 +89,10 @@ const struct chemistry_data *chem_data) {
   kernel_deval(ui, &wi, &wi_dx);
 
   /* Compute contribution to the smooth metallicity */
-  for(size_t i=0; i < chemistry_element_count; i++) {
-    chi->smoothed_metal_mass_fraction[i] += mj * chj->smoothed_metal_mass_fraction[i] * wi;
+  for(int i=0; i < chemistry_element_count; i++) {
+    chi->smoothed_metal_mass_fraction[i] += mj * chj->metal_mass_fraction[i] * wi;
   }
+
 }
 
 

src/chemistry/gear/chemistry_io.h

@@ -33,15 +33,11 @@
  */
 int chemistry_read_particles(struct part* parts, struct io_props* list) {
 
-  for(size_t i=0; i < chemistry_element_count; i++) {
-    /* List what we want to read */
-    char buffer[20];
-    strcpy(buffer, chemistry_get_element_name(i));
-    list[i] =
-      io_make_input_field(buffer, FLOAT, 1, OPTIONAL, UNIT_CONV_NO_UNITS,
-			  parts, chemistry_data.metal_mass_fraction[i]);
-  }
-
+  /* List what we want to read */
+  list[0] =
+    io_make_input_field("ElementAbundance", FLOAT, chemistry_element_count, OPTIONAL,
+			UNIT_CONV_NO_UNITS,
+			parts, chemistry_data.metal_mass_fraction);
   return 1;
 }
 
@@ -59,7 +55,12 @@ int chemistry_write_particles(const struct part* parts, struct io_props* list) {
   list[0] = io_make_output_field("SmoothedElementAbundance", FLOAT, chemistry_element_count,
 				 UNIT_CONV_NO_UNITS,
 				 parts, chemistry_data.smoothed_metal_mass_fraction);
-  return 1;
+
+  list[1] = io_make_output_field("ElementAbundance", FLOAT, chemistry_element_count,
+				 UNIT_CONV_NO_UNITS,
+				 parts, chemistry_data.metal_mass_fraction);
+
+  return 2;
 }
 
 /**

src/runner.c

@@ -701,7 +701,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
           /* Finish the density calculation */
           hydro_end_density(p);
-	  chemistry_end(p, e->chemistry);
+	  chemistry_end_density(p, e->chemistry);
 
           /* Compute one step of the Newton-Raphson scheme */
           const float n_sum = p->density.wcount * h_old_dim;