From a01df06babcdf16e97917043a30139c799b64810 Mon Sep 17 00:00:00 2001 From: Jonathan Davies <j.j.davies@ljmu.ac.uk> Date: Fri, 10 Jan 2025 12:32:01 +0000 Subject: [PATCH] Smoothing lengths for sinks --- configure.ac | 17 + doc/RTD/source/SubgridModels/Basic/index.rst | 19 + .../SubgridModels/GEAR/sinks/params.rst | 16 +- .../GEAR/params.yml | 3 +- .../IsolatedGalaxy_sink/isolated_galaxy.yml | 3 +- .../SinkParticles/HomogeneousBox/params.yml | 3 +- .../HomogeneousBoxSinkParticles/params.yml | 3 +- .../SinkParticles/PlummerSphere/params.yml | 3 +- examples/SinkParticles/SingleSink/README | 31 + .../SinkParticles/SingleSink/make_sink_ic.py | 234 +++++++ .../SingleSink/make_test_plots.py | 123 ++++ examples/SinkParticles/SingleSink/params.yml | 61 ++ examples/parameter_example.yml | 6 +- src/Makefile.am | 4 +- src/cell.h | 8 +- src/cell_convert_part.c | 2 +- src/cell_drift.c | 20 +- src/cell_pack.c | 4 +- src/cell_sinks.h | 13 +- src/cell_unskip.c | 4 +- src/debug.c | 22 +- src/engine.c | 26 +- src/part.h | 2 + src/proxy.c | 4 +- src/random.h | 1 + src/runner_doiact_functions_hydro.h | 90 +-- src/runner_doiact_functions_sinks.h | 339 ++++++++- src/runner_doiact_sinks.h | 33 + src/runner_ghost.c | 376 +++++++++- src/sink.c | 341 ++++++++++ src/sink.h | 9 + src/sink/Basic/sink.h | 644 ++++++++++++++++++ src/sink/Basic/sink_debug.h | 29 + src/sink/Basic/sink_iact.h | 373 ++++++++++ src/sink/Basic/sink_io.h | 219 ++++++ src/sink/Basic/sink_part.h | 168 +++++ src/sink/Basic/sink_properties.h | 148 ++++ src/sink/Basic/sink_struct.h | 115 ++++ src/sink/Default/sink.h | 54 +- src/sink/Default/sink_iact.h | 42 +- src/sink/Default/sink_io.h | 10 +- src/sink/Default/sink_part.h | 37 +- src/sink/Default/sink_properties.h | 81 ++- src/sink/GEAR/sink.h | 69 +- src/sink/GEAR/sink_iact.h | 93 ++- src/sink/GEAR/sink_io.h | 20 +- src/sink/GEAR/sink_part.h | 36 +- src/sink/GEAR/sink_properties.h | 84 ++- src/sink_debug.h | 2 + src/sink_iact.h | 2 + src/sink_io.h | 2 + src/sink_properties.h | 2 + src/sink_struct.h | 2 + src/space_regrid.c | 13 +- src/space_split.c | 17 +- src/tools.c | 12 +- swift.c | 2 +- 57 files changed, 3794 insertions(+), 302 deletions(-) create mode 100644 examples/SinkParticles/SingleSink/README create mode 100755 examples/SinkParticles/SingleSink/make_sink_ic.py create mode 100644 examples/SinkParticles/SingleSink/make_test_plots.py create mode 100755 examples/SinkParticles/SingleSink/params.yml create mode 100644 src/sink.c create mode 100644 src/sink/Basic/sink.h create mode 100644 src/sink/Basic/sink_debug.h create mode 100644 src/sink/Basic/sink_iact.h create mode 100644 src/sink/Basic/sink_io.h create mode 100644 src/sink/Basic/sink_part.h create mode 100644 src/sink/Basic/sink_properties.h create mode 100644 src/sink/Basic/sink_struct.h diff --git a/configure.ac b/configure.ac index 80c3a8d8d4..944b5a97aa 100644 --- a/configure.ac +++ b/configure.ac @@ -433,6 +433,20 @@ elif test "$stars_density_checks" != "no"; then AC_DEFINE_UNQUOTED([SWIFT_STARS_DENSITY_CHECKS], [$enableval] ,[Enable stars density brute-force checks]) fi +# Check if sink density checks are on for some particles. +AC_ARG_ENABLE([sink-density-checks], + [AS_HELP_STRING([--enable-sink-density-checks], + [Activate expensive brute-force sink density checks for a fraction 1/N of all particles @<:@N@:>@] + )], + [sink_density_checks="$enableval"], + [sink_density_checks="no"] +) +if test "$sink_density_checks" = "yes"; then + AC_MSG_ERROR(Need to specify the fraction of particles to check when using --enable-sink-density-checks!) +elif test "$sink_density_checks" != "no"; then + AC_DEFINE_UNQUOTED([SWIFT_SINK_DENSITY_CHECKS], [$enableval] ,[Enable sink density brute-force checks]) +fi + # Check if ghost statistics are enabled AC_ARG_ENABLE([ghost-statistics], [AS_HELP_STRING([--enable-ghost-statistics], @@ -2865,6 +2879,9 @@ case "$with_sink" in none) AC_DEFINE([SINK_NONE], [1], [No sink particle model]) ;; + Basic) + AC_DEFINE([SINK_BASIC], [1], [Simple, self-contained sink model with only Bondi-Hoyle accretion and no SF.]) + ;; GEAR) AC_DEFINE([SINK_GEAR], [1], [GEAR sink particle model]) ;; diff --git a/doc/RTD/source/SubgridModels/Basic/index.rst b/doc/RTD/source/SubgridModels/Basic/index.rst index ea193ef9eb..710da13715 100644 --- a/doc/RTD/source/SubgridModels/Basic/index.rst +++ b/doc/RTD/source/SubgridModels/Basic/index.rst @@ -5,6 +5,25 @@ Basic model (others) ==================== +Sinks: Simple Bondi-Hoyle accretion +~~~~~~~~~~~~~~~~~~~~~~~~ + +The ``Basic`` sink model provides a foundation on which new sink implementations could be built. It includes a prescription for Bondi-Hoyle gas accretion, and a method for sink-sink mergers that is a slightly simplified version of the implementation used in GEAR. + +No other physics is implemented for this model. Sinks cannot form - to use this model, sink particles must already be present in the initial conditions. They also cannot spawn stars from the gas they accrete. + +Bondi-Hoyle accretion can be done in one of two ways: + + * Gas particles within the sink's kernel are stochastically swallowed entirely, with a probability set by the Bondi-Hoyle rate. Specifically, the probability is set by the current difference between the sink's subgrid mass (determined by the accretion rate) and its dynamical mass (which tracks the number of particles/sinks actually swallowed). This mode is equivalent to the EAGLE black hole accretion model. + * Gas particles within the sink's kernel are "nibbled" down to some minimal mass, which can be specified by the user. This method is equivalent to the black hole accretion model of Bahe et al. 2022. + +This model has only two parameters that must be specified in your parameter ``yml`` file: + + * ``BasicSink:use_nibbling``: determines whether accretion is done by "nibbling" or by swallowing outright. + * ``BasicSink:min_gas_mass_for_nibbling_Msun``: if using "nibbling", the minimum mass to which gas particles can be nibbled. A good default is half the original particle mass. + +For an even more bare-bones starting point, the ``Default`` sink model contains no physics at all, and is a totally blank canvas on which to build your sink model. + Cooling: Analytic models ~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst b/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst index 9f8409de7e..e542d867ce 100644 --- a/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst +++ b/doc/RTD/source/SubgridModels/GEAR/sinks/params.rst @@ -10,12 +10,17 @@ Model parameters The parameters of the GEAR sink model are grouped into the ``GEARSink`` section of the parameter file. -The first two parameters are: +The first three parameters are: -* The sink cut-off radius for gas and sink accretion: ``cut_off_radius``, +* Whether we are using a fixed cut-off radius for gas and sink accretion: ``use_fixed_cut_off_radius``, +* The sink cut-off radius: ``cut_off_radius``, * The sink inner accretion radius fraction in terms of the cut-off radius: ``f_acc``. -The ``f_acc`` parameter is optional. Its default value is :math:`0.1`. Its value must respect :math:`0 \leq f_\text{acc} \leq 1` . It describes the inner radius :math:`f_{\text{acc}} \cdot r_{\text{cut-off}}` in which gas particles are swallowed without any further checks, as explained below. +The ``use_fixed_cut_off_radius`` is mandatory and should be set to 0 or 1. If set to 1, the GEAR model will use a fixed cutoff radius equal to the value of ``cut_off_radius``. If not, the cutoff radius is allowed to vary according to the local gas density. In the code, the cutoff radius is always equal to the sink's smoothing length multiplied by a constant factor ``kernel_gamma`` - setting a fixed cutoff will fix the smoothing length at the appropriate value for all sinks. + +The ``cut_off_radius`` parameter is optional, and is ignored if ``use_fixed_cut_off_radius`` is 0. If a fixed cut-off is used, every sink's smoothing length will be permanently set to this number divided by ``kernel_gamma`` on formation. + +The ``f_acc`` parameter is also optional. Its default value is :math:`0.1`. Its value must respect :math:`0 \leq f_\text{acc} \leq 1` . It describes the inner radius :math:`f_{\text{acc}} \cdot r_{\text{cut-off}}` in which gas particles are swallowed without any further checks, as explained below. The next three mandatory parameters are: @@ -63,7 +68,8 @@ The full section is: .. code:: YAML GEARSink: - cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). This parameter should be adapted with the resolution. + use_fixed_cut_off_radius: 1 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma. + cut_off_radius: 1e-3 # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0. f_acc: 0.1 # (Optional) Fraction of the cut_off_radius that determines if a gas particle should be swallowed wihtout additional check. It has to respect 0 <= f_acc <= 1. (Default: 0.1) temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. density_threshold_Hpcm3: 1e3 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. @@ -104,6 +110,8 @@ This problem can be mitigated by choosing a higher value of ``stellar_particle_m *If you do not want to change your parameters*, you can increase the ``sort_stack_size`` variable at the beginning ``runner_sort.c``. The default value is 10 in powers of 2 (so the stack size is 1024 particles). Increase it to the desired value. Be careful to not overestimate this. +Note that if a cutoff radius is not specified, and the radius is instead left to vary with the local gas density, the smoothing length criterion is always satisfied. + Guide to choose the the accretion radius or the density threshold ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml index 69ce818141..e5a251acf6 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_multi_component/GEAR/params.yml @@ -78,7 +78,8 @@ GEARStarFormation: min_mass_frac: 0.5 GEARSink: - cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. + use_fixed_cut_off_radius: 1 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma. + cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0. f_acc: 0.1 temperature_threshold_K: 1e4 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. density_threshold_Hpcm3: 1e0 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml index ce46711254..c35f58e4a6 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml @@ -88,7 +88,8 @@ DefaultSink: cut_off_radius: 0.1 # Cut off radius of all the sinks in internal units. GEARSink: - cut_off_radius: 1e-2 # Cut off radius of all the sinks in internal units. + use_fixed_cut_off_radius: 1 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma. + cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0. f_acc: 1e-2 temperature_threshold_K: 5e3 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. density_threshold_Hpcm3: 1e0 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. diff --git a/examples/SinkParticles/HomogeneousBox/params.yml b/examples/SinkParticles/HomogeneousBox/params.yml index a805594d05..ac851c5e3e 100755 --- a/examples/SinkParticles/HomogeneousBox/params.yml +++ b/examples/SinkParticles/HomogeneousBox/params.yml @@ -85,7 +85,8 @@ GEARFeedback: GEARSink: - cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. + use_fixed_cut_off_radius: 1 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma. + cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0. f_acc: 1e-2 temperature_threshold_K: 1000 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. density_threshold_Hpcm3: 1e1 # Minimum gas density (in g/cm3) required to form a sink particle. diff --git a/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml b/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml index 90bd18a419..3e7dfb0759 100755 --- a/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml +++ b/examples/SinkParticles/HomogeneousBoxSinkParticles/params.yml @@ -83,7 +83,8 @@ GEARFeedback: elements: [Fe, Mg, O, C, Al, Ca, Ba, Zn, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). GEARSink: - cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. + use_fixed_cut_off_radius: 1 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma. + cut_off_radius: 5e-3 # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0. f_acc: 0.1 temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_g_per_cm3 <= density <= maximal_density_threshold_g_per_cm3. density_threshold_Hpcm3: 1e1 # Minimum gas density (in g/cm3) required to form a sink particle. diff --git a/examples/SinkParticles/PlummerSphere/params.yml b/examples/SinkParticles/PlummerSphere/params.yml index 436773e49e..262c7dba97 100755 --- a/examples/SinkParticles/PlummerSphere/params.yml +++ b/examples/SinkParticles/PlummerSphere/params.yml @@ -84,7 +84,8 @@ GEARFeedback: elements: [Fe, Mg, O, C, Al, Ca, Ba, Zn, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). GEARSink: - cut_off_radius: 1e-2 # Cut off radius of all the sinks in internal units. + use_fixed_cut_off_radius: 1 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma. + cut_off_radius: 1e-2 # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0. f_acc: 0.4 temperature_threshold_K: 3e4 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. density_threshold_Hpcm3: 1e0 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. diff --git a/examples/SinkParticles/SingleSink/README b/examples/SinkParticles/SingleSink/README new file mode 100644 index 0000000000..1498298fce --- /dev/null +++ b/examples/SinkParticles/SingleSink/README @@ -0,0 +1,31 @@ +# Intro +This example is a non cosmological homogeneous box containing gas and a single sink in the centre. It's designed for testing that your sink model accretes the amount of gas that you're expecting. It's mainly designed to test the Basic model, which requires no extra subgrid physics to function. + +# ICs +The included python script `make_sink_ic.py` creates the initial conditions. You can set the particle masses of the gas and sink, the gas density and velocity dispersion (which sets the internal energy), and set an overall "level", which determines how many particles to simulate (N_particle = 2^(3*level)). Run `python make_sink_ic.py --help` to get the list of parameters. + +Running the script with no options produces a "level 6" box at "m5" resolution, with a sink mass 50x that of the gas mass, and a density and velocity dispersion of 0.1 atoms/cc and 10 km/s respectively, representative of the ISM. The ICs have the default name `ics.hdf5`. + +# Configure +To run this example with GEAR model, + +./configure --with-sink=Basic --with-kernel=wendland-C2 + +and then + +make -j + +You can also add `--enable-sink-density-checks=<CHECK_FREQUENCY>` to run the brute-force density checks for the sink, to make sure things are being calculated properly. + +# Run +We can now run with + +swift --hydro --self-gravity --sinks --threads=<NUM_THREADS> params.yml + +By default the simulation will run for 500 Myr, and output a snapshot every 10 Myr. You can change these options in your paramfile. + +# Testing the output + +Included in this directory is a very simple script for making some test plots. You can run it with + +python make_test_plots.py <snapshot-directory-name> \ No newline at end of file diff --git a/examples/SinkParticles/SingleSink/make_sink_ic.py b/examples/SinkParticles/SingleSink/make_sink_ic.py new file mode 100755 index 0000000000..be55d45f26 --- /dev/null +++ b/examples/SinkParticles/SingleSink/make_sink_ic.py @@ -0,0 +1,234 @@ +################################################################################ +# This file is part of SWIFT. +# Copyright (c) 2022 Yves Revaz (yves.revaz@epfl.ch) +# Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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 +from optparse import OptionParser +from astropy import units +from astropy import constants + + +def parse_options(): + + usage = "usage: %prog [options] file" + parser = OptionParser(usage=usage) + + parser.add_option( + "--rho", + action="store", + dest="rho", + type="float", + default=0.1, + help="Mean gas density in atom/cm3" + ) + + parser.add_option( + "--sigma", + action="store", + dest="sigma", + type="float", + default=10., + help="Velocity dispersion of the gas in km/s (sets internal energy)." + ) + + parser.add_option( + "--gas-mass", + action="store", + dest="gas_mass", + type="float", + default=1e5, + help="Gas particle mass in solar masses" + ) + + parser.add_option( + "--sink-mass", + action="store", + dest="sink_mass", + type="float", + default=5e6, + help="Sink particle mass in solar masses" + ) + + parser.add_option( + "--sink-velocity", + action="store", + dest="sink_mach", + type="float", + default=0., + help="Sink velocity as a multiple of the sound speed (i.e. Mach number)" + ) + + parser.add_option( + "--level", + action="store", + dest="level", + type="int", + default=6, + help="Resolution level: N = (2**l)**3" + ) + + parser.add_option( + "-o", + action="store", + dest="outputfilename", + type="string", + default="ics.hdf5", + help="output filename", + ) + + (options, args) = parser.parse_args() + + files = args + + return files, options + + +######################################## +# main +######################################## + +files, opt = parse_options() + +# define standard units +UnitMass_in_cgs = 1.988409870698051e43 # 10^10 M_sun in grams +UnitLength_in_cgs = 3.0856775814913673e21 # kpc in centimeters +UnitVelocity_in_cgs = 1e5 # km/s in centimeters per second +UnitCurrent_in_cgs = 1 # Amperes +UnitTemp_in_cgs = 1 # Kelvin +UnitTime_in_cgs = UnitLength_in_cgs / UnitVelocity_in_cgs + +UnitMass = UnitMass_in_cgs * units.g +UnitLength = UnitLength_in_cgs * units.cm +UnitTime = UnitTime_in_cgs * units.s +UnitVelocity = UnitVelocity_in_cgs * units.cm / units.s + +np.random.seed(1) + +# Number of particles +N = (2 ** opt.level) ** 3 # number of particles + +# Mean density +rho = opt.rho # atom/cc +rho = rho * constants.m_p / units.cm ** 3 + +# Gas particle mass +m = opt.gas_mass * units.Msun # in solar mass + +# Sink particle mass +sm = opt.sink_mass * units.Msun # in solar mass + +# Gas mass in the box +M = N * m + +# Size of the box +L = (M / rho) ** (1 / 3.0) + +# Gravitational constant +G = constants.G + +print("Number of particles : {}".format(N)) + +# Convert to code units +m = m.to(UnitMass).value +sm = sm.to(UnitMass).value +L = L.to(UnitLength).value +rho = rho.to(UnitMass / UnitLength ** 3).value +sigma = (opt.sigma * units.km/units.s).to(UnitVelocity).value + +# Generate the particles +pos = np.random.random([N, 3]) * np.array([L, L, L]) +vel = np.zeros([N, 3]) +mass = np.ones(N) * m +u = np.ones(N) * sigma**2 +ids = np.arange(N) +h = np.ones(N) * 3 * L / N ** (1 / 3.0) +rho = np.ones(N) * rho + +print("Inter-particle distance (code unit) : {}".format(L / N ** (1 / 3.0))) + +# Generate the sink +# Always put it 10% of the way through the box to give it room to move (if it's going to) + +NSINK = 1 + +sink_pos = np.ones([NSINK, 3]) +sink_pos[:,0] = L/10. +sink_pos[:,1] = L/2. +sink_pos[:,2] = L/2. + +sink_mass = np.array([sm,]) +sink_ids = np.array([2*ids[-1]]) +sink_h = np.array([3 * L / N ** (1 / 3.0),]) + +gas_cs = np.sqrt(sigma**2 * 5./3. * ((5./3.)-1)) + +sink_vel = np.zeros([NSINK, 3]) +sink_vel[:,0] += gas_cs * opt.sink_mach + +print(f"Sink velocity: {gas_cs * opt.sink_mach}") + +if gas_cs * opt.sink_mach > 0: + sink_time_in_box = L*0.9 / (gas_cs * opt.sink_mach) + print(f"Sink will leave box (neglecting dynamical friction) after time: {sink_time_in_box}") + + +# File +fileOutput = h5py.File(opt.outputfilename, "w") + +# Header +grp = fileOutput.create_group("/Header") +grp.attrs["BoxSize"] = [L, L, L] +grp.attrs["NumPart_Total"] = [N, 0, 0, NSINK, 0, 0] +grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] +grp.attrs["NumPart_ThisFile"] = [N, 0, 0, NSINK, 0, 0] +grp.attrs["Time"] = 0.0 +grp.attrs["NumFileOutputsPerSnapshot"] = 1 +grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0] +grp.attrs["Dimension"] = 3 + +# Units +grp = fileOutput.create_group("/Units") +grp.attrs["Unit length in cgs (U_L)"] = UnitLength_in_cgs +grp.attrs["Unit mass in cgs (U_M)"] = UnitMass_in_cgs +grp.attrs["Unit time in cgs (U_t)"] = UnitTime_in_cgs +grp.attrs["Unit current in cgs (U_I)"] = UnitCurrent_in_cgs +grp.attrs["Unit temperature in cgs (U_T)"] = UnitTemp_in_cgs + +# Particle groups +grp = fileOutput.create_group("/PartType0") +grp.create_dataset("Coordinates", data=pos, dtype="d") +grp.create_dataset("Velocities", data=vel, dtype="f") +grp.create_dataset("Masses", data=mass, 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("Densities", data=rho, dtype="f") + +grp = fileOutput.create_group("/PartType3") +grp.create_dataset("Coordinates", data=sink_pos, dtype="d") +grp.create_dataset("Velocities", data=sink_vel, dtype="f") +grp.create_dataset("Masses", data=sink_mass, dtype="f") +grp.create_dataset("SmoothingLength", data=sink_h, dtype="f") +grp.create_dataset("ParticleIDs", data=sink_ids, dtype="L") + +fileOutput.close() + +print(f"{opt.outputfilename} saved.") \ No newline at end of file diff --git a/examples/SinkParticles/SingleSink/make_test_plots.py b/examples/SinkParticles/SingleSink/make_test_plots.py new file mode 100644 index 0000000000..d4d4c52ba6 --- /dev/null +++ b/examples/SinkParticles/SingleSink/make_test_plots.py @@ -0,0 +1,123 @@ +################################################################################ +# This file is part of SWIFT. +# Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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/>. +# +################################################################################ + +''' +Make some very simple plots showing the evolution of the single sink. + +Run with python make_test_plots.py <snapshot-directory-name> + +''' + + +import numpy as np +from sys import argv +from glob import glob +import matplotlib.pyplot as plt +import unyt +import h5py + +params = {"text.usetex": True, 'axes.labelsize': 16, 'xtick.labelsize': 13, 'ytick.labelsize': 13, 'lines.linewidth' : 2, 'axes.titlesize' : 16, 'font.family' : 'serif'} +plt.rcParams.update(params) + +# Units +G_cgs = 6.6743e-8 * unyt.cm**3 * unyt.g**-1 * unyt.s**-2 +unit_mass_cgs = 1.988409870698051e+43 * unyt.g +unit_density_cgs = 6.767905323247329e-22 * unyt.Unit('g/cm**3') +unit_velocity_cgs = (1.*unyt.Unit("km/s")).to('cm/s') + + +# Basic Bondi-Hoyle prediction from the sink's starting mass in a constant medium +def simple_bondi_hoyle(t, m, rho, v, cs): + + m_sink = np.zeros(len(t),dtype=np.float64) * unyt.g + + m_sink[0] = m[0] * unit_mass_cgs + + rho_0 = rho[0] * unit_density_cgs + v_0 = v[0] * unit_velocity_cgs + cs_0 = cs[0] * unit_velocity_cgs + + timestep = ((t[1] - t[0]) * unyt.Gyr).to('s') + + for i in range(len(times)-1): + + numerator = 4. * np.pi * G_cgs**2 * m_sink[i]**2 * rho_0 + denominator = np.power(v_0**2 + cs_0**2, 3./2.) + + m_sink[i+1] = m_sink[i] + (numerator/denominator) * timestep + + return m_sink.to('Msun') + + +snapshots = glob(f'./{argv[1]}/*.hdf5') + +times = np.empty(len(snapshots),dtype=np.float32) +mass_evo = np.empty(len(snapshots),dtype=np.float32) +subgrid_mass_evo = np.empty(len(snapshots),dtype=np.float32) +rho_evo = np.empty(len(snapshots),dtype=np.float32) +v_evo = np.empty(len(snapshots),dtype=np.float32) +cs_evo = np.empty(len(snapshots),dtype=np.float32) +hsml_evo = np.empty(len(snapshots),dtype=np.float32) + + +for s, snap in enumerate(snapshots): + + with h5py.File(snap) as f: + times[s] = f['Header'].attrs['Time'][0] + mass_evo[s] = f['PartType3/Masses'][0] + subgrid_mass_evo[s] = f['PartType3/SubgridMasses'][0] + rho_evo[s] = f['PartType3/GasDensities'][0] + cs_evo[s] = f['PartType3/GasSoundSpeeds'][0] + hsml_evo[s] = f['PartType3/SmoothingLengths'][0] + + v = f['PartType3/GasVelocities'][0] - f['PartType3/Velocities'][0] + v_evo[s] = np.sqrt(v[0]**2+v[1]**2+v[2]**2) + + +# Run the Bondi-Hoyle prediction +m_sink_bondi_prediction = simple_bondi_hoyle(times,mass_evo,rho_evo,v_evo,cs_evo) + +# Normalise time to go up to 1 +t = times/times[-1] + +# Time evolution of the sink mass, target mass, and Bondi-Hoyle prediction +fig, ax = plt.subplots(1,figsize=(8,6)) +ax.plot(t,np.log10(m_sink_bondi_prediction),label=r'Simple Bondi-Hoyle, constant $\rho$, $v$, $c_{\rm s}$') +ax.plot(t,np.log10(subgrid_mass_evo * 1e10),label=r'Subgrid mass') +ax.plot(t,np.log10(mass_evo * 1e10),label=r'Dynamical mass') + +ax.set_xlabel(r"$t/t_{\rm final}$") +ax.set_ylabel(r"$\log_{10}(M_{\rm sink})$") +ax.legend(fontsize=14) +ax.set_ylim(6.68,7.) +fig.savefig('mass_evolution.png',bbox_inches='tight') + +# Time evolution of the total mass accreted relative to the target mass +fig, ax = plt.subplots(1,figsize=(8,6)) +ax.plot(t,mass_evo/subgrid_mass_evo) +ax.set_xlabel(r"$t/t_{\rm final}$") +ax.set_ylabel(r"$m_{\rm accr}^{\rm gas}/m_{\rm target}^{\rm gas}$") +fig.savefig('target_mass_ratio.png',bbox_inches='tight') + +# Time evolution of the smoothing length +fig, ax = plt.subplots(1,figsize=(8,6)) +ax.plot(t,hsml_evo) +ax.set_xlabel(r"$t/t_{\rm final}$") +ax.set_ylabel(r"$h_{\rm sink}$ [kpc]") +fig.savefig('smoothing_length.png',bbox_inches='tight') diff --git a/examples/SinkParticles/SingleSink/params.yml b/examples/SinkParticles/SingleSink/params.yml new file mode 100755 index 0000000000..eafa709d36 --- /dev/null +++ b/examples/SinkParticles/SingleSink/params.yml @@ -0,0 +1,61 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.988409870698051e+43 # 10^10 solar masses + UnitLength_in_cgs: 3.0856775814913673e+21 # 1 kpc + UnitVelocity_in_cgs: 1e5 # km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters for the self-gravity scheme +Gravity: + MAC: adaptive # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'. + epsilon_fmm: 0.001 # Tolerance parameter for the adaptive multipole acceptance criterion. + theta_cr: 0.7 # Opening angle for the purely gemoetric criterion. + eta: 0.025 # Constant dimensionless multiplier for time integration. + theta: 0.7 # Opening angle (Multipole acceptance criterion). + max_physical_baryon_softening: 0.35 # Physical softening length (in internal units) + mesh_side_length: 32 + +# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 0.5 # 500 Myr # 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-2 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + subdir: snap + basename: snapshot # Common part of the name of output files + time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units) + delta_time: 1e-2 # Time difference between consecutive outputs (in internal units) + +Scheduler: + cell_extra_gparts: 10000 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation. + cell_extra_sinks: 10000 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation. + cell_extra_sparts: 10000 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation. + max_top_level_cells: 3 # + dependency_graph_frequency: 0 # (Optional) Dumping frequency of the dependency graph. By default, writes only at the first step. + + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-3 # Time between statistics output + time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units) + +# Parameters related to the initial conditions +InitialConditions: + file_name: ics.hdf5 + periodic: 1 # Are we running with periodic ICs? + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 57Ngbs with the Wendland C2 kernel). + # resolution_eta: 1.55 # Double the number of neighbours + CFL_condition: 0.2 # Courant-Friedrich-Levy condition for time integration. + h_max: 250. + minimal_temperature: 1 + +BasicSink: + use_nibbling: 1 + min_gas_mass_for_nibbling_Msun: 5e4 # half the gas particle mass in msun diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 36e23ed82b..9af0363637 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -876,9 +876,13 @@ XrayEmissivity: DefaultSink: cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). This parameter should be adapted with the resolution. +BasicSink: + use_nibbling: 1 # Use nibbling to accrete gas? If zero, do swallowing instead. + min_gas_mass_for_nibbling_Msun: 5e4 # The mass in Msun that particles cannot be nibbled below. A good default is half the initial mass. + # GEAR sink particles GEARSink: - cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). This parameter should be adapted with the resolution. + cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). In GEAR, if this is specified, the cutoff radius is fixed, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma. f_acc: 0.1 # (Optional) Fraction of the cut_off_radius that determines if a gas particle should be swallowed wihtout additional check. It has to respect 0 <= f_acc <= 1. (Default: 0.1) temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3. density_threshold_Hpcm3: 1e3 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle. diff --git a/src/Makefile.am b/src/Makefile.am index 3d2a06d80f..cbf6d4f785 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -179,7 +179,7 @@ AM_SOURCES += output_options.c line_of_sight.c restart.c parser.c xmf.c AM_SOURCES += kernel_hydro.c tools.c map.c part.c partition.c clocks.c AM_SOURCES += physical_constants.c units.c potential.c hydro_properties.c AM_SOURCES += threadpool.c cooling.c star_formation.c -AM_SOURCES += hydro.c stars.c +AM_SOURCES += hydro.c stars.c sink.c AM_SOURCES += statistics.c profiler.c csds.c part_type.c AM_SOURCES += gravity_properties.c gravity.c multipole.c AM_SOURCES += collectgroup.c hydro_space.c equation_of_state.c io_compression.c @@ -495,6 +495,8 @@ nobase_noinst_HEADERS += sink/Default/sink_iact.h sink/Default/sink_struct.h sin nobase_noinst_HEADERS += sink/GEAR/sink.h sink/GEAR/sink_io.h sink/GEAR/sink_part.h nobase_noinst_HEADERS += sink/GEAR/sink_properties.h sink/GEAR/sink_setters.h sink/GEAR/sink_getters.h nobase_noinst_HEADERS += sink/GEAR/sink_iact.h sink/GEAR/sink_struct.h sink/GEAR/sink_debug.h +nobase_noinst_HEADERS += sink/Basic/sink.h sink/Basic/sink_io.h sink/Basic/sink_part.h sink/Basic/sink_properties.h +nobase_noinst_HEADERS += sink/Basic/sink_iact.h sink/Basic/sink_struct.h sink/Basic/sink_debug.h nobase_noinst_HEADERS += neutrino.h neutrino_properties.h neutrino_io.h nobase_noinst_HEADERS += neutrino/Default/neutrino.h neutrino/Default/relativity.h neutrino/Default/fermi_dirac.h nobase_noinst_HEADERS += neutrino/Default/neutrino_properties.h neutrino/Default/neutrino_io.h diff --git a/src/cell.h b/src/cell.h index 9cb04b32c4..7ad7b17e42 100644 --- a/src/cell.h +++ b/src/cell.h @@ -202,7 +202,7 @@ struct pcell { int count; /*! Maximal cut off radius. */ - float r_cut_max; + float h_max; /*! Minimal integer end-of-timestep in this cell for sinks tasks */ integertime_t ti_end_min; @@ -865,7 +865,7 @@ cell_can_recurse_in_pair_sinks_task(const struct cell *ci, /* smaller than the sub-cell sizes ? */ /* Note: We use the _old values as these might have been updated by a drift */ return ci->split && cj->split && - ((ci->sinks.r_cut_max_old + ci->sinks.dx_max_part_old) < + ((kernel_gamma * ci->sinks.h_max_old + ci->sinks.dx_max_part_old) < 0.5f * ci->dmin) && ((kernel_gamma * cj->hydro.h_max_old + cj->hydro.dx_max_part_old) < 0.5f * cj->dmin); @@ -918,7 +918,7 @@ __attribute__((always_inline)) INLINE static int cell_can_recurse_in_self_sinks_task(const struct cell *c) { /* Is the cell split and not smaller than the smoothing length? */ - return c->split && (c->sinks.r_cut_max_old < 0.5f * c->dmin) && + return c->split && (kernel_gamma * c->sinks.h_max_old < 0.5f * c->dmin) && (kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin); } @@ -1097,7 +1097,7 @@ cell_need_rebuild_for_sinks_pair(const struct cell *ci, const struct cell *cj) { /* Is the cut-off radius plus the max distance the parts in both cells have */ /* moved larger than the cell size ? */ /* Note ci->dmin == cj->dmin */ - if (max(ci->sinks.r_cut_max, kernel_gamma * cj->hydro.h_max) + + if (max(kernel_gamma * ci->sinks.h_max, kernel_gamma * cj->hydro.h_max) + ci->sinks.dx_max_part + cj->hydro.dx_max_part > cj->dmin) { return 1; diff --git a/src/cell_convert_part.c b/src/cell_convert_part.c index 1d9154f832..1bee5146c7 100644 --- a/src/cell_convert_part.c +++ b/src/cell_convert_part.c @@ -1140,7 +1140,7 @@ struct spart *cell_spawn_new_spart_from_sink(struct engine *e, struct cell *c, #endif /* Set a smoothing length */ - sp->h = s->r_cut; + sp->h = s->h; /* Here comes the Sun! */ return sp; diff --git a/src/cell_drift.c b/src/cell_drift.c index f394645555..f4f8d146b8 100644 --- a/src/cell_drift.c +++ b/src/cell_drift.c @@ -991,8 +991,8 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) { struct sink *const sinks = c->sinks.parts; float dx_max = 0.f, dx2_max = 0.f; - float cell_r_max = 0.f; - float cell_r_max_active = 0.f; + float cell_h_max = 0.f; + float cell_h_max_active = 0.f; /* Drift irrespective of cell flags? */ force = (force || cell_get_flag(c, cell_flag_do_sink_drift)); @@ -1032,14 +1032,14 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) { /* Update */ dx_max = max(dx_max, cp->sinks.dx_max_part); - cell_r_max = max(cell_r_max, cp->sinks.r_cut_max); - cell_r_max_active = max(cell_r_max_active, cp->sinks.r_cut_max_active); + cell_h_max = max(cell_h_max, cp->sinks.h_max); + cell_h_max_active = max(cell_h_max_active, cp->sinks.h_max_active); } } /* Store the values */ - c->sinks.r_cut_max = cell_r_max; - c->sinks.r_cut_max_active = cell_r_max_active; + c->sinks.h_max = cell_h_max; + c->sinks.h_max_active = cell_h_max_active; c->sinks.dx_max_part = dx_max; /* Update the time of the last drift */ @@ -1118,7 +1118,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) { dx2_max = max(dx2_max, dx2); /* Maximal smoothing length */ - cell_r_max = max(cell_r_max, sink->r_cut); + cell_h_max = max(cell_h_max, sink->h); /* Mark the particle has not being swallowed */ sink_mark_sink_as_not_swallowed(&sink->merger_data); @@ -1127,7 +1127,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) { if (sink_is_active(sink, e)) { sink_init_sink(sink); - cell_r_max_active = max(cell_r_max_active, sink->r_cut); + cell_h_max_active = max(cell_h_max_active, sink->h); } } @@ -1135,8 +1135,8 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) { dx_max = sqrtf(dx2_max); /* Store the values */ - c->sinks.r_cut_max = cell_r_max; - c->sinks.r_cut_max_active = cell_r_max_active; + c->sinks.h_max = cell_h_max; + c->sinks.h_max_active = cell_h_max_active; c->sinks.dx_max_part = dx_max; /* Update the time of the last drift */ diff --git a/src/cell_pack.c b/src/cell_pack.c index 83f8c4cab9..d3bcc175f5 100644 --- a/src/cell_pack.c +++ b/src/cell_pack.c @@ -44,7 +44,7 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc, pc->hydro.h_max = c->hydro.h_max; pc->stars.h_max = c->stars.h_max; pc->black_holes.h_max = c->black_holes.h_max; - pc->sinks.r_cut_max = c->sinks.r_cut_max; + pc->sinks.h_max = c->sinks.h_max; pc->hydro.ti_end_min = c->hydro.ti_end_min; pc->grav.ti_end_min = c->grav.ti_end_min; @@ -246,7 +246,7 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c, c->hydro.h_max = pc->hydro.h_max; c->stars.h_max = pc->stars.h_max; c->black_holes.h_max = pc->black_holes.h_max; - c->sinks.r_cut_max = pc->sinks.r_cut_max; + c->sinks.h_max = pc->sinks.h_max; c->hydro.ti_end_min = pc->hydro.ti_end_min; c->grav.ti_end_min = pc->grav.ti_end_min; diff --git a/src/cell_sinks.h b/src/cell_sinks.h index 9043ad1915..2a8bb3c6a1 100644 --- a/src/cell_sinks.h +++ b/src/cell_sinks.h @@ -91,11 +91,12 @@ struct cell_sinks { /*! Nr of #sink this cell can hold after addition of new one. */ int count_total; - /*! Max cut off radius of active particles in this cell. */ - float r_cut_max_active; + /*! Max smoothing length of active particles in this cell. + */ + float h_max_active; - /*! Values of r_cut_max before the drifts, used for sub-cell tasks. */ - float r_cut_max_old; + /*! Values of h_max before the drifts, used for sub-cell tasks. */ + float h_max_old; /*! Maximum part movement in this cell since last construction. */ float dx_max_part; @@ -117,8 +118,8 @@ struct cell_sinks { /*! Spin lock for various uses (#sink case). */ swift_lock_type lock; - /*! Max cut off radius in this cell. */ - float r_cut_max; + /*! Max smoothing length in this cell. */ + float h_max; /*! Number of #sink updated in this cell. */ int updated; diff --git a/src/cell_unskip.c b/src/cell_unskip.c index adfd6eb6ef..f7b93c3176 100644 --- a/src/cell_unskip.c +++ b/src/cell_unskip.c @@ -1151,13 +1151,13 @@ void cell_activate_subcell_sinks_tasks(struct cell *ci, struct cell *cj, /* Store the current dx_max and h_max values. */ ci->sinks.dx_max_part_old = ci->sinks.dx_max_part; - ci->sinks.r_cut_max_old = ci->sinks.r_cut_max; + ci->sinks.h_max_old = ci->sinks.h_max; ci->hydro.dx_max_part_old = ci->hydro.dx_max_part; ci->hydro.h_max_old = ci->hydro.h_max; if (cj != NULL) { cj->sinks.dx_max_part_old = cj->sinks.dx_max_part; - cj->sinks.r_cut_max_old = cj->sinks.r_cut_max; + cj->sinks.h_max_old = cj->sinks.h_max; cj->hydro.dx_max_part_old = cj->hydro.dx_max_part; cj->hydro.h_max_old = cj->hydro.h_max; } diff --git a/src/debug.c b/src/debug.c index a82363365c..0b1538b5c4 100644 --- a/src/debug.c +++ b/src/debug.c @@ -242,8 +242,8 @@ int checkSpacehmax(struct space *s) { float cell_sinks_h_max = 0.0f; for (int k = 0; k < s->nr_cells; k++) { if (s->cells_top[k].nodeID == s->e->nodeID && - s->cells_top[k].sinks.r_cut_max > cell_sinks_h_max) { - cell_sinks_h_max = s->cells_top[k].sinks.r_cut_max; + s->cells_top[k].sinks.h_max > cell_sinks_h_max) { + cell_sinks_h_max = s->cells_top[k].sinks.h_max; } } @@ -266,8 +266,8 @@ int checkSpacehmax(struct space *s) { /* Now all the sinks. */ float sink_h_max = 0.0f; for (size_t k = 0; k < s->nr_sinks; k++) { - if (s->sinks[k].r_cut > sink_h_max) { - sink_h_max = s->sinks[k].r_cut; + if (s->sinks[k].h > sink_h_max) { + sink_h_max = s->sinks[k].h; } } @@ -315,17 +315,17 @@ int checkSpacehmax(struct space *s) { /* sink */ for (int k = 0; k < s->nr_cells; k++) { if (s->cells_top[k].nodeID == s->e->nodeID) { - if (s->cells_top[k].sinks.r_cut_max > sink_h_max) { + if (s->cells_top[k].sinks.h_max > sink_h_max) { message("cell %d is inconsistent (%f > %f)", k, - s->cells_top[k].sinks.r_cut_max, sink_h_max); + s->cells_top[k].sinks.h_max, sink_h_max); } } } for (size_t k = 0; k < s->nr_sinks; k++) { - if (s->sinks[k].r_cut > cell_sinks_h_max) { + if (s->sinks[k].h > cell_sinks_h_max) { message("spart %lld is inconsistent (%f > %f)", s->sinks[k].id, - s->sinks[k].r_cut, cell_sinks_h_max); + s->sinks[k].h, cell_sinks_h_max); } } @@ -436,7 +436,7 @@ int checkCellhdxmax(const struct cell *c, int *depth) { sp->x_diff[1] * sp->x_diff[1] + sp->x_diff[2] * sp->x_diff[2]; - sinks_h_max = max(sinks_h_max, sp->r_cut); + sinks_h_max = max(sinks_h_max, sp->h); sinks_dx_max = max(sinks_dx_max, sqrtf(dx2)); } @@ -476,9 +476,9 @@ int checkCellhdxmax(const struct cell *c, int *depth) { result = 0; } - if (c->sinks.r_cut_max != sinks_h_max) { + if (c->sinks.h_max != sinks_h_max) { message("%d Inconsistent sinks_h_max: cell %f != parts %f", *depth, - c->sinks.r_cut_max, sinks_h_max); + c->sinks.h_max, sinks_h_max); message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]); result = 0; } diff --git a/src/engine.c b/src/engine.c index f5e6cbd5af..6e596aedf9 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2360,6 +2360,15 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs, stars_exact_density_check(e->s, e, /*rel_tol=*/1e-3); #endif +#ifdef SWIFT_SINK_DENSITY_CHECKS + /* Run the brute-force sink calculation for some sinks */ + if (e->policy & engine_policy_sinks) sink_exact_density_compute(e->s, e); + + /* Check the accuracy of the sink calculation */ + if (e->policy & engine_policy_sinks) + sink_exact_density_check(e->s, e, /*rel_tol=*/1e-3); +#endif + #ifdef SWIFT_GRAVITY_FORCE_CHECKS /* Check the accuracy of the gravity calculation */ if (e->policy & engine_policy_self_gravity) @@ -2472,12 +2481,12 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs, for (int i = 0; i < s->nr_cells; i++) { struct cell *c = &s->cells_top[i]; if (c->nodeID == engine_rank && c->sinks.count > 0) { - float sink_h_max = c->sinks.parts[0].r_cut; + float sink_h_max = c->sinks.parts[0].h; for (int k = 1; k < c->sinks.count; k++) { - if (c->sinks.parts[k].r_cut > sink_h_max) - sink_h_max = c->sinks.parts[k].r_cut; + if (c->sinks.parts[k].h > sink_h_max) + sink_h_max = c->sinks.parts[k].h; } - c->sinks.r_cut_max = max(sink_h_max, c->sinks.r_cut_max); + c->sinks.h_max = max(sink_h_max, c->sinks.h_max); } } } @@ -2906,6 +2915,15 @@ int engine_step(struct engine *e) { stars_exact_density_check(e->s, e, /*rel_tol=*/1e-2); #endif +#ifdef SWIFT_SINK_DENSITY_CHECKS + /* Run the brute-force sink calculation for some sinks */ + if (e->policy & engine_policy_sinks) sink_exact_density_compute(e->s, e); + + /* Check the accuracy of the sink calculation */ + if (e->policy & engine_policy_sinks) + sink_exact_density_check(e->s, e, /*rel_tol=*/1e-2); +#endif + #ifdef SWIFT_GRAVITY_FORCE_CHECKS /* Check if we want to run force checks this timestep. */ if (e->policy & engine_policy_self_gravity) { diff --git a/src/part.h b/src/part.h index c0e87bfde7..bb574651b5 100644 --- a/src/part.h +++ b/src/part.h @@ -132,6 +132,8 @@ struct threadpool; /* Import the right sink particle definition */ #if defined(SINK_NONE) #include "./sink/Default/sink_part.h" +#elif defined(SINK_BASIC) +#include "./sink/Basic/sink_part.h" #elif defined(SINK_GEAR) #include "./sink/GEAR/sink_part.h" #else diff --git a/src/proxy.c b/src/proxy.c index 3830162ada..03608afaf1 100644 --- a/src/proxy.c +++ b/src/proxy.c @@ -796,7 +796,7 @@ void proxy_parts_exchange_first(struct proxy *p) { for (int k = 0; k < p->nr_sinks_out; k++) message("sending sinks %lli, x=[%.3e %.3e %.3e], h=%.3e, to node %i.", p->sinks_out[k].id, p->sinks_out[k].x[0], p->sinks_out[k].x[1], - p->sinks_out[k].x[2], p->sinks_out[k].r_cut, p->nodeID); + p->sinks_out[k].x[2], p->sinks_out[k].h, p->nodeID); #endif /* SWIFT_DEBUG_CHECKS */ } @@ -953,7 +953,7 @@ void proxy_parts_exchange_second(struct proxy *p) { for (int k = 0; k < p->nr_sinks_in; k++) message("receiving sinks %lli, x=[%.3e %.3e %.3e], h=%.3e, from node %i.", p->sinks_in[k].id, p->sinks_in[k].x[0], p->sinks_in[k].x[1], - p->sinks_in[k].x[2], p->sinks_in[k].r_cut, p->nodeID); + p->sinks_in[k].x[2], p->sinks_in[k].h, p->nodeID); #endif /* SWIFT_DEBUG_CHECKS */ } diff --git a/src/random.h b/src/random.h index 9b5d75ddc1..bce43ccca2 100644 --- a/src/random.h +++ b/src/random.h @@ -62,6 +62,7 @@ enum random_number_type { random_number_BH_reposition = 59969537LL, random_number_BH_spin = 193877777LL, random_number_BH_kick = 303595777LL, + random_number_sink_swallow = 7337737777LL, random_number_snapshot_sampling = 6561001LL, random_number_stellar_winds = 5947309451LL, random_number_HII_regions = 8134165677LL, diff --git a/src/runner_doiact_functions_hydro.h b/src/runner_doiact_functions_hydro.h index 1bcf1af207..d8bf99b2f6 100644 --- a/src/runner_doiact_functions_hydro.h +++ b/src/runner_doiact_functions_hydro.h @@ -123,8 +123,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -145,8 +144,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -260,8 +258,7 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci, runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -277,8 +274,7 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -298,8 +294,7 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -402,8 +397,7 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) { runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -419,8 +413,7 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -440,8 +433,7 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -543,8 +535,7 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) { runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -560,8 +551,7 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -581,8 +571,7 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -683,8 +672,7 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, pj->h, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, pj->h, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, pj->h, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, pj->h, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, pj->h, pi, pj, a, H); @@ -792,8 +780,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -858,8 +845,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -1025,8 +1011,7 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -1208,8 +1193,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -1308,8 +1292,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -1629,8 +1612,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -1712,8 +1694,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -1728,8 +1709,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -1844,8 +1824,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -1929,8 +1908,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -1945,8 +1923,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid, runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -2151,8 +2128,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -2213,8 +2189,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -2230,8 +2205,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -2251,8 +2225,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -2395,8 +2368,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H); @@ -2452,8 +2424,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { runner_iact_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H); @@ -2468,8 +2439,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H); #endif #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE) runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H); diff --git a/src/runner_doiact_functions_sinks.h b/src/runner_doiact_functions_sinks.h index 6ff95ba1de..87b993882a 100644 --- a/src/runner_doiact_functions_sinks.h +++ b/src/runner_doiact_functions_sinks.h @@ -59,8 +59,8 @@ void DOSELF1_SINKS(struct runner *r, struct cell *c, int timer) { /* Skip inactive particles */ if (!sink_is_active(si, e)) continue; - const float ri = si->r_cut; - const float ri2 = ri * ri; + const float hi = si->h; + const float hig2 = hi * hi * kernel_gamma2; const float six[3] = {(float)(si->x[0] - c->loc[0]), (float)(si->x[1] - c->loc[1]), (float)(si->x[2] - c->loc[2])}; @@ -90,8 +90,8 @@ void DOSELF1_SINKS(struct runner *r, struct cell *c, int timer) { error("Particle pj not drifted to current time"); #endif - if (r2 < ri2) { - IACT_SINKS_GAS(r2, dx, ri, hj, si, pj, with_cosmology, cosmo, + if (r2 < hig2) { + IACT_SINKS_GAS(r2, dx, hi, hj, si, pj, with_cosmology, cosmo, e->gravity_properties, e->sink_properties, e->ti_current, e->time); } @@ -112,8 +112,8 @@ void DOSELF1_SINKS(struct runner *r, struct cell *c, int timer) { /* Skip inactive particles */ if (!sink_is_active(si, e)) continue; - const float ri = si->r_cut; - const float ri2 = ri * ri; + const float hi = si->h; + const float hig2 = hi * hi * kernel_gamma2; const float six[3] = {(float)(si->x[0] - c->loc[0]), (float)(si->x[1] - c->loc[1]), (float)(si->x[2] - c->loc[2])}; @@ -126,8 +126,8 @@ void DOSELF1_SINKS(struct runner *r, struct cell *c, int timer) { /* Get a pointer to the jth particle. */ struct sink *restrict sj = &sinks[sjd]; - const float rj = sj->r_cut; - const float rj2 = rj * rj; + const float hj = sj->h; + const float hjg2 = hj * hj * kernel_gamma2; /* Early abort? */ if (sink_is_inhibited(sj, e)) continue; @@ -142,13 +142,13 @@ void DOSELF1_SINKS(struct runner *r, struct cell *c, int timer) { #ifdef SWIFT_DEBUG_CHECKS /* Check that particles have been drifted to the current time */ if (si->ti_drift != e->ti_current) - error("Particle bi not drifted to current time"); + error("Particle si not drifted to current time"); if (sj->ti_drift != e->ti_current) error("Particle bj not drifted to current time"); #endif - if (r2 < ri2 || r2 < rj2) { - IACT_SINKS_SINK(r2, dx, ri, rj, si, sj, with_cosmology, cosmo, + if (r2 < hig2 || r2 < hjg2) { + IACT_SINKS_SINK(r2, dx, hi, hj, si, sj, with_cosmology, cosmo, e->gravity_properties, e->sink_properties, e->ti_current, e->time); } @@ -210,8 +210,8 @@ void DO_NONSYM_PAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, /* Skip inactive particles */ if (!sink_is_active(si, e)) continue; - const float ri = si->r_cut; - const float ri2 = ri * ri; + const float hi = si->h; + const float hig2 = hi * hi * kernel_gamma2; const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])), (float)(si->x[1] - (cj->loc[1] + shift[1])), (float)(si->x[2] - (cj->loc[2] + shift[2]))}; @@ -241,8 +241,8 @@ void DO_NONSYM_PAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, error("Particle pj not drifted to current time"); #endif - if (r2 < ri2) { - IACT_SINKS_GAS(r2, dx, ri, hj, si, pj, with_cosmology, cosmo, + if (r2 < hig2) { + IACT_SINKS_GAS(r2, dx, hi, hj, si, pj, with_cosmology, cosmo, e->gravity_properties, e->sink_properties, e->ti_current, e->time); } @@ -266,8 +266,8 @@ void DO_NONSYM_PAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, /* Skip inactive particles */ if (!sink_is_active(si, e)) continue; - const float ri = si->r_cut; - const float ri2 = ri * ri; + const float hi = si->h; + const float hig2 = hi * hi * kernel_gamma2; const float six[3] = {(float)(si->x[0] - (cj->loc[0] + shift[0])), (float)(si->x[1] - (cj->loc[1] + shift[1])), (float)(si->x[2] - (cj->loc[2] + shift[2]))}; @@ -277,8 +277,8 @@ void DO_NONSYM_PAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, /* Get a pointer to the jth particle. */ struct sink *restrict sj = &sinks_j[sjd]; - const float rj = sj->r_cut; - const float rj2 = rj * rj; + const float hj = sj->h; + const float hjg2 = hj * hj * kernel_gamma2; /* Skip inhibited particles. */ if (sink_is_inhibited(sj, e)) continue; @@ -298,8 +298,8 @@ void DO_NONSYM_PAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, error("Particle sj not drifted to current time"); #endif - if (r2 < ri2 || r2 < rj2) { - IACT_SINKS_SINK(r2, dx, ri, rj, si, sj, with_cosmology, cosmo, + if (r2 < hig2 || r2 < hjg2) { + IACT_SINKS_SINK(r2, dx, hi, hj, si, sj, with_cosmology, cosmo, e->gravity_properties, e->sink_properties, e->ti_current, e->time); } @@ -337,6 +337,299 @@ void DOPAIR1_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, if (timer) TIMER_TOC(TIMER_DOPAIR_SINKS); } +/** + * @brief Compute the interactions between a cell pair, but only for the + * given indices in ci. + * + * Version using a brute-force algorithm. + * + * @param r The #runner. + * @param ci The first #cell. + * @param sinks_i The #sink to interact with @c cj. + * @param ind The list of indices of particles in @c ci to interact with. + * @param scount The number of particles in @c ind. + * @param cj The second #cell. + * @param shift The shift vector to apply to the particles in ci. + */ +void DOPAIR1_SUBSET_SINKS_NAIVE(struct runner *r, struct cell *restrict ci, + struct sink *restrict sinks_i, + int *restrict ind, const int scount, + struct cell *restrict cj, const double *shift) { + +#ifdef SWIFT_DEBUG_CHECKS + if (ci->nodeID != engine_rank) error("Should be run on a different node"); +#endif + + const struct engine *e = r->e; + const struct cosmology *cosmo = e->cosmology; + const int with_cosmology = e->policy & engine_policy_cosmology; + const int count_j = cj->hydro.count; + struct part *restrict parts_j = cj->hydro.parts; + + /* Early abort? */ + if (count_j == 0) return; + + /* Loop over the parts_i. */ + for (int sid = 0; sid < scount; sid++) { + + /* Get a hold of the ith part in ci. */ + struct sink *restrict si = &sinks_i[ind[sid]]; + + const double six = si->x[0] - (shift[0]); + const double siy = si->x[1] - (shift[1]); + const double siz = si->x[2] - (shift[2]); + const float hi = si->h; + const float hig2 = hi * hi * kernel_gamma2; + +#ifdef SWIFT_DEBUG_CHECKS + if (!sink_is_active(si, e)) + error("Trying to correct smoothing length of inactive particle !"); +#endif + + /* Loop over the parts in cj. */ + for (int pjd = 0; pjd < count_j; pjd++) { + + /* Get a pointer to the jth particle. */ + struct part *restrict pj = &parts_j[pjd]; + + /* Skip inhibited particles */ + if (part_is_inhibited(pj, e)) continue; + + const double pjx = pj->x[0]; + const double pjy = pj->x[1]; + const double pjz = pj->x[2]; + const float hj = pj->h; + + /* Compute the pairwise distance. */ + const float dx[3] = {(float)(six - pjx), (float)(siy - pjy), + (float)(siz - pjz)}; + const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; + +#ifdef SWIFT_DEBUG_CHECKS + /* Check that particles have been drifted to the current time */ + if (pj->ti_drift != e->ti_current) + error("Particle pj not drifted to current time"); +#endif + /* Hit or miss? */ + if (r2 < hig2) { + IACT_SINKS_GAS(r2, dx, hi, hj, si, pj, with_cosmology, cosmo, + e->gravity_properties, e->sink_properties, e->ti_current, + e->time); + } + } /* loop over the parts in cj. */ + } /* loop over the parts in ci. */ +} + +/** + * @brief Compute the interactions between a cell pair, but only for the + * given indices in ci. + * + * @param r The #runner. + * @param ci The first #cell. + * @param sinks The #sink to interact. + * @param ind The list of indices of particles in @c ci to interact with. + * @param scount The number of particles in @c ind. + */ +void DOSELF1_SUBSET_SINKS(struct runner *r, struct cell *restrict ci, + struct sink *restrict sinks, int *restrict ind, + const int scount) { + +#ifdef SWIFT_DEBUG_CHECKS + if (ci->nodeID != engine_rank) error("Should be run on a different node"); +#endif + + const struct engine *e = r->e; + const struct cosmology *cosmo = e->cosmology; + const int with_cosmology = e->policy & engine_policy_cosmology; + const int count_i = ci->hydro.count; + struct part *restrict parts_j = ci->hydro.parts; + + /* Early abort? */ + if (count_i == 0) return; + + /* Loop over the parts in ci. */ + for (int sid = 0; sid < scount; sid++) { + + /* Get a hold of the ith part in ci. */ + struct sink *si = &sinks[ind[sid]]; + const float six[3] = {(float)(si->x[0] - ci->loc[0]), + (float)(si->x[1] - ci->loc[1]), + (float)(si->x[2] - ci->loc[2])}; + const float hi = si->h; + const float hig2 = hi * hi * kernel_gamma2; + +#ifdef SWIFT_DEBUG_CHECKS + if (!sink_is_active(si, e)) error("Inactive particle in subset function!"); +#endif + + /* Loop over the parts in cj. */ + for (int pjd = 0; pjd < count_i; pjd++) { + + /* Get a pointer to the jth particle. */ + struct part *restrict pj = &parts_j[pjd]; + + /* Early abort? */ + if (part_is_inhibited(pj, e)) continue; + + /* Compute the pairwise distance. */ + const float pjx[3] = {(float)(pj->x[0] - ci->loc[0]), + (float)(pj->x[1] - ci->loc[1]), + (float)(pj->x[2] - ci->loc[2])}; + const float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]}; + const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; + +#ifdef SWIFT_DEBUG_CHECKS + /* Check that particles have been drifted to the current time */ + if (pj->ti_drift != e->ti_current) + error("Particle pj not drifted to current time"); +#endif + + /* Hit or miss? */ + if (r2 < hig2) { + IACT_SINKS_GAS(r2, dx, hi, pj->h, si, pj, with_cosmology, cosmo, + e->gravity_properties, e->sink_properties, e->ti_current, + e->time); + } + } /* loop over the parts in cj. */ + } /* loop over the parts in ci. */ +} + +/** + * @brief Determine which version of DOSELF1_SUBSET_SINKS needs to be called + * depending on the optimisation level. + * + * @param r The #runner. + * @param ci The first #cell. + * @param sinks The #sink to interact. + * @param ind The list of indices of particles in @c ci to interact with. + * @param scount The number of particles in @c ind. + */ +void DOSELF1_SUBSET_BRANCH_SINKS(struct runner *r, struct cell *restrict ci, + struct sink *restrict sinks, int *restrict ind, + const int scount) { + + DOSELF1_SUBSET_SINKS(r, ci, sinks, ind, scount); +} + +/** + * @brief Determine which version of DOPAIR1_SUBSET_SINKS needs to be called + * depending on the orientation of the cells or whether DOPAIR1_SUBSET_SINKS + * needs to be called at all. + * + * @param r The #runner. + * @param ci The first #cell. + * @param sinks_i The #sink to interact with @c cj. + * @param ind The list of indices of particles in @c ci to interact with. + * @param scount The number of particles in @c ind. + * @param cj The second #cell. + */ +void DOPAIR1_SUBSET_BRANCH_SINKS(struct runner *r, struct cell *restrict ci, + struct sink *restrict sinks_i, + int *restrict ind, int const scount, + struct cell *restrict cj) { + + const struct engine *e = r->e; + + /* Anything to do here? */ + if (cj->hydro.count == 0) return; + + /* Get the relative distance between the pairs, wrapping. */ + double shift[3] = {0.0, 0.0, 0.0}; + for (int k = 0; k < 3; k++) { + if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2) + shift[k] = e->s->dim[k]; + else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2) + shift[k] = -e->s->dim[k]; + } + + DOPAIR1_SUBSET_SINKS_NAIVE(r, ci, sinks_i, ind, scount, cj, shift); +} + +void DOSUB_SUBSET_SINKS(struct runner *r, struct cell *ci, struct sink *sinks, + int *ind, const int scount, struct cell *cj, + int gettimer) { + + const struct engine *e = r->e; + struct space *s = e->s; + + /* Should we even bother? */ + if (!cell_is_active_sinks(ci, e) && + (cj == NULL || !cell_is_active_sinks(cj, e))) + return; + + /* Find out in which sub-cell of ci the parts are. */ + struct cell *sub = NULL; + if (ci->split) { + for (int k = 0; k < 8; k++) { + if (ci->progeny[k] != NULL) { + if (&sinks[ind[0]] >= &ci->progeny[k]->sinks.parts[0] && + &sinks[ind[0]] < + &ci->progeny[k]->sinks.parts[ci->progeny[k]->sinks.count]) { + sub = ci->progeny[k]; + break; + } + } + } + } + + /* Is this a single cell? */ + if (cj == NULL) { + + /* Recurse? */ + if (cell_can_recurse_in_self_sinks_task(ci)) { + + /* Loop over all progeny. */ + DOSUB_SUBSET_SINKS(r, sub, sinks, ind, scount, NULL, 0); + for (int j = 0; j < 8; j++) + if (ci->progeny[j] != sub && ci->progeny[j] != NULL) + DOSUB_SUBSET_SINKS(r, sub, sinks, ind, scount, ci->progeny[j], 0); + + } + + /* Otherwise, compute self-interaction. */ + else + DOSELF1_SUBSET_BRANCH_SINKS(r, ci, sinks, ind, scount); + } /* self-interaction. */ + + /* Otherwise, it's a pair interaction. */ + else { + + /* Recurse? */ + if (cell_can_recurse_in_pair_sinks_task(ci, cj) && + cell_can_recurse_in_pair_sinks_task(cj, ci)) { + + /* Get the type of pair and flip ci/cj if needed. */ + double shift[3] = {0.0, 0.0, 0.0}; + const int sid = space_getsid_and_swap_cells(s, &ci, &cj, shift); + + struct cell_split_pair *csp = &cell_split_pairs[sid]; + for (int k = 0; k < csp->count; k++) { + const int pid = csp->pairs[k].pid; + const int pjd = csp->pairs[k].pjd; + if (ci->progeny[pid] == sub && cj->progeny[pjd] != NULL) + DOSUB_SUBSET_SINKS(r, ci->progeny[pid], sinks, ind, scount, + cj->progeny[pjd], 0); + if (ci->progeny[pid] != NULL && cj->progeny[pjd] == sub) + DOSUB_SUBSET_SINKS(r, cj->progeny[pjd], sinks, ind, scount, + ci->progeny[pid], 0); + } + } + + /* Otherwise, compute the pair directly. */ + else if (cell_is_active_sinks(ci, e) && cj->hydro.count > 0) { + + /* Do any of the cells need to be drifted first? */ + if (cell_is_active_sinks(ci, e)) { + if (!cell_are_sink_drifted(ci, e)) error("Cell should be drifted!"); + if (!cell_are_part_drifted(cj, e)) error("Cell should be drifted!"); + } + + DOPAIR1_SUBSET_BRANCH_SINKS(r, ci, sinks, ind, scount, cj); + } + + } /* otherwise, pair interaction. */ +} + /** * @brief Wrapper to runner_doself_sinks_swallow * @@ -359,8 +652,8 @@ void DOSELF1_BRANCH_SINKS(struct runner *r, struct cell *c) { if (!cell_is_active_sinks(c, e)) return; /* Did we mess up the recursion? */ - if (c->sinks.r_cut_max_old > c->dmin) - error("Cell smaller than the cut off radius"); + if (c->sinks.h_max_old * kernel_gamma > c->dmin) + error("Cell smaller than the cut off radius or smoothing length"); DOSELF1_SINKS(r, c, 1); } diff --git a/src/runner_doiact_sinks.h b/src/runner_doiact_sinks.h index 8195924b5e..679db023aa 100644 --- a/src/runner_doiact_sinks.h +++ b/src/runner_doiact_sinks.h @@ -38,6 +38,27 @@ #define _DOPAIR1_SINKS_NAIVE(f) PASTE(runner_dopair_sinks_naive, f) #define DOPAIR1_SINKS_NAIVE _DOPAIR1_SINKS_NAIVE(FUNCTION) +#define _DOPAIR1_SUBSET_SINKS(f) PASTE(runner_dopair_subset_sinks, f) +#define DOPAIR1_SUBSET_SINKS _DOPAIR1_SUBSET_SINKS(FUNCTION) + +#define _DOPAIR1_SUBSET_SINKS_NAIVE(f) \ + PASTE(runner_dopair_subset_sinks_naive, f) +#define DOPAIR1_SUBSET_SINKS_NAIVE _DOPAIR1_SUBSET_SINKS_NAIVE(FUNCTION) + +#define _DOSELF1_SUBSET_SINKS(f) PASTE(runner_doself_subset_sinks, f) +#define DOSELF1_SUBSET_SINKS _DOSELF1_SUBSET_SINKS(FUNCTION) + +#define _DOSELF1_SUBSET_BRANCH_SINKS(f) \ + PASTE(runner_doself_subset_branch_sinks, f) +#define DOSELF1_SUBSET_BRANCH_SINKS _DOSELF1_SUBSET_BRANCH_SINKS(FUNCTION) + +#define _DOPAIR1_SUBSET_BRANCH_SINKS(f) \ + PASTE(runner_dopair_subset_branch_sinks, f) +#define DOPAIR1_SUBSET_BRANCH_SINKS _DOPAIR1_SUBSET_BRANCH_SINKS(FUNCTION) + +#define _DOSUB_SUBSET_SINKS(f) PASTE(runner_dosub_subset_sinks, f) +#define DOSUB_SUBSET_SINKS _DOSUB_SUBSET_SINKS(FUNCTION) + #define _DOSELF1_BRANCH_SINKS(f) PASTE(runner_doself_branch_sinks, f) #define DOSELF1_BRANCH_SINKS _DOSELF1_BRANCH_SINKS(FUNCTION) @@ -74,3 +95,15 @@ void DOPAIR1_BRANCH_SINKS(struct runner *r, struct cell *ci, struct cell *cj); void DOSUB_SELF1_SINKS(struct runner *r, struct cell *ci, int gettimer); void DOSUB_PAIR1_SINKS(struct runner *r, struct cell *ci, struct cell *cj, int gettimer); + +void DOSELF1_SUBSET_BRANCH_SINKS(struct runner *r, struct cell *restrict ci, + struct sink *restrict sinks, int *restrict ind, + const int scount); +void DOPAIR1_SUBSET_BRANCH_SINKS(struct runner *r, struct cell *restrict ci, + struct sink *restrict sinks_i, + int *restrict ind, int const scount, + struct cell *restrict cj); + +void DOSUB_SUBSET_SINKS(struct runner *r, struct cell *ci, struct sink *sinks, + int *ind, const int scount, struct cell *cj, + int gettimer); \ No newline at end of file diff --git a/src/runner_ghost.c b/src/runner_ghost.c index 38aaf9fd52..9929331a3e 100644 --- a/src/runner_ghost.c +++ b/src/runner_ghost.c @@ -63,6 +63,13 @@ #undef FUNCTION_TASK_LOOP #undef FUNCTION +/* Import the sink density loop functions. */ +#define FUNCTION density +#define FUNCTION_TASK_LOOP TASK_LOOP_DENSITY +#include "runner_doiact_sinks.h" +#undef FUNCTION_TASK_LOOP +#undef FUNCTION + /** * @brief Intermediate task after the density to check that the smoothing * lengths are correct. @@ -1710,8 +1717,9 @@ void runner_do_rt_ghost2(struct runner *r, struct cell *c, int timer) { } /** - * @brief Intermediate task after the density to finish density calculation - * and calculate accretion rates for the particle swallowing step + * @brief Intermediate task after the density to check that the smoothing + * lengths are correct, finish density calculations + * and calculate accretion rates for the particle swallowing step * * @param r The runner thread. * @param c The cell. @@ -1724,6 +1732,16 @@ void runner_do_sinks_density_ghost(struct runner *r, struct cell *c, const struct engine *e = r->e; const struct cosmology *cosmo = e->cosmology; const int with_cosmology = e->policy & engine_policy_cosmology; + const float sinks_h_max = e->hydro_properties->h_max; + const float sinks_h_min = e->hydro_properties->h_min; + const float eps = e->sink_properties->h_tolerance; + const float sinks_eta_dim = pow_dimension(e->sink_properties->eta_neighbours); + const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations; + int redo = 0, scount = 0; + + /* Running value of the maximal smoothing length */ + float h_max = c->sinks.h_max; + float h_max_active = c->sinks.h_max_active; TIMER_TIC; @@ -1736,58 +1754,356 @@ void runner_do_sinks_density_ghost(struct runner *r, struct cell *c, for (int k = 0; k < 8; k++) { if (c->progeny[k] != NULL) { runner_do_sinks_density_ghost(r, c->progeny[k], 0); + + /* Update h_max */ + h_max = max(h_max, c->progeny[k]->sinks.h_max); + h_max_active = max(h_max_active, c->progeny[k]->sinks.h_max_active); } } } else { /* Init the list of active particles that have to be updated. */ int *sid = NULL; + float *h_0 = NULL; + float *left = NULL; + float *right = NULL; if ((sid = (int *)malloc(sizeof(int) * c->sinks.count)) == NULL) error("Can't allocate memory for sid."); - - int scount = 0; + if ((h_0 = (float *)malloc(sizeof(float) * c->sinks.count)) == NULL) + error("Can't allocate memory for h_0."); + if ((left = (float *)malloc(sizeof(float) * c->sinks.count)) == NULL) + error("Can't allocate memory for left."); + if ((right = (float *)malloc(sizeof(float) * c->sinks.count)) == NULL) + error("Can't allocate memory for right."); for (int k = 0; k < c->sinks.count; k++) if (sink_is_active(&sinks[k], e)) { sid[scount] = k; + h_0[scount] = sinks[k].h; + left[scount] = 0.f; + right[scount] = sinks_h_max; ++scount; } - /* Loop over the remaining active parts in this cell. */ - for (int i = 0; i < scount; i++) { + if (e->sink_properties->use_fixed_r_cut) { + /* If we're using a fixed cutoff rather than a smoothing length, just + * finish up the density task and leave sp->h untouched. */ - /* Get a direct pointer on the part. */ - struct sink *sp = &sinks[sid[i]]; + /* Loop over the active sinks in this cell. */ + for (int i = 0; i < scount; i++) { + + /* Get a direct pointer on the part. */ + struct sink *sp = &sinks[sid[i]]; #ifdef SWIFT_DEBUG_CHECKS - /* Is this part within the timestep? */ - if (!sink_is_active(sp, e)) error("Ghost applied to inactive particle"); + /* Is this part within the timestep? */ + if (!sink_is_active(sp, e)) error("Ghost applied to inactive particle"); #endif - /* Finish the density calculation */ - sink_end_density(sp, cosmo); + /* Finish the density calculation */ + sink_end_density(sp, cosmo); - if (sp->num_ngbs == 0) { - sinks_sink_has_no_neighbours(sp, cosmo); + /* Set these variables to the fixed cutoff radius for the rest of the + * ghost task */ + h_max = sp->h; + h_max_active = sp->h; } - /* Get particle time-step */ - double dt; - if (with_cosmology) { - const integertime_t ti_step = get_integer_timestep(sp->time_bin); - const integertime_t ti_begin = - get_integer_time_begin(e->ti_current - 1, sp->time_bin); - - dt = cosmology_get_delta_time(e->cosmology, ti_begin, - ti_begin + ti_step); - } else { - dt = get_timestep(sp->time_bin, e->time_base); + } else { + /* Otherwise we need to iterate to update the smoothing lengths */ + + /* While there are particles that need to be updated... */ + for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter; + num_reruns++) { + + /* Reset the redo-count. */ + redo = 0; + + /* Loop over the remaining active parts in this cell. */ + for (int i = 0; i < scount; i++) { + + /* Get a direct pointer on the part. */ + struct sink *sp = &sinks[sid[i]]; + +#ifdef SWIFT_DEBUG_CHECKS + /* Is this part within the timestep? */ + if (!sink_is_active(sp, e)) + error("Ghost applied to inactive particle"); +#endif + + /* Get some useful values */ + const float h_init = h_0[i]; + const float h_old = sp->h; + const float h_old_dim = pow_dimension(h_old); + const float h_old_dim_minus_one = pow_dimension_minus_one(h_old); + + float h_new; + int has_no_neighbours = 0; + + if (sp->density.wcount < + 1.e-5 * kernel_root) { /* No neighbours case */ + + /* Flag that there were no neighbours */ + has_no_neighbours = 1; + + /* Double h and try again */ + h_new = 2.f * h_old; + + } else { + + /* Finish the density calculation */ + sink_end_density(sp, cosmo); + + /* Compute one step of the Newton-Raphson scheme */ + const float n_sum = sp->density.wcount * h_old_dim; + const float n_target = sinks_eta_dim; + const float f = n_sum - n_target; + const float f_prime = + sp->density.wcount_dh * h_old_dim + + hydro_dimension * sp->density.wcount * h_old_dim_minus_one; + + /* Improve the bisection bounds */ + if (n_sum < n_target) + left[i] = max(left[i], h_old); + else if (n_sum > n_target) + right[i] = min(right[i], h_old); + +#ifdef SWIFT_DEBUG_CHECKS + /* Check the validity of the left and right bounds */ + if (left[i] > right[i]) + error("Invalid left (%e) and right (%e)", left[i], right[i]); +#endif + + /* Skip if h is already h_max and we don't have enough neighbours + */ + /* Same if we are below h_min */ + if (((sp->h >= sinks_h_max) && (f < 0.f)) || + ((sp->h <= sinks_h_min) && (f > 0.f))) { + + /* Ok, we are done with this particle */ + continue; + } + + /* Normal case: Use Newton-Raphson to get a better value of h */ + + /* Avoid floating point exception from f_prime = 0 */ + h_new = h_old - f / (f_prime + FLT_MIN); + + /* Be verbose about the particles that struggle to converge */ + if (num_reruns > max_smoothing_iter - 10) { + + message( + "Smoothing length convergence problem: iter=%d p->id=%lld " + "h_init=%12.8e h_old=%12.8e h_new=%12.8e f=%f f_prime=%f " + "n_sum=%12.8e n_target=%12.8e left=%12.8e right=%12.8e", + num_reruns, sp->id, h_init, h_old, h_new, f, f_prime, n_sum, + n_target, left[i], right[i]); + } + + /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */ + h_new = min(h_new, 2.f * h_old); + h_new = max(h_new, 0.5f * h_old); + + /* Verify that we are actually progrssing towards the answer */ + h_new = max(h_new, left[i]); + h_new = min(h_new, right[i]); + } + + /* Check whether the particle has an inappropriate smoothing length + */ + if (fabsf(h_new - h_old) > eps * h_old) { + + /* Ok, correct then */ + + /* Case where we have been oscillating around the solution */ + if ((h_new == left[i] && h_old == right[i]) || + (h_old == left[i] && h_new == right[i])) { + + /* Bisect the remaining interval */ + sp->h = pow_inv_dimension( + 0.5f * (pow_dimension(left[i]) + pow_dimension(right[i]))); + + } else { + + /* Normal case */ + sp->h = h_new; + } + + /* If below the absolute maximum, try again */ + if (sp->h < sinks_h_max && sp->h > sinks_h_min) { + + /* Flag for another round of fun */ + sid[redo] = sid[i]; + h_0[redo] = h_0[i]; + left[redo] = left[i]; + right[redo] = right[i]; + redo += 1; + + /* Re-initialise everything */ + sink_init_sink(sp); + + /* Off we go ! */ + continue; + + } else if (sp->h <= sinks_h_min) { + + /* Ok, this particle is a lost cause... */ + sp->h = sinks_h_min; + + } else if (sp->h >= sinks_h_max) { + + /* Ok, this particle is a lost cause... */ + sp->h = sinks_h_max; + + /* Do some damage control if no neighbours at all were found */ + if (has_no_neighbours) { + sinks_sink_has_no_neighbours(sp, cosmo); + } + + } else { + error( + "Fundamental problem with the smoothing length iteration " + "logic."); + } + } + + /* We now have a particle whose smoothing length has converged */ + + /* Check if h_max has increased */ + h_max = max(h_max, sp->h); + h_max_active = max(h_max_active, sp->h); + } + + /* We now need to treat the particles whose smoothing length had not + * converged again */ + + /* Re-set the counter for the next loop (potentially). */ + scount = redo; + if (scount > 0) { + + /* Climb up the cell hierarchy. */ + for (struct cell *finger = c; finger != NULL; + finger = finger->parent) { + + /* Run through this cell's density interactions. */ + for (struct link *l = finger->sinks.density; l != NULL; + l = l->next) { + +#ifdef SWIFT_DEBUG_CHECKS + if (l->t->ti_run < r->e->ti_current) + error("Density task should have been run."); +#endif + + /* Self-interaction? */ + if (l->t->type == task_type_self) + runner_doself_subset_branch_sinks_density(r, finger, sinks, sid, + scount); + + /* Otherwise, pair interaction? */ + else if (l->t->type == task_type_pair) { + + /* Left or right? */ + if (l->t->ci == finger) + runner_dopair_subset_branch_sinks_density( + r, finger, sinks, sid, scount, l->t->cj); + else + runner_dopair_subset_branch_sinks_density( + r, finger, sinks, sid, scount, l->t->ci); + } + + /* Otherwise, sub-self interaction? */ + else if (l->t->type == task_type_sub_self) + runner_dosub_subset_sinks_density(r, finger, sinks, sid, scount, + NULL, 1); + + /* Otherwise, sub-pair interaction? */ + else if (l->t->type == task_type_sub_pair) { + + /* Left or right? */ + if (l->t->ci == finger) + runner_dosub_subset_sinks_density(r, finger, sinks, sid, + scount, l->t->cj, 1); + else + runner_dosub_subset_sinks_density(r, finger, sinks, sid, + scount, l->t->ci, 1); + } + } + } + } + } + + if (scount) { + warning( + "Smoothing length failed to converge for the following sink " + "particles:"); + for (int i = 0; i < scount; i++) { + struct sink *sp = &sinks[sid[i]]; + warning("ID: %lld, h: %g, wcount: %g", sp->id, sp->h, + sp->density.wcount); + } + + error("Smoothing length failed to converge on %i particles.", scount); } - /* Calculate the accretion rate and accreted mass this timestep, for use - * in swallow loop */ - sink_prepare_swallow(sp, e->sink_properties, e->physical_constants, - e->cosmology, e->cooling_func, e->entropy_floor, - e->time, with_cosmology, dt, e->ti_current); + /* Be clean */ + free(left); + free(right); + free(sid); + free(h_0); + } + + /* We need one more quick loop over the sinks to run prepare_swallow */ + for (int i = 0; i < c->sinks.count; i++) { + + /* Get a direct pointer on the part. */ + struct sink *sp = &sinks[i]; + + if (sink_is_active(sp, e)) { + + /* Get particle time-step */ + double dt; + if (with_cosmology) { + const integertime_t ti_step = get_integer_timestep(sp->time_bin); + const integertime_t ti_begin = + get_integer_time_begin(e->ti_current - 1, sp->time_bin); + + dt = cosmology_get_delta_time(e->cosmology, ti_begin, + ti_begin + ti_step); + } else { + dt = get_timestep(sp->time_bin, e->time_base); + } + + /* Calculate the accretion rate and accreted mass this timestep, for use + * in swallow loop */ + sink_prepare_swallow(sp, e->sink_properties, e->physical_constants, + e->cosmology, e->cooling_func, e->entropy_floor, + e->time, with_cosmology, dt, e->ti_current); + } + } + } + + /* Update h_max */ + c->sinks.h_max = h_max; + c->sinks.h_max_active = h_max_active; + +#ifdef SWIFT_DEBUG_CHECKS + for (int i = 0; i < c->sinks.count; ++i) { + const struct sink *sp = &c->sinks.parts[i]; + const float h = c->sinks.parts[i].h; + if (sink_is_inhibited(sp, e)) continue; + + if (h > c->sinks.h_max) + error("Particle has h larger than h_max (id=%lld)", sp->id); + if (sink_is_active(sp, e) && h > c->sinks.h_max_active) + error("Active particle has h larger than h_max_active (id=%lld)", sp->id); + } +#endif + + /* The ghost may not always be at the top level. + * Therefore we need to update h_max between the super- and top-levels */ + if (c->sinks.density_ghost) { + for (struct cell *tmp = c->parent; tmp != NULL; tmp = tmp->parent) { + atomic_max_f(&tmp->sinks.h_max, h_max); + atomic_max_f(&tmp->sinks.h_max_active, h_max_active); } } diff --git a/src/sink.c b/src/sink.c new file mode 100644 index 0000000000..c3adf40237 --- /dev/null +++ b/src/sink.c @@ -0,0 +1,341 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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/>. + * + ******************************************************************************/ + +/* Config parameters. */ +#include <config.h> + +/* Some standard headers. */ +#include <float.h> +#include <stdio.h> +#include <stdlib.h> +#include <unistd.h> + +/* Local headers. */ +#include "active.h" +#include "error.h" +#include "sink_properties.h" +#include "version.h" + +struct exact_density_data { + const struct engine *e; + const struct space *s; + int counter_global; +}; + +/** + * @brief Mapper function for the exact sink checks. + * + * @brief map_data The #sinks. + * @brief nr_sinks The number of star particles. + * @brief extra_data Pointers to the structure containing global interaction + * counters. + */ +void sink_exact_density_compute_mapper(void *map_data, int nr_sinks, + void *extra_data) { +#ifdef SWIFT_SINK_DENSITY_CHECKS + + /* Unpack the data */ + struct sink *restrict sinks = (struct sink *)map_data; + struct exact_density_data *data = (struct exact_density_data *)extra_data; + const struct space *s = data->s; + const struct engine *e = data->e; + const int periodic = s->periodic; + const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; + int counter = 0; + + for (int i = 0; i < nr_sinks; ++i) { + + struct sink *si = &sinks[i]; + const long long id = si->id; + + /* Is the particle active and part of the subset to be tested ? */ + if (id % SWIFT_SINK_DENSITY_CHECKS == 0 && sink_is_starting(si, e)) { + + /* Get some information about the particle */ + const double pix[3] = {si->x[0], si->x[1], si->x[2]}; + const double hi = si->h; + const float hi_inv = 1.f / hi; + const float hig2 = hi * hi * kernel_gamma2; + + /* Be ready for the calculation */ + int N_density_exact = 0; + double rho_exact = 0.; + double n_exact = 0.; + + /* Interact it with all other particles in the space.*/ + for (int j = 0; j < (int)s->nr_parts; ++j) { + + const struct part *pj = &s->parts[j]; + + /* Compute the pairwise distance. */ + double dx = pj->x[0] - pix[0]; + double dy = pj->x[1] - pix[1]; + double dz = pj->x[2] - pix[2]; + + /* Now apply periodic BC */ + if (periodic) { + dx = nearest(dx, dim[0]); + dy = nearest(dy, dim[1]); + dz = nearest(dz, dim[2]); + } + + const double r2 = dx * dx + dy * dy + dz * dz; + + /* Interact loop of type 1? */ + if (r2 < hig2) { + + const float mj = pj->mass; + + float wi, wi_dx; + + /* Kernel function */ + const float r = sqrtf(r2); + const float ui = r * hi_inv; + kernel_deval(ui, &wi, &wi_dx); + + /* Flag that we found an inhibited neighbour */ + if (part_is_inhibited(pj, e)) { + si->inhibited_check_exact = 1; + } else { + + /* Density */ + rho_exact += mj * wi; + + /* Number density */ + n_exact += wi; + + /* Number of neighbours */ + N_density_exact++; + } + } + } + + /* Store the exact answer */ + si->N_check_density_exact = N_density_exact; + si->rho_check_exact = rho_exact * pow_dimension(hi_inv); + si->n_check_exact = n_exact * pow_dimension(hi_inv); + + counter++; + } + } + atomic_add(&data->counter_global, counter); + +#else + error("Sink checking function called without the corresponding flag."); +#endif +} + +/** + * @brief Compute the exact interactions for a selection of star particles + * by running a brute force loop over all the particles in the simulation. + * + * Will be incorrect over MPI. + * + * @param s The #space. + * @param e The #engine. + */ +void sink_exact_density_compute(struct space *s, const struct engine *e) { + +#ifdef SWIFT_SINK_DENSITY_CHECKS + + const ticks tic = getticks(); + + struct exact_density_data data; + data.e = e; + data.s = s; + data.counter_global = 0; + + threadpool_map(&s->e->threadpool, sink_exact_density_compute_mapper, s->sinks, + s->nr_sinks, sizeof(struct sink), 0, &data); + + if (e->verbose) + message("Computed exact densities for %d sinks (took %.3f %s). ", + data.counter_global, clocks_from_ticks(getticks() - tic), + clocks_getunit()); +#else + error("Sink checking function called without the corresponding flag."); +#endif +} + +/** + * @brief Check the star particles' density and force calculations against the + * values obtained via the brute-force summation. + * + * @param s The #space. + * @param e The #engine. + * @param rel_tol Relative tolerance for the checks + */ +void sink_exact_density_check(struct space *s, const struct engine *e, + const double rel_tol) { + +#ifdef SWIFT_SINK_DENSITY_CHECKS + + const ticks tic = getticks(); + + const struct sink *sinks = s->sinks; + const size_t nr_sinks = s->nr_sinks; + + const double eta = e->sink_properties->eta_neighbours; + const double N_ngb_target = + (4. / 3.) * M_PI * pow_dimension(kernel_gamma * eta); + const double N_ngb_max = + N_ngb_target + 2. * e->sink_properties->delta_neighbours; + const double N_ngb_min = + N_ngb_target - 2. * e->sink_properties->delta_neighbours; + + /* File name */ + char file_name_swift[100]; + sprintf(file_name_swift, "sink_checks_swift_step%.4d.dat", e->step); + + /* Creare files and write header */ + FILE *file_swift = fopen(file_name_swift, "w"); + if (file_swift == NULL) error("Could not create file '%s'.", file_name_swift); + fprintf(file_swift, "# Sink accuracy test - SWIFT DENSITIES\n"); + fprintf(file_swift, "# N= %d\n", SWIFT_SINK_DENSITY_CHECKS); + fprintf(file_swift, "# periodic= %d\n", s->periodic); + fprintf(file_swift, "# N_ngb_target= %f +/- %f\n", N_ngb_target, + e->sink_properties->delta_neighbours); + fprintf(file_swift, "# Git Branch: %s\n", git_branch()); + fprintf(file_swift, "# Git Revision: %s\n", git_revision()); + fprintf(file_swift, "# %16s %16s %16s %16s %16s %7s %7s %16s %16s %16s\n", + "id", "pos[0]", "pos[1]", "pos[2]", "h", "Nd", "Nf", "rho", "n_rho", + "N_ngb"); + + /* Output particle SWIFT densities */ + for (size_t i = 0; i < nr_sinks; ++i) { + + const struct sink *si = &sinks[i]; + const long long id = si->id; + + const double N_ngb = (4. / 3.) * M_PI * kernel_gamma * kernel_gamma * + kernel_gamma * si->h * si->h * si->h * si->n_check; + + if (id % SWIFT_SINK_DENSITY_CHECKS == 0 && sink_is_starting(si, e)) { + + fprintf( + file_swift, + "%18lld %16.8e %16.8e %16.8e %16.8e %7d %7d %16.8e %16.8e %16.8e\n", + id, si->x[0], si->x[1], si->x[2], si->h, si->N_check_density, 0, + si->rho_check, si->n_check, N_ngb); + } + } + + if (e->verbose) + message("Written SWIFT densities in file '%s'.", file_name_swift); + + /* Be nice */ + fclose(file_swift); + + /* File name */ + char file_name_exact[100]; + sprintf(file_name_exact, "sink_checks_exact_step%.4d.dat", e->step); + + /* Creare files and write header */ + FILE *file_exact = fopen(file_name_exact, "w"); + if (file_exact == NULL) error("Could not create file '%s'.", file_name_exact); + fprintf(file_exact, "# Sink accuracy test - EXACT DENSITIES\n"); + fprintf(file_exact, "# N= %d\n", SWIFT_SINK_DENSITY_CHECKS); + fprintf(file_exact, "# periodic= %d\n", s->periodic); + fprintf(file_exact, "# N_ngb_target= %f +/- %f\n", N_ngb_target, + e->sink_properties->delta_neighbours); + fprintf(file_exact, "# Git Branch: %s\n", git_branch()); + fprintf(file_exact, "# Git Revision: %s\n", git_revision()); + fprintf(file_exact, "# %16s %16s %16s %16s %16s %7s %7s %16s %16s %16s\n", + "id", "pos[0]", "pos[1]", "pos[2]", "h", "Nd", "Nf", "rho_exact", + "n_rho_exact", "N_ngb"); + + int wrong_rho = 0; + int wrong_n_ngb = 0; + int counter = 0; + + /* Output particle SWIFT densities */ + for (size_t i = 0; i < nr_sinks; ++i) { + + const struct sink *si = &sinks[i]; + const long long id = si->id; + const int found_inhibited = si->inhibited_check_exact; + + const double N_ngb = (4. / 3.) * M_PI * kernel_gamma * kernel_gamma * + kernel_gamma * si->h * si->h * si->h * + si->n_check_exact; + + if (id % SWIFT_SINK_DENSITY_CHECKS == 0 && sink_is_starting(si, e)) { + + counter++; + + fprintf( + file_exact, + "%18lld %16.8e %16.8e %16.8e %16.8e %7d %7d %16.8e %16.8e %16.8e\n", + id, si->x[0], si->x[1], si->x[2], si->h, si->N_check_density_exact, 0, + si->rho_check_exact, si->n_check_exact, N_ngb); + + /* Check that we did not go above the threshold. + * Note that we ignore particles that saw an inhibted particle as a + * neighbour as we don't know whether that neighbour became inhibited in + * that step or not. */ + if (!found_inhibited && + si->N_check_density_exact != si->N_check_density && + (fabsf(si->rho_check / si->rho_check_exact - 1.f) > rel_tol || + fabsf(si->rho_check_exact / si->rho_check - 1.f) > rel_tol)) { + message("RHO: id=%lld swift=%e exact=%e N_swift=%d N_true=%d", id, + si->rho_check, si->rho_check_exact, si->N_check_density, + si->N_check_density_exact); + wrong_rho++; + } + + if (!found_inhibited && (N_ngb > N_ngb_max || N_ngb < N_ngb_min)) { + + message("N_NGB: id=%lld exact=%f N_true=%d N_swift=%d", id, N_ngb, + si->N_check_density_exact, si->N_check_density); + + wrong_n_ngb++; + } + } + } + + if (e->verbose) + message("Written exact densities in file '%s'.", file_name_exact); + + /* Be nice */ + fclose(file_exact); + + if (wrong_rho) + error( + "Density difference larger than the allowed tolerance for %d " + "sink particles! (out of %d particles)", + wrong_rho, counter); + else + message("Verified %d sink particles", counter); + + /* if (wrong_n_ngb) */ + /* error( */ + /* "N_ngb difference larger than the allowed tolerance for %d " */ + /* "star particles! (out of %d particles)", */ + /* wrong_n_ngb, counter); */ + /* else */ + /* message("Verified %d star particles", counter); */ + + if (e->verbose) + message("Writting brute-force density files took %.3f %s. ", + clocks_from_ticks(getticks() - tic), clocks_getunit()); + +#else + error("Sink checking function called without the corresponding flag."); +#endif +} diff --git a/src/sink.h b/src/sink.h index 128eaff03e..e891283826 100644 --- a/src/sink.h +++ b/src/sink.h @@ -25,10 +25,19 @@ /* Select the correct sink model */ #if defined(SINK_NONE) #include "./sink/Default/sink.h" +#elif defined(SINK_BASIC) +#include "./sink/Basic/sink.h" #elif defined(SINK_GEAR) #include "./sink/GEAR/sink.h" #else #error "Invalid choice of sink model" #endif +struct engine; +struct space; + +void sink_exact_density_compute(struct space *s, const struct engine *e); +void sink_exact_density_check(struct space *s, const struct engine *e, + const double rel_tol); + #endif /* SWIFT_SINK_H */ diff --git a/src/sink/Basic/sink.h b/src/sink/Basic/sink.h new file mode 100644 index 0000000000..238a2bf477 --- /dev/null +++ b/src/sink/Basic/sink.h @@ -0,0 +1,644 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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_BASIC_SINK_H +#define SWIFT_BASIC_SINK_H + +#include <float.h> + +/* Put pragma if gsl around here */ +#ifdef HAVE_LIBGSL +#include <gsl/gsl_cdf.h> +#endif + +/* Local includes */ +#include "active.h" +#include "chemistry.h" +#include "cooling.h" +#include "feedback.h" +#include "minmax.h" +#include "random.h" +#include "sink_part.h" +#include "sink_properties.h" +#include "star_formation.h" + +/** + * @brief Computes the time-step of a given sink particle. + * + * @param sp Pointer to the sink-particle data. + */ +__attribute__((always_inline)) INLINE static float sink_compute_timestep( + const struct sink* const sink, const struct sink_props* sink_properties, + const int with_cosmology, const struct cosmology* cosmo, + const struct gravity_props* grav_props, const double time, + const double time_base) { + + return FLT_MAX; +} + +/** + * @brief Initialises the sink-particles for the first time + * + * This function is called only once just after the ICs have been + * read in to do some conversions. + * + * @param sp The #sink particle to act upon. + * @param sink_props The properties of the sink particles scheme. + * @param e The #engine + */ +__attribute__((always_inline)) INLINE static void sink_first_init_sink( + struct sink* sp, const struct sink_props* sink_props, + const struct engine* e) { + + sp->time_bin = 0; + + sp->number_of_gas_swallows = 0; + sp->number_of_direct_gas_swallows = 0; + sp->number_of_sink_swallows = 0; + sp->number_of_direct_sink_swallows = 0; + sp->swallowed_angular_momentum[0] = 0.f; + sp->swallowed_angular_momentum[1] = 0.f; + sp->swallowed_angular_momentum[2] = 0.f; + + /* Initially set the subgrid mass equal to the dynamical mass read from the + * ICs. */ + sp->subgrid_mass = sp->mass; + + sink_mark_sink_as_not_swallowed(&sp->merger_data); +} + +/** + * @brief Initialisation of particle data before the hydro density loop. + * Note: during initalisation (space_init) + * + * @param p The #part to act upon. + * @param sink_props The properties of the sink particles scheme. + */ +__attribute__((always_inline)) INLINE static void sink_init_part( + struct part* restrict p, const struct sink_props* sink_props) {} + +/** + * @brief Initialisation of sink particle data before sink loops. + * Note: during initalisation (space_init_sinks) + * + * @param sp The #sink particle to act upon. + */ +__attribute__((always_inline)) INLINE static void sink_init_sink( + struct sink* sp) { + + sp->density.wcount = 0.f; + sp->density.wcount_dh = 0.f; + sp->rho_gas = 0.f; + sp->sound_speed_gas = 0.f; + sp->velocity_gas[0] = 0.f; + sp->velocity_gas[1] = 0.f; + sp->velocity_gas[2] = 0.f; + sp->ngb_mass = 0.f; + sp->num_ngbs = 0; + sp->accretion_rate = 0.f; + sp->mass_at_start_of_step = sp->mass; /* sp->mass may grow in nibbling mode */ + +#ifdef DEBUG_INTERACTIONS_SINKS + for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_SINKS; ++i) + sp->ids_ngbs_accretion[i] = -1; + sp->num_ngb_accretion = 0; + + for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_SINKS; ++i) + sp->ids_ngbs_merger[i] = -1; + sp->num_ngb_merger = 0; + + for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_SINKS; ++i) + sp->ids_ngbs_formation[i] = -1; + sp->num_ngb_formation = 0; +#endif + +#ifdef SWIFT_SINK_DENSITY_CHECKS + sp->N_check_density = 0; + sp->N_check_density_exact = 0; + sp->rho_check = 0.f; + sp->rho_check_exact = 0.f; + sp->n_check = 0.f; + sp->n_check_exact = 0.f; + sp->inhibited_check_exact = 0; +#endif +} + +/** + * @brief Predict additional particle fields forward in time when drifting + * + * @param sp The #sink. + * @param dt_drift The drift time-step for positions. + */ +__attribute__((always_inline)) INLINE static void sink_predict_extra( + struct sink* restrict sp, float dt_drift) {} + +/** + * @brief Sets the values to be predicted in the drifts to their values at a + * kick time + * + * @param sp The #sink particle. + */ +__attribute__((always_inline)) INLINE static void sink_reset_predicted_values( + struct sink* restrict sp) {} + +/** + * @brief Kick the additional variables + * + * @param sp The #sink particle to act upon + * @param dt The time-step for this kick + */ +__attribute__((always_inline)) INLINE static void sink_kick_extra( + struct sink* sp, float dt) {} + +/** + * @brief Finishes the calculation of density on sinks + * + * @param si The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void sink_end_density( + struct sink* si, const struct cosmology* cosmo) { + + /* Some smoothing length multiples. */ + const float h = si->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */ + + /* Finish the calculation by inserting the missing h-factors */ + si->density.wcount *= h_inv_dim; + si->density.wcount_dh *= h_inv_dim_plus_one; + + /* Finish the density calculation */ + si->rho_gas *= h_inv_dim; + + /* For the following, we also have to undo the mass smoothing + * (N.B.: bp->velocity_gas is in BH frame, in internal units). */ + const float rho_inv = 1.f / si->rho_gas; + si->sound_speed_gas *= h_inv_dim * rho_inv; + si->velocity_gas[0] *= h_inv_dim * rho_inv; + si->velocity_gas[1] *= h_inv_dim * rho_inv; + si->velocity_gas[2] *= h_inv_dim * rho_inv; + +#ifdef SWIFT_SINK_DENSITY_CHECKS + si->rho_check *= h_inv_dim; + si->n_check *= h_inv_dim; +#endif +} + +/** + * @brief Sets all particle fields to sensible values when the #sink has 0 + * ngbs. + * + * @param sp The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void sinks_sink_has_no_neighbours( + struct sink* restrict sp, const struct cosmology* cosmo) { + + warning( + "Sink particle with ID %lld treated as having no neighbours (h: %g, " + "wcount: %g).", + sp->id, sp->h, sp->density.wcount); + + /* Some smoothing length multiples. */ + const float h = sp->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + + /* Re-set problematic values */ + sp->density.wcount = kernel_root * h_inv_dim; + sp->density.wcount_dh = 0.f; +} + +/** + * @brief Compute the accretion rate of the sink and any quantities + * required swallowing based on an accretion rate + * + * Adapted from black_holes_prepare_feedback + * + * @param si The sink particle. + * @param props The properties of the sink scheme. + * @param constants The physical constants (in internal units). + * @param cosmo The cosmological model. + * @param cooling Properties of the cooling model. + * @param floor_props Properties of the entropy floor. + * @param time Time since the start of the simulation (non-cosmo mode). + * @param with_cosmology Are we running with cosmology? + * @param dt The time-step size (in physical internal units). + * @param ti_begin Integer time value at the beginning of timestep + */ +__attribute__((always_inline)) INLINE static void sink_prepare_swallow( + struct sink* restrict si, const struct sink_props* props, + const struct phys_const* constants, const struct cosmology* cosmo, + const struct cooling_function_data* cooling, + const struct entropy_floor_properties* floor_props, const double time, + const int with_cosmology, const double dt, const integertime_t ti_begin) { + + if (dt == 0. || si->rho_gas == 0.) return; + + /* Gather some physical constants (all in internal units) */ + const double G = constants->const_newton_G; + + /* (Subgrid) mass of the sink (internal units) */ + const double sink_mass = si->subgrid_mass; + + /* We can now compute the accretion rate (internal units) */ + /* Standard approach: compute accretion rate for all gas simultaneously. + * + * Convert the quantities we gathered to physical frame (all internal + * units). Note: velocities are already in black hole frame. */ + const double gas_rho_phys = si->rho_gas * cosmo->a3_inv; + const double gas_c_phys = si->sound_speed_gas * cosmo->a_factor_sound_speed; + const double gas_c_phys2 = gas_c_phys * gas_c_phys; + const double gas_v_phys[3] = {si->velocity_gas[0] * cosmo->a_inv, + si->velocity_gas[1] * cosmo->a_inv, + si->velocity_gas[2] * cosmo->a_inv}; + const double gas_v_norm2 = gas_v_phys[0] * gas_v_phys[0] + + gas_v_phys[1] * gas_v_phys[1] + + gas_v_phys[2] * gas_v_phys[2]; + + const double denominator2 = gas_v_norm2 + gas_c_phys2; + +#ifdef SWIFT_DEBUG_CHECKS + /* Make sure that the denominator is strictly positive */ + if (denominator2 <= 0) + error( + "Invalid denominator for sink particle %lld in Bondi rate " + "calculation.", + si->id); +#endif + + const double denominator_inv = 1. / sqrt(denominator2); + + double accr_rate = 4. * M_PI * G * G * sink_mass * sink_mass * gas_rho_phys * + denominator_inv * denominator_inv * denominator_inv; + + /* Integrate forward in time */ + si->subgrid_mass += accr_rate * dt; +} + +/** + * @brief Calculate if the gas has the potential of becoming + * a sink. + * + * Return 0 if no sink formation should occur. + * Note: called in runner_do_sink_formation + * + * @param p the gas particles. + * @param xp the additional properties of the gas particles. + * @param sink_props the sink properties to use. + * @param phys_const the physical constants in internal units. + * @param cosmo the cosmological parameters and properties. + * @param hydro_props The properties of the hydro scheme. + * @param us The internal system of units. + * @param cooling The cooling data struct. + * @param entropy_floor The entropy_floor properties. + * + */ +INLINE static int sink_is_forming( + const struct part* restrict p, const struct xpart* restrict xp, + const struct sink_props* sink_props, const struct phys_const* phys_const, + const struct cosmology* cosmo, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cooling_function_data* restrict cooling, + const struct entropy_floor_properties* restrict entropy_floor) { + + /* Sink formation is not implemented in this model. */ + return 0; +} + +/** + * @brief Decides whether a particle should be converted into a + * sink or not. + * + * No SF should occur, so return 0. + * Note: called in runner_do_sink_formation + * + * @param p The #part. + * @param xp The #xpart. + * @param sink_props The properties of the sink model. + * @param e The #engine (for random numbers). + * @param dt_sink The time-step of this particle + * @return 1 if a conversion should be done, 0 otherwise. + */ +INLINE static int sink_should_convert_to_sink( + const struct part* p, const struct xpart* xp, + const struct sink_props* sink_props, const struct engine* e, + const double dt_sink) { + + return 0; +} + +/** + * @brief Copies the properties of the gas particle over to the + * sink particle. + * + * @param p The gas particles. + * @param xp The additional properties of the gas particles. + * @param sink the new created #sink particle. + * @param e The #engine. + * @param sink_props The sink properties to use. + * @param cosmo the cosmological parameters and properties. + * @param with_cosmology if we run with cosmology. + * @param phys_const The physical constants in internal units. + * @param hydro_props The hydro properties to use. + * @param us The internal unit system. + * @param cooling The cooling function to use. + */ +INLINE static void sink_copy_properties( + const struct part* p, const struct xpart* xp, struct sink* sink, + const struct engine* e, const struct sink_props* sink_props, + const struct cosmology* cosmo, const int with_cosmology, + const struct phys_const* phys_const, + const struct hydro_props* restrict hydro_props, + const struct unit_system* restrict us, + const struct cooling_function_data* restrict cooling) {} + +/** + * @brief Update the properties of a sink particles by swallowing + * a gas particle. + * + * @param sp The #sink to update. + * @param p The #part that is swallowed. + * @param xp The #xpart that is swallowed. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void sink_swallow_part( + struct sink* sp, const struct part* p, const struct xpart* xp, + const struct cosmology* cosmo) { + + /* Get the current dynamical masses */ + const float gas_mass = hydro_get_mass(p); + const float sink_mass = sp->mass; + + /* Increase the dynamical mass of the sink. */ + sp->mass += gas_mass; + sp->gpart->mass += gas_mass; + + /* Comoving and physical distance between the particles */ + const float dx[3] = {sp->x[0] - p->x[0], sp->x[1] - p->x[1], + sp->x[2] - p->x[2]}; + const float dx_physical[3] = {dx[0] * cosmo->a, dx[1] * cosmo->a, + dx[2] * cosmo->a}; + + /* Relative velocity between the sink and the part */ + const float dv[3] = {sp->v[0] - p->v[0], sp->v[1] - p->v[1], + sp->v[2] - p->v[2]}; + + const float a = cosmo->a; + const float H = cosmo->H; + const float a2H = a * a * H; + + /* Calculate the velocity with the Hubble flow */ + const float v_plus_H_flow[3] = {a2H * dx[0] + dv[0], a2H * dx[1] + dv[1], + a2H * dx[2] + dv[2]}; + + /* Compute the physical relative velocity between the particles */ + const float dv_physical[3] = {v_plus_H_flow[0] * cosmo->a_inv, + v_plus_H_flow[1] * cosmo->a_inv, + v_plus_H_flow[2] * cosmo->a_inv}; + + /* Collect the swallowed angular momentum */ + sp->swallowed_angular_momentum[0] += + gas_mass * + (dx_physical[1] * dv_physical[2] - dx_physical[2] * dv_physical[1]); + sp->swallowed_angular_momentum[1] += + gas_mass * + (dx_physical[2] * dv_physical[0] - dx_physical[0] * dv_physical[2]); + sp->swallowed_angular_momentum[2] += + gas_mass * + (dx_physical[0] * dv_physical[1] - dx_physical[1] * dv_physical[0]); + + /* Update the sink momentum */ + const float sink_mom[3] = {sink_mass * sp->v[0] + gas_mass * p->v[0], + sink_mass * sp->v[1] + gas_mass * p->v[1], + sink_mass * sp->v[2] + gas_mass * p->v[2]}; + + sp->v[0] = sink_mom[0] / sp->mass; + sp->v[1] = sink_mom[1] / sp->mass; + sp->v[2] = sink_mom[2] / sp->mass; + sp->gpart->v_full[0] = sp->v[0]; + sp->gpart->v_full[1] = sp->v[1]; + sp->gpart->v_full[2] = sp->v[2]; + + /* This sink swallowed a gas particle */ + sp->number_of_gas_swallows++; + sp->number_of_direct_gas_swallows++; + +#ifdef SWIFT_DEBUG_CHECKS + const float dr = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]); + message( + "sink %lld swallow gas particle %lld. " + "(Mass = %e, " + "Delta_v = [%f, %f, %f] U_V, " + "Delta_x = [%f, %f, %f] U_L, " + "Delta_v_rad = %f)", + sp->id, p->id, sp->mass, -dv[0], -dv[1], -dv[2], -dx[0], -dx[1], -dx[2], + (dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]) / dr); +#endif +} + +/** + * @brief Update the properties of a sink particles by swallowing + * a sink particle. + * + * @param spi The #sink to update. + * @param spj The #sink that is swallowed. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void sink_swallow_sink( + struct sink* spi, const struct sink* spj, const struct cosmology* cosmo) { + + /* Get the current dynamical masses */ + const float spi_dyn_mass = spi->mass; + const float spj_dyn_mass = spj->mass; + + /* Increase the masses of the sink. */ + spi->mass += spj->mass; + spi->gpart->mass += spj->mass; + spi->subgrid_mass += spj->subgrid_mass; + + /* Collect the swallowed angular momentum */ + spi->swallowed_angular_momentum[0] += spj->swallowed_angular_momentum[0]; + spi->swallowed_angular_momentum[1] += spj->swallowed_angular_momentum[1]; + spi->swallowed_angular_momentum[2] += spj->swallowed_angular_momentum[2]; + + /* Update the sink momentum */ + const float sink_mom[3] = { + spi_dyn_mass * spi->v[0] + spj_dyn_mass * spj->v[0], + spi_dyn_mass * spi->v[1] + spj_dyn_mass * spj->v[1], + spi_dyn_mass * spi->v[2] + spj_dyn_mass * spj->v[2]}; + + spi->v[0] = sink_mom[0] / spi->mass; + spi->v[1] = sink_mom[1] / spi->mass; + spi->v[2] = sink_mom[2] / spi->mass; + spi->gpart->v_full[0] = spi->v[0]; + spi->gpart->v_full[1] = spi->v[1]; + spi->gpart->v_full[2] = spi->v[2]; + + /* This sink swallowed a sink particle */ + spi->number_of_sink_swallows++; + spi->number_of_direct_sink_swallows++; + + /* Add all other swallowed particles swallowed by the swallowed sink */ + spi->number_of_sink_swallows += spj->number_of_sink_swallows; + spi->number_of_gas_swallows += spj->number_of_gas_swallows; + + message("sink %lld swallow sink particle %lld. New mass: %e.", spi->id, + spj->id, spi->mass); +} + +/** + * @brief Should the sink spawn a star particle? + * + * @param sink the sink particle. + * @param e The #engine + * @param sink_props The sink properties to use. + * @param cosmo The cosmological parameters and properties. + * @param with_cosmology If we run with cosmology. + * @param phys_const The physical constants in internal units. + * @param us The internal unit system. + */ +INLINE static int sink_spawn_star(struct sink* sink, const struct engine* e, + const struct sink_props* sink_props, + const struct cosmology* cosmo, + const int with_cosmology, + const struct phys_const* phys_const, + const struct unit_system* restrict us) { + + /* Star formation from sinks is disabled in this model. */ + return 0; +} + +/** + * @brief Copy the properties of the sink particle towards the new star. Also, + * give the stars some properties such as position and velocity. + * + * This function also needs to update the sink particle. + * + * @param sink The #sink particle. + * @param sp The star particle. + * @param e The #engine + * @param sink_props The sink properties to use. + * @param cosmo The cosmological parameters and properties. + * @param with_cosmology If we run with cosmology. + * @param phys_const The physical constants in internal units. + * @param us The internal unit system. + */ +INLINE static void sink_copy_properties_to_star( + struct sink* sink, struct spart* sp, const struct engine* e, + const struct sink_props* sink_props, const struct cosmology* cosmo, + const int with_cosmology, const struct phys_const* phys_const, + const struct unit_system* restrict us) {} + +/** + * @brief Update the #sink particle properties before spawning a star. + * + * In GEAR, we check if the sink had an IMF change from pop III to pop II + * during the last gas/sink accretion loops. If so, we draw a new target mass + * with the correct IMF so that stars have the right metallicities. + * + * @param sink The #sink particle. + * @param e The #engine + * @param sink_props The sink properties to use. + * @param phys_const The physical constants in internal units. + */ +INLINE static void sink_update_sink_properties_before_star_formation( + struct sink* sink, const struct engine* e, + const struct sink_props* sink_props, const struct phys_const* phys_const) {} + +/** + * @brief Update the #sink particle properties right after spawning a star. + * + * In GEAR: Important properties that are updated are the sink mass and the + * sink->target_mass_Msun to draw the next star mass. + * + * @param sink The #sink particle that spawed stars. + * @param sp The #spart particle spawned. + * @param e The #engine + * @param sink_props the sink properties to use. + * @param phys_const the physical constants in internal units. + * @param star_counter The star loop counter. + */ +INLINE static void sink_update_sink_properties_during_star_formation( + struct sink* sink, const struct spart* sp, const struct engine* e, + const struct sink_props* sink_props, const struct phys_const* phys_const, + int star_counter) {} + +/** + * @brief Update the #sink particle properties after star formation. + * + * In GEAR, this is unused. + * + * @param sink The #sink particle. + * @param e The #engine + * @param sink_props The sink properties to use. + * @param phys_const The physical constants in internal units. + */ +INLINE static void sink_update_sink_properties_after_star_formation( + struct sink* sink, const struct engine* e, + const struct sink_props* sink_props, const struct phys_const* phys_const) {} + +/** + * @brief Store the gravitational potential of a particle by copying it from + * its #gpart friend. + * + * @param p_data The sink data of a gas particle. + * @param gp The part's #gpart. + */ +__attribute__((always_inline)) INLINE static void sink_store_potential_in_part( + struct sink_part_data* p_data, const struct gpart* gp) {} + +/** + * @brief Compute all quantities required for the formation of a sink such as + * kinetic energy, potential energy, etc. This function works on the + * neighbouring gas particles. + * + * @param e The #engine. + * @param p The #part for which we compute the quantities. + * @param xp The #xpart data of the particle #p. + * @param pi A neighbouring #part of #p. + * @param xpi The #xpart data of the particle #pi. + * @param cosmo The cosmological parameters and properties. + * @param sink_props The sink properties to use. + */ +INLINE static void sink_prepare_part_sink_formation_gas_criteria( + struct engine* e, struct part* restrict p, struct xpart* restrict xp, + struct part* restrict pi, struct xpart* restrict xpi, + const struct cosmology* cosmo, const struct sink_props* sink_props) {} + +/** + * @brief Compute all quantities required for the formation of a sink. This + * function works on the neighbouring sink particles. + * + * @param e The #engine. + * @param p The #part for which we compute the quantities. + * @param xp The #xpart data of the particle #p. + * @param si A neighbouring #sink of #p. + * @param cosmo The cosmological parameters and properties. + * @param sink_props The sink properties to use. + */ +INLINE static void sink_prepare_part_sink_formation_sink_criteria( + struct engine* e, struct part* restrict p, struct xpart* restrict xp, + struct sink* restrict si, const int with_cosmology, + const struct cosmology* cosmo, const struct sink_props* sink_props, + const double time) {} + +#endif /* SWIFT_BASIC_SINK_H */ diff --git a/src/sink/Basic/sink_debug.h b/src/sink/Basic/sink_debug.h new file mode 100644 index 0000000000..441af55c3e --- /dev/null +++ b/src/sink/Basic/sink_debug.h @@ -0,0 +1,29 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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_SINK_BASIC_DEBUG_H +#define SWIFT_SINK_BASIC_DEBUG_H + +__attribute__((always_inline)) INLINE static void sink_debug_particle( + const struct part* p, const struct xpart* xp) { + + warning("[PID%lld] sink_part_data:", p->id); + warning("[PID%lld] swallow_id = %lld", p->id, p->sink_data.swallow_id); +} + +#endif /* SWIFT_SINK_BASIC_DEBUG_H */ diff --git a/src/sink/Basic/sink_iact.h b/src/sink/Basic/sink_iact.h new file mode 100644 index 0000000000..f6cdf999dd --- /dev/null +++ b/src/sink/Basic/sink_iact.h @@ -0,0 +1,373 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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_BASIC_SINKS_IACT_H +#define SWIFT_BASIC_SINKS_IACT_H + +/* Local includes */ +#include "gravity.h" +#include "gravity_iact.h" +#include "random.h" +#include "sink_properties.h" + +/** + * @brief Gas particle interactions relevant for sinks, to run in hydro density + * interaction + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi First particle. + * @param pj Second particle. + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_sink( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, struct part *restrict pj, const float a, + const float H) {} + +/** + * @brief Gas particle interactions relevant for sinks, to run in hydro density + * interaction (non symmetric version) + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param pi First particle. + * @param pj Second particle (not updated). + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( + const float r2, const float dx[3], const float hi, const float hj, + struct part *restrict pi, const struct part *restrict pj, const float a, + const float H) {} + +/** + * @brief Density interaction between sinks and gas (non-symmetric). + * + * The gas particle cannot be touched. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing length or cut off radius of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param si First particle (sink). + * @param pj Second particle (gas, not updated). + * @param with_cosmology Are we doing a cosmological run? + * @param cosmo The cosmological model. + * @param grav_props The properties of the gravity scheme (softening, G, ...). + * @param sink_props the sink properties to use. + * @param ti_current Current integer time value (for random numbers). + * @param time current physical time in the simulation + */ +__attribute__((always_inline)) INLINE static void +runner_iact_nonsym_sinks_gas_density( + const float r2, const float dx[3], const float hi, const float hj, + struct sink *si, const struct part *pj, const int with_cosmology, + const struct cosmology *cosmo, const struct gravity_props *grav_props, + const struct sink_props *sink_props, const integertime_t ti_current, + const double time) { + + float wi, wi_dx; + + /* Get r. */ + const float r = sqrtf(r2); + + /* Compute the kernel function */ + const float hi_inv = 1.0f / hi; + const float ui = r * hi_inv; + kernel_deval(ui, &wi, &wi_dx); + + /* Compute contribution to the number of neighbours */ + si->density.wcount += wi; + si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx); + + /* Contribution to the number of neighbours */ + si->num_ngbs++; + + /* Neighbour gas mass */ + const float mj = hydro_get_mass(pj); + + /* Contribution to the BH gas density */ + si->rho_gas += mj * wi; + + /* Contribution to the total neighbour mass */ + si->ngb_mass += mj; + + /* Neighbour's sound speed */ + const float cj = hydro_get_comoving_soundspeed(pj); + + /* Contribution to the smoothed sound speed */ + si->sound_speed_gas += mj * cj * wi; + + /* Neighbour's (drifted) velocity in the frame of the black hole + * (we do include a Hubble term) */ + const float dv[3] = {pj->v[0] - si->v[0], pj->v[1] - si->v[1], + pj->v[2] - si->v[2]}; + + /* Contribution to the smoothed velocity (gas w.r.t. black hole) */ + si->velocity_gas[0] += mj * dv[0] * wi; + si->velocity_gas[1] += mj * dv[1] * wi; + si->velocity_gas[2] += mj * dv[2] * wi; + +#ifdef SWIFT_SINK_DENSITY_CHECKS + si->rho_check += mj * wi; + si->n_check += wi; + si->N_check_density++; +#endif +} + +/** + * @brief Compute sink-sink swallow interaction (non-symmetric). + * + * Note: Energies are computed with physical quantities, not the comoving ones. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing length or cut off radius of particle i. + * @param hj Comoving smoothing length or cut off radius of particle j. + * @param si First sink particle. + * @param sj Second sink particle. + * @param with_cosmology if we run with cosmology. + * @param cosmo The cosmological parameters and properties. + * @param grav_props The gravity scheme parameters and properties. + * @param sink_props the sink properties to use. + * @param ti_current Current integer time value (for random numbers). + * @param time current physical time in the simulation + */ +__attribute__((always_inline)) INLINE static void +runner_iact_nonsym_sinks_sink_swallow( + const float r2, const float dx[3], const float hi, const float hj, + struct sink *restrict si, struct sink *restrict sj, + const int with_cosmology, const struct cosmology *cosmo, + const struct gravity_props *grav_props, + const struct sink_props *sink_properties, const integertime_t ti_current, + const double time) { + + /* Simpler version of GEAR as a placeholder. Sinks bound to each other are + * merged. */ + + /* Relative velocity between the sinks */ + const float dv[3] = {sj->v[0] - si->v[0], sj->v[1] - si->v[1], + sj->v[2] - si->v[2]}; + + const float a = cosmo->a; + const float H = cosmo->H; + const float a2H = a * a * H; + + /* Calculate the velocity with the Hubble flow */ + const float v_plus_H_flow[3] = {a2H * dx[0] + dv[0], a2H * dx[1] + dv[1], + a2H * dx[2] + dv[2]}; + + /* Binding energy check */ + /* Compute the physical relative velocity between the particles */ + const float dv_physical[3] = {v_plus_H_flow[0] * cosmo->a_inv, + v_plus_H_flow[1] * cosmo->a_inv, + v_plus_H_flow[2] * cosmo->a_inv}; + + const float dv_physical_squared = dv_physical[0] * dv_physical[0] + + dv_physical[1] * dv_physical[1] + + dv_physical[2] * dv_physical[2]; + + /* Kinetic energy per unit mass of the gas */ + const float E_kin_rel = 0.5f * dv_physical_squared; + + /* Compute the Newtonian or softened potential the sink exherts onto the + gas particle */ + /* TODO: needs updating for MPI safety. We don't have access to foreign gparts + * here. */ + const float eps = gravity_get_softening(si->gpart, grav_props); + const float eps2 = eps * eps; + const float eps_inv = 1.f / eps; + const float eps_inv3 = eps_inv * eps_inv * eps_inv; + const float si_mass = si->mass; + const float sj_mass = sj->mass; + + float dummy, pot_ij, pot_ji; + runner_iact_grav_pp_full(r2, eps2, eps_inv, eps_inv3, si_mass, &dummy, + &pot_ij); + runner_iact_grav_pp_full(r2, eps2, eps_inv, eps_inv3, sj_mass, &dummy, + &pot_ji); + + /* Compute the physical potential energies per unit mass : + E_pot_phys = G*pot_grav*a^(-1) + c(a). + The normalization is c(a) = 0. */ + const float E_pot_ij = grav_props->G_Newton * pot_ij * cosmo->a_inv; + const float E_pot_ji = grav_props->G_Newton * pot_ji * cosmo->a_inv; + + /* Mechanical energy per unit mass of the pair i-j and j-i */ + const float E_mec_si = E_kin_rel + E_pot_ij; + const float E_mec_sj = E_kin_rel + E_pot_ji; + + /* Now, check if one is bound to the other */ + if ((E_mec_si > 0) || (E_mec_sj > 0)) { + return; + } + + /* The sink with the smaller mass will be merged onto the one with the + * larger mass. + * To avoid rounding issues, we additionally check for IDs if the sink + * have the exact same mass. */ + if ((sj->mass < si->mass) || (sj->mass == si->mass && sj->id < si->id)) { + /* This particle is swallowed by the sink with the largest mass of all the + * candidates wanting to swallow it (we use IDs to break ties)*/ + if ((sj->merger_data.swallow_mass < si->mass) || + (sj->merger_data.swallow_mass == si->mass && + sj->merger_data.swallow_id < si->id)) { + sj->merger_data.swallow_id = si->id; + sj->merger_data.swallow_mass = si->mass; + } + } + +#ifdef DEBUG_INTERACTIONS_SINKS + /* Update ngb counters */ + if (si->num_ngb_formation < MAX_NUM_OF_NEIGHBOURS_SINKS) + si->ids_ngbs_formation[si->num_ngb_formation] = pj->id; + + /* Update ngb counters */ + ++si->num_ngb_formation; +#endif +} + +/** + * @brief Compute sink-gas swallow interaction (non-symmetric). + * + * Note: Energies are computed with physical quantities, not the comoving ones. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing length or cut off radius of particle i. + * @param hj Comoving smoothing-length of particle j. + * @param si First sink particle. + * @param pj Second particle. + * @param ti_current Current integer time value (for random numbers). + * @param time current physical time in the simulation + */ +__attribute__((always_inline)) INLINE static void +runner_iact_nonsym_sinks_gas_swallow( + const float r2, const float dx[3], const float hi, const float hj, + struct sink *restrict si, struct part *restrict pj, + const int with_cosmology, const struct cosmology *cosmo, + const struct gravity_props *grav_props, + const struct sink_props *sink_properties, const integertime_t ti_current, + const double time) { + + float wi; + + /* Get r. */ + const float r = sqrtf(r2); + + /* Compute the kernel function */ + const float hi_inv = 1.0f / hi; + const float hi_inv_dim = pow_dimension(hi_inv); + const float ui = r * hi_inv; + kernel_eval(ui, &wi); + + /* Check if the sink needs to be fed. If not, we're done here */ + const float sink_mass_deficit = si->subgrid_mass - si->mass_at_start_of_step; + if (sink_mass_deficit <= 0) return; + + if (sink_properties->use_nibbling) { + + /* If we do nibbling, things are quite straightforward. We transfer + * the mass and all associated quantities right here. */ + + const float si_mass_orig = si->mass; + const float pj_mass_orig = hydro_get_mass(pj); + + /* Don't nibble from particles that are too small already */ + if (pj_mass_orig < sink_properties->min_gas_mass_for_nibbling) return; + + /* Next line is equivalent to w_ij * m_j / Sum_j (w_ij * m_j) */ + const float particle_weight = hi_inv_dim * wi * pj_mass_orig / si->rho_gas; + float nibble_mass = sink_mass_deficit * particle_weight; + + /* Need to check whether nibbling would push gas mass below minimum + * allowed mass */ + float new_gas_mass = pj_mass_orig - nibble_mass; + if (new_gas_mass < sink_properties->min_gas_mass_for_nibbling) { + new_gas_mass = sink_properties->min_gas_mass_for_nibbling; + nibble_mass = pj_mass_orig - sink_properties->min_gas_mass_for_nibbling; + } + + /* Transfer (dynamical) mass from the gas particle to the sink */ + si->mass += nibble_mass; + hydro_set_mass(pj, new_gas_mass); + + /* Add the angular momentum of the accreted gas to the sink total. + * Note no change to gas here. The cosmological conversion factors for + * velocity (a^-1) and distance (a) cancel out, so the angular momentum + * is already in physical units. */ + const float dv[3] = {si->v[0] - pj->v[0], si->v[1] - pj->v[1], + si->v[2] - pj->v[2]}; + si->swallowed_angular_momentum[0] += + nibble_mass * (dx[1] * dv[2] - dx[2] * dv[1]); + si->swallowed_angular_momentum[1] += + nibble_mass * (dx[2] * dv[0] - dx[0] * dv[2]); + si->swallowed_angular_momentum[2] += + nibble_mass * (dx[0] * dv[1] - dx[1] * dv[0]); + + /* Update the BH momentum and velocity. Again, no change to gas here. */ + const float si_mom[3] = {si_mass_orig * si->v[0] + nibble_mass * pj->v[0], + si_mass_orig * si->v[1] + nibble_mass * pj->v[1], + si_mass_orig * si->v[2] + nibble_mass * pj->v[2]}; + + si->v[0] = si_mom[0] / si->mass; + si->v[1] = si_mom[1] / si->mass; + si->v[2] = si_mom[2] / si->mass; + + } else { /* ends nibbling section, below comes swallowing */ + + /* Probability to swallow this particle + * Recall that in SWIFT the SPH kernel is recovered by computing + * kernel_eval() and muliplying by (1/h^d) */ + + const float prob = + (si->subgrid_mass - si->mass) * hi_inv_dim * wi / si->rho_gas; + + /* Draw a random number (Note mixing both IDs) */ + const float rand = random_unit_interval_two_IDs(si->id, pj->id, ti_current, + random_number_sink_swallow); + + /* Are we lucky? */ + if (rand < prob) { + + /* This particle is swallowed by the BH with the largest ID of all the + * candidates wanting to swallow it */ + if (pj->sink_data.swallow_id < si->id) { + + message("Sink %lld wants to swallow gas particle %lld", si->id, pj->id); + + pj->sink_data.swallow_id = si->id; + + } else { + + message( + "Sink %lld wants to swallow gas particle %lld BUT CANNOT (old " + "swallow id=%lld)", + si->id, pj->id, pj->sink_data.swallow_id); + } + } + } /* ends section for swallowing */ +} + +#endif diff --git a/src/sink/Basic/sink_io.h b/src/sink/Basic/sink_io.h new file mode 100644 index 0000000000..0e0cabc695 --- /dev/null +++ b/src/sink/Basic/sink_io.h @@ -0,0 +1,219 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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_BASIC_SINK_IO_H +#define SWIFT_BASIC_SINK_IO_H + +#include "io_properties.h" +#include "sink_part.h" + +/** + * @brief Specifies which sink-particle fields to read from a dataset + * + * @param sinks The sink-particle array. + * @param list The list of i/o properties to read. + * @param num_fields The number of i/o fields to read. + */ +INLINE static void sink_read_particles(struct sink* sinks, + struct io_props* list, int* num_fields) { + + /* Say how much we want to read */ + *num_fields = 5; + + /* List what we want to read */ + list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, + UNIT_CONV_LENGTH, sinks, x); + list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY, + UNIT_CONV_SPEED, sinks, v); + list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS, + sinks, mass); + list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY, + UNIT_CONV_NO_UNITS, sinks, id); + list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL, + UNIT_CONV_LENGTH, sinks, h); +} + +INLINE static void convert_sink_pos(const struct engine* e, + const struct sink* sp, double* ret) { + + const struct space* s = e->s; + if (s->periodic) { + ret[0] = box_wrap(sp->x[0], 0.0, s->dim[0]); + ret[1] = box_wrap(sp->x[1], 0.0, s->dim[1]); + ret[2] = box_wrap(sp->x[2], 0.0, s->dim[2]); + } else { + ret[0] = sp->x[0]; + ret[1] = sp->x[1]; + ret[2] = sp->x[2]; + } +} + +INLINE static void convert_sink_vel(const struct engine* e, + const struct sink* sp, float* ret) { + + const int with_cosmology = (e->policy & engine_policy_cosmology); + const struct cosmology* cosmo = e->cosmology; + const integertime_t ti_current = e->ti_current; + const double time_base = e->time_base; + + const integertime_t ti_beg = get_integer_time_begin(ti_current, sp->time_bin); + const integertime_t ti_end = get_integer_time_end(ti_current, sp->time_bin); + + /* Get time-step since the last kick */ + float dt_kick_grav; + if (with_cosmology) { + dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current); + dt_kick_grav -= + cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2); + } else { + dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base; + } + + /* Extrapolate the velocites to the current time */ + const struct gpart* gp = sp->gpart; + ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav; + ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav; + ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav; + + /* Conversion from internal units to peculiar velocities */ + ret[0] *= cosmo->a_inv; + ret[1] *= cosmo->a_inv; + ret[2] *= cosmo->a_inv; +} + +INLINE static void convert_sink_gas_vel(const struct engine* e, + const struct sink* sink, float* ret) { + const struct cosmology* cosmo = e->cosmology; + ret[0] = sink->velocity_gas[0] * cosmo->a_inv; + ret[1] = sink->velocity_gas[1] * cosmo->a_inv; + ret[2] = sink->velocity_gas[2] * cosmo->a_inv; +} + +INLINE static void convert_sink_gas_sound_speed(const struct engine* e, + const struct sink* sink, + double* ret) { + const struct cosmology* cosmo = e->cosmology; + ret[0] = sink->sound_speed_gas * cosmo->a_factor_sound_speed; +} + +INLINE static void convert_sink_swallowed_angular_momentum( + const struct engine* e, const struct sink* sink, float* ret) { + ret[0] = sink->swallowed_angular_momentum[0]; + ret[1] = sink->swallowed_angular_momentum[1]; + ret[2] = sink->swallowed_angular_momentum[2]; +} + +/** + * @brief Specifies which sink-particle fields to write to a dataset + * + * @param sinks The sink-particle array. + * @param list The list of i/o properties to write. + * @param num_fields The number of i/o fields to write. + * @param with_cosmology Are we running a cosmological simulation? + */ +INLINE static void sink_write_particles(const struct sink* sinks, + struct io_props* list, int* num_fields, + int with_cosmology) { + + /* Say how much we want to write */ + *num_fields = 12; + + /* List what we want to write */ + list[0] = io_make_output_field_convert_sink( + "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, sinks, convert_sink_pos, + "Co-moving position of the particles"); + + list[1] = io_make_output_field_convert_sink( + "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, sinks, convert_sink_vel, + "Peculiar velocities of the particles. This is a * dx/dt where x is the " + "co-moving position of the particles."); + + list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, sinks, + mass, "Dynamical masses of the sinks."); + + list[3] = io_make_physical_output_field( + "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, id, + /*can convert to comoving=*/0, "Unique ID of the particles"); + + list[4] = io_make_output_field( + "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sinks, h, + "Co-moving smoothing lengths (FWHM of the kernel) of the particles"); + + list[5] = io_make_physical_output_field( + "NumberOfSinkSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, + number_of_sink_swallows, /*can convert to comoving=*/0, + "Total number of sink merger events"); + + list[6] = io_make_physical_output_field( + "NumberOfGasSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, + number_of_gas_swallows, /*can convert to comoving=*/0, + "Total number of gas merger events"); + + /* Note: Since the swallowed momentum is computed with the physical velocity, + i.e. including the Hubble flow term, it is not convertible to comoving + frame. */ + list[7] = io_make_physical_output_field_convert_sink( + "SwallowedAngularMomentum", FLOAT, 3, UNIT_CONV_ANGULAR_MOMENTUM, 0.f, + sinks, + /*can convert to comoving=*/0, convert_sink_swallowed_angular_momentum, + "Physical swallowed angular momentum of the particles"); + + list[8] = io_make_output_field( + "SubgridMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f, sinks, subgrid_mass, + "Subgrid mass of the sink. Summed on sink-sink mergers"); + + list[9] = io_make_physical_output_field( + "GasDensities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, sinks, rho_gas, + /*can convert to comoving=*/1, "Gas density at the location of the sink"); + + list[10] = io_make_output_field_convert_sink( + "GasVelocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, sinks, + convert_sink_gas_vel, + "Gas velocity at the location of the sink. Velocities are peculiar, i.e. " + "a * dx/dt where x is the " + "co-moving position of the gas."); + + list[11] = io_make_output_field_convert_sink( + "GasSoundSpeeds", DOUBLE, 1, UNIT_CONV_SPEED, 0.f, sinks, + convert_sink_gas_sound_speed, + "Gas sound speed at the location of the sink (physical units)"); + +#ifdef DEBUG_INTERACTIONS_SINKS + + list += *num_fields; + *num_fields += 4; + + list[0] = + io_make_output_field("Num_ngb_formation", INT, 1, UNIT_CONV_NO_UNITS, 0.f, + sinks, num_ngb_formation, "Number of neighbors"); + list[1] = + io_make_output_field("Ids_ngb_formation", LONGLONG, + MAX_NUM_OF_NEIGHBOURS_SINKS, UNIT_CONV_NO_UNITS, 0.f, + sinks, ids_ngbs_formation, "IDs of the neighbors"); + + list[2] = + io_make_output_field("Num_ngb_merger", INT, 1, UNIT_CONV_NO_UNITS, 0.f, + sinks, num_ngb_merger, "Number of merger"); + list[3] = io_make_output_field( + "Ids_ngb_merger", LONGLONG, MAX_NUM_OF_NEIGHBOURS_SINKS, + UNIT_CONV_NO_UNITS, 0.f, sinks, ids_ngbs_merger, "IDs of the neighbors"); + +#endif +} + +#endif /* SWIFT_BASIC_SINK_IO_H */ diff --git a/src/sink/Basic/sink_part.h b/src/sink/Basic/sink_part.h new file mode 100644 index 0000000000..9e1ccbbe85 --- /dev/null +++ b/src/sink/Basic/sink_part.h @@ -0,0 +1,168 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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_BASIC_SINK_PART_H +#define SWIFT_BASIC_SINK_PART_H + +#include "timeline.h" + +#define sink_need_unique_id 1 + +/** + * @brief Particle fields for the sink particles. + * + * All quantities related to gravity are stored in the associate #gpart. + */ +struct sink { + + /*! Particle ID. */ + long long id; + + /*! Pointer to corresponding gravity part. */ + struct gpart* gpart; + + /*! Particle position. */ + double x[3]; + + /* Offset between current position and position at last tree rebuild. */ + float x_diff[3]; + + /*! Particle velocity. */ + float v[3]; + + /* Particle smoothing length */ + float h; + + struct { + + /* Number of neighbours. */ + float wcount; + + /* Number of neighbours spatial derivative. */ + float wcount_dh; + + } density; + + /*! Sink particle mass */ + float mass; + + /*! Total mass of the gas neighbours. */ + float ngb_mass; + + /*! Integer number of neighbours */ + int num_ngbs; + + /*! Instantaneous accretion rate */ + float accretion_rate; + + /*! Subgrid mass of the sink */ + float subgrid_mass; + + /*! Sink mass at the start of each step, prior to any nibbling */ + float mass_at_start_of_step; + + /*! Density of the gas surrounding the black hole. */ + float rho_gas; + + /*! Smoothed sound speed of the gas surrounding the black hole. */ + float sound_speed_gas; + + /*! Smoothed velocity of the gas surrounding the black hole, + * in the frame of the black hole (internal units) */ + float velocity_gas[3]; + + /*! Particle time bin */ + timebin_t time_bin; + + /*! Total (physical) angular momentum accumulated by swallowing particles */ + float swallowed_angular_momentum[3]; + + /*! Total number of sink merger events (including sink swallowed + * by merged-in sinks) */ + int number_of_sink_swallows; + + /*! Total number of sink merger events (excluding sink swallowed + * by merged-in sinks) */ + int number_of_direct_sink_swallows; + + /*! Total number of gas particles swallowed (including particles swallowed + * by merged-in sinks) */ + int number_of_gas_swallows; + + /*! Total number of gas particles swallowed (excluding particles swallowed + * by merged-in sinks) */ + int number_of_direct_gas_swallows; + + /*! sink merger information (e.g. merging ID) */ + struct sink_sink_data merger_data; + +#ifdef SWIFT_DEBUG_CHECKS + + /* Time of the last drift */ + integertime_t ti_drift; + + /* Time of the last kick */ + integertime_t ti_kick; + +#endif + +#ifdef DEBUG_INTERACTIONS_SINKS + /*! Number of interactions in merger SELF and PAIR */ + int num_ngb_merger; + + /*! List of interacting particles in merger SELF and PAIR */ + long long ids_ngbs_merger[MAX_NUM_OF_NEIGHBOURS_SINKS]; + + /*! Number of interactions in compute formation SELF and PAIR */ + int num_ngb_formation; + + /*! List of interacting particles in compute formation SELF and PAIR */ + long long ids_ngbs_formation[MAX_NUM_OF_NEIGHBOURS_SINKS]; + + /*! Number of interactions in compute formation SELF and PAIR */ + int num_ngb_accretion; + + /*! List of interacting particles in compute formation SELF and PAIR */ + long long ids_ngbs_accretion[MAX_NUM_OF_NEIGHBOURS_SINKS]; +#endif + +#ifdef SWIFT_SINK_DENSITY_CHECKS + + /* Integer number of neighbours in the density loop */ + int N_check_density; + + /* Exact integer number of neighbours in the density loop */ + int N_check_density_exact; + + /*! Has this particle interacted with any unhibited neighbour? */ + char inhibited_check_exact; + + float n_check; + + float n_check_exact; + + float rho_check; + + /*! Exact value of the density field obtained via brute-force loop */ + float rho_check_exact; + +#endif + +} SWIFT_STRUCT_ALIGN; + +#endif /* SWIFT_BASIC_SINK_PART_H */ diff --git a/src/sink/Basic/sink_properties.h b/src/sink/Basic/sink_properties.h new file mode 100644 index 0000000000..12a4d21c40 --- /dev/null +++ b/src/sink/Basic/sink_properties.h @@ -0,0 +1,148 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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_BASIC_SINK_PROPERTIES_H +#define SWIFT_BASIC_SINK_PROPERTIES_H + +/* Local header */ +#include "feedback_properties.h" +#include "parser.h" + +/** + * @brief Properties of sink in the Default model. + */ +struct sink_props { + + /* ----- Basic neighbour search properties ------ */ + + /*! Resolution parameter */ + float eta_neighbours; + + /*! Target weightd number of neighbours (for info only)*/ + float target_neighbours; + + /*! Smoothing length tolerance */ + float h_tolerance; + + /*! Tolerance on neighbour number (for info only)*/ + float delta_neighbours; + + /*! Maximal number of iterations to converge h */ + int max_smoothing_iterations; + + /*! Maximal change of h over one time-step */ + float log_max_h_change; + + /*! Are we using a fixed cutoff radius? (all smoothing length calculations are + * disabled if so) */ + char use_fixed_r_cut; + + /* Use nibbling rather than swallowing for gas? */ + float use_nibbling; + + /* Gas mass below which sinks will not nibble. */ + float min_gas_mass_for_nibbling; +}; + +/** + * @brief Initialise the sink properties from the parameter file. + * + * @param sp The #sink_props. + * @param phys_const The physical constants in the internal unit system. + * @param us The internal unit system. + * @param params The parsed parameters. + * @param cosmo The cosmological model. + * @param with_feedback Are we running with feedback? + */ +INLINE static void sink_props_init( + struct sink_props *sp, struct feedback_props *fp, + const struct phys_const *phys_const, const struct unit_system *us, + struct swift_params *params, const struct hydro_props *hydro_props, + const struct cosmology *cosmo, const int with_feedback) { + + /* We don't use a fixed cutoff radius in this model. This property must always + * be present, as we use it to skip smoothing length iterations in + * runner_ghost if a fixed cutoff is being used. */ + sp->use_fixed_r_cut = 0; + + /* Read in the basic neighbour search properties or default to the hydro + ones if the user did not provide any different values */ + + /* Kernel properties */ + sp->eta_neighbours = parser_get_opt_param_float( + params, "BasicSink:resolution_eta", hydro_props->eta_neighbours); + + /* Tolerance for the smoothing length Newton-Raphson scheme */ + sp->h_tolerance = parser_get_opt_param_float(params, "BasicSink:h_tolerance", + hydro_props->h_tolerance); + + /* Get derived properties */ + sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm; + const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance); + sp->delta_neighbours = + (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) * + kernel_norm; + + /* Number of iterations to converge h */ + sp->max_smoothing_iterations = + parser_get_opt_param_int(params, "BasicSink:max_ghost_iterations", + hydro_props->max_smoothing_iterations); + + /* Time integration properties */ + const float max_volume_change = + parser_get_opt_param_float(params, "BasicSink:max_volume_change", -1); + if (max_volume_change == -1) + sp->log_max_h_change = hydro_props->log_max_h_change; + else + sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv)); + + sp->use_nibbling = parser_get_param_int(params, "BasicSink:use_nibbling"); + if (sp->use_nibbling) { + sp->min_gas_mass_for_nibbling = parser_get_param_float( + params, "BasicSink:min_gas_mass_for_nibbling_Msun"); + sp->min_gas_mass_for_nibbling *= phys_const->const_solar_mass; + } +} + +/** + * @brief Write a sink_props struct to the given FILE as a stream of + * bytes. + * + * @param props the sink properties struct + * @param stream the file stream + */ +INLINE static void sink_struct_dump(const struct sink_props *props, + FILE *stream) { + restart_write_blocks((void *)props, sizeof(struct sink_props), 1, stream, + "sink props", "Sink props"); +} + +/** + * @brief Restore a sink_props struct from the given FILE as a stream of + * bytes. + * + * @param props the sink properties struct + * @param stream the file stream + */ +INLINE static void sink_struct_restore(const struct sink_props *props, + FILE *stream) { + restart_read_blocks((void *)props, sizeof(struct sink_props), 1, stream, NULL, + "Sink props"); +} + +#endif /* SWIFT_BASIC_SINK_PROPERTIES_H */ diff --git a/src/sink/Basic/sink_struct.h b/src/sink/Basic/sink_struct.h new file mode 100644 index 0000000000..a255368243 --- /dev/null +++ b/src/sink/Basic/sink_struct.h @@ -0,0 +1,115 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2024 Jonathan Davies (j.j.davies@ljmu.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_BASIC_STRUCT_DEFAULT_H +#define SWIFT_BASIC_STRUCT_DEFAULT_H + +/** + * @brief Sink-related fields carried by each *gas* particle. + */ +struct sink_part_data { + + /*! ID of the sink that will swallow this #part. */ + long long swallow_id; +}; + +/** + * @brief Sink-related fields carried by each *sink* particle. + */ +struct sink_sink_data { + + /*! ID of the sink that will swallow this #sink. */ + long long swallow_id; + + /*! Mass of the sink that will swallow this #sink. */ + float swallow_mass; +}; + +/** + * @brief Return the ID of the sink that should swallow this #part. + * + * @param s_data The #part's #sink_part_data structure. + */ +__attribute__((always_inline)) INLINE static long long sink_get_part_swallow_id( + struct sink_part_data* s_data) { + + return s_data->swallow_id; +} + +/** + * @brief Update a given #part's sink data field to mark the particle has + * not yet been swallowed. + * + * @param s_data The #part's #sink_part_data structure. + */ +__attribute__((always_inline)) INLINE static void +sink_mark_part_as_not_swallowed(struct sink_part_data* s_data) { + + s_data->swallow_id = -1; +} + +/** + * @brief Update a given #part's sink data field to mark the particle has + * having been been swallowed. + * + * @param p_data The #part's #sink_part_data structure. + */ +__attribute__((always_inline)) INLINE static void sink_mark_part_as_swallowed( + struct sink_part_data* s_data) { + + s_data->swallow_id = -2; +} + +/** + * @brief Update a given #sink's sink data field to mark the particle has + * not yet been swallowed. + * + * @param s_data The #sink's #sink_sink_data structure. + */ +__attribute__((always_inline)) INLINE static void +sink_mark_sink_as_not_swallowed(struct sink_sink_data* s_data) { + + s_data->swallow_id = -1; + s_data->swallow_mass = 0.f; +} + +/** + * @brief Update a given #sink's sink data field to mark the particle has + * having been been swallowed. + * + * @param s_data The #sink's #bsink_sink_data structure. + */ +__attribute__((always_inline)) INLINE static void sink_mark_sink_as_merged( + struct sink_sink_data* s_data) { + + s_data->swallow_id = -2; + s_data->swallow_mass = -1.f; +} + +/** + * @brief Return the ID of the sink that should swallow this #sink. + * + * @param s_data The #sink's #sink_sink_data structure. + */ +__attribute__((always_inline)) INLINE static long long sink_get_sink_swallow_id( + struct sink_sink_data* s_data) { + + return s_data->swallow_id; +} + +#endif /* SWIFT_BASIC_STRUCT_DEFAULT_H */ diff --git a/src/sink/Default/sink.h b/src/sink/Default/sink.h index e0d3f3b4a3..1001987ae5 100644 --- a/src/sink/Default/sink.h +++ b/src/sink/Default/sink.h @@ -81,7 +81,18 @@ __attribute__((always_inline)) INLINE static void sink_init_part( __attribute__((always_inline)) INLINE static void sink_init_sink( struct sink* sp) { - sp->num_ngbs = 0; + sp->density.wcount = 0.f; + sp->density.wcount_dh = 0.f; + +#ifdef SWIFT_SINK_DENSITY_CHECKS + sp->N_check_density = 0; + sp->N_check_density_exact = 0; + sp->rho_check = 0.f; + sp->rho_check_exact = 0.f; + sp->n_check = 0.f; + sp->n_check_exact = 0.f; + sp->inhibited_check_exact = 0; +#endif } /** @@ -118,7 +129,23 @@ __attribute__((always_inline)) INLINE static void sink_kick_extra( * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void sink_end_density( - struct sink* si, const struct cosmology* cosmo) {} + struct sink* si, const struct cosmology* cosmo) { + + /* Some smoothing length multiples. */ + const float h = si->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */ + + /* Finish the calculation by inserting the missing h-factors */ + si->density.wcount *= h_inv_dim; + si->density.wcount_dh *= h_inv_dim_plus_one; + +#ifdef SWIFT_SINK_DENSITY_CHECKS + si->rho_check *= h_inv_dim; + si->n_check *= h_inv_dim; +#endif +} /** * @brief Sets all particle fields to sensible values when the #sink has 0 @@ -128,7 +155,22 @@ __attribute__((always_inline)) INLINE static void sink_end_density( * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void sinks_sink_has_no_neighbours( - struct sink* restrict sp, const struct cosmology* cosmo) {} + struct sink* restrict sp, const struct cosmology* cosmo) { + + warning( + "Sink particle with ID %lld treated as having no neighbours (h: %g, " + "wcount: %g).", + sp->id, sp->h, sp->density.wcount); + + /* Some smoothing length multiples. */ + const float h = sp->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + + /* Re-set problematic values */ + sp->density.wcount = kernel_root * h_inv_dim; + sp->density.wcount_dh = 0.f; +} /** * @brief Compute the accretion rate of the sink and any quantities @@ -231,7 +273,11 @@ INLINE static void sink_copy_properties( const struct phys_const* phys_const, const struct hydro_props* restrict hydro_props, const struct unit_system* restrict us, - const struct cooling_function_data* restrict cooling) {} + const struct cooling_function_data* restrict cooling) { + + /* Set a smoothing length */ + sink->h = p->h; +} /** * @brief Update the properties of a sink particles by swallowing diff --git a/src/sink/Default/sink_iact.h b/src/sink/Default/sink_iact.h index bae7b6129f..653ae2fe05 100644 --- a/src/sink/Default/sink_iact.h +++ b/src/sink/Default/sink_iact.h @@ -31,12 +31,11 @@ * @param pj Second particle. * @param a Current scale factor. * @param H Current Hubble parameter. - * @param cut_off_radius Sink cut off radius. */ __attribute__((always_inline)) INLINE static void runner_iact_sink( const float r2, const float dx[3], const float hi, const float hj, struct part *restrict pi, struct part *restrict pj, const float a, - const float H, const float cut_off_radius) {} + const float H) {} /** * @brief do sink computation after the runner_iact_density (non symmetric @@ -50,21 +49,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_sink( * @param pj Second particle (not updated). * @param a Current scale factor. * @param H Current Hubble parameter. - * @param cut_off_radius Sink cut off radius. - * @param ti_current Current integer time value (for random numbers). - * @param time current physical time in the simulation */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( const float r2, const float dx[3], const float hi, const float hj, struct part *restrict pi, const struct part *restrict pj, const float a, - const float H, const float cut_off_radius) {} + const float H) {} /** * @brief Density interaction between two particles (non-symmetric). * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). - * @param ri Comoving cut off radius of particle i. + * @param hi Comoving smoothing length of particle i. * @param hj Comoving smoothing-length of particle j. * @param si First particle (sink). * @param pj Second particle (gas, not updated). @@ -77,19 +73,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sinks_gas_density( - const float r2, const float dx[3], const float ri, const float hj, + const float r2, const float dx[3], const float hi, const float hj, struct sink *si, const struct part *pj, const int with_cosmology, const struct cosmology *cosmo, const struct gravity_props *grav_props, const struct sink_props *sink_props, const integertime_t ti_current, const double time) { + float wi, wi_dx; + /* Get r. */ const float r = sqrtf(r2); - if (r < ri) { - /* Contribution to the number of neighbours in cutoff radius */ - si->num_ngbs++; - } + /* Compute the kernel function */ + const float hi_inv = 1.0f / hi; + const float ui = r * hi_inv; + kernel_deval(ui, &wi, &wi_dx); + + /* Compute contribution to the number of neighbours */ + si->density.wcount += wi; + si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx); + +#ifdef SWIFT_SINK_DENSITY_CHECKS + si->rho_check += hydro_get_mass(pj) * wi; + si->n_check += wi; + si->N_check_density++; +#endif } /** @@ -97,8 +105,8 @@ runner_iact_nonsym_sinks_gas_density( * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). - * @param ri Comoving cut off radius of particle i. - * @param rj Comoving cut off radius of particle j. + * @param hi Comoving smoothing length of particle i. + * @param hj Comoving smoothing length of particle j. * @param si First sink particle. * @param sj Second sink particle. * @param with_cosmology if we run with cosmology. @@ -110,7 +118,7 @@ runner_iact_nonsym_sinks_gas_density( */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sinks_sink_swallow( - const float r2, const float dx[3], const float ri, const float rj, + const float r2, const float dx[3], const float hi, const float hj, struct sink *restrict si, struct sink *restrict sj, const int with_cosmology, const struct cosmology *cosmo, const struct gravity_props *grav_props, @@ -122,7 +130,7 @@ runner_iact_nonsym_sinks_sink_swallow( * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). - * @param ri Comoving cut off radius of particle i. + * @param hi Comoving smoothing length of particle i. * @param hj Comoving smoothing-length of particle j. * @param si First sink particle. * @param pj Second particle. @@ -135,7 +143,7 @@ runner_iact_nonsym_sinks_sink_swallow( */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sinks_gas_swallow( - const float r2, const float dx[3], const float ri, const float hj, + const float r2, const float dx[3], const float hi, const float hj, struct sink *restrict si, struct part *restrict pj, const int with_cosmology, const struct cosmology *cosmo, const struct gravity_props *grav_props, diff --git a/src/sink/Default/sink_io.h b/src/sink/Default/sink_io.h index 734dbb3309..c4fd307376 100644 --- a/src/sink/Default/sink_io.h +++ b/src/sink/Default/sink_io.h @@ -33,7 +33,7 @@ INLINE static void sink_read_particles(struct sink* sinks, struct io_props* list, int* num_fields) { /* Say how much we want to read */ - *num_fields = 4; + *num_fields = 5; /* List what we want to read */ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, @@ -44,6 +44,8 @@ INLINE static void sink_read_particles(struct sink* sinks, sinks, mass); list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY, UNIT_CONV_NO_UNITS, sinks, id); + list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL, + UNIT_CONV_LENGTH, sinks, h); } INLINE static void convert_sink_pos(const struct engine* e, @@ -107,7 +109,7 @@ INLINE static void sink_write_particles(const struct sink* sinks, int with_cosmology) { /* Say how much we want to write */ - *num_fields = 4; + *num_fields = 5; /* List what we want to write */ list[0] = io_make_output_field_convert_sink( @@ -126,6 +128,10 @@ INLINE static void sink_write_particles(const struct sink* sinks, "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, id, /*can convert to comoving=*/0, "Unique ID of the particles"); + list[4] = io_make_output_field( + "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sinks, h, + "Co-moving smoothing lengths (FWHM of the kernel) of the particles"); + #ifdef DEBUG_INTERACTIONS_SINKS list += *num_fields; diff --git a/src/sink/Default/sink_part.h b/src/sink/Default/sink_part.h index 79ce4ed81a..893a01c514 100644 --- a/src/sink/Default/sink_part.h +++ b/src/sink/Default/sink_part.h @@ -45,8 +45,18 @@ struct sink { /*! Particle velocity. */ float v[3]; - /*! Cut off radius. */ - float r_cut; + /* Particle smoothing length */ + float h; + + struct { + + /* Number of neighbours. */ + float wcount; + + /* Number of neighbours spatial derivative. */ + float wcount_dh; + + } density; /*! Sink particle mass */ float mass; @@ -96,6 +106,29 @@ struct sink { /*! List of interacting particles in compute formation SELF and PAIR */ long long ids_ngbs_accretion[MAX_NUM_OF_NEIGHBOURS_SINKS]; #endif + +#ifdef SWIFT_SINK_DENSITY_CHECKS + + /* Integer number of neighbours in the density loop */ + int N_check_density; + + /* Exact integer number of neighbours in the density loop */ + int N_check_density_exact; + + /*! Has this particle interacted with any unhibited neighbour? */ + char inhibited_check_exact; + + float n_check; + + float n_check_exact; + + float rho_check; + + /*! Exact value of the density field obtained via brute-force loop */ + float rho_check_exact; + +#endif + } SWIFT_STRUCT_ALIGN; #endif /* SWIFT_DEFAULT_SINK_PART_H */ diff --git a/src/sink/Default/sink_properties.h b/src/sink/Default/sink_properties.h index 08b458c775..f9bd54e749 100644 --- a/src/sink/Default/sink_properties.h +++ b/src/sink/Default/sink_properties.h @@ -19,13 +19,38 @@ #ifndef SWIFT_DEFAULT_SINK_PROPERTIES_H #define SWIFT_DEFAULT_SINK_PROPERTIES_H +/* Local header */ +#include "feedback_properties.h" +#include "parser.h" + /** * @brief Properties of sink in the Default model. */ struct sink_props { - /*! Cut off radius */ - float cut_off_radius; + /* ----- Basic neighbour search properties ------ */ + + /*! Resolution parameter */ + float eta_neighbours; + + /*! Target weightd number of neighbours (for info only)*/ + float target_neighbours; + + /*! Smoothing length tolerance */ + float h_tolerance; + + /*! Tolerance on neighbour number (for info only)*/ + float delta_neighbours; + + /*! Maximal number of iterations to converge h */ + int max_smoothing_iterations; + + /*! Maximal change of h over one time-step */ + float log_max_h_change; + + /*! Are we using a fixed cutoff radius? (all smoothing length calculations are + * disabled if so) */ + char use_fixed_r_cut; }; /** @@ -38,16 +63,48 @@ struct sink_props { * @param cosmo The cosmological model. * @param with_feedback Are we running with feedback? */ -INLINE static void sink_props_init(struct sink_props *sp, - struct feedback_props *fp, - const struct phys_const *phys_const, - const struct unit_system *us, - struct swift_params *params, - const struct cosmology *cosmo, - const int with_feedback) { - - sp->cut_off_radius = - parser_get_param_float(params, "DefaultSink:cut_off_radius"); +INLINE static void sink_props_init( + struct sink_props *sp, struct feedback_props *fp, + const struct phys_const *phys_const, const struct unit_system *us, + struct swift_params *params, const struct hydro_props *hydro_props, + const struct cosmology *cosmo, const int with_feedback) { + + /* Read in the basic neighbour search properties or default to the hydro + ones if the user did not provide any different values */ + + /* Kernel properties */ + sp->eta_neighbours = parser_get_opt_param_float( + params, "Sinks:resolution_eta", hydro_props->eta_neighbours); + + /* Tolerance for the smoothing length Newton-Raphson scheme */ + sp->h_tolerance = parser_get_opt_param_float(params, "Sinks:h_tolerance", + hydro_props->h_tolerance); + + /* Get derived properties */ + sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm; + const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance); + sp->delta_neighbours = + (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) * + kernel_norm; + + /* Number of iterations to converge h */ + sp->max_smoothing_iterations = + parser_get_opt_param_int(params, "Sinks:max_ghost_iterations", + hydro_props->max_smoothing_iterations); + + /* Time integration properties */ + const float max_volume_change = + parser_get_opt_param_float(params, "Sinks:max_volume_change", -1); + if (max_volume_change == -1) + sp->log_max_h_change = hydro_props->log_max_h_change; + else + sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv)); + + /* This property determines whether we're using a fixed cutoff radius (and + * smoothing length) for the sink. It must always be present, as we use it + * to skip smoothing length iterations in runner_ghost if a fixed cutoff is + * being used. */ + sp->use_fixed_r_cut = 0; } /** diff --git a/src/sink/GEAR/sink.h b/src/sink/GEAR/sink.h index 627713379e..f2b5c00099 100644 --- a/src/sink/GEAR/sink.h +++ b/src/sink/GEAR/sink.h @@ -77,9 +77,7 @@ __attribute__((always_inline)) INLINE static float sink_compute_timestep( sink->to_collect.sound_speed_gas * cosmo->a_factor_sound_speed; const double gas_c_phys2 = gas_c_phys * gas_c_phys; const float denominator = sqrtf(gas_c_phys2 + gas_v_norm2); - const float h_min = - cosmo->a * - min(sink->r_cut, sink->to_collect.minimal_h_gas * kernel_gamma); + const float h_min = cosmo->a * min(sink->h, sink->to_collect.minimal_h_gas); float dt_cfl = 0.0; /* This case can happen if the sink is just born. */ @@ -167,7 +165,6 @@ __attribute__((always_inline)) INLINE static void sink_first_init_sink( struct sink* sp, const struct sink_props* sink_props, const struct engine* e) { - sp->r_cut = sink_props->cut_off_radius; sp->time_bin = 0; sp->number_of_gas_swallows = 0; @@ -255,6 +252,9 @@ __attribute__((always_inline)) INLINE static void sink_init_part( __attribute__((always_inline)) INLINE static void sink_init_sink( struct sink* sp) { + sp->density.wcount = 0.f; + sp->density.wcount_dh = 0.f; + /* Reset to the mass of the sink */ sp->mass_tot_before_star_spawning = sp->mass; @@ -284,6 +284,16 @@ __attribute__((always_inline)) INLINE static void sink_init_sink( sp->ids_ngbs_formation[i] = -1; sp->num_ngb_formation = 0; #endif + +#ifdef SWIFT_SINK_DENSITY_CHECKS + sp->N_check_density = 0; + sp->N_check_density_exact = 0; + sp->rho_check = 0.f; + sp->rho_check_exact = 0.f; + sp->n_check = 0.f; + sp->n_check_exact = 0.f; + sp->inhibited_check_exact = 0; +#endif } /** @@ -322,9 +332,10 @@ __attribute__((always_inline)) INLINE static void sink_kick_extra( __attribute__((always_inline)) INLINE static void sink_end_density( struct sink* si, const struct cosmology* cosmo) { - const float h = si->r_cut / kernel_gamma; - const float h_inv = 1.0f / h; /* 1/h */ - const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + const float h = si->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */ /* --- Finish the calculation by inserting the missing h factors --- */ si->to_collect.rho_gas *= h_inv_dim; @@ -336,6 +347,15 @@ __attribute__((always_inline)) INLINE static void sink_end_density( si->to_collect.velocity_gas[0] *= h_inv_dim * rho_inv; si->to_collect.velocity_gas[1] *= h_inv_dim * rho_inv; si->to_collect.velocity_gas[2] *= h_inv_dim * rho_inv; + + /* Finish the calculation by inserting the missing h-factors */ + si->density.wcount *= h_inv_dim; + si->density.wcount_dh *= h_inv_dim_plus_one; + +#ifdef SWIFT_SINK_DENSITY_CHECKS + si->rho_check *= h_inv_dim; + si->n_check *= h_inv_dim; +#endif } /** @@ -349,15 +369,22 @@ __attribute__((always_inline)) INLINE static void sinks_sink_has_no_neighbours( struct sink* restrict sp, const struct cosmology* cosmo) { warning( - "Sink particle with ID %lld treated as having no neighbours (r_cut: %g, " + "Sink particle with ID %lld treated as having no neighbours (h: %g, " "numb_ngbs: %i).", - sp->id, sp->r_cut, sp->num_ngbs); + sp->id, sp->h, sp->num_ngbs); + + /* Some smoothing length multiples. */ + const float h = sp->h; + const float h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ /* Reset problematic values */ sp->to_collect.velocity_gas[0] = sp->v[0]; sp->to_collect.velocity_gas[1] = sp->v[1]; sp->to_collect.velocity_gas[2] = sp->v[2]; - sp->to_collect.minimal_h_gas = sp->r_cut / kernel_gamma; + sp->to_collect.minimal_h_gas = h; + sp->density.wcount = kernel_root * h_inv_dim; + sp->density.wcount_dh = 0.f; } /** @@ -548,7 +575,11 @@ INLINE static void sink_copy_properties( sink_init_sink(sink); /* Set a smoothing length */ - sink->r_cut = sink_props->cut_off_radius; + if (sink_props->use_fixed_r_cut) { + sink->h = sink_props->cut_off_radius / kernel_gamma; + } else { + sink->h = p->h; + } /* Flag it as not swallowed */ sink_mark_sink_as_not_swallowed(&sink->merger_data); @@ -776,7 +807,7 @@ INLINE static int sink_spawn_star(struct sink* sink, const struct engine* e, * @brief Give the #spart a new position. * * In GEAR: Positions are set by randomly sampling coordinates in an homogeneous - * sphere centered on the #sink with radius the sink's r_cut. + * sphere centered on the #sink with radius the sink's r_cut. * * @param e The #engine. * @param si The #sink generating a star. @@ -798,8 +829,9 @@ INLINE static void sink_star_formation_give_new_position(const struct engine* e, const double phi = 2 * M_PI * random_unit_interval(sp->id, e->ti_current, (enum random_number_type)3); - const double r = si->r_cut * random_unit_interval(sp->id, e->ti_current, - (enum random_number_type)4); + const float rmax = si->h * kernel_gamma; + const double r = rmax * random_unit_interval(sp->id, e->ti_current, + (enum random_number_type)4); const double cos_theta = 1.0 - 2.0 * random_unit_interval(sp->id, e->ti_current, (enum random_number_type)5); @@ -837,8 +869,8 @@ INLINE static void sink_star_formation_give_new_velocity( and subtracted from the sink. */ double v_given[3] = {0.0, 0.0, 0.0}; const double G_newton = e->physical_constants->const_newton_G; - const double sigma_2 = - G_newton * si->mass_tot_before_star_spawning / si->r_cut; + const float rmax = si->h * kernel_gamma; + const double sigma_2 = G_newton * si->mass_tot_before_star_spawning / rmax; const double sigma = sink_props->star_spawning_sigma_factor * sqrt(sigma_2); for (int i = 0; i < 3; ++i) { @@ -900,7 +932,7 @@ INLINE static void sink_copy_properties_to_star( sink_star_formation_give_new_velocity(e, sink, sp, sink_props); /* Sph smoothing length */ - sp->h = sink->r_cut; + sp->h = sink->h; /* Feedback related initialisation */ /* ------------------------------- */ @@ -1205,7 +1237,8 @@ INLINE static void sink_prepare_part_sink_formation_sink_criteria( const float r_acc_p = sink_props->cut_off_radius * cosmo->a; /* Physical accretion radius of sink si */ - const float r_acc_si = si->r_cut * cosmo->a; + const float rmax = si->h * kernel_gamma; + const float r_acc_si = rmax * cosmo->a; /* Comoving distance of particl p */ const float px[3] = {(float)(p->x[0]), (float)(p->x[1]), (float)(p->x[2])}; diff --git a/src/sink/GEAR/sink_iact.h b/src/sink/GEAR/sink_iact.h index b2c0dd86f1..d52d6b94f9 100644 --- a/src/sink/GEAR/sink_iact.h +++ b/src/sink/GEAR/sink_iact.h @@ -42,21 +42,21 @@ * @param pj Second particle. * @param a Current scale factor. * @param H Current Hubble parameter. - * @param cut_off_radius Sink cut off radius. */ __attribute__((always_inline)) INLINE static void runner_iact_sink( const float r2, const float dx[3], const float hi, const float hj, struct part *restrict pi, struct part *restrict pj, const float a, - const float H, const float cut_off_radius) { + const float H) { - /* In order to prevent the formation of two sink particles at a distance - * smaller than the sink cutoff radius, we keep only gas particles with - * the smallest potential. */ + /* In order to prevent the formation of two sink particles too close together, + * we keep only gas particles with the smallest potential. The distance at + * which to prevent sink formation is the cutoff radius if this is fixed, or + * it is the variable smoothing length times gamma. */ const float r = sqrtf(r2); + const float rmax = max(hi, hj) * kernel_gamma; - /* if the distance is less than the cut off radius */ - if (r < cut_off_radius) { + if (r < rmax) { float potential_i = pi->sink_data.potential; float potential_j = pj->sink_data.potential; @@ -88,20 +88,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_sink( * @param pj Second particle (not updated). * @param a Current scale factor. * @param H Current Hubble parameter. - * @param cut_off_radius Sink cut off radius. */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( const float r2, const float dx[3], const float hi, const float hj, struct part *restrict pi, const struct part *restrict pj, const float a, - const float H, const float cut_off_radius) { + const float H) { - /* In order to prevent the formation of two sink particles at a distance - * smaller than the sink cutoff radius, we keep only gas particles with - * the smallest potential. */ + /* In order to prevent the formation of two sink particles too close together, + * we keep only gas particles with the smallest potential. The distance at + * which to prevent sink formation is the cutoff radius if this is fixed, or + * it is the variable smoothing length times gamma. */ const float r = sqrtf(r2); + const float rmax = max(hi, hj) * kernel_gamma; - if (r < cut_off_radius) { + if (r < rmax) { float potential_i = pi->sink_data.potential; float potential_j = pj->sink_data.potential; @@ -111,7 +112,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( } } -/* @brief Density interaction between two particles (non-symmetric). +/** + * @brief Density interaction between two particles (non-symmetric). + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing length of particle i. * @param hj Comoving smoothing-length of particle j. * @param si First particle (sink). * @param pj Second particle (gas, not updated). @@ -124,7 +130,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sink( */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sinks_gas_density( - const float r2, const float dx[3], const float ri, const float hj, + const float r2, const float dx[3], const float hi, const float hj, struct sink *si, const struct part *pj, const int with_cosmology, const struct cosmology *cosmo, const struct gravity_props *grav_props, const struct sink_props *sink_props, const integertime_t ti_current, @@ -137,7 +143,6 @@ runner_iact_nonsym_sinks_gas_density( /* Compute the kernel function */ const float r = sqrtf(r2); - const float hi = ri / kernel_gamma; const float hi_inv = 1.0f / hi; const float ui = r * hi_inv; kernel_deval(ui, &wi, &wi_dx); @@ -165,6 +170,16 @@ runner_iact_nonsym_sinks_gas_density( si->to_collect.velocity_gas[0] += mj * dv[0] * wi; si->to_collect.velocity_gas[1] += mj * dv[1] * wi; si->to_collect.velocity_gas[2] += mj * dv[2] * wi; + + /* Compute contribution to the number of neighbours */ + si->density.wcount += wi; + si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx); + +#ifdef SWIFT_SINK_DENSITY_CHECKS + si->rho_check += hydro_get_mass(pj) * wi; + si->n_check += wi; + si->N_check_density++; +#endif } /** @@ -175,14 +190,14 @@ runner_iact_nonsym_sinks_gas_density( * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). - * @param ri Comoving cut off radius of particle i. - * @param rj Comoving cut off radius of particle j. + * @param hi Comoving smoothing length of particle i. + * @param hj Comoving smoothing length of particle j. * @param si First sink particle. * @param sj Second sink particle. */ __attribute__((always_inline)) INLINE static void sink_collect_properties_from_sink(const float r2, const float dx[3], - const float ri, const float rj, + const float hi, const float hj, struct sink *restrict si, struct sink *restrict sj, const struct gravity_props *grav_props) { @@ -229,8 +244,8 @@ sink_collect_properties_from_sink(const float r2, const float dx[3], * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). - * @param ri Comoving cut off radius of particle i. - * @param rj Comoving cut off radius of particle j. + * @param hi Comoving smoothing length of particle i. + * @param hj Comoving smoothing length of particle j. * @param si First sink particle. * @param sj Second sink particle. * @param with_cosmology if we run with cosmology. @@ -242,15 +257,18 @@ sink_collect_properties_from_sink(const float r2, const float dx[3], */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sinks_sink_swallow( - const float r2, const float dx[3], const float ri, const float rj, + const float r2, const float dx[3], const float hi, const float hj, struct sink *restrict si, struct sink *restrict sj, const int with_cosmology, const struct cosmology *cosmo, const struct gravity_props *grav_props, const struct sink_props *sink_properties, const integertime_t ti_current, const double time) { + /* Convert the smoothing length back into a cutoff radius */ + const float hig = hi * kernel_gamma; + const float r = sqrtf(r2); - const float f_acc_r_acc_i = sink_properties->f_acc * ri; + const float f_acc_r_acc_i = sink_properties->f_acc * hig; /* Determine if the sink is dead, i.e. if its age is bigger than the age_threshold_unlimited */ @@ -264,7 +282,7 @@ runner_iact_nonsym_sinks_sink_swallow( they are both dead, we do not want to restrict the timesteps for 2-body encounters since they won't merge. */ if (!si_is_dead || !sj_is_dead) { - sink_collect_properties_from_sink(r2, dx, ri, rj, si, sj, grav_props); + sink_collect_properties_from_sink(r2, dx, hi, hj, si, sj, grav_props); } /* If si is dead, do not swallow sj. However, sj can swallow si if it alive. @@ -318,9 +336,9 @@ runner_iact_nonsym_sinks_sink_swallow( /* Momentum check------------------------------------------------------- */ float L2_j = 0.0; /* Relative momentum of the sink j */ float L2_kepler = 0.0; /* Keplerian angular momentum squared */ - sink_compute_angular_momenta_criterion(dx, v_plus_H_flow, r, si->r_cut, - si->mass, cosmo, grav_props, - &L2_kepler, &L2_j); + sink_compute_angular_momenta_criterion( + dx, v_plus_H_flow, r, si->h * kernel_gamma, si->mass, cosmo, grav_props, + &L2_kepler, &L2_j); /* To be accreted, the sink momentum should lower than the keplerian orbit * momentum. */ @@ -416,7 +434,7 @@ runner_iact_nonsym_sinks_sink_swallow( * * @param r2 Comoving square distance between the two particles. * @param dx Comoving vector separating both particles (pi - pj). - * @param ri Comoving cut off radius of particle i. + * @param hi Comoving smoothing length of particle i. * @param hj Comoving smoothing-length of particle j. * @param si First sink particle. * @param pj Second particle. @@ -429,15 +447,18 @@ runner_iact_nonsym_sinks_sink_swallow( */ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_sinks_gas_swallow( - const float r2, const float dx[3], const float ri, const float hj, + const float r2, const float dx[3], const float hi, const float hj, struct sink *restrict si, struct part *restrict pj, const int with_cosmology, const struct cosmology *cosmo, const struct gravity_props *grav_props, const struct sink_props *sink_properties, const integertime_t ti_current, const double time) { + /* Convert the smoothing length back into a cutoff radius */ + const float hig = hi * kernel_gamma; + const float r = sqrtf(r2); - const float f_acc_r_acc = sink_properties->f_acc * ri; + const float f_acc_r_acc = sink_properties->f_acc * hig; /* Determine if the sink is dead, i.e. if its age is bigger than the age_threshold_unlimited */ @@ -460,7 +481,7 @@ runner_iact_nonsym_sinks_gas_swallow( } /* f_acc*r_acc <= r <= r_acc, we perform other checks */ - } else if ((r >= f_acc_r_acc) && (r < ri)) { + } else if ((r >= f_acc_r_acc) && (r < hig)) { /* Relative velocity between the sinks */ const float dv[3] = {pj->v[0] - si->v[0], pj->v[1] - si->v[1], @@ -486,9 +507,9 @@ runner_iact_nonsym_sinks_gas_swallow( /* Momentum check------------------------------------------------------- */ float L2_gas_j = 0.0; /* Relative momentum of the gas */ float L2_kepler = 0.0; /* Keplerian angular momentum squared */ - sink_compute_angular_momenta_criterion(dx, v_plus_H_flow, r, si->r_cut, - si->mass, cosmo, grav_props, - &L2_kepler, &L2_gas_j); + sink_compute_angular_momenta_criterion( + dx, v_plus_H_flow, r, si->h * kernel_gamma, si->mass, cosmo, grav_props, + &L2_kepler, &L2_gas_j); /* To be accreted, the gas momentum should lower than the keplerian orbit * momentum. */ @@ -528,10 +549,10 @@ runner_iact_nonsym_sinks_gas_swallow( if (E_mec_sink_part >= 0) return; /* To be accreted, the gas smoothing length must be smaller than the sink - accretion radius. This is similar to AMR codes requesting the maximum + smoothing length. This is similar to AMR codes requesting the maximum refinement level close to the sink. */ if (sink_properties->sink_formation_smoothing_length_criterion && - (pj->h * kernel_gamma >= si->r_cut)) + (pj->h >= si->h)) return; /* Most bound pair check------------------------------------------------ */ diff --git a/src/sink/GEAR/sink_io.h b/src/sink/GEAR/sink_io.h index 293d58d501..fe86573da2 100644 --- a/src/sink/GEAR/sink_io.h +++ b/src/sink/GEAR/sink_io.h @@ -34,7 +34,7 @@ INLINE static void sink_read_particles(struct sink* sinks, struct io_props* list, int* num_fields) { /* Say how much we want to read */ - *num_fields = 5; + *num_fields = 6; /* List what we want to read */ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, @@ -45,7 +45,9 @@ INLINE static void sink_read_particles(struct sink* sinks, sinks, mass); list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY, UNIT_CONV_NO_UNITS, sinks, id); - list[4] = io_make_input_field("BirthTime", FLOAT, 1, OPTIONAL, UNIT_CONV_MASS, + list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL, + UNIT_CONV_LENGTH, sinks, h); + list[5] = io_make_input_field("BirthTime", FLOAT, 1, OPTIONAL, UNIT_CONV_MASS, sinks, birth_time); } @@ -144,21 +146,25 @@ INLINE static void sink_write_particles(const struct sink* sinks, "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, id, /*can convert to comoving=*/0, "Unique ID of the particles"); - list[4] = io_make_physical_output_field( + list[4] = io_make_output_field( + "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sinks, h, + "Co-moving smoothing lengths (FWHM of the kernel) of the particles"); + + list[5] = io_make_physical_output_field( "NumberOfSinkSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, number_of_sink_swallows, /*can convert to comoving=*/0, "Total number of sink merger events"); - list[5] = io_make_physical_output_field( + list[6] = io_make_physical_output_field( "NumberOfGasSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, number_of_gas_swallows, /*can convert to comoving=*/0, "Total number of gas merger events"); - list[6] = io_make_output_field_convert_sink( + list[7] = io_make_output_field_convert_sink( "TargetMass", FLOAT, 1, UNIT_CONV_MASS, 0.f, sinks, convert_sink_target_mass, "Sink target mass to spawn star particles"); - list[7] = io_make_physical_output_field( + list[8] = io_make_physical_output_field( "Nstars", INT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, n_stars, /*can convert to comoving=*/0, "Number of stars spawned by the sink particles"); @@ -166,7 +172,7 @@ INLINE static void sink_write_particles(const struct sink* sinks, /* Note: Since the swallowed momentum is computed with the physical velocity, i.e. including the Hubble flow term, it is not convertible to comoving frame. */ - list[8] = io_make_physical_output_field_convert_sink( + list[9] = io_make_physical_output_field_convert_sink( "SwallowedAngularMomentum", FLOAT, 3, UNIT_CONV_ANGULAR_MOMENTUM, 0.f, sinks, /*can convert to comoving=*/0, convert_sink_swallowed_angular_momentum, diff --git a/src/sink/GEAR/sink_part.h b/src/sink/GEAR/sink_part.h index c35fc7de57..1f53903977 100644 --- a/src/sink/GEAR/sink_part.h +++ b/src/sink/GEAR/sink_part.h @@ -46,8 +46,18 @@ struct sink { /*! Particle velocity. */ float v[3]; - /*! Cut off radius. */ - float r_cut; + /* Particle smoothing length, or r_cut/kernel_gamma if using a fixed cutoff*/ + float h; + + struct { + + /* Number of neighbours. */ + float wcount; + + /* Number of neighbours spatial derivative. */ + float wcount_dh; + + } density; /*! Sink particle mass */ float mass; @@ -173,6 +183,28 @@ struct sink { /*! List of interacting particles in compute formation SELF and PAIR */ long long ids_ngbs_accretion[MAX_NUM_OF_NEIGHBOURS_SINKS]; #endif + +#ifdef SWIFT_SINK_DENSITY_CHECKS + + /* Integer number of neighbours in the density loop */ + int N_check_density; + + /* Exact integer number of neighbours in the density loop */ + int N_check_density_exact; + + /*! Has this particle interacted with any unhibited neighbour? */ + char inhibited_check_exact; + + float n_check; + + float n_check_exact; + + float rho_check; + + /*! Exact value of the density field obtained via brute-force loop */ + float rho_check_exact; + +#endif } SWIFT_STRUCT_ALIGN; #endif /* SWIFT_GEAR_SINK_PART_H */ diff --git a/src/sink/GEAR/sink_properties.h b/src/sink/GEAR/sink_properties.h index 5b60b61d36..4134912e5f 100644 --- a/src/sink/GEAR/sink_properties.h +++ b/src/sink/GEAR/sink_properties.h @@ -41,6 +41,30 @@ */ struct sink_props { + /* ----- Basic neighbour search properties ------ */ + + /*! Resolution parameter */ + float eta_neighbours; + + /*! Target weightd number of neighbours (for info only)*/ + float target_neighbours; + + /*! Smoothing length tolerance */ + float h_tolerance; + + /*! Tolerance on neighbour number (for info only)*/ + float delta_neighbours; + + /*! Maximal number of iterations to converge h */ + int max_smoothing_iterations; + + /*! Maximal change of h over one time-step */ + float log_max_h_change; + + /*! Are we using a fixed cutoff radius? (all smoothing length calculations are + * disabled if so) */ + char use_fixed_r_cut; + /*! Cut off radius */ float cut_off_radius; @@ -200,13 +224,42 @@ INLINE static void sink_props_init_probabilities( * @param cosmo The cosmological model. * @param with_feedback Are we running with feedback? */ -INLINE static void sink_props_init(struct sink_props *sp, - struct feedback_props *fp, - const struct phys_const *phys_const, - const struct unit_system *us, - struct swift_params *params, - const struct cosmology *cosmo, - const int with_feedback) { +INLINE static void sink_props_init( + struct sink_props *sp, struct feedback_props *fp, + const struct phys_const *phys_const, const struct unit_system *us, + struct swift_params *params, const struct hydro_props *hydro_props, + const struct cosmology *cosmo, const int with_feedback) { + + /* Read in the basic neighbour search properties or default to the hydro + ones if the user did not provide any different values */ + + /* Kernel properties */ + sp->eta_neighbours = parser_get_opt_param_float( + params, "Sinks:resolution_eta", hydro_props->eta_neighbours); + + /* Tolerance for the smoothing length Newton-Raphson scheme */ + sp->h_tolerance = parser_get_opt_param_float(params, "Sinks:h_tolerance", + hydro_props->h_tolerance); + + /* Get derived properties */ + sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm; + const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance); + sp->delta_neighbours = + (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) * + kernel_norm; + + /* Number of iterations to converge h */ + sp->max_smoothing_iterations = + parser_get_opt_param_int(params, "Sinks:max_ghost_iterations", + hydro_props->max_smoothing_iterations); + + /* Time integration properties */ + const float max_volume_change = + parser_get_opt_param_float(params, "Sinks:max_volume_change", -1); + if (max_volume_change == -1) + sp->log_max_h_change = hydro_props->log_max_h_change; + else + sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv)); /* If we do not run with feedback, abort and print an error */ if (!with_feedback) @@ -214,9 +267,20 @@ INLINE static void sink_props_init(struct sink_props *sp, "ERROR: Running with sink but without feedback. GEAR sink model needs " "to be run with --sink and --feedback"); - /* Read the parameters from the parameter file */ - sp->cut_off_radius = - parser_get_param_float(params, "GEARSink:cut_off_radius"); + /* This property is used in all models to flag if we're using a fixed cutoff. + * If it is set to 1, we use a fixed r_cut (read in below) and don't + * (re)calculate h. If not, the code will use the variable h*kernel_gamma as a + * cutoff radius. + */ + sp->use_fixed_r_cut = + parser_get_param_char(params, "GEARSink:use_fixed_cut_off_radius"); + + /* The property cut_off_radius is now only used in the GEAR model. + * It is ignored if use_fixed_r_cut is 0. */ + if (sp->use_fixed_r_cut) { + sp->cut_off_radius = + parser_get_param_float(params, "GEARSink:cut_off_radius"); + } sp->f_acc = parser_get_opt_param_float(params, "GEARSink:f_acc", sink_gear_f_acc_default); diff --git a/src/sink_debug.h b/src/sink_debug.h index 2346711e73..44b1937d43 100644 --- a/src/sink_debug.h +++ b/src/sink_debug.h @@ -25,6 +25,8 @@ /* Import the debug routines of the right sink definition */ #if defined(SINK_NONE) #include "./sink/Default/sink_debug.h" +#elif defined(SINK_BASIC) +#include "./sink/Basic/sink_debug.h" #elif defined(SINK_GEAR) #include "./sink/GEAR/sink_debug.h" #else diff --git a/src/sink_iact.h b/src/sink_iact.h index d47b9971b7..4d5b6467f8 100644 --- a/src/sink_iact.h +++ b/src/sink_iact.h @@ -25,6 +25,8 @@ /* Select the correct sink model */ #if defined(SINK_NONE) #include "./sink/Default/sink_iact.h" +#elif defined(SINK_BASIC) +#include "./sink/Basic/sink_iact.h" #elif defined(SINK_GEAR) #include "./sink/GEAR/sink_iact.h" #else diff --git a/src/sink_io.h b/src/sink_io.h index f5b0be176c..685824f93b 100644 --- a/src/sink_io.h +++ b/src/sink_io.h @@ -24,6 +24,8 @@ /* Load the correct sink type */ #if defined(SINK_NONE) #include "./sink/Default/sink_io.h" +#elif defined(SINK_BASIC) +#include "./sink/Basic/sink_io.h" #elif defined(SINK_GEAR) #include "./sink/GEAR/sink_io.h" #else diff --git a/src/sink_properties.h b/src/sink_properties.h index 2ec9bd9d00..7700fa5395 100644 --- a/src/sink_properties.h +++ b/src/sink_properties.h @@ -25,6 +25,8 @@ /* Select the correct sink model */ #if defined(SINK_NONE) #include "./sink/Default/sink_properties.h" +#elif defined(SINK_BASIC) +#include "./sink/Basic/sink_properties.h" #elif defined(SINK_GEAR) #include "./sink/GEAR/sink_properties.h" #else diff --git a/src/sink_struct.h b/src/sink_struct.h index 00e6380298..7164c1cb20 100644 --- a/src/sink_struct.h +++ b/src/sink_struct.h @@ -33,6 +33,8 @@ /* Import the right black holes definition */ #if defined(SINK_NONE) #include "./sink/Default/sink_struct.h" +#elif defined(SINK_BASIC) +#include "./sink/Basic/sink_struct.h" #elif defined(SINK_GEAR) #include "./sink/GEAR/sink_struct.h" #else diff --git a/src/space_regrid.c b/src/space_regrid.c index 1ebc1074d0..d11ec4b65d 100644 --- a/src/space_regrid.c +++ b/src/space_regrid.c @@ -64,9 +64,8 @@ void space_regrid(struct space *s, int verbose) { if (c->black_holes.h_max > h_max) { h_max = c->black_holes.h_max; } - /* Note: For sinks, r_cut <=> h*kernel_gamma */ - if (c->sinks.r_cut_max > h_max * kernel_gamma) { - h_max = c->sinks.r_cut_max / kernel_gamma; + if (c->sinks.h_max > h_max) { + h_max = c->sinks.h_max; } } @@ -83,10 +82,8 @@ void space_regrid(struct space *s, int verbose) { if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) { h_max = c->black_holes.h_max; } - /* Note: For sinks, r_cut <=> h*kernel_gamma */ - if (c->nodeID == engine_rank && - c->sinks.r_cut_max > h_max * kernel_gamma) { - h_max = c->sinks.r_cut_max / kernel_gamma; + if (c->nodeID == engine_rank && c->sinks.h_max > h_max) { + h_max = c->sinks.h_max; } } @@ -102,7 +99,7 @@ void space_regrid(struct space *s, int verbose) { if (s->bparts[k].h > h_max) h_max = s->bparts[k].h; } for (size_t k = 0; k < nr_sinks; k++) { - if (s->sinks[k].r_cut > h_max) h_max = s->sinks[k].r_cut / kernel_gamma; + if (s->sinks[k].h > h_max) h_max = s->sinks[k].h; } } } diff --git a/src/space_split.c b/src/space_split.c index 43e020d5e5..d7495d1ac9 100644 --- a/src/space_split.c +++ b/src/space_split.c @@ -240,8 +240,8 @@ void space_split_recursive(struct space *s, struct cell *c, cp->stars.h_max_active = 0.f; cp->stars.dx_max_part = 0.f; cp->stars.dx_max_sort = 0.f; - cp->sinks.r_cut_max = 0.f; - cp->sinks.r_cut_max_active = 0.f; + cp->sinks.h_max = 0.f; + cp->sinks.h_max_active = 0.f; cp->sinks.dx_max_part = 0.f; cp->black_holes.h_max = 0.f; cp->black_holes.h_max_active = 0.f; @@ -305,9 +305,8 @@ void space_split_recursive(struct space *s, struct cell *c, black_holes_h_max = max(black_holes_h_max, cp->black_holes.h_max); black_holes_h_max_active = max(black_holes_h_max_active, cp->black_holes.h_max_active); - sinks_h_max = max(sinks_h_max, cp->sinks.r_cut_max); - sinks_h_max_active = - max(sinks_h_max_active, cp->sinks.r_cut_max_active); + sinks_h_max = max(sinks_h_max, cp->sinks.h_max); + sinks_h_max_active = max(sinks_h_max_active, cp->sinks.h_max_active); ti_hydro_end_min = min(ti_hydro_end_min, cp->hydro.ti_end_min); ti_hydro_beg_max = max(ti_hydro_beg_max, cp->hydro.ti_beg_max); @@ -579,10 +578,10 @@ void space_split_recursive(struct space *s, struct cell *c, ti_sinks_end_min = min(ti_sinks_end_min, ti_end); ti_sinks_beg_max = max(ti_sinks_beg_max, ti_beg); - sinks_h_max = max(sinks_h_max, sinks[k].r_cut); + sinks_h_max = max(sinks_h_max, sinks[k].h); if (sink_is_active(&sinks[k], e)) - sinks_h_max_active = max(sinks_h_max_active, sinks[k].r_cut); + sinks_h_max_active = max(sinks_h_max_active, sinks[k].h); /* Reset x_diff */ sinks[k].x_diff[0] = 0.f; @@ -666,8 +665,8 @@ void space_split_recursive(struct space *s, struct cell *c, c->stars.h_max_active = stars_h_max_active; c->sinks.ti_end_min = ti_sinks_end_min; c->sinks.ti_beg_max = ti_sinks_beg_max; - c->sinks.r_cut_max = sinks_h_max; - c->sinks.r_cut_max_active = sinks_h_max_active; + c->sinks.h_max = sinks_h_max; + c->sinks.h_max_active = sinks_h_max_active; c->black_holes.ti_end_min = ti_black_holes_end_min; c->black_holes.ti_beg_max = ti_black_holes_beg_max; c->black_holes.h_max = black_holes_h_max; diff --git a/src/tools.c b/src/tools.c index 03c4e29325..65d7d976cc 100644 --- a/src/tools.c +++ b/src/tools.c @@ -243,8 +243,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) { runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hi, pj->h, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dx, hi, pj->h, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dx, hi, pj->h, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hi, pj->h, pi, pj, a, H); } } } @@ -279,8 +278,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) { runner_iact_nonsym_chemistry(r2, dx, hj, pi->h, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dx, hj, pi->h, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dx, hj, pi->h, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dx, hj, pi->h, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dx, hj, pi->h, pj, pi, a, H); } } } @@ -559,8 +557,7 @@ void self_all_density(struct runner *r, struct cell *ci) { runner_iact_nonsym_chemistry(r2, dxi, hi, hj, pi, pj, a, H); runner_iact_nonsym_pressure_floor(r2, dxi, hi, hj, pi, pj, a, H); runner_iact_nonsym_star_formation(r2, dxi, hi, hj, pi, pj, a, H); - runner_iact_nonsym_sink(r2, dxi, hi, hj, pi, pj, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dxi, hi, hj, pi, pj, a, H); } /* Hit or miss? */ @@ -575,8 +572,7 @@ void self_all_density(struct runner *r, struct cell *ci) { runner_iact_nonsym_chemistry(r2, dxi, hj, hi, pj, pi, a, H); runner_iact_nonsym_pressure_floor(r2, dxi, hj, hi, pj, pi, a, H); runner_iact_nonsym_star_formation(r2, dxi, hj, hi, pj, pi, a, H); - runner_iact_nonsym_sink(r2, dxi, hj, hi, pj, pi, a, H, - e->sink_properties->cut_off_radius); + runner_iact_nonsym_sink(r2, dxi, hj, hi, pj, pi, a, H); } } } diff --git a/swift.c b/swift.c index 41c96e0afb..17890b5313 100644 --- a/swift.c +++ b/swift.c @@ -1174,7 +1174,7 @@ int main(int argc, char *argv[]) { /* Initialise the sink properties */ if (with_sinks) { sink_props_init(&sink_properties, &feedback_properties, &prog_const, &us, - params, &cosmo, with_feedback); + params, &hydro_properties, &cosmo, with_feedback); } else bzero(&sink_properties, sizeof(struct sink_props)); -- GitLab