Merge branch 'master' into intrinsic-vectorisation

Conflicts:
	src/hydro/Default/hydro_iact.h
	src/hydro/Gadget2/hydro.h
	src/hydro/Gadget2/hydro_debug.h
	src/hydro/Gadget2/hydro_iact.h
	src/hydro/Gadget2/hydro_part.h
	tests/testInteractions.c
	tests/testSymInt.c

.gitignore

@@ -28,10 +28,13 @@ examples/*.xmf
 examples/used_parameters.yml
 examples/energy.txt
 examples/*/*.xmf
-examples/*/*.dat
 examples/*/*.hdf5
+examples/*/*.txt
 examples/*/used_parameters.yml
-examples/*/energy.txt
+examples/*/*/*.xmf
+examples/*/*/*.hdf5
+examples/*/*/*.txt
+examples/*/*/used_parameters.yml
 
 tests/testPair
 tests/brute_force_standard.dat

examples/ExternalPointMass/externalPointMass.yml

@@ -43,4 +43,5 @@ PointMass:
   position_y:      50.
   position_z:      50.	
   mass:            1e10     # mass of external point mass in internal units
+  timestep_mult:   0.03     # controls time step
 

examples/IsothermalPotential/GravityOnly/README

+;
+; this probelm generates a set of gravity particles in an isothermal
+; potential and follows their orbits. Tests verify consdevation of
+; energy and angular momentum
+;
+;

examples/IsothermalPotential/GravityOnly/isothermal.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  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 governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   8.    # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            Isothermal # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          0.02       # Time difference between consecutive outputs (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_smoothing_length:  40.      # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  Isothermal.hdf5       # The file to read
+  shift_x:    100.                  # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    100.
+  shift_z:    100.
+ 
+# External potential parameters
+IsothermalPotential:
+  position_x:      100.     # location of centre of isothermal potential in internal units
+  position_y:      100.
+  position_z:      100.	
+  vrot:            200.     # rotation speed of isothermal potential in internal units
+  timestep_mult:   0.03     # controls time step
+

examples/IsothermalPotential/GravityOnly/makeIC.py

+###############################################################################
+ # 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 N particles in a spherical distribution centred on [0,0,0], to be moved in an isothermal potential
+# usage: python makeIC.py 1000 0 : generate 1000 particles on circular orbits
+#        python makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
+# all particles move in the xy plane, and start at y=0
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS  = 6.672e-8
+SOLAR_MASS_IN_CGS   = 1.9885e33
+PARSEC_IN_CGS       = 3.0856776e18
+PROTON_MASS_IN_CGS  = 1.6726231e24
+YEAR_IN_CGS         = 3.154e+7
+
+# 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
+
+
+# rotation speed of isothermal potential [km/s]
+vrot_kms = 200.
+
+# 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
+vrot   = vrot_kms * 1e5 / const_unit_velocity_in_cgs
+
+
+# Parameters
+periodic= 1            # 1 For periodic box
+boxSize = 400.         #  [kpc]
+Radius  = 100.          # maximum radius of particles [kpc]
+G       = const_G 
+
+N       = int(sys.argv[1])  # Number of particles
+icirc   = int(sys.argv[2])  # if = 0, all particles are on circular orbits, if = 1, Lz/Lcirc uniform in ]0,1[
+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 = "Isothermal.hdf5" 
+
+
+#---------------------------------------------------
+numPart        = N
+mass           = 1
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs 
+grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# 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
+
+# set seed for random number
+numpy.random.seed(1234)
+
+#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  = vrot
+v      = numpy.zeros((numPart, 3))
+omega  = speed / radius
+period = 2.*math.pi/omega
+print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+omegav = omega
+if (icirc != 0):
+    omegav = omega * numpy.random.rand(N)
+
+v[:,0] = -omegav * r[:,1]
+v[:,1] =  omegav * 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, dtype='L')
+ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ds[()] = ids
+
+ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = r
+
+
+file.close()

examples/IsothermalPotential/GravityOnly/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Isothermal.hdf5 ]
+then
+    echo "Generating initial conditions for the isothermal potential box example..."
+    python makeIC.py 1000 1
+fi
+
+../../swift -g -t 2 isothermal.yml

examples/IsothermalPotential/GravityOnly/test.pro

+;
+;  test energy / angular momentum conservation of test problem
+;
+
+iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
+
+; set physical constants
+@physunits
+
+indir    = './'
+basefile = 'Isothermal_'
+
+; set properties of potential
+uL   = 1e3 * phys.pc             ; unit of length
+uM   = phys.msun                 ; unit of mass
+uV   = 1d5                       ; unit of velocity
+vrot = 200.                      ; km/s
+r200 = 100.                      ; virial radius
+
+; derived units
+constG   = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre  = [100.,100.,100.] * 1d3 * pc / uL
+
+;
+infile = indir + basefile + '*'
+spawn,'ls -1 '+infile,res
+nfiles = n_elements(res)
+
+
+
+; choose: calculate change of energy and Lz, comparing first and last
+; snapshots for all particles, or do so for a subset
+
+; compare all
+ifile   = 0
+inf     = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+id      = h5rd(inf,'PartType1/ParticleIDs')
+nfollow = n_elements(id)
+
+; follow a subset
+nfollow  = 500                    ; number of particles to follow
+
+;
+if (iplot eq 1) then begin
+   nskip = 1
+   nsave = nfiles
+endif else begin
+   nskip = nfiles - 2
+   nsave = 2
+endelse
+
+;
+lout     = fltarr(nfollow, nsave) ; Lz
+xout     = fltarr(nfollow, nsave) ; x
+yout     = fltarr(nfollow, nsave) ; y
+zout     = fltarr(nfollow, nsave) ; z
+eout     = fltarr(nfollow, nsave) ; energies
+ekin     = fltarr(nfollow, nsave)
+epot     = fltarr(nfollow, nsave)
+tout     = fltarr(nsave)
+
+
+
+ifile  = 0
+isave = 0
+for ifile=0,nfiles-1,nskip 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
+   xout[*,isave] = p[0,0:nfollow-1]-pcentre[0]
+   yout[*,isave] = p[1,0:nfollow-1]-pcentre[1]
+   zout[*,isave] = p[2,0:nfollow-1]-pcentre[2]
+   Lz  = (p[0,*]-pcentre[0]) * v[1,*] - (p[1,*]-pcentre[1]) * v[0,*]
+   dr = sqrt(dr)
+;   print,'time = ',time,p[0,0],v[0,0],id[0]
+   ek   = 0.5 * dv
+;   ep   = - constG * mextern / dr
+   ep   = -vrot*vrot * (1 + alog(r200/dr))
+   ener = ek + ep
+   tout(isave) = time
+   lout[*,isave] = lz[0:nfollow-1]
+   eout(*,isave) = ener[0:nfollow-1]
+   ekin(*,isave) = ek[0:nfollow-1]
+   epot(*,isave) = ep[0:nfollow-1]
+
+;  write some output
+;   print,' time= ',time,' e= ',eout[0],' Lz= ',lz[0],format='(%a %f %a
+;   %f)'
+   print,format='('' time= '',f7.1,'' E= '',f9.2,'' Lz= '',e9.2)', time,eout[0],lz[0]
+   isave = isave + 1
+   
+endfor
+x0 = reform(xout[0,*])
+y0 = reform(xout[1,*])
+z0 = reform(xout[2,*])
+
+; calculate relative energy change
+de    = 0.0 * eout
+dl    = 0.0 * lout
+nsave = isave
+for ifile=1, nsave-1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
+for ifile=1, nsave-1 do dl[*,ifile] = (lout[*,ifile] - lout[*,0])/lout[*,0]
+
+
+; calculate statistics of energy changes
+print,' relatve energy change: (per cent) ',minmax(de) * 100.
+print,' relative Lz    change: (per cent) ',minmax(dl) * 100.
+
+; plot enery and Lz conservation for some particles
+if(iplot eq 1) then begin
+; plot results on energy conservation for some particles
+   nplot = min(10, nfollow)
+   wset,0
+   xr = [min(tout), max(tout)]
+   yr = [-2,2]*1d-2             ; in percent
+   plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='time',ytitle='dE/E, dL/L (%)'
+   for i=0,nplot-1 do oplot,tout,de[i,*]
+   for i=0,nplot-1 do oplot,tout,dl[i,*],color=red
+   legend,['dE/E','dL/L'],linestyle=[0,0],color=[black,red],box=0,/bottom,/left
+   screen_to_png,'e-time.png'
+
+;  plot orbits of those particles
+   wset,2
+   xr = [-100,100]
+   yr = xr
+   plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/iso,/nodata,xtitle='x',ytitle='y'
+   color = floor(findgen(nplot)*255/float(nplot))
+   for i=0,nplot-1 do oplot,xout[i,*],yout[i,*],color=color(i)
+   screen_to_png,'orbit.png'
+
+; plot radial position of these particles
+   wset,4
+   xr = [min(tout), max(tout)]
+   yr = [0,80]
+   plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='r'
+   color = floor(findgen(nplot)*255/float(nplot))
+for i=0,nplot-1 do begin dr = sqrt(reform(xout[i,*])^2 + reform(yout[i,*])^2) &  oplot,tout,dr,color=color[i] & endfor
+   screen_to_png,'r-time.png'
+
+; make histogram of energy changes at end
+   wset,6
+   ohist,de,x,y,-0.05,0.05,0.001
+   plot,x,y,psym=10,xtitle='de (%)'
+   screen_to_png,'de-hist.png'
+
+
+endif
+
+end
+
+

examples/parameter_example.yml

@@ -54,7 +54,7 @@ InitialConditions:
   shift_y:    0.
   shift_z:    0.
 
-# Parameters govering domain decomposition
+# Parameters governing domain decomposition
 DomainDecomposition:
   initial_type:       m     # (Optional) The initial strategy ("g", "m", "w", or "v").
   initial_grid_x:    10     # (Optional) Grid size if the "g" strategy is chosen.
@@ -64,9 +64,17 @@ DomainDecomposition:
  
 # Parameters related to external potentials
   
-# Point mass external potential
+# Point mass external potentials
 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
+  timestep_mult:   0.03     # Pre-factor for the time-step condition
+
+IsothermalPotential:
+  position_x:      100.     # Location of centre of isothermal potential in internal units
+  position_y:      100.
+  position_z:      100.
+  vrot:            200.     # Rotation speed of isothermal potential in internal units
+  timestep_mult:   0.03     # Pre-factor for the time-step condition

src/cell.c

@@ -110,13 +110,13 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
       temp->loc[0] = c->loc[0];
       temp->loc[1] = c->loc[1];
       temp->loc[2] = c->loc[2];
-      temp->h[0] = c->h[0] / 2;
-      temp->h[1] = c->h[1] / 2;
-      temp->h[2] = c->h[2] / 2;
+      temp->width[0] = c->width[0] / 2;
+      temp->width[1] = c->width[1] / 2;
+      temp->width[2] = c->width[2] / 2;
       temp->dmin = c->dmin / 2;
-      if (k & 4) temp->loc[0] += temp->h[0];
-      if (k & 2) temp->loc[1] += temp->h[1];
-      if (k & 1) temp->loc[2] += temp->h[2];
+      if (k & 4) temp->loc[0] += temp->width[0];
+      if (k & 2) temp->loc[1] += temp->width[1];
+      if (k & 1) temp->loc[2] += temp->width[2];
       temp->depth = c->depth + 1;
       temp->split = 0;
       temp->dx_max = 0.f;
@@ -425,7 +425,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) {
   double pivot[3];
 
   /* Init the pivots. */
-  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->h[k] / 2;
+  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->width[k] / 2;
 
   /* Split along the x-axis. */
   i = 0;

src/cell.h

@@ -68,7 +68,7 @@ struct cell {
   double loc[3];
 
   /* The cell dimensions. */
-  double h[3];
+  double width[3];
 
   /* Max radii in this cell. */
   double h_max;

src/const.h

@@ -59,8 +59,9 @@
 #define GADGET2_SPH
 //#define DEFAULT_SPH
 
-/* Gravity properties */
+/* External potential properties */
 #define EXTERNAL_POTENTIAL_POINTMASS
+//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 
 /* Are we debugging ? */
 //#define SWIFT_DEBUG_CHECKS

src/drift.h

@@ -77,7 +77,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
   p->v[2] += p->a_hydro[2] * dt;
 
   /* Predict smoothing length */
-  const float w1 = p->h_dt * h_inv * dt;
+  const float w1 = p->force.h_dt * h_inv * dt;
   if (fabsf(w1) < 0.2f)
     p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
   else

src/engine.c

@@ -265,7 +265,7 @@ void engine_redistribute(struct engine *e) {
   struct cell *cells = s->cells;
   const int nr_cells = s->nr_cells;
   const int *cdim = s->cdim;
-  const double ih[3] = {s->ih[0], s->ih[1], s->ih[2]};
+  const double iwidth[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]};
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
   struct part *parts = s->parts;
   struct xpart *xparts = s->xparts;
@@ -299,8 +299,9 @@ void engine_redistribute(struct engine *e) {
       else if (parts[k].x[j] >= dim[j])
         parts[k].x[j] -= dim[j];
     }
-    const int cid = cell_getid(cdim, parts[k].x[0] * ih[0],
-                               parts[k].x[1] * ih[1], parts[k].x[2] * ih[2]);
+    const int cid =
+        cell_getid(cdim, parts[k].x[0] * iwidth[0], parts[k].x[1] * iwidth[1],
+                   parts[k].x[2] * iwidth[2]);
 #ifdef SWIFT_DEBUG_CHECKS
     if (cid < 0 || cid >= s->nr_cells)
       error("Bad cell id %i for part %zi at [%.3e,%.3e,%.3e].", cid, k,
@@ -354,8 +355,9 @@ void engine_redistribute(struct engine *e) {
       else if (gparts[k].x[j] >= dim[j])
         gparts[k].x[j] -= dim[j];
     }
-    const int cid = cell_getid(cdim, gparts[k].x[0] * ih[0],
-                               gparts[k].x[1] * ih[1], gparts[k].x[2] * ih[2]);
+    const int cid =
+        cell_getid(cdim, gparts[k].x[0] * iwidth[0], gparts[k].x[1] * iwidth[1],
+                   gparts[k].x[2] * iwidth[2]);
 #ifdef SWIFT_DEBUG_CHECKS
     if (cid < 0 || cid >= s->nr_cells)
       error("Bad cell id %i for part %zi at [%.3e,%.3e,%.3e].", cid, k,
@@ -540,8 +542,9 @@ void engine_redistribute(struct engine *e) {
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that all parts are in the right place. */
   for (int k = 0; k < nr_parts; k++) {
-    int cid = cell_getid(cdim, parts_new[k].x[0] * ih[0],
-                         parts_new[k].x[1] * ih[1], parts_new[k].x[2] * ih[2]);
+    int cid = cell_getid(cdim, parts_new[k].x[0] * iwidth[0],
+                         parts_new[k].x[1] * iwidth[1],
+                         parts_new[k].x[2] * iwidth[2]);
     if (cells[cid].nodeID != nodeID)
       error("Received particle (%i) that does not belong here (nodeID=%i).", k,
             cells[cid].nodeID);

src/gravity/Default/gravity.h

@@ -42,6 +42,10 @@ gravity_compute_timestep_external(const struct external_potential* potential,
   dt =
       fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
 #endif
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+  dt = fminf(dt, external_gravity_isothermalpotential_timestep(potential,
+                                                               phys_const, gp));
+#endif
 
   return dt;
 }
@@ -117,6 +121,9 @@ __attribute__((always_inline)) INLINE static void external_gravity(
 #ifdef EXTERNAL_POTENTIAL_POINTMASS
   external_gravity_pointmass(potential, phys_const, gp);
 #endif
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+  external_gravity_isothermalpotential(potential, phys_const, gp);
+#endif
 }