From 730192c653705c215e8dc622c87bf1776e4252aa Mon Sep 17 00:00:00 2001
From: Stefan Arridge <stefan.arridge@durham.ac.uk>
Date: Fri, 16 Sep 2016 14:45:26 +0100
Subject: [PATCH] Added space struct as an argument to potential_init. Can use
to set centre of potential to be the box centre
---
examples/HydrostaticHalo/hydrostatic.yml | 2 +-
examples/HydrostaticHalo/makeIC.py | 7 ++++++-
examples/HydrostaticHalo/run.sh | 2 ++
examples/main.c | 2 +-
src/const.h | 4 ++--
src/potential.c | 3 ++-
src/potential.h | 1 +
src/potential/disc_patch/potential.h | 2 ++
src/potential/isothermal/potential.h | 8 +++++---
src/potential/none/potential.h | 2 ++
src/potential/point_mass/potential.h | 2 ++
11 files changed, 26 insertions(+), 9 deletions(-)
diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index ad7b4abe67..177f950ea2 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -9,7 +9,7 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 0.1 # The end time of the simulation (in internal units).
+ time_end: 1.0 # 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-2 # The maximal time-step size of the simulation (in internal units).
diff --git a/examples/HydrostaticHalo/makeIC.py b/examples/HydrostaticHalo/makeIC.py
index 49bd73f4fa..09e7734b0e 100644
--- a/examples/HydrostaticHalo/makeIC.py
+++ b/examples/HydrostaticHalo/makeIC.py
@@ -73,7 +73,7 @@ const_G = ((CONST_G_CGS*const_unit_mass_in_cgs*const_unit_time_in
print 'G=', const_G
# Parameters
-periodic= 0 # 1 For periodic box
+periodic= 1 # 1 For periodic box
boxSize = 4.
G = const_G
N = int(sys.argv[1]) # Number of particles
@@ -128,6 +128,11 @@ coords += np.full((N,3),boxSize/2.)
print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
+
+print np.mean(coords[:,0])
+print np.mean(coords[:,1])
+print np.mean(coords[:,2])
+
ds = grp.create_dataset('Coordinates', (N, 3), 'd')
ds[()] = coords
coords = np.zeros(1)
diff --git a/examples/HydrostaticHalo/run.sh b/examples/HydrostaticHalo/run.sh
index 7b462f447b..b37c3b4532 100755
--- a/examples/HydrostaticHalo/run.sh
+++ b/examples/HydrostaticHalo/run.sh
@@ -9,3 +9,5 @@ python makeIC.py 1000
python radial_profile.py 10
python internal_energy_profile.py 10
+
+python test_energy_conservation.py
diff --git a/examples/main.c b/examples/main.c
index 895bb9003e..d7080ee1f1 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -435,7 +435,7 @@ int main(int argc, char *argv[]) {
/* Initialise the external potential properties */
struct external_potential potential;
if (with_external_gravity)
- potential_init(params, &prog_const, &us, &potential);
+ potential_init(params, &prog_const, &us, &s, &potential);
if (with_external_gravity && myrank == 0) potential_print(&potential);
/* Initialise the cooling function properties */
diff --git a/src/const.h b/src/const.h
index d905f636dc..1855d04de7 100644
--- a/src/const.h
+++ b/src/const.h
@@ -92,9 +92,9 @@
/* External gravity properties */
-#define EXTERNAL_POTENTIAL_NONE
+//#define EXTERNAL_POTENTIAL_NONE
//#define EXTERNAL_POTENTIAL_POINTMASS
-//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//#define EXTERNAL_POTENTIAL_DISC_PATCH
diff --git a/src/potential.c b/src/potential.c
index 5433a05e3e..534540ea5b 100644
--- a/src/potential.c
+++ b/src/potential.c
@@ -36,9 +36,10 @@
void potential_init(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_init_backend(parameter_file, phys_const, us, potential);
+ potential_init_backend(parameter_file, phys_const, us, s, potential);
}
/**
diff --git a/src/potential.h b/src/potential.h
index 77bd41794a..de23882770 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -47,6 +47,7 @@
void potential_init(const struct swift_params* parameter_file,
const struct phys_const* phys_const,
const struct UnitSystem* us,
+ const struct space* s,
struct external_potential* potential);
void potential_print(const struct external_potential* potential);
diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h
index 21d168818e..116b9dab36 100644
--- a/src/potential/disc_patch/potential.h
+++ b/src/potential/disc_patch/potential.h
@@ -34,6 +34,7 @@
#include "part.h"
#include "physical_constants.h"
#include "units.h"
+#include "space.h"
/**
* @brief External Potential Properties - Disc patch case
@@ -161,6 +162,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
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->surface_density = parser_get_param_double(
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
index a993c09a97..b1f0ef95b2 100644
--- a/src/potential/isothermal/potential.h
+++ b/src/potential/isothermal/potential.h
@@ -32,6 +32,7 @@
#include "part.h"
#include "physical_constants.h"
#include "units.h"
+#include "space.h"
/**
* @brief External Potential Properties - Isothermal sphere case
@@ -126,13 +127,14 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
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 =
+ potential->x = s->dim[0]/2. +
parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
- potential->y =
+ potential->y = s->dim[1]/2. +
parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
- potential->z =
+ potential->z = s->dim[2]/2. +
parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
potential->vrot =
parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
diff --git a/src/potential/none/potential.h b/src/potential/none/potential.h
index 8b1c3e8415..3f075f8fdc 100644
--- a/src/potential/none/potential.h
+++ b/src/potential/none/potential.h
@@ -31,6 +31,7 @@
#include "part.h"
#include "physical_constants.h"
#include "units.h"
+#include "space.h"
/**
* @brief External Potential Properties
@@ -81,6 +82,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
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) {}
/**
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index f718af2e2c..28f9e81b43 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -32,6 +32,7 @@
#include "part.h"
#include "physical_constants.h"
#include "units.h"
+#include "space.h"
/**
* @brief External Potential Properties - Point mass case
@@ -127,6 +128,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
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 =
--
GitLab