diff --git a/configure.ac b/configure.ac
index a02dcc57c720f1a9a792b160485caee728a91b98..38b99e39b0383ef0ef50198de16cd61f3aa80ed6 100644
--- a/configure.ac
+++ b/configure.ac
@@ -766,9 +766,6 @@ case "$with_potential" in
isothermal)
AC_DEFINE([EXTERNAL_POTENTIAL_ISOTHERMAL], [1], [Isothermal external potential])
;;
- softened-isothermal)
- AC_DEFINE([EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL], [1], [Softened isothermal external potential])
- ;;
disc-patch)
AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
;;
diff --git a/examples/CoolingHalo/cooling_halo.yml b/examples/CoolingHalo/cooling_halo.yml
index c06b099eb0dd06d39040e0ecc8e8f1320a89ac6b..e8978ad6c96017d9b5fbe35346555e6b59bc7e7d 100644
--- a/examples/CoolingHalo/cooling_halo.yml
+++ b/examples/CoolingHalo/cooling_halo.yml
@@ -37,7 +37,7 @@ InitialConditions:
shift_z: 0.
# External potential parameters
-SoftenedIsothermalPotential:
+IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0.
position_z: 0.
diff --git a/examples/CoolingHaloWithSpin/cooling_halo.yml b/examples/CoolingHaloWithSpin/cooling_halo.yml
index 684dd11fcf7adc9477d199e599dfb5b76faa91f6..3c2266f6fc28a2989c2b2a1909adddfb8e41abf6 100644
--- a/examples/CoolingHaloWithSpin/cooling_halo.yml
+++ b/examples/CoolingHaloWithSpin/cooling_halo.yml
@@ -34,7 +34,7 @@ InitialConditions:
file_name: CoolingHalo.hdf5 # The file to read
# External potential parameters
-SoftenedIsothermalPotential:
+IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0.
position_z: 0.
diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index 39a91a4ec475a70ef4e61b9cdc59b8221a74093e..d20d6018f323de0628a0500d8ba767018711fd0a 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -34,7 +34,7 @@ InitialConditions:
file_name: Hydrostatic.hdf5 # The file to read
# External potential parameters
-SoftenedIsothermalPotential:
+IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0.
position_z: 0.
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 899bfb02243c69d3c0d0a94cb05c4f63cb62df32..9681e880e9f4b44d4f77195eecf5b8642266ce52 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -81,15 +81,7 @@ IsothermalPotential:
position_z: 100.
vrot: 200. # Rotation speed of isothermal potential (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
-
-# External potential parameters
-SoftenedIsothermalPotential:
- position_x: 0. # Location of centre of isothermal potential with respect to centre of the box (internal units)
- position_y: 0.
- position_z: 0.
- vrot: 200. # rotation speed of isothermal potential (internal units)
epsilon: 0.1 # Softening size (internal units)
- timestep_mult: 0.03 # controls time step
# Disk-patch potential parameters
DiscPatchPotential:
diff --git a/src/Makefile.am b/src/Makefile.am
index bce5e2e36518afeb148240e1067b3f4e8c817820..32262b8acfb0b7f2765b875fc6098cc8799a0307 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -79,7 +79,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
riemann/riemann_exact.h riemann/riemann_vacuum.h \
potential/none/potential.h potential/point_mass/potential.h \
potential/isothermal/potential.h potential/disc_patch/potential.h \
- potential/softened_isothermal/potential.h \
cooling/none/cooling.h cooling/none/cooling_struct.h \
cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \
cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \
diff --git a/src/potential.h b/src/potential.h
index c462806e206e0e0455bf7094708ab003b7ca9682..116ea8302e7f706cdb861540a89d562174d73408 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -34,8 +34,6 @@
#include "./potential/point_mass/potential.h"
#elif defined(EXTERNAL_POTENTIAL_ISOTHERMAL)
#include "./potential/isothermal/potential.h"
-#elif defined(EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL)
-#include "./potential/softened_isothermal/potential.h"
#elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
#include "./potential/disc_patch/potential.h"
#else
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
index a582dce17daba0ac9705ef4ae1fc6be9db19315a..78e074eec363164b43668f2d716da3aeaacecb9b 100644
--- a/src/potential/isothermal/potential.h
+++ b/src/potential/isothermal/potential.h
@@ -1,7 +1,8 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Stefan Arridge (stefan.arridge@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
@@ -35,7 +36,8 @@
#include "units.h"
/**
- * @brief External Potential Properties - Isothermal sphere case
+ * @brief External Potential Properties - Isothermal sphere case with
+ * central softening
*/
struct external_potential {
@@ -45,9 +47,14 @@ struct external_potential {
/*! Rotation velocity */
double vrot;
- /*! Square of vrot divided by G \f$ \frac{v_{rot}^2}{G} \f$ */
+ /*! Square of vrot, the circular velocity which defines the isothermal
+ * potential devided by Newton's constant */
double vrot2_over_G;
+ /*! Square of the softening length. Acceleration tends to zero within this
+ * distance from the origin */
+ double epsilon2;
+
/*! Time-step condition pre-factor */
double timestep_mult;
};
@@ -70,17 +77,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
- const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
+ const float r2_plus_epsilon2_inv =
+ 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
const float drdv =
dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
const double vrot = potential->vrot;
- const float dota_x =
- vrot * vrot * rinv2 * (g->v_full[0] - 2.f * drdv * dx * rinv2);
- const float dota_y =
- vrot * vrot * rinv2 * (g->v_full[1] - 2.f * drdv * dy * rinv2);
- const float dota_z =
- vrot * vrot * rinv2 * (g->v_full[2] - 2.f * drdv * dz * rinv2);
+ const float dota_x = vrot * vrot * r2_plus_epsilon2_inv *
+ (g->v_full[0] - 2.f * drdv * dx * r2_plus_epsilon2_inv);
+ const float dota_y = vrot * vrot * r2_plus_epsilon2_inv *
+ (g->v_full[1] - 2.f * drdv * dy * r2_plus_epsilon2_inv);
+ const float dota_z = vrot * vrot * r2_plus_epsilon2_inv *
+ (g->v_full[2] - 2.f * drdv * dz * r2_plus_epsilon2_inv);
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];
@@ -94,6 +102,7 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
* Note that the accelerations are multiplied by Newton's G constant
* later on.
*
+ * a = v_rot^2 * (x,y,z) / (r^2 + epsilon^2)
* @param time The current time.
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
@@ -106,10 +115,10 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
+ const float r2_plus_epsilon2_inv =
+ 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
- const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
-
- const double term = -potential->vrot2_over_G * rinv2;
+ const double term = -potential->vrot2_over_G * r2_plus_epsilon2_inv;
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
@@ -134,9 +143,8 @@ external_gravity_get_potential_energy(
const float dz = g->x[2] - potential->z;
return 0.5f * potential->vrot * potential->vrot *
- logf(dx * dx + dy * dy * dz * dz);
+ logf(dx * dx + dy * dy * dz * dz + potential->epsilon2);
}
-
/**
* @brief Initialises the external potential properties in the internal system
* of units.
@@ -164,9 +172,11 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
potential->timestep_mult = parser_get_param_float(
parameter_file, "IsothermalPotential:timestep_mult");
-
+ const double epsilon =
+ parser_get_param_double(parameter_file, "IsothermalPotential:epsilon");
potential->vrot2_over_G =
potential->vrot * potential->vrot / phys_const->const_newton_G;
+ potential->epsilon2 = epsilon * epsilon;
}
/**
@@ -180,9 +190,9 @@ static INLINE void potential_print_backend(
message(
"External potential is 'Isothermal' with properties are (x,y,z) = (%e, "
"%e, %e), vrot = %e "
- "timestep multiplier = %e.",
+ "timestep multiplier = %e, epsilon = %e",
potential->x, potential->y, potential->z, potential->vrot,
- potential->timestep_mult);
+ potential->timestep_mult, sqrtf(potential->epsilon2));
}
-#endif /* SWIFT_POTENTIAL_ISOTHERMAL_H */
+#endif /* SWIFT_ISOTHERMAL_H */
diff --git a/src/potential/softened_isothermal/potential.h b/src/potential/softened_isothermal/potential.h
deleted file mode 100644
index 24e59b12a5745728fb1189fbbfbc7cc3c06fbfa6..0000000000000000000000000000000000000000
--- a/src/potential/softened_isothermal/potential.h
+++ /dev/null
@@ -1,196 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Stefan Arridge (stefan.arridge@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_POTENTIAL_SOFTENED_ISOTHERMAL_H
-#define SWIFT_POTENTIAL_SOFTENED_ISOTHERMAL_H
-
-/* Config parameters. */
-#include "../config.h"
-
-/* Some standard headers. */
-#include <math.h>
-
-/* Local includes. */
-#include "error.h"
-#include "parser.h"
-#include "part.h"
-#include "physical_constants.h"
-#include "space.h"
-#include "units.h"
-
-/**
- * @brief External Potential Properties - Softened Isothermal sphere case
- */
-struct external_potential {
-
- /*! Position of the centre of potential */
- double x, y, z;
-
- /*! Rotation velocity */
- double vrot;
-
- /*! Square of vrot, the circular velocity which defines the isothermal
- * potential */
- double vrot2_over_G;
-
- /*! Square of the softening length. Acceleration tends to zero within this
- * distance from the origin */
- double epsilon2;
-
- /*! Time-step condition pre-factor */
- double timestep_mult;
-};
-
-/**
- * @brief Computes the time-step due to the acceleration from an isothermal
- * potential.
- *
- * @param time The current time.
- * @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_timestep(
- double time, const struct external_potential* restrict potential,
- const struct phys_const* restrict phys_const,
- const struct gpart* restrict g) {
-
- const float dx = g->x[0] - potential->x;
- const float dy = g->x[1] - potential->y;
- const float dz = g->x[2] - potential->z;
-
- const float r2_plus_epsilon2_inv =
- 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
- const float drdv =
- dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
- const double vrot = potential->vrot;
-
- const float dota_x = vrot * vrot * r2_plus_epsilon2_inv *
- (g->v_full[0] - 2.f * drdv * dx * r2_plus_epsilon2_inv);
- const float dota_y = vrot * vrot * r2_plus_epsilon2_inv *
- (g->v_full[1] - 2.f * drdv * dy * r2_plus_epsilon2_inv);
- const float dota_z = vrot * vrot * r2_plus_epsilon2_inv *
- (g->v_full[2] - 2.f * drdv * dz * r2_plus_epsilon2_inv);
- 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->timestep_mult * sqrtf(a_2 / dota_2);
-}
-
-/**
- * @brief Computes the gravitational acceleration from an isothermal potential.
- *
- * Note that the accelerations are multiplied by Newton's G constant
- * later on.
- *
- * a = v_rot^2 * (x,y,z) / (r^2 + epsilon^2)
- * @param time The current time.
- * @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_acceleration(
- double time, const struct external_potential* potential,
- const struct phys_const* const phys_const, struct gpart* g) {
-
- const float dx = g->x[0] - potential->x;
- const float dy = g->x[1] - potential->y;
- const float dz = g->x[2] - potential->z;
- const float r2_plus_epsilon2_inv =
- 1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
-
- const double term = -potential->vrot2_over_G * r2_plus_epsilon2_inv;
-
- g->a_grav[0] += term * dx;
- g->a_grav[1] += term * dy;
- g->a_grav[2] += term * dz;
-}
-
-/**
- * @brief Computes the gravitational potential energy of a particle in an
- * isothermal potential.
- *
- * @param potential The #external_potential used in the run.
- * @param phys_const Physical constants in internal units.
- * @param g Pointer to the particle data.
- */
-__attribute__((always_inline)) INLINE static float
-external_gravity_get_potential_energy(
- const struct external_potential* potential,
- const struct phys_const* const phys_const, const struct gpart* g) {
-
- const float dx = g->x[0] - potential->x;
- const float dy = g->x[1] - potential->y;
- const float dz = g->x[2] - potential->z;
-
- return 0.5f * potential->vrot * potential->vrot *
- logf(dx * dx + dy * dy * dz * dz + potential->epsilon2);
-}
-/**
- * @brief Initialises the external potential properties in the internal system
- * of units.
- *
- * @param parameter_file The parsed parameter file
- * @param phys_const Physical constants in internal units
- * @param us The current internal system of units
- * @param potential The external potential properties to initialize
- */
-static INLINE void potential_init_backend(
- const struct swift_params* parameter_file,
- const struct phys_const* phys_const, const struct UnitSystem* us,
- const struct space* s, struct external_potential* potential) {
-
- potential->x = s->dim[0] / 2. +
- parser_get_param_double(
- parameter_file, "SoftenedIsothermalPotential:position_x");
- potential->y = s->dim[1] / 2. +
- parser_get_param_double(
- parameter_file, "SoftenedIsothermalPotential:position_y");
- potential->z = s->dim[2] / 2. +
- parser_get_param_double(
- parameter_file, "SoftenedIsothermalPotential:position_z");
- potential->vrot = parser_get_param_double(parameter_file,
- "SoftenedIsothermalPotential:vrot");
- potential->timestep_mult = parser_get_param_float(
- parameter_file, "SoftenedIsothermalPotential:timestep_mult");
- const double epsilon = parser_get_param_float(
- parameter_file, "SoftenedIsothermalPotential:epsilon");
- potential->vrot2_over_G =
- potential->vrot * potential->vrot / phys_const->const_newton_G;
- potential->epsilon2 = epsilon * epsilon;
-}
-
-/**
- * @brief Prints the properties of the external potential to stdout.
- *
- * @param potential The external potential properties.
- */
-static INLINE void potential_print_backend(
- const struct external_potential* potential) {
-
- message(
- "External potential is 'Isothermal' with properties are (x,y,z) = (%e, "
- "%e, %e), vrot = %e "
- "timestep multiplier = %e, epsilon = %e",
- potential->x, potential->y, potential->z, potential->vrot,
- potential->timestep_mult, sqrtf(potential->epsilon2));
-}
-
-#endif /* SWIFT_POTENTIAL_ISOTHERMAL_H */