Merge branch 'isothermal_potential' into 'master'

Isothermal potential

added some more external potentials to gravity, plus a test for energy and angular momentum conservation in the IsothermalPotential/GravityOnly example.

See merge request !205

.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/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/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
 }
 
 /**

src/potentials.c

@@ -46,8 +46,25 @@ void potential_init(const struct swift_params* parameter_file,
       parser_get_param_double(parameter_file, "PointMass:position_z");
   potential->point_mass.mass =
       parser_get_param_double(parameter_file, "PointMass:mass");
+  potential->point_mass.timestep_mult =
+      parser_get_param_double(parameter_file, "PointMass:timestep_mult");
 
 #endif /* EXTERNAL_POTENTIAL_POINTMASS */
+
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+
+  potential->isothermal_potential.x =
+      parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
+  potential->isothermal_potential.y =
+      parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
+  potential->isothermal_potential.z =
+      parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
+  potential->isothermal_potential.vrot =
+      parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
+  potential->isothermal_potential.timestep_mult = parser_get_param_double(
+      parameter_file, "IsothermalPotential:timestep_mult");
+
+#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
 }
 
 /**
@@ -58,8 +75,23 @@ void potential_init(const struct swift_params* parameter_file,
 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);
+
+  message(
+      "Point mass properties are (x,y,z) = (%e, %e, %e), M = %e timestep "
+      "multiplier = %e",
+      potential->point_mass.x, potential->point_mass.y, potential->point_mass.z,
+      potential->point_mass.mass, potential->point_mass.timestep_mult);
+
 #endif /* EXTERNAL_POTENTIAL_POINTMASS */
+
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+
+  message(
+      "Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e "
+      "timestep multiplier= %e",
+      potential->isothermal_potential.x, potential->isothermal_potential.y,
+      potential->isothermal_potential.z, potential->isothermal_potential.vrot,
+      potential->isothermal_potential.timestep_mult);
+
+#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
 }

src/potentials.h

@@ -42,15 +42,85 @@ struct external_potential {
   struct {
     double x, y, z;
     double mass;
+    double timestep_mult;
   } point_mass;
 #endif
+
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+  struct {
+    double x, y, z;
+    double vrot;
+    double timestep_mult;
+  } isothermal_potential;
+#endif
 };
 
+/* Include external isothermal potential */
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+
+/**
+ * @brief Computes the time-step due to the acceleration from a point mass
+ *
+ * @param potential The #external_potential used in the run.
+ * @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_isothermalpotential_timestep(
+    const struct external_potential* potential,
+    const struct phys_const* const phys_const, const struct gpart* const g) {
+
+  const float dx = g->x[0] - potential->isothermal_potential.x;
+  const float dy = g->x[1] - potential->isothermal_potential.y;
+  const float dz = g->x[2] - potential->isothermal_potential.z;
+  const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
+  const float drdv =
+      dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
+  const double vrot = potential->isothermal_potential.vrot;
+
+  const float dota_x =
+      vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
+  const float dota_y =
+      vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
+  const float dota_z =
+      vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
+  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 potential->isothermal_potential.timestep_mult * sqrtf(a_2 / dota_2);
+}
+
+/**
+ * @brief Computes the gravitational acceleration of a particle due to a point
+ * mass
+ *
+ * @param potential The #external_potential used in the run.
+ * @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_isothermalpotential(const struct external_potential* potential,
+                                     const struct phys_const* const phys_const,
+                                     struct gpart* g) {
+
+  const float dx = g->x[0] - potential->isothermal_potential.x;
+  const float dy = g->x[1] - potential->isothermal_potential.y;
+  const float dz = g->x[2] - potential->isothermal_potential.z;
+  const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
+
+  const double vrot = potential->isothermal_potential.vrot;
+  g->a_grav[0] += -vrot * vrot * rinv2 * dx;
+  g->a_grav[1] += -vrot * vrot * rinv2 * dy;
+  g->a_grav[2] += -vrot * vrot * rinv2 * dz;
+  // error(" %f %f %f %f", vrot, rinv2, dx, g->a_grav[0]);
+}
+
+#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+
 /* 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
  *
@@ -81,7 +151,7 @@ external_gravity_pointmass_timestep(const struct external_potential* potential,
   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);
+  return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2);
 }
 
 /**
@@ -109,10 +179,10 @@ __attribute__((always_inline)) INLINE static void external_gravity_pointmass(
   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);