Commit c0d32233 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Removed the softened isothermal potential. The isothermal potential is now...

Removed the softened isothermal potential. The isothermal potential is now softened by default. Set epsilon to 0 to restore the old behaviour.
parent e8c2ac7a
...@@ -766,9 +766,6 @@ case "$with_potential" in ...@@ -766,9 +766,6 @@ case "$with_potential" in
isothermal) isothermal)
AC_DEFINE([EXTERNAL_POTENTIAL_ISOTHERMAL], [1], [Isothermal external potential]) AC_DEFINE([EXTERNAL_POTENTIAL_ISOTHERMAL], [1], [Isothermal external potential])
;; ;;
softened-isothermal)
AC_DEFINE([EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL], [1], [Softened isothermal external potential])
;;
disc-patch) disc-patch)
AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential]) AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
;; ;;
......
...@@ -37,7 +37,7 @@ InitialConditions: ...@@ -37,7 +37,7 @@ InitialConditions:
shift_z: 0. shift_z: 0.
# External potential parameters # External potential parameters
SoftenedIsothermalPotential: IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0. position_y: 0.
position_z: 0. position_z: 0.
......
...@@ -34,7 +34,7 @@ InitialConditions: ...@@ -34,7 +34,7 @@ InitialConditions:
file_name: CoolingHalo.hdf5 # The file to read file_name: CoolingHalo.hdf5 # The file to read
# External potential parameters # External potential parameters
SoftenedIsothermalPotential: IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0. position_y: 0.
position_z: 0. position_z: 0.
......
...@@ -34,7 +34,7 @@ InitialConditions: ...@@ -34,7 +34,7 @@ InitialConditions:
file_name: Hydrostatic.hdf5 # The file to read file_name: Hydrostatic.hdf5 # The file to read
# External potential parameters # External potential parameters
SoftenedIsothermalPotential: IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0. position_y: 0.
position_z: 0. position_z: 0.
......
...@@ -81,15 +81,7 @@ IsothermalPotential: ...@@ -81,15 +81,7 @@ IsothermalPotential:
position_z: 100. position_z: 100.
vrot: 200. # Rotation speed of isothermal potential (internal units) vrot: 200. # Rotation speed of isothermal potential (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition 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) epsilon: 0.1 # Softening size (internal units)
timestep_mult: 0.03 # controls time step
# Disk-patch potential parameters # Disk-patch potential parameters
DiscPatchPotential: DiscPatchPotential:
......
...@@ -79,7 +79,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h ...@@ -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 \ riemann/riemann_exact.h riemann/riemann_vacuum.h \
potential/none/potential.h potential/point_mass/potential.h \ potential/none/potential.h potential/point_mass/potential.h \
potential/isothermal/potential.h potential/disc_patch/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/none/cooling.h cooling/none/cooling_struct.h \
cooling/const_du/cooling.h cooling/const_du/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 \ cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \
......
...@@ -34,8 +34,6 @@ ...@@ -34,8 +34,6 @@
#include "./potential/point_mass/potential.h" #include "./potential/point_mass/potential.h"
#elif defined(EXTERNAL_POTENTIAL_ISOTHERMAL) #elif defined(EXTERNAL_POTENTIAL_ISOTHERMAL)
#include "./potential/isothermal/potential.h" #include "./potential/isothermal/potential.h"
#elif defined(EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL)
#include "./potential/softened_isothermal/potential.h"
#elif defined(EXTERNAL_POTENTIAL_DISC_PATCH) #elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
#include "./potential/disc_patch/potential.h" #include "./potential/disc_patch/potential.h"
#else #else
......
/******************************************************************************* /*******************************************************************************
* This file is part of SWIFT. * This file is part of SWIFT.
* Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk) * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@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 * 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 * it under the terms of the GNU Lesser General Public License as published
...@@ -35,7 +36,8 @@ ...@@ -35,7 +36,8 @@
#include "units.h" #include "units.h"
/** /**
* @brief External Potential Properties - Isothermal sphere case * @brief External Potential Properties - Isothermal sphere case with
* central softening
*/ */
struct external_potential { struct external_potential {
...@@ -45,9 +47,14 @@ struct external_potential { ...@@ -45,9 +47,14 @@ struct external_potential {
/*! Rotation velocity */ /*! Rotation velocity */
double vrot; 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; 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 */ /*! Time-step condition pre-factor */
double timestep_mult; double timestep_mult;
}; };
...@@ -70,17 +77,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( ...@@ -70,17 +77,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const float dy = g->x[1] - potential->y; const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z; 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 = const float drdv =
dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]); dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
const double vrot = potential->vrot; const double vrot = potential->vrot;
const float dota_x = const float dota_x = vrot * vrot * r2_plus_epsilon2_inv *
vrot * vrot * rinv2 * (g->v_full[0] - 2.f * drdv * dx * rinv2); (g->v_full[0] - 2.f * drdv * dx * r2_plus_epsilon2_inv);
const float dota_y = const float dota_y = vrot * vrot * r2_plus_epsilon2_inv *
vrot * vrot * rinv2 * (g->v_full[1] - 2.f * drdv * dy * rinv2); (g->v_full[1] - 2.f * drdv * dy * r2_plus_epsilon2_inv);
const float dota_z = const float dota_z = vrot * vrot * r2_plus_epsilon2_inv *
vrot * vrot * rinv2 * (g->v_full[2] - 2.f * drdv * dz * rinv2); (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 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] + 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]; g->a_grav[2] * g->a_grav[2];
...@@ -94,6 +102,7 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( ...@@ -94,6 +102,7 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
* Note that the accelerations are multiplied by Newton's G constant * Note that the accelerations are multiplied by Newton's G constant
* later on. * later on.
* *
* a = v_rot^2 * (x,y,z) / (r^2 + epsilon^2)
* @param time The current time. * @param time The current time.
* @param potential The #external_potential used in the run. * @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units. * @param phys_const The physical constants in internal units.
...@@ -106,10 +115,10 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( ...@@ -106,10 +115,10 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
const float dx = g->x[0] - potential->x; const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y; const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z; 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 * r2_plus_epsilon2_inv;
const double term = -potential->vrot2_over_G * rinv2;
g->a_grav[0] += term * dx; g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy; g->a_grav[1] += term * dy;
...@@ -134,9 +143,8 @@ external_gravity_get_potential_energy( ...@@ -134,9 +143,8 @@ external_gravity_get_potential_energy(
const float dz = g->x[2] - potential->z; const float dz = g->x[2] - potential->z;
return 0.5f * potential->vrot * potential->vrot * 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 * @brief Initialises the external potential properties in the internal system
* of units. * of units.
...@@ -164,9 +172,11 @@ static INLINE void potential_init_backend( ...@@ -164,9 +172,11 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "IsothermalPotential:vrot"); parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
potential->timestep_mult = parser_get_param_float( potential->timestep_mult = parser_get_param_float(
parameter_file, "IsothermalPotential:timestep_mult"); parameter_file, "IsothermalPotential:timestep_mult");
const double epsilon =
parser_get_param_double(parameter_file, "IsothermalPotential:epsilon");
potential->vrot2_over_G = potential->vrot2_over_G =
potential->vrot * potential->vrot / phys_const->const_newton_G; potential->vrot * potential->vrot / phys_const->const_newton_G;
potential->epsilon2 = epsilon * epsilon;
} }
/** /**
...@@ -180,9 +190,9 @@ static INLINE void potential_print_backend( ...@@ -180,9 +190,9 @@ static INLINE void potential_print_backend(
message( message(
"External potential is 'Isothermal' with properties are (x,y,z) = (%e, " "External potential is 'Isothermal' with properties are (x,y,z) = (%e, "
"%e, %e), vrot = %e " "%e, %e), vrot = %e "
"timestep multiplier = %e.", "timestep multiplier = %e, epsilon = %e",
potential->x, potential->y, potential->z, potential->vrot, 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 */
/*******************************************************************************
* 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 */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment