Commit 1823a387 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'smooth_metal' into 'master'

Smoothed metallicity

See merge request !507
parents 3cc0dba0 0a29be20
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
###############################################################################
# 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 Smoothed Metallicity test in a periodic cubic box
# Parameters
gamma = 5./3. # Gas adiabatic index
rho0 = 1. # Background density
P0 = 1.e-6 # Background pressure
Nelem = 9 # Gear: 9, EAGLE: 9
low_metal = -6 # Low iron fraction
high_metal = -5 # high iron fraction
sigma_metal = 0.1 # relative standard deviation for the metallicities
fileName = "smoothed_metallicity.hdf5"
# shift all metals in order to obtain nicer plots
low_metal = [low_metal] * Nelem + np.linspace(0, 3, Nelem)
high_metal = [high_metal] * Nelem + np.linspace(0, 3, Nelem)
# ---------------------------------------------------
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)
Z[select, :] = low_metal * (1 + np.random.normal(loc=0., scale=sigma_metal,
size=(nber, Nelem)))
nber = numPart - nber
Z[np.logical_not(select), :] = high_metal * (1 + np.random.normal(
loc=0., scale=sigma_metal, 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()
#!/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
low_metal = -6 # low metal abundance
high_metal = -5 # High metal abundance
sigma_metal = 0.1 # relative standard deviation for Z
Nelem = 9
# shift all metals in order to obtain nicer plots
low_metal = [low_metal] * Nelem + np.linspace(0, 3, Nelem)
high_metal = [high_metal] * Nelem + np.linspace(0, 3, Nelem)
# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------
# 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
smooth_metal = sim["/PartType0/SmoothedElementAbundance"][:, :]
metal = sim["/PartType0/ElementAbundance"][:, :]
h = sim["/PartType0/SmoothingLength"][:]
h = np.mean(h)
if (Nelem != metal.shape[1]):
print("Unexpected number of element, please check makeIC.py"
" and plotSolution.py")
exit(1)
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[i]
# compute high metallicity side
s = d > h
a[s, i] = high_metal[i]
# compute non constant parts
m = (high_metal[i] - low_metal[i]) / (2.0 * h)
c = (high_metal[i] + low_metal[i]) / 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 --------------------------------
plt.subplot(221)
for e in range(Nelem):
plt.plot(metal[:, e], smooth_metal[:, e], '.', ms=0.5, alpha=0.2)
xmin, xmax = metal.min(), metal.max()
ymin, ymax = smooth_metal.min(), smooth_metal.max()
x = max(xmin, ymin)
y = min(xmax, ymax)
plt.plot([x, y], [x, y], "--k", lw=1.0)
plt.xlabel("${\\rm{Metallicity}}~Z_\\textrm{part}$", labelpad=0)
plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
# Metallicity --------------------------------
e = 0
plt.subplot(223)
plt.plot(d, smooth_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{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
plt.xlim(-0.5, 0.5)
plt.ylim(low_metal[e]-1, high_metal[e]+1)
# Information -------------------------------------
plt.subplot(222, 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.tight_layout()
plt.savefig("SmoothedMetallicity.png", dpi=200)
#!/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 -n 1 -s -t 4 smoothed_metallicity.yml 2>&1 | tee output.log
# Plot the solution
python plotSolution.py 1
# 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
......@@ -122,12 +122,15 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
chemistry/none/chemistry.h \
chemistry/none/chemistry_io.h \
chemistry/none/chemistry_struct.h \
chemistry/none/chemistry_iact.h \
chemistry/gear/chemistry.h \
chemistry/gear/chemistry_io.h \
chemistry/gear/chemistry_struct.h \
chemistry/gear/chemistry_iact.h \
chemistry/EAGLE/chemistry.h \
chemistry/EAGLE/chemistry_io.h \
chemistry/EAGLE/chemistry_struct.h
chemistry/EAGLE/chemistry_struct.h\
chemistry/EAGLE/chemistry_iact.h
# Sources and flags for regular library
......
......@@ -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);
}
}
......
......@@ -31,10 +31,13 @@
/* Import the right chemistry definition */
#if defined(CHEMISTRY_NONE)
#include "./chemistry/none/chemistry.h"
#include "./chemistry/none/chemistry_iact.h"
#elif defined(CHEMISTRY_GEAR)
#include "./chemistry/gear/chemistry.h"
#include "./chemistry/gear/chemistry_iact.h"
#elif defined(CHEMISTRY_EAGLE)
#include "./chemistry/EAGLE/chemistry.h"
#include "./chemistry/EAGLE/chemistry_iact.h"
#else
#error "Invalid choice of chemistry function."
#endif
......
......@@ -50,6 +50,32 @@ chemistry_get_element_name(enum chemistry_element elem) {
return chemistry_element_names[elem];
}
/**
* @brief Prepares a particle for the smooth metal calculation.
*
* Zeroes all the relevant arrays in preparation for the sums taking place in
* the various smooth metallicity tasks
*
* @param p The particle to act upon
* @param cd #chemistry_data containing chemistry informations.
*/
__attribute__((always_inline)) INLINE static void chemistry_init_part(
struct part* restrict p, const struct chemistry_data* cd) {}
/**
* @brief Finishes the smooth metal calculation.
*
* Multiplies the smoothed metallicity and number of neighbours by the
* appropiate constants and add the self-contribution term.
*
* This function requires the #hydro_end_density to have been called.
*
* @param p The particle to act upon.
* @param cd #chemistry_data containing chemistry informations.
*/
__attribute__((always_inline)) INLINE static void chemistry_end_density(
struct part* restrict p, const struct chemistry_data* cd) {}
/**
* @brief Sets the chemistry properties of the (x-)particles to a valid start
* state.
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2018 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/>.
*
******************************************************************************/
#ifndef SWIFT_EAGLE_CHEMISTRY_IACT_H
#define SWIFT_EAGLE_CHEMISTRY_IACT_H
/**
* @file EAGLE/chemistry_iact.h
* @brief Smooth metal interaction functions following the EAGLE model.
*/
#include "chemistry_struct.h"
/**
* @brief Do chemistry computation after the runner_iact_density (symmetric
* version)
*
* @param r2 Distance squared between particles
* @param dx Distance between particles
* @param hi Smoothing length of i
* @param hj Smoothing length of j
* @param pi #part i
* @param pj #part j
*/
__attribute__((always_inline)) INLINE static void runner_iact_chemistry(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
}
/**
* @brief Do chemistry computation after the runner_iact_density (non-symmetric
* version)
*
* @param r2 Distance squared between particles
* @param dx Distance between particles
* @param hi Smoothing length of i
* @param hj Smoothing length of j
* @param pi #part i
* @param pj #part j
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
float r2, float *dx, float hi, float hj, struct part *pi,
const struct part *pj) {}
#endif /* SWIFT_EAGLE_CHEMISTRY_IACT_H */
......@@ -38,18 +38,17 @@
#include "units.h"
/**
* @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.
* @brief Return a string containing the name of a given #chemistry_element.
*/
__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
const struct part* restrict p, struct xpart* restrict xp,
const struct chemistry_data* data) {}
__attribute__((always_inline)) INLINE static const char*
chemistry_get_element_name(enum chemistry_element elem) {
static const char* chemistry_element_names[chemistry_element_count] = {
"Oxygen", "Magnesium", "Sulfur", "Iron", "Zinc",
"Strontium", "Yttrium", "Barium", "Europium"};
return chemistry_element_names[elem];
}
/**
* @brief Initialises the chemistry properties.
......@@ -72,7 +71,74 @@ 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'.");
}
/**
* @brief Prepares a particle for the smooth metal calculation.
*
* Zeroes all the relevant arrays in preparation for the sums taking place in
* the various smooth metallicity tasks
*
* @param p The particle to act upon
* @param cd #chemistry_data containing chemistry informations.
*/
__attribute__((always_inline)) INLINE static void chemistry_init_part(
struct part* restrict p, const struct chemistry_data* cd) {
struct chemistry_part_data* cpd = &p->chemistry_data;
for (int i = 0; i < chemistry_element_count; i++) {
cpd->smoothed_metal_mass_fraction[i] = 0.f;
}
}
/**
* @brief Finishes the smooth metal calculation.
*
* Multiplies the smoothed metallicity and number of neighbours by the
* appropiate constants and add the self-contribution term.
*
* This function requires the #hydro_end_density to have been called.
*
* @param p The particle to act upon.
* @param cd #chemistry_data containing chemistry informations.
*/
__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 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 (int i = 0; i < chemistry_element_count; i++) {
/* Final operation on the density (add self-contribution). */
cpd->smoothed_metal_mass_fraction[i] +=