diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
new file mode 100644
index 0000000000000000000000000000000000000000..cf96535b64c078d063b6f773823ca085876c2ecc
--- /dev/null
+++ b/examples/ExternalPointMass/externalPointMass.yml
@@ -0,0 +1,55 @@
+
+# Define the system of units to use internally.
+UnitSystem:
+ UnitMass_in_cgs: 1.9885e33 # Grams
+ UnitLength_in_cgs: 3.0856776e21 # Centimeters
+ UnitVelocity_in_cgs: 1e5 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters for the task scheduling
+Scheduler:
+ nr_queues: 0 # The number of task queues to use. Use 0 to let the system decide.
+ cell_max_size: 8000000 # Maximal number of interactions per task (this is the default value).
+ cell_sub_size: 8000000 # Maximal number of interactions per sub-task (this is the default value).
+ cell_split_size: 400 # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1. # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
+ max_smoothing_length: 10. # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: Sphere.hdf5 # The file to read
+ h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
+ shift_x: 50. # A shift to apply to all particles read from the ICs (in internal units).
+ shift_y: 50.
+ shift_z: 50.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+ initial_type: m # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+ initial_grid_x: 10 # Grid size if the 'g' strategy is chosen.
+ initial_grid_y: 10
+ initial_grid_z: 10
+ repartition_type: b # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+
+
+# External potential parameters
+PointMass:
+ position_x: 50. # location of external point mass in internal units
+ position_y: 50.
+ position_z: 50.
+ mass: 1e10 # mass of external point mass in internal units
+
diff --git a/examples/ExternalPointMass/makeIC.py b/examples/ExternalPointMass/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..b1b032e47383d6e757306f09063a8931635ea8e9
--- /dev/null
+++ b/examples/ExternalPointMass/makeIC.py
@@ -0,0 +1,142 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ # Tom Theuns (tom.theuns@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+import sys
+import numpy
+import math
+import random
+
+# Generates a random distriution of particles, for motion in an external potnetial centred at (0,0,0)
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS = 6.672e-8
+SOLAR_MASS_IN_CGS = 1.9885e33
+PARSEC_IN_CGS = 3.0856776e18
+
+# choice of units
+const_unit_length_in_cgs = (1000*PARSEC_IN_CGS)
+const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
+const_unit_velocity_in_cgs = (1e5)
+
+print "UnitMass_in_cgs: ", const_unit_mass_in_cgs
+print "UnitLength_in_cgs: ", const_unit_length_in_cgs
+print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
+
+# derived units
+const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+const_G = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
+print 'G=', const_G
+
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 100. #
+Radius = boxSize / 4. # maximum radius of particles
+G = const_G
+Mass = 1e10
+
+N = int(sys.argv[1]) # Number of particles
+L = N**(1./3.)
+
+# these are not used but necessary for I/O
+rho = 2. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
+fileName = "Sphere.hdf5"
+
+
+#---------------------------------------------------
+numPart = N
+mass = 1
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#Particle group
+#grp0 = file.create_group("/PartType0")
+grp1 = file.create_group("/PartType1")
+#generate particle positions
+radius = Radius * (numpy.random.rand(N))**(1./3.)
+ctheta = -1. + 2 * numpy.random.rand(N)
+stheta = numpy.sqrt(1.-ctheta**2)
+phi = 2 * math.pi * numpy.random.rand(N)
+r = numpy.zeros((numPart, 3))
+# r[:,0] = radius * stheta * numpy.cos(phi)
+# r[:,1] = radius * stheta * numpy.sin(phi)
+# r[:,2] = radius * ctheta
+r[:,0] = radius
+#
+speed = numpy.sqrt(G * Mass / radius)
+v = numpy.zeros((numPart, 3))
+omega = speed / radius
+period = 2.*math.pi/omega
+print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+v[:,0] = -omega * r[:,1]
+v[:,1] = omega * r[:,0]
+
+ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = numpy.zeros(1)
+
+m = numpy.full((numPart, ), mass)
+ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ds[()] = m
+m = numpy.zeros(1)
+
+h = numpy.full((numPart, ), 1.1255 * boxSize / L)
+ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
+ds[()] = h
+h = numpy.zeros(1)
+
+u = numpy.full((numPart, ), internalEnergy)
+ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
+ds[()] = u
+u = numpy.zeros(1)
+
+
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
+ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ds[()] = ids
+
+ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = r
+
+file.close()
diff --git a/examples/ExternalPointMass/run.sh b/examples/ExternalPointMass/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..45f0e697b3ebfd5d510ee1f9f895ed3d7e4b1caf
--- /dev/null
+++ b/examples/ExternalPointMass/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e uniformBox.hdf5 ]
+then
+ echo "Generating initial conditions for the point mass potential box example..."
+ python makeIC.py 10000
+fi
+
+../swift -g -t 2 externalPointMass.yml
diff --git a/examples/ExternalPointMass/test.pro b/examples/ExternalPointMass/test.pro
new file mode 100644
index 0000000000000000000000000000000000000000..21c10e9d27daa45b085c6a659ba3cf7260f017fb
--- /dev/null
+++ b/examples/ExternalPointMass/test.pro
@@ -0,0 +1,65 @@
+;
+; test energy / angular momentum conservation of test problem
+;
+@physunits
+
+indir = '/gpfs/data/tt/Codes/Swift-git/swiftsim/examples/'
+basefile = 'output_'
+nfiles = 657
+nfollow = 100 ; number of particles to follow
+eout = fltarr(nfollow, nfiles)
+ekin = fltarr(nfollow, nfiles)
+epot = fltarr(nfollow, nfiles)
+tout = fltarr(nfiles)
+; set properties of potential
+uL = 1e3 * phys.pc ; unit of length
+uM = phys.msun ; unit of mass
+uV = 1d5 ; unit of velocity
+
+; derived units
+constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre = [50.,50.,50.] * 1d3 * pc / uL
+mextern = 1d10 * msun / uM
+;
+;
+;
+ifile = 0
+for ifile=0,nfiles-1 do begin
+;for ifile=0,3 do begin
+ inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+ time = h5ra(inf, 'Header','Time')
+ p = h5rd(inf,'PartType1/Coordinates')
+ v = h5rd(inf,'PartType1/Velocities')
+ id = h5rd(inf,'PartType1/ParticleIDs')
+ indx = sort(id)
+;
+ id = id[indx]
+ for ic=0,2 do begin
+ tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
+ tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
+ endfor
+; calculate energy
+ dd = size(p,/dimen) & npart = dd[1]
+ ener = fltarr(npart)
+ dr = fltarr(npart) & dv = dr
+ for ic=0,2 do dr[*] = dr[*] + (p[ic,*]-pcentre[ic])^2
+ for ic=0,2 do dv[*] = dv[*] + v[ic,*]^2
+ dr = sqrt(dr)
+; print,'time = ',time,p[0,0],v[0,0],id[0]
+ ek = 0.5 * dv
+ ep = - constG * mextern / dr
+ ener = ek + ep
+ tout(ifile) = time
+ eout(*,ifile) = ener[0:nfollow-1]
+ ekin(*,ifile) = ek[0:nfollow-1]
+ epot(*,ifile) = ep[0:nfollow-1]
+endfor
+
+; calculate relative energy change
+de = 0.0 * eout
+for ifile=1, nfiles -1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
+
+
+end
+
+
diff --git a/examples/main.c b/examples/main.c
index 5a20125fe84bac5b793a2e1e2694cdfd76042051..d4538473031bc410606efc79ab6451d6273f61ad 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -3,6 +3,9 @@
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * Angus Lepper (angus.lepper@ed.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -268,17 +271,6 @@ int main(int argc, char *argv[]) {
MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
- /* Initialize unit system */
- struct UnitSystem us;
- units_init(&us, params);
- if (myrank == 0) {
- message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
- message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
- message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
- message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
- message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
- }
-
/* Prepare the domain decomposition scheme */
#ifdef WITH_MPI
struct partition initial_partition;
@@ -296,6 +288,25 @@ int main(int argc, char *argv[]) {
}
#endif
+ /* Initialize unit system and constants */
+ struct UnitSystem us;
+ struct phys_const prog_const;
+ units_init(&us, params);
+ phys_const_init(&us, &prog_const);
+ if (myrank == 0 && verbose > 0) {
+ message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
+ message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
+ message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
+ message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
+ message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
+ phys_const_print(&prog_const);
+ }
+
+ /* Initialise the external potential properties */
+ struct external_potential potential;
+ if (with_external_gravity) potential_init(params, &us, &potential);
+ if (with_external_gravity && myrank == 0) potential_print(&potential);
+
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
@@ -333,6 +344,12 @@ int main(int argc, char *argv[]) {
for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
Ngpart = 0;
}
+ if (!with_hydro) {
+ free(parts);
+ parts = NULL;
+ for (size_t k = 0; k < Ngpart; ++k) if(gparts[k].id > 0) error("Linking problem");
+ Ngas = 0;
+ }
/* Get the total number of particles across all nodes. */
long long N_total[2] = {0, 0};
@@ -396,7 +413,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, engine_policies,
- talking);
+ talking, &prog_const, &potential);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index e3cd5b209c9d36f9774364a661b77f3d649c398e..228e748f3b2899639303597c92ea29ad09987313 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -45,3 +45,11 @@ DomainDecomposition:
initial_grid_z: 10
repartition_type: b # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+# Parameters related to external potentials
+
+# Point mass external potential
+PointMass:
+ position_x: 50. # location of external point mass in internal units
+ position_y: 50.
+ position_z: 50.
+ mass: 1e10 # mass of external point mass in internal units
diff --git a/src/Makefile.am b/src/Makefile.am
index a96f35b3cf0d8a23aec4f8c0f8d16bec8638cbcd..65f093b014ed16577c4f0f36bf4636da223e6929 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -35,13 +35,15 @@ endif
# List required headers
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
- common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h
+ common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
+ physical_constants.h physical_constants_cgs.h potentials.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
- kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c
+ kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c \
+ physical_constants.c potentials.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/cell.c b/src/cell.c
index 8e47b415ac5b43558fea9a557cc352cc6dea2018..793a2e10c23cb0dd355a181e2ad658ce19dc0014 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1,6 +1,10 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
diff --git a/src/cell.h b/src/cell.h
index a471eac44bfd3533c4220ab8c5ff2ddec724e87f..fbdaa928bda51aa76c149df1249d982fbeef64c3 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -1,6 +1,10 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -122,6 +126,9 @@ struct cell {
/* Tasks for gravity tree. */
struct task *grav_up, *grav_down;
+ /* Task for external gravity */
+ struct task *grav_external;
+
/* Number of tasks that are associated with this cell. */
int nr_tasks;
diff --git a/src/const.h b/src/const.h
index 6a52ec4796a4904629a57ffa8b32a3107bde263e..dd3e34016c1953efffca67e9c7ee1b42f455dec1 100644
--- a/src/const.h
+++ b/src/const.h
@@ -70,4 +70,7 @@
#define GADGET2_SPH
//#define DEFAULT_SPH
+/* Gravity properties */
+#define EXTERNAL_POTENTIAL_POINTMASS
+
#endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index 7287f47b40b06a691c778a5f8d3b6c0f04721564..f022f95a2a51930a24ffa9835c5d8e46aa0088cd 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1,7 +1,11 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * Angus Lepper (angus.lepper@ed.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -101,9 +105,12 @@ void engine_make_ghost_tasks(struct engine *e, struct cell *c,
struct cell *super) {
struct scheduler *s = &e->sched;
+ const int is_with_external_gravity =
+ (e->policy & engine_policy_external_gravity) ==
+ engine_policy_external_gravity;
/* Am I the super-cell? */
- if (super == NULL && c->nr_tasks > 0) {
+ if (super == NULL && (c->count > 0 || c->gcount > 0)) {
/* Remember me. */
super = c;
@@ -111,18 +118,32 @@ void engine_make_ghost_tasks(struct engine *e, struct cell *c,
/* Local tasks only... */
if (c->nodeID == e->nodeID) {
- /* Generate the ghost task. */
- c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
- c, NULL, 0);
- /* Add the drift task. */
- c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
- c, NULL, 0);
/* Add the init task. */
c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0, c,
NULL, 0);
+
+ /* Add the drift task. */
+ c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
+ c, NULL, 0);
+
/* Add the kick task. */
c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0, c,
NULL, 0);
+
+ if (c->count > 0) {
+
+ /* Generate the ghost task. */
+ c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0,
+ 0, c, NULL, 0);
+ }
+
+ if (c->gcount > 0) {
+
+ /* Add the external gravity tasks */
+ if (is_with_external_gravity)
+ c->grav_external = scheduler_addtask(
+ s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0);
+ }
}
}
@@ -211,27 +232,29 @@ void engine_redistribute(struct engine *e) {
space_parts_sort(s, dest, s->nr_parts, 0, nr_nodes - 1, e->verbose);
/* We need to re-link the gpart partners of parts. */
- int current_dest = dest[0];
- size_t count_this_dest = 0;
- for (size_t k = 0; k < s->nr_parts; ++k) {
- if (s->parts[k].gpart != NULL) {
-
- /* As the addresses will be invalidated by the communications, we will */
- /* instead store the absolute index from the start of the sub-array */
- /* of particles to be sent to a given node. */
- /* Recall that gparts without partners have a negative id. */
- /* We will restore the pointers on the receiving node later on. */
- if (dest[k] != current_dest) {
- current_dest = dest[k];
- count_this_dest = 0;
- }
+ if (s->nr_parts > 0) {
+ int current_dest = dest[0];
+ size_t count_this_dest = 0;
+ for (size_t k = 0; k < s->nr_parts; ++k) {
+ if (s->parts[k].gpart != NULL) {
+
+ /* As the addresses will be invalidated by the communications, we will */
+ /* instead store the absolute index from the start of the sub-array */
+ /* of particles to be sent to a given node. */
+ /* Recall that gparts without partners have a negative id. */
+ /* We will restore the pointers on the receiving node later on. */
+ if (dest[k] != current_dest) {
+ current_dest = dest[k];
+ count_this_dest = 0;
+ }
- /* Debug */
- /* if(s->parts[k].gpart->id < 0) */
- /* error("Trying to link a partnerless gpart !"); */
+ /* Debug */
+ /* if(s->parts[k].gpart->id < 0) */
+ /* error("Trying to link a partnerless gpart !"); */
- s->parts[k].gpart->id = count_this_dest;
- count_this_dest++;
+ s->parts[k].gpart->id = count_this_dest;
+ count_this_dest++;
+ }
}
}
@@ -819,7 +842,10 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
ticks tic = getticks();
/* Re-set the proxies. */
- for (int k = 0; k < e->nr_proxies; k++) e->proxies[k].nr_parts_out = 0;
+ for (int k = 0; k < e->nr_proxies; k++) {
+ e->proxies[k].nr_parts_out = 0;
+ e->proxies[k].nr_gparts_out = 0;
+ }
/* Put the parts and gparts into the corresponding proxies. */
for (size_t k = 0; k < *Npart; k++) {
@@ -839,7 +865,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
/* Re-link the associated gpart with the buffer offset of the part. */
if (s->parts[offset_parts + k].gpart != NULL) {
- s->parts[offset_parts + k].gpart->id = e->proxies[pid].nr_parts_in;
+ s->parts[offset_parts + k].gpart->id = e->proxies[pid].nr_parts_out;
}
/* Load the part and xpart into the proxy. */
@@ -898,6 +924,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
count_parts_in, count_gparts_in);
}
if (offset_parts + count_parts_in > s->size_parts) {
+ message("re-allocating parts array.");
s->size_parts = (offset_parts + count_parts_in) * engine_parts_size_grow;
struct part *parts_new = NULL;
struct xpart *xparts_new = NULL;
@@ -912,8 +939,14 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
free(s->xparts);
s->parts = parts_new;
s->xparts = xparts_new;
+ for (size_t k = 0; k < offset_parts; k++) {
+ if (s->parts[k].gpart != NULL) {
+ s->parts[k].gpart->part = &s->parts[k];
+ }
+ }
}
if (offset_gparts + count_gparts_in > s->size_gparts) {
+ message("re-allocating gparts array.");
s->size_gparts = (offset_gparts + count_gparts_in) * engine_parts_size_grow;
struct gpart *gparts_new = NULL;
if (posix_memalign((void **)&gparts_new, gpart_align,
@@ -922,6 +955,11 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
memcpy(gparts_new, s->gparts, sizeof(struct gpart) * offset_gparts);
free(s->gparts);
s->gparts = gparts_new;
+ for (size_t k = 0; k < offset_gparts; k++) {
+ if (s->gparts[k].id > 0) {
+ s->gparts[k].part->gpart = &s->gparts[k];
+ }
+ }
}
/* Collect the requests for the particle data from the proxies. */
@@ -976,7 +1014,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
reqs_in[pid + 1] == MPI_REQUEST_NULL &&
reqs_in[pid + 2] == MPI_REQUEST_NULL) {
/* Copy the particle data to the part/xpart/gpart arrays. */
- struct proxy *p = &e->proxies[pid >> 1];
+ struct proxy *p = &e->proxies[pid / 3];
memcpy(&s->parts[offset_parts + count_parts], p->parts_in,
sizeof(struct part) * p->nr_parts_in);
memcpy(&s->xparts[offset_parts + count_parts], p->xparts_in,
@@ -993,7 +1031,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
for (int k = 0; k < p->nr_gparts_in; k++) {
struct gpart *gp = &s->gparts[offset_gparts + count_gparts + k];
if (gp->id >= 0) {
- struct part *p = &s->parts[offset_gparts + count_parts + gp->id];
+ struct part *p = &s->parts[offset_parts + count_parts + gp->id];
gp->part = p;
p->gpart = gp;
}
@@ -1004,7 +1042,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
count_gparts += p->nr_gparts_in;
}
}
-
+
/* Wait for all the sends to have finished too. */
if (nr_out > 0)
if (MPI_Waitall(3 * e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) !=
@@ -1276,6 +1314,12 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
/* /\* Kick tasks should rely on the grav_down tasks of their cell. *\/ */
/* else if (t->type == task_type_kick && t->ci->grav_down != NULL) */
/* scheduler_addunlock(sched, t->ci->grav_down, t); */
+
+ /* External gravity tasks should depend on init and unlock the kick */
+ else if (t->type == task_type_grav_external) {
+ scheduler_addunlock(sched, t->ci->init, t);
+ scheduler_addunlock(sched, t, t->ci->kick);
+ }
}
}
@@ -1377,7 +1421,7 @@ void engine_maketasks(struct engine *e) {
engine_make_hydroloop_tasks(e);
/* Add the gravity mm tasks. */
- if ((e->policy & engine_policy_self_gravity) == engine_policy_self_gravity)
+ if (e->policy & engine_policy_self_gravity)
engine_make_gravityinteraction_tasks(e);
/* Split the tasks. */
@@ -1393,7 +1437,7 @@ void engine_maketasks(struct engine *e) {
e->nr_links = 0;
/* Add the gravity up/down tasks at the top-level cells and push them down. */
- if ((e->policy & engine_policy_self_gravity) == engine_policy_self_gravity)
+ if (e->policy & engine_policy_self_gravity)
engine_make_gravityrecursive_tasks(e);
/* Count the number of tasks associated with each cell and
@@ -1412,7 +1456,7 @@ void engine_maketasks(struct engine *e) {
engine_make_extra_hydroloop_tasks(e);
/* Add the communication tasks if MPI is being used. */
- if ((e->policy & engine_policy_mpi) == engine_policy_mpi) {
+ if (e->policy & engine_policy_mpi) {
/* Loop over the proxies. */
for (int pid = 0; pid < e->nr_proxies; pid++) {
@@ -1465,7 +1509,7 @@ int engine_marktasks(struct engine *e) {
const ticks tic = getticks();
/* Much less to do here if we're on a fixed time-step. */
- if ((e->policy & engine_policy_fixdt) == engine_policy_fixdt) {
+ if (e->policy & engine_policy_fixdt) {
/* Run through the tasks and mark as skip or not. */
for (int k = 0; k < nr_tasks; k++) {
@@ -1754,7 +1798,7 @@ void engine_collect_kick(struct cell *c) {
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
/* Only do something is the cell is non-empty */
- if (c->count != 0) {
+ if (c->count != 0 || c->gcount != 0) {
/* If this cell is not split, I'm in trouble. */
if (!c->split) error("Cell has no super-cell.");
@@ -1853,13 +1897,11 @@ void engine_init_particles(struct engine *e) {
/* Make sure all particles are ready to go */
/* i.e. clean-up any stupid state in the ICs */
- if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
+ if (e->policy & engine_policy_hydro) {
space_map_cells_pre(s, 1, cell_init_parts, NULL);
}
- if (((e->policy & engine_policy_self_gravity) ==
- engine_policy_self_gravity) ||
- ((e->policy & engine_policy_external_gravity) ==
- engine_policy_external_gravity)) {
+ if ((e->policy & engine_policy_self_gravity) ||
+ (e->policy & engine_policy_external_gravity)) {
space_map_cells_pre(s, 1, cell_init_gparts, NULL);
}
@@ -1867,22 +1909,17 @@ void engine_init_particles(struct engine *e) {
engine_marktasks(e);
- // printParticle(e->s->parts, 1000, e->s->nr_parts);
- // printParticle(e->s->parts, 515050, e->s->nr_parts);
-
- // message("\n0th DENSITY CALC\n");
-
/* Build the masks corresponding to the policy */
unsigned int mask = 0;
unsigned int submask = 0;
/* We always have sort tasks */
mask |= 1 << task_type_sort;
+ mask |= 1 << task_type_init;
/* Add the tasks corresponding to hydro operations to the masks */
- if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
+ if (e->policy & engine_policy_hydro) {
- mask |= 1 << task_type_init;
mask |= 1 << task_type_self;
mask |= 1 << task_type_pair;
mask |= 1 << task_type_sub;
@@ -1892,20 +1929,19 @@ void engine_init_particles(struct engine *e) {
}
/* Add the tasks corresponding to self-gravity to the masks */
- if ((e->policy & engine_policy_self_gravity) == engine_policy_self_gravity) {
+ if (e->policy & engine_policy_self_gravity) {
/* Nothing here for now */
}
- /* Add the tasks corresponding to self-gravity to the masks */
- if ((e->policy & engine_policy_external_gravity) ==
- engine_policy_external_gravity) {
+ /* Add the tasks corresponding to external gravity to the masks */
+ if (e->policy & engine_policy_external_gravity) {
- /* Nothing here for now */
+ mask |= 1 << task_type_grav_external;
}
/* Add MPI tasks if need be */
- if ((e->policy & engine_policy_mpi) == engine_policy_mpi) {
+ if (e->policy & engine_policy_mpi) {
mask |= 1 << task_type_send;
mask |= 1 << task_type_recv;
@@ -1916,14 +1952,9 @@ void engine_init_particles(struct engine *e) {
engine_launch(e, e->nr_threads, mask, submask);
TIMER_TOC(timer_runners);
- // message("\n0th ENTROPY CONVERSION\n")
-
/* Apply some conversions (e.g. internal energy -> entropy) */
space_map_cells_pre(s, 1, cell_convert_hydro, NULL);
- // printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts);
- // printParticle(e->s->parts, e->s->xparts,515050, e->s->nr_parts);
-
/* Ready to go */
e->step = -1;
}
@@ -2005,8 +2036,6 @@ void engine_step(struct engine *e) {
}
#endif
- // message("\nDRIFT\n");
-
/* Move forward in time */
e->ti_old = e->ti_current;
e->ti_current = ti_end_min;
@@ -2018,11 +2047,6 @@ void engine_step(struct engine *e) {
/* Drift everybody */
engine_launch(e, e->nr_threads, 1 << task_type_drift, 0);
- // printParticle(e->s->parts, e->s->xparts, 1000, e->s->nr_parts);
- // printParticle(e->s->parts, e->s->xparts, 515050, e->s->nr_parts);
-
- // if(e->step == 2) exit(0);
-
if (e->nodeID == 0) {
/* Print some information to the screen */
@@ -2037,8 +2061,6 @@ void engine_step(struct engine *e) {
fflush(e->file_stats);
}
- // message("\nACCELERATION AND KICK\n");
-
/* Re-distribute the particles amongst the nodes? */
if (e->forcerepart != REPART_NONE) engine_repartition(e);
@@ -2051,11 +2073,11 @@ void engine_step(struct engine *e) {
/* We always have sort tasks and kick tasks */
mask |= 1 << task_type_sort;
mask |= 1 << task_type_kick;
+ mask |= 1 << task_type_init;
/* Add the tasks corresponding to hydro operations to the masks */
- if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
+ if (e->policy & engine_policy_hydro) {
- mask |= 1 << task_type_init;
mask |= 1 << task_type_self;
mask |= 1 << task_type_pair;
mask |= 1 << task_type_sub;
@@ -2066,20 +2088,18 @@ void engine_step(struct engine *e) {
}
/* Add the tasks corresponding to self-gravity to the masks */
- if ((e->policy & engine_policy_self_gravity) == engine_policy_self_gravity) {
+ if (e->policy & engine_policy_self_gravity) {
/* Nothing here for now */
}
- /* Add the tasks corresponding to self-gravity to the masks */
- if ((e->policy & engine_policy_external_gravity) ==
- engine_policy_external_gravity) {
-
- /* Nothing here for now */
+ /* Add the tasks corresponding to external gravity to the masks */
+ if (e->policy & engine_policy_external_gravity) {
+ mask |= 1 << task_type_grav_external;
}
/* Add MPI tasks if need be */
- if ((e->policy & engine_policy_mpi) == engine_policy_mpi) {
+ if (e->policy & engine_policy_mpi) {
mask |= 1 << task_type_send;
mask |= 1 << task_type_recv;
@@ -2095,8 +2115,6 @@ void engine_step(struct engine *e) {
clocks_gettime(&time2);
e->wallclock_time = (float)clocks_diff(&time1, &time2);
- // printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts);
- // printParticle(e->s->parts, e->s->xparts,515050, e->s->nr_parts);
}
/**
@@ -2326,11 +2344,15 @@ static bool hyperthreads_present(void) {
* @param nr_threads The number of threads per MPI rank.
* @param policy The queuing policy to use.
* @param verbose Is this #engine talkative ?
+ * @param physical_constants The #phys_const used for this run.
+ * @param potential The properties of the external potential.
*/
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
- int nr_threads, int policy, int verbose) {
+ int nr_threads, int policy, int verbose,
+ const struct phys_const *physical_constants,
+ const struct external_potential *potential) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
@@ -2361,6 +2383,8 @@ void engine_init(struct engine *e, struct space *s,
e->verbose = verbose;
e->count_step = 0;
e->wallclock_time = 0.f;
+ e->physical_constants = physical_constants;
+ e->external_potential = potential;
engine_rank = nodeID;
/* Make the space link back to the engine. */
@@ -2466,7 +2490,7 @@ void engine_init(struct engine *e, struct space *s,
engine_print_policy(e);
/* Print information about the hydro scheme */
- if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
+ if (e->policy & engine_policy_hydro) {
if (e->nodeID == 0) message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
if (e->nodeID == 0)
message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours.",
@@ -2492,7 +2516,7 @@ void engine_init(struct engine *e, struct space *s,
e->ti_current = 0;
/* Fixed time-step case */
- if ((e->policy & engine_policy_fixdt) == engine_policy_fixdt) {
+ if (e->policy & engine_policy_fixdt) {
e->dt_min = e->dt_max;
/* Find timestep on the timeline */
@@ -2572,7 +2596,7 @@ void engine_init(struct engine *e, struct space *s,
if (pthread_create(&e->runners[k].thread, NULL, &runner_main,
&e->runners[k]) != 0)
error("Failed to create runner thread.");
- if ((e->policy & engine_policy_setaffinity) == engine_policy_setaffinity) {
+ if (e->policy & engine_policy_setaffinity) {
#if defined(HAVE_SETAFFINITY)
/* Set a reasonable queue ID. */
@@ -2598,8 +2622,8 @@ void engine_init(struct engine *e, struct space *s,
e->runners[k].cpuid = k;
e->runners[k].qid = k * nr_queues / e->nr_threads;
}
- // message( "runner %i on cpuid=%i with qid=%i." , e->runners[k].id ,
- // e->runners[k].cpuid , e->runners[k].qid );
+ /* message( "runner %i on cpuid=%i with qid=%i." , e->runners[k].id , */
+ /* e->runners[k].cpuid , e->runners[k].qid ); */
}
/* Wait for the runner threads to be in place. */
diff --git a/src/engine.h b/src/engine.h
index c8b9d7a46f8c788237e3b9ba7cb9d7d0a6e1370d..82e76481f8d1b17c6f23ee70b15b2013bb0fce64 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -1,6 +1,11 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * Angus Lepper (angus.lepper@ed.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -40,6 +45,8 @@
#include "task.h"
#include "parser.h"
#include "partition.h"
+#include "physical_constants.h"
+#include "potentials.h"
/* Some constants. */
enum engine_policy {
@@ -164,13 +171,21 @@ struct engine {
/* Are we talkative ? */
int verbose;
+
+ /* Physical constants definition */
+ const struct phys_const *physical_constants;
+
+ /* Properties of external gravitational potential */
+ const struct external_potential *external_potential;
};
/* Function prototypes. */
void engine_barrier(struct engine *e, int tid);
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
- int nr_threads, int policy, int verbose);
+ int nr_threads, int policy, int verbose,
+ const struct phys_const *physical_constants,
+ const struct external_potential *potential);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
void engine_prepare(struct engine *e);
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 92a9f64c1f84a9e949f4c0e9485f892b5c808cdc..4dd8d5b7c7194ae072f90e7703907500c208b9b9 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2016 Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -18,19 +19,50 @@
******************************************************************************/
#include <float.h>
+#include "potentials.h"
/**
- * @brief Computes the gravity time-step of a given particle
+ * @brief Computes the gravity time-step of a given particle due to an external
+ *potential.
*
- * @param gp Pointer to the g-particle data
+ * This function only branches towards the potential chosen by the user.
*
+ * @param potential The properties of the external potential.
+ * @param phys_const The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
*/
+__attribute__((always_inline))
+ INLINE static float gravity_compute_timestep_external(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const gp) {
+
+ float dt = FLT_MAX;
+
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+ dt =
+ fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
+#endif
+
+ return dt;
+}
+/**
+ * @brief Computes the gravity time-step of a given particle due to self-gravity
+ *
+ * This function only branches towards the potential chosen by the user.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
+ */
__attribute__((always_inline))
- INLINE static float gravity_compute_timestep(struct gpart* gp) {
+ INLINE static float gravity_compute_timestep_self(
+ const struct phys_const* const phys_const,
+ const struct gpart* const gp) {
- /* Currently no limit is imposed */
- return FLT_MAX;
+ float dt = FLT_MAX;
+
+ return dt;
}
/**
@@ -71,6 +103,24 @@ __attribute__((always_inline))
__attribute__((always_inline))
INLINE static void gravity_end_force(struct gpart* gp) {}
+/**
+ * @brief Computes the gravitational acceleration induced by external potentials
+ *
+ * This function only branches towards the potential chosen by the user.
+ *
+ * @param potential The properties of the external potential.
+ * @param phys_const The physical constants in internal units.
+ * @param gp The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* gp) {
+
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+ external_gravity_pointmass(potential, phys_const, gp);
+#endif
+}
+
/**
* @brief Kick the additional variables
*
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 654745bfeb70dddba772af9e23797713376377a7..72ae5fd55f399448d78c6f2d1533c7f007f16dd1 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -24,5 +24,6 @@ __attribute__((always_inline))
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"mass=%.3e t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
- p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin, p->ti_end);
+ p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin,
+ p->ti_end);
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index ec7e1e7e9750d957547e04e95410017fd5b0d453..a54fc1d18e7e84a1ca7a7e7913117ed207549271 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -43,6 +43,10 @@ struct gpart {
/* Particle time of end of time-step. */
int ti_end;
+ /* /\* current time of x, and of v_full *\/ */
+ /* float tx; */
+ /* float tv; */
+
/* Anonymous union for id/part. */
union {
diff --git a/src/partition.c b/src/partition.c
index c4145422de72cd03d5e76c3e5b2a9ac217420d47..70249a679edfe72bda97eddfc918a8526659d995 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -1044,7 +1044,6 @@ static int check_complete(struct space *s, int verbose, int nregions) {
return (!failed);
}
-
/**
* @brief Partition a space of cells based on another space of cells.
*
@@ -1062,7 +1061,8 @@ static int check_complete(struct space *s, int verbose, int nregions) {
*
* @param oldh the cell dimensions of old space.
* @param oldcdim number of cells per dimension in old space.
- * @param oldnodeIDs the nodeIDs of cells in the old space, indexed by old cellid.
+ * @param oldnodeIDs the nodeIDs of cells in the old space, indexed by old
+ *cellid.
* @param s the space to be partitioned.
*
* @return 1 if the new space contains nodeIDs from all nodes, 0 otherwise.
diff --git a/src/physical_constants.c b/src/physical_constants.c
new file mode 100644
index 0000000000000000000000000000000000000000..d00d63df1c1e6a4821ce4ba50dfef9c4e0def9d2
--- /dev/null
+++ b/src/physical_constants.c
@@ -0,0 +1,119 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "physical_constants.h"
+
+/* Local headers. */
+#include "error.h"
+#include "physical_constants_cgs.h"
+
+/**
+ * @brief Converts physical constants to the internal unit system
+ *
+ * @param us The current internal system of units.
+ * @param internal_const The physical constants to initialize.
+ */
+void phys_const_init(struct UnitSystem* us, struct phys_const* internal_const) {
+
+ /* Units are declared as {U_M, U_L, U_t, U_I, U_T} */
+
+ const float dimension_G[5] = {-1, 3, -2, 0, 0};
+ internal_const->const_newton_G =
+ const_newton_G_cgs / units_general_conversion_factor(us, dimension_G);
+
+ const float dimension_c[5] = {0, 1, -1, 0, 0};
+ internal_const->const_speed_light_c =
+ const_speed_light_c_cgs /
+ units_general_conversion_factor(us, dimension_c);
+
+ const float dimension_h[5] = {1, -2, -1, 0, 0};
+ internal_const->const_planck_h =
+ const_planck_h_cgs / units_general_conversion_factor(us, dimension_h);
+ internal_const->const_planck_hbar =
+ const_planck_hbar_cgs / units_general_conversion_factor(us, dimension_h);
+
+ const float dimension_k[5] = {1, 2, -2, 0, -1};
+ internal_const->const_boltzmann_k =
+ const_boltzmann_k_cgs / units_general_conversion_factor(us, dimension_k);
+
+ const float dimension_thomson[5] = {0, 2, 0, 0, 0};
+ internal_const->const_thomson_cross_section =
+ const_thomson_cross_section_cgs /
+ units_general_conversion_factor(us, dimension_thomson);
+
+ const float dimension_ev[5] = {1, 2, -2, 0, 0};
+ internal_const->const_electron_volt =
+ const_electron_volt_cgs /
+ units_general_conversion_factor(us, dimension_ev);
+
+ const float dimension_charge[5] = {0, 0, -1, 1, 0};
+ internal_const->const_electron_charge =
+ const_electron_charge_cgs /
+ units_general_conversion_factor(us, dimension_charge);
+
+ const float dimension_mass[5] = {1, 0, 0, 0, 0};
+ internal_const->const_electron_mass =
+ const_electron_mass_cgs /
+ units_general_conversion_factor(us, dimension_mass);
+ internal_const->const_proton_mass =
+ const_proton_mass_cgs /
+ units_general_conversion_factor(us, dimension_mass);
+ internal_const->const_solar_mass =
+ const_solar_mass_cgs /
+ units_general_conversion_factor(us, dimension_mass);
+ internal_const->const_earth_mass =
+ const_earth_mass_cgs /
+ units_general_conversion_factor(us, dimension_mass);
+
+ const float dimension_time[5] = {0, 0, 1, 0, 0};
+ internal_const->const_year =
+ const_year_cgs / units_general_conversion_factor(us, dimension_time);
+
+ const float dimension_length[5] = {0, 1, 0, 0, 0};
+ internal_const->const_astronomical_unit =
+ const_astronomical_unit_cgs /
+ units_general_conversion_factor(us, dimension_length);
+ internal_const->const_parsec =
+ const_parsec_cgs / units_general_conversion_factor(us, dimension_length);
+ internal_const->const_light_year =
+ const_light_year_cgs /
+ units_general_conversion_factor(us, dimension_length);
+}
+
+void phys_const_print(struct phys_const* internal_const) {
+
+ message("%25s = %e", "Gravitational constant",
+ internal_const->const_newton_G);
+ message("%25s = %e", "Speed of light", internal_const->const_speed_light_c);
+ message("%25s = %e", "Planck constant", internal_const->const_planck_h);
+ message("%25s = %e", "Boltzmann constant", internal_const->const_boltzmann_k);
+ message("%25s = %e", "Thomson cross-section",
+ internal_const->const_thomson_cross_section);
+ message("%25s = %e", "Electron-Volt", internal_const->const_electron_volt);
+ message("%25s = %e", "Year", internal_const->const_year);
+ message("%25s = %e", "Astronomical Unit",
+ internal_const->const_astronomical_unit);
+ message("%25s = %e", "Parsec", internal_const->const_parsec);
+ message("%25s = %e", "Solar mass", internal_const->const_solar_mass);
+}
diff --git a/src/physical_constants.h b/src/physical_constants.h
new file mode 100644
index 0000000000000000000000000000000000000000..aac0fa0840e08a1140bbab67e944f7f134f4222b
--- /dev/null
+++ b/src/physical_constants.h
@@ -0,0 +1,85 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PHYSICAL_CONSTANTS_H
+#define SWIFT_PHYSICAL_CONSTANTS_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "units.h"
+
+/* physical constants in in defined programme units */
+struct phys_const {
+
+ /* Newton's gravitationl constant */
+ double const_newton_G;
+
+ /* Speed of light in vacuum */
+ double const_speed_light_c;
+
+ /* Planck's constant */
+ double const_planck_h;
+
+ /* Planck's reduced constant */
+ double const_planck_hbar;
+
+ /* Boltzmann's constant */
+ double const_boltzmann_k;
+
+ /* Thomson cross-section */
+ double const_thomson_cross_section;
+
+ /* Charge of the electron */
+ double const_electron_charge;
+
+ /* Electron-Volt */
+ double const_electron_volt;
+
+ /* Mass of the electron */
+ double const_electron_mass;
+
+ /* Mass of the proton */
+ double const_proton_mass;
+
+ /* (Tropical) Year */
+ double const_year;
+
+ /* Astronomical unit */
+ double const_astronomical_unit;
+
+ /* Parsec */
+ double const_parsec;
+
+ /* Light-year */
+ double const_light_year;
+
+ /* Mass of the Sun */
+ double const_solar_mass;
+
+ /* Mass of the Earth */
+ double const_earth_mass;
+};
+
+void phys_const_init(struct UnitSystem* us, struct phys_const* internal_const);
+
+void phys_const_print(struct phys_const* internal_const);
+
+#endif /* SWIFT_PHYSICAL_CONSTANTS_H */
diff --git a/src/physical_constants_cgs.h b/src/physical_constants_cgs.h
new file mode 100644
index 0000000000000000000000000000000000000000..9bd435bde9f0f80c06f0de5fee980f080c2f57d3
--- /dev/null
+++ b/src/physical_constants_cgs.h
@@ -0,0 +1,81 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PHYSICAL_CONSTANTS_CGS_H
+#define SWIFT_PHYSICAL_CONSTANTS_CGS_H
+
+/* The constants declared in this file should _NOT_ be used directly */
+/* Users should use the converted values in the phys_const structure */
+/* where all the constants are defined in the system of units specified */
+/* in the parameter file. */
+
+/* All values are taken from K.A. Olive et al. (Particle Data Group), Chin. */
+/* Phys. C, 38, 090001 (2014) and 2015 update. */
+/* http://pdg.lbl.gov/2015/reviews/rpp2015-rev-phys-constants.pdf */
+/* http://pdg.lbl.gov/2015/reviews/rpp2015-rev-astrophysical-constants.pdf */
+
+/* Newton's gravitation constant */
+const double const_newton_G_cgs = 6.67408e-8; /* g^-1 cm^3 s^-2 */
+
+/* Speed of light in vacuum */
+const double const_speed_light_c_cgs = 2.99792458e10; /* cm s^-1 */
+
+/* Planck's constant */
+const double const_planck_h_cgs = 6.626070040e-27; /* g cm^-2 s^-1 */
+
+/* Planck's reduced constant */
+const double const_planck_hbar_cgs = 1.054571800e-27; /* g cm^-2 s^-1 */
+
+/* Boltzmann's constant */
+const double const_boltzmann_k_cgs = 1.38064852e-16; /* g cm^2 s^-2 K^-1 */
+
+/* Thomson cross-section */
+const double const_thomson_cross_section_cgs = 6.6524587158e-25; /* cm^2 */
+
+/* Elementary charge */
+const double const_electron_charge_cgs = 1.6021766208e-19; /* A s^-1 */
+
+/* Electron-Volt */
+const double const_electron_volt_cgs = 1.6021766208e-12; /* g cm^2 s^-2 */
+
+/* Mass of the electron */
+const double const_electron_mass_cgs = 9.10938356e-28; /* g */
+
+/* Mass of the proton */
+const double const_proton_mass_cgs = 1.672621898e-24; /* g */
+
+/* Tropical year */
+const double const_year_cgs = 3.15569252e7; /* s */
+
+/* Astronomical unit */
+const double const_astronomical_unit_cgs = 1.49597870700e13; /* cm */
+
+/* Parsec */
+const double const_parsec_cgs = 3.08567758149e18; /* cm */
+
+/* Light-year */
+const double const_light_year_cgs = 9.46053e17; /* cm */
+
+/* Mass of the Sun */
+const double const_solar_mass_cgs = 1.9885e33; /* g */
+
+/* Mass of the Earth */
+const double const_earth_mass_cgs = 5.9726e27; /* g */
+
+#endif /* SWIFT_PHYSICAL_CONSTANTS_CGS_H */
diff --git a/src/potentials.c b/src/potentials.c
new file mode 100644
index 0000000000000000000000000000000000000000..98e57d7959963e07cfbbd2fdb60a0df504e6a487
--- /dev/null
+++ b/src/potentials.c
@@ -0,0 +1,65 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "potentials.h"
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param us The current internal system of units
+ * @param potential The external potential properties to initialize
+ */
+void potential_init(const struct swift_params* parameter_file,
+ struct UnitSystem* us,
+ struct external_potential* potential) {
+
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+
+ potential->point_mass.x =
+ parser_get_param_double(parameter_file, "PointMass:position_x");
+ potential->point_mass.y =
+ parser_get_param_double(parameter_file, "PointMass:position_y");
+ potential->point_mass.z =
+ parser_get_param_double(parameter_file, "PointMass:position_z");
+ potential->point_mass.mass =
+ parser_get_param_double(parameter_file, "PointMass:mass");
+
+#endif /* EXTERNAL_POTENTIAL_POINTMASS */
+}
+
+/**
+ * @brief Prints the properties of the external potential to stdout.
+ *
+ * @param potential The external potential properties.
+ */
+void potential_print(const struct external_potential* potential) {
+
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+ message("Point mass properties are (x,y,z) = (%e, %e, %e), M = %e",
+ potential->point_mass.x, potential->point_mass.y,
+ potential->point_mass.z, potential->point_mass.mass);
+#endif /* EXTERNAL_POTENTIAL_POINTMASS */
+}
diff --git a/src/potentials.h b/src/potentials.h
new file mode 100644
index 0000000000000000000000000000000000000000..9c687fa4c4a64bcc0ab8cac1a5e495487fa938f1
--- /dev/null
+++ b/src/potentials.h
@@ -0,0 +1,122 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_POTENTIALS_H
+#define SWIFT_POTENTIALS_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/* External Potential Properties */
+struct external_potential {
+
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+ struct {
+ double x, y, z;
+ double mass;
+ } point_mass;
+#endif
+};
+
+/* Include exteral pointmass potential */
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+
+#define EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR 0.03f
+
+/**
+ * @brief Computes the time-step due to the acceleration from a point mass
+ *
+ * @param potential The properties of the externa potential.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline))
+ INLINE static float external_gravity_pointmass_timestep(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const g) {
+
+ const float G_newton = phys_const->const_newton_G;
+ const float dx = g->x[0] - potential->point_mass.x;
+ const float dy = g->x[1] - potential->point_mass.y;
+ const float dz = g->x[2] - potential->point_mass.z;
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
+ const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
+ (g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
+ (g->x[2] - potential->point_mass.z) * (g->v_full[2]);
+ const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
+ const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
+ const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
+ const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
+ const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
+ g->a_grav[2] * g->a_grav[2];
+
+ return EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR * sqrtf(a_2 / dota_2);
+}
+
+/**
+ * @brief Computes the gravitational acceleration of a particle due to a point
+ * mass
+ *
+ * @param potential The proerties of the external potential.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_pointmass(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* g) {
+
+ const float G_newton = phys_const->const_newton_G;
+ const float dx = g->x[0] - potential->point_mass.x;
+ const float dy = g->x[1] - potential->point_mass.y;
+ const float dz = g->x[2] - potential->point_mass.z;
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
+
+ g->a_grav[0] +=
+ -G_newton * potential->point_mass.mass * dx * rinv * rinv * rinv;
+ g->a_grav[1] +=
+ -G_newton * potential->point_mass.mass * dy * rinv * rinv * rinv;
+ g->a_grav[2] +=
+ -G_newton * potential->point_mass.mass * dz * rinv * rinv * rinv;
+}
+#endif /* EXTERNAL_POTENTIAL_POINTMASS */
+
+/* Now, some generic functions, defined in the source file */
+
+void potential_init(const struct swift_params* parameter_file,
+ struct UnitSystem* us,
+ struct external_potential* potential);
+
+void potential_print(const struct external_potential* potential);
+
+#endif /* SWIFT_POTENTIALS_H */
diff --git a/src/runner.c b/src/runner.c
index 25ffc9f2c2b9fff4596d680ef505e714c7595c07..2fd272039d04f2c42db3d39b2acdbe9009b19125 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1,6 +1,10 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -68,6 +72,21 @@ const float runner_shift[13 * 3] = {
const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+/* #define MY_CELL 428428428 */
+/* #define DX 0.1 */
+/* #define NX 1000 */
+/* #define CELL_ID \ */
+/* ((int)(c->loc[0] / DX) * NX *NX + (int)(c->loc[1] / DX) * NX + \ */
+/* (int)(c->loc[2] / DX)) */
+/* #define OUT \ */
+/* if (CELL_ID == MY_CELL) { \ */
+/* message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount,
+ * r->e->time); \ */
+/* } */
+/* #define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time); */
+/* #define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d */
+/* %f\n",c->loc[0],c->loc[1],c->loc[2], CELL_ID, c->count, r->e->time); */
+
/* Import the density loop functions. */
#define FUNCTION density
#include "runner_doiact.h"
@@ -80,6 +99,95 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
/* Import the gravity loop functions. */
#include "runner_doiact_grav.h"
+/**
+ * @brief Calculate gravity acceleration from external potential
+ *
+ * @param r runner task
+ * @param c cell
+ * @param timer 1 if the time is to be recorded.
+ */
+void runner_dograv_external(struct runner *r, struct cell *c, int timer) {
+
+ struct gpart *g, *gparts = c->gparts;
+ int i, k, gcount = c->gcount;
+ const int ti_current = r->e->ti_current;
+ const struct external_potential *potential = r->e->external_potential;
+ const struct phys_const *constants = r->e->physical_constants;
+
+ TIMER_TIC;
+
+ /* Recurse? */
+ if (c->split) {
+ for (k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) runner_dograv_external(r, c->progeny[k], 0);
+ return;
+ }
+
+#ifdef TASK_VERBOSE
+ OUT;
+#endif
+
+ /* Loop over the parts in this cell. */
+ for (i = 0; i < gcount; i++) {
+
+ /* Get a direct pointer on the part. */
+ g = &gparts[i];
+
+ /* Is this part within the time step? */
+ if (g->ti_end <= ti_current) {
+
+ external_gravity(potential, constants, g);
+
+ /* /\* check for energy and angular momentum conservation - begin by */
+ /* * synchronizing velocity*\/ */
+ /* const float dx = g->x[0] - r->e->potential->point_mass.x; */
+ /* const float dy = g->x[1] - r->e->potential->point_mass.y; */
+ /* const float dz = g->x[2] - r->e->potential->point_mass.z; */
+ /* const float dr = sqrtf((dx * dx) + (dy * dy) + (dz * dz)); */
+ /* const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz); */
+
+ /* const int current_dti = g->ti_end - g->ti_begin; */
+ /* const float dt = 0.5f * current_dti * r->e->timeBase; */
+ /* const float vx = g->v_full[0] + dt * g->a_grav[0]; */
+ /* const float vy = g->v_full[1] + dt * g->a_grav[1]; */
+ /* const float vz = g->v_full[2] + dt * g->a_grav[2]; */
+
+ /* /\* E/L *\/ */
+ /* float L[3], E; */
+ /* E = 0.5 * ((vx * vx) + (vy * vy) + (vz * vz)) - */
+ /* r->e->physical_constants->newton_gravity * */
+ /* r->e->potential->point_mass.mass * rinv; */
+ /* L[0] = dy * vz - dz * vy; */
+ /* L[1] = dz * vx - dx * vz; */
+ /* L[2] = dx * vy - dy * vx; */
+ /* if (abs(g->id) == 1) { */
+ /* float v2 = vx * vx + vy * vy + vz * vz; */
+ /* float fg = r->e->physical_constants->newton_gravity * */
+ /* r->e->potential->point_mass.mass * rinv; */
+ /* float fga = sqrtf((g->a_grav[0] * g->a_grav[0]) + */
+ /* (g->a_grav[1] * g->a_grav[1]) + */
+ /* (g->a_grav[2] * g->a_grav[2])) * */
+ /* dr; */
+ /* // message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]=
+ * %f */
+ /* // x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2],
+ * g->x[0], */
+ /* // g->x[1], vx, vy); */
+ /* message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time,
+ */
+ /* g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1],
+ */
+ /* vx, vy); */
+ /* /\* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity,
+ */
+ /* * r->e->potential->point_mass.mass); *\/ */
+ /* /\* exit(-1); *\/ */
+ /* } */
+ }
+ }
+ if (timer) TIMER_TOC(timer_dograv_external);
+}
+
/**
* @brief Sort the entries in ascending order using QuickSort.
*
@@ -466,7 +574,6 @@ void runner_dogsort(struct runner *r, struct cell *c, int flags, int clock) {
* @param c The cell.
* @param timer 1 if the time is to be recorded.
*/
-
void runner_doinit(struct runner *r, struct cell *c, int timer) {
struct part *const parts = c->parts;
@@ -682,6 +789,9 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
TIMER_TIC
+#ifdef TASK_VERBOSE
+ OUT;
+#endif
/* No children? */
if (!c->split) {
@@ -808,6 +918,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
+ const struct external_potential *potential = r->e->external_potential;
+ const struct phys_const *constants = r->e->physical_constants;
const int is_fixdt =
(r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
@@ -819,6 +931,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
TIMER_TIC
+#ifdef TASK_VERBOSE
+ OUT;
+#endif
+
/* No children? */
if (!c->split) {
@@ -829,70 +945,78 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
- if (gp->id < 0 && (is_fixdt || gp->ti_end <= ti_current)) {
+ if (gp->id < 0) {
- /* First, finish the force calculation */
- gravity_end_force(gp);
+ if (is_fixdt || gp->ti_end <= ti_current) {
- /* Now we are ready to compute the next time-step size */
- int new_dti;
+ /* First, finish the force calculation */
+ gravity_end_force(gp);
- if (is_fixdt) {
+ /* Now we are ready to compute the next time-step size */
+ int new_dti;
- /* Now we have a time step, proceed with the kick */
- new_dti = global_dt_max * timeBase_inv;
+ if (is_fixdt) {
- } else {
+ /* Now we have a time step, proceed with the kick */
+ new_dti = global_dt_max * timeBase_inv;
- /* Compute the next timestep (gravity condition) */
- float new_dt = gravity_compute_timestep(gp);
+ } else {
- /* Limit timestep within the allowed range */
- new_dt = fminf(new_dt, global_dt_max);
- new_dt = fmaxf(new_dt, global_dt_min);
+ /* Compute the next timestep (gravity condition) */
+ const float new_dt_external =
+ gravity_compute_timestep_external(potential, constants, gp);
+ const float new_dt_self =
+ gravity_compute_timestep_self(constants, gp);
- /* Convert to integer time */
- new_dti = new_dt * timeBase_inv;
+ float new_dt = fminf(new_dt_external, new_dt_self);
- /* Recover the current timestep */
- const int current_dti = gp->ti_end - gp->ti_begin;
+ /* Limit timestep within the allowed range */
+ new_dt = fminf(new_dt, global_dt_max);
+ new_dt = fmaxf(new_dt, global_dt_min);
- /* Limit timestep increase */
- if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+ /* Convert to integer time */
+ new_dti = new_dt * timeBase_inv;
- /* Put this timestep on the time line */
- int dti_timeline = max_nr_timesteps;
- while (new_dti < dti_timeline) dti_timeline /= 2;
+ /* Recover the current timestep */
+ const int current_dti = gp->ti_end - gp->ti_begin;
- /* Now we have a time step, proceed with the kick */
- new_dti = dti_timeline;
- }
+ /* Limit timestep increase */
+ if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
- /* Compute the time step for this kick */
- const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
- const int ti_end = gp->ti_end + new_dti / 2;
- const double dt = (ti_end - ti_start) * timeBase;
- const double half_dt = (ti_end - gp->ti_end) * timeBase;
+ /* Put this timestep on the time line */
+ int dti_timeline = max_nr_timesteps;
+ while (new_dti < dti_timeline) dti_timeline /= 2;
- /* Move particle forward in time */
- gp->ti_begin = gp->ti_end;
- gp->ti_end = gp->ti_begin + new_dti;
+ /* Now we have a time step, proceed with the kick */
+ new_dti = dti_timeline;
+ }
- /* Kick particles in momentum space */
- gp->v_full[0] += gp->a_grav[0] * dt;
- gp->v_full[1] += gp->a_grav[1] * dt;
- gp->v_full[2] += gp->a_grav[2] * dt;
+ /* Compute the time step for this kick */
+ const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
+ const int ti_end = gp->ti_end + new_dti / 2;
+ const double dt = (ti_end - ti_start) * timeBase;
+ const double half_dt = (ti_end - gp->ti_end) * timeBase;
- /* Extra kick work */
- gravity_kick_extra(gp, dt, half_dt);
+ /* Move particle forward in time */
+ gp->ti_begin = gp->ti_end;
+ gp->ti_end = gp->ti_begin + new_dti;
- /* Number of updated g-particles */
- g_updated++;
- }
+ /* Kick particles in momentum space */
+ gp->v_full[0] += gp->a_grav[0] * dt;
+ gp->v_full[1] += gp->a_grav[1] * dt;
+ gp->v_full[2] += gp->a_grav[2] * dt;
- /* Minimal time for next end of time-step */
- ti_end_min = min(gp->ti_end, ti_end_min);
- ti_end_max = max(gp->ti_end, ti_end_max);
+ /* Extra kick work */
+ gravity_kick_extra(gp, dt, half_dt);
+
+ /* Number of updated g-particles */
+ g_updated++;
+ }
+
+ /* Minimal time for next end of time-step */
+ ti_end_min = min(gp->ti_end, ti_end_min);
+ ti_end_max = max(gp->ti_end, ti_end_max);
+ }
}
/* Now do the hydro ones... */
@@ -929,9 +1053,17 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Compute the next timestep (gravity condition) */
float new_dt_grav = FLT_MAX;
- if (p->gpart != NULL)
- new_dt_grav = gravity_compute_timestep(p->gpart);
+ if (p->gpart != NULL) {
+ const float new_dt_external = gravity_compute_timestep_external(
+ potential, constants, p->gpart);
+ const float new_dt_self =
+ gravity_compute_timestep_self(constants, p->gpart);
+
+ new_dt_grav = fminf(new_dt_external, new_dt_self);
+ }
+
+ /* Final time-step is minimum of hydro and gravity */
float new_dt = fminf(new_dt_hydro, new_dt_grav);
/* Limit change in h */
@@ -1038,7 +1170,6 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
ti_end_min = min(p->ti_end, ti_end_min);
ti_end_max = max(p->ti_end, ti_end_max);
}
-
}
/* Otherwise, aggregate data from children. */
@@ -1199,6 +1330,9 @@ void *runner_main(void *data) {
case task_type_grav_down:
runner_dograv_down(r, t->ci);
break;
+ case task_type_grav_external:
+ runner_dograv_external(r, t->ci, 1);
+ break;
case task_type_part_sort:
space_do_parts_sort();
break;
diff --git a/src/runner.h b/src/runner.h
index ff2f93db6eae9d9e85cd9c3de6f398ed0f64c681..f052bb57cdfc92a632cc4143d74250e7b6f6e847 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -1,6 +1,10 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
diff --git a/src/space.c b/src/space.c
index ee2b8013cdfa7a885e9e9d61bf78e20aaca31366..6d4c256d502e121448238cadbe179d0000f00e88 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1,22 +1,25 @@
/*******************************************************************************
-* This file is part of SWIFT.
-* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
-* 2016 Peter W. Draper (p.w.draper@durham.ac.uk)
-*
-* This program is free software: you can redistribute it and/or modify
-* it under the terms of the GNU Lesser General Public License as published
-* by the Free Software Foundation, either version 3 of the License, or
-* (at your option) any later version.
-*
-* This program is distributed in the hope that it will be useful,
-* but WITHOUT ANY WARRANTY; without even the implied warranty of
-* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-* GNU General Public License for more details.
-*
-* You should have received a copy of the GNU Lesser General Public License
-* along with this program. If not, see <http://www.gnu.org/licenses/>.
-*
-******************************************************************************/
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
/* Config parameters. */
#include "../config.h"
@@ -164,15 +167,17 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
/* Run through the parts and get the current h_max. */
// tic = getticks();
float h_max = s->cell_min / kernel_gamma / space_stretch;
- if (s->cells != NULL) {
- for (int k = 0; k < s->nr_cells; k++) {
- if (s->cells[k].h_max > h_max) h_max = s->cells[k].h_max;
- }
- } else {
- for (size_t k = 0; k < nr_parts; k++) {
- if (s->parts[k].h > h_max) h_max = s->parts[k].h;
+ if (nr_parts > 0) {
+ if (s->cells != NULL) {
+ for (int k = 0; k < s->nr_cells; k++) {
+ if (s->cells[k].h_max > h_max) h_max = s->cells[k].h_max;
+ }
+ } else {
+ for (size_t k = 0; k < nr_parts; k++) {
+ if (s->parts[k].h > h_max) h_max = s->parts[k].h;
+ }
+ s->h_max = h_max;
}
- s->h_max = h_max;
}
/* If we are running in parallel, make sure everybody agrees on
@@ -372,7 +377,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
const ticks tic = getticks();
/* Be verbose about this. */
- // message( "re)building space..." ); fflush(stdout);
+ // message("re)building space..."); fflush(stdout);
/* Re-grid if necessary, or just re-set the cell data. */
space_regrid(s, cell_max, verbose);
@@ -1280,9 +1285,11 @@ void space_do_split(struct space *s, struct cell *c) {
if (s->nr_parts > 0)
c->owner =
((c->parts - s->parts) % s->nr_parts) * s->nr_queues / s->nr_parts;
- else
+ else if (s->nr_gparts > 0)
c->owner =
((c->gparts - s->gparts) % s->nr_gparts) * s->nr_queues / s->nr_gparts;
+ else
+ c->owner = 0; /* Ok, there is really nothing on this rank... */
}
/**
@@ -1475,10 +1482,12 @@ void space_init(struct space *s, const struct swift_params *params,
}
/* Allocate the extra parts array. */
- if (posix_memalign((void *)&s->xparts, xpart_align,
- Npart * sizeof(struct xpart)) != 0)
- error("Failed to allocate xparts.");
- bzero(s->xparts, Npart * sizeof(struct xpart));
+ if (Npart > 0) {
+ if (posix_memalign((void *)&s->xparts, xpart_align,
+ Npart * sizeof(struct xpart)) != 0)
+ error("Failed to allocate xparts.");
+ bzero(s->xparts, Npart * sizeof(struct xpart));
+ }
/* Init the space lock. */
if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
diff --git a/src/space.h b/src/space.h
index 88e2f6f52774651217c4ff24e25f549d8ae1e347..ad7b9e77d484c997527d32c934bf1b85d7d54e41 100644
--- a/src/space.h
+++ b/src/space.h
@@ -1,6 +1,10 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
diff --git a/src/swift.h b/src/swift.h
index e568a28c888295affc9ec45b6d059d34f5b4bf04..4b570447345465684fc51bcd4ca96d07a4b7f518 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -40,6 +40,8 @@
#include "parser.h"
#include "part.h"
#include "partition.h"
+#include "physical_constants.h"
+#include "potentials.h"
#include "queue.h"
#include "runner.h"
#include "scheduler.h"
diff --git a/src/task.c b/src/task.c
index 5f1475a46e4626e1f51db673d73fd84f86e6edb6..31c2bfcb3462b9f91547ebd631f645e190c07143 100644
--- a/src/task.c
+++ b/src/task.c
@@ -1,6 +1,10 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -43,10 +47,10 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
- "none", "sort", "self", "pair", "sub",
- "init", "ghost", "drift", "kick", "send",
- "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
- "part_sort", "gpart_sort", "split_cell", "rewait"};
+ "none", "sort", "self", "pair", "sub",
+ "init", "ghost", "drift", "kick", "send",
+ "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
+ "grav_external", "part_sort", "gpart_sort", "split_cell", "rewait"};
const char *subtaskID_names[task_type_count] = {"none", "density",
"force", "grav"};
diff --git a/src/task.h b/src/task.h
index 9c0ba6087d772d7362a98bc40a838c6fa3713166..25cc886f4b38456a0431fb6c7d0b7b1864053dd9 100644
--- a/src/task.h
+++ b/src/task.h
@@ -2,6 +2,9 @@
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -45,6 +48,7 @@ enum task_types {
task_type_grav_mm,
task_type_grav_up,
task_type_grav_down,
+ task_type_grav_external,
task_type_part_sort,
task_type_gpart_sort,
task_type_split_cell,
diff --git a/src/timers.c b/src/timers.c
index 2501d347c8cea608650ece4c2883dab85ceee058..b621d27c90902f06c3760cbef6a88237a2b3b95b 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -1,6 +1,9 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
diff --git a/src/timers.h b/src/timers.h
index de2d9edb9ed54717472e5ae1222dfb33235c3e95..bf288d73929bd664992ef3da6d531d77aa85ac3b 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -1,6 +1,9 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ * Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -37,6 +40,7 @@ enum {
timer_dopair_density,
timer_dopair_force,
timer_dopair_grav,
+ timer_dograv_external,
timer_dosub_density,
timer_dosub_force,
timer_dosub_grav,
diff --git a/src/units.c b/src/units.c
index 184dbe8a0df000008dba1d7003558d83b1f08cad..c0276fd4fa74136f46f2f4b7d719c124c92fa48c 100644
--- a/src/units.c
+++ b/src/units.c
@@ -39,7 +39,6 @@
/* Includes. */
#include "const.h"
#include "error.h"
-#include "units.h"
/**
* @brief Initialises the UnitSystem structure with the constants given in
@@ -331,7 +330,7 @@ void units_conversion_string(char* buffer, const struct UnitSystem* us,
* the desired quantity. See conversionFactor() for a working example
*/
double units_general_conversion_factor(const struct UnitSystem* us,
- float baseUnitsExponants[5]) {
+ const float baseUnitsExponants[5]) {
double factor = 1.;
int i;
@@ -349,7 +348,7 @@ double units_general_conversion_factor(const struct UnitSystem* us,
* the desired quantity. See conversionFactor() for a working example
*/
float units_general_h_factor(const struct UnitSystem* us,
- float baseUnitsExponants[5]) {
+ const float baseUnitsExponants[5]) {
float factor_exp = 0.f;
factor_exp += -baseUnitsExponants[UNIT_MASS];
@@ -367,7 +366,7 @@ float units_general_h_factor(const struct UnitSystem* us,
* the desired quantity. See conversionFactor() for a working example
*/
float units_general_a_factor(const struct UnitSystem* us,
- float baseUnitsExponants[5]) {
+ const float baseUnitsExponants[5]) {
float factor_exp = 0.f;
factor_exp += baseUnitsExponants[UNIT_LENGTH];
@@ -385,7 +384,7 @@ float units_general_a_factor(const struct UnitSystem* us,
* the desired quantity. See conversionFactor() for a working example
*/
void units_general_conversion_string(char* buffer, const struct UnitSystem* us,
- float baseUnitsExponants[5]) {
+ const float baseUnitsExponants[5]) {
char temp[14];
double a_exp = units_general_a_factor(us, baseUnitsExponants);
double h_exp = units_general_h_factor(us, baseUnitsExponants);
diff --git a/src/units.h b/src/units.h
index 3e349dc16787cd4052a3e9205b21dce3c3732448..197d3e1f1dac2e01015d44758dd079ac21c6a0b7 100644
--- a/src/units.h
+++ b/src/units.h
@@ -97,19 +97,19 @@ double units_get_base_unit(const struct UnitSystem*, enum BaseUnits);
const char* units_get_base_unit_symbol(enum BaseUnits);
const char* units_get_base_unit_CGS_symbol(enum BaseUnits);
double units_general_conversion_factor(const struct UnitSystem* us,
- float baseUnitsExponants[5]);
+ const float baseUnitsExponants[5]);
double units_conversion_factor(const struct UnitSystem* us,
enum UnitConversionFactor unit);
float units_general_h_factor(const struct UnitSystem* us,
- float baseUnitsExponants[5]);
+ const float baseUnitsExponants[5]);
float units_h_factor(const struct UnitSystem* us,
enum UnitConversionFactor unit);
float units_general_a_factor(const struct UnitSystem* us,
- float baseUnitsExponants[5]);
+ const float baseUnitsExponants[5]);
float units_a_factor(const struct UnitSystem* us,
enum UnitConversionFactor unit);
void units_general_conversion_string(char* buffer, const struct UnitSystem* us,
- float baseUnitsExponants[5]);
+ const float baseUnitsExponants[5]);
void units_conversion_string(char* buffer, const struct UnitSystem* us,
enum UnitConversionFactor unit);