From fe50afb3a36d2453eefe29142f265f8733c2ab84 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Thu, 18 Aug 2016 23:42:39 +0100
Subject: [PATCH] Disk-patch now without drag force. Disk growth time-scale is
now a parameter.
---
examples/Disk-Patch/GravityOnly/README | 13 +++----
.../Disk-Patch/GravityOnly/disk-patch.yml | 26 +++-----------
examples/Disk-Patch/HydroStatic/README | 5 ++-
.../Disk-Patch/HydroStatic/disk-patch-icc.yml | 31 ++++-------------
.../Disk-Patch/HydroStatic/disk-patch.yml | 34 +++++--------------
examples/parameter_example.yml | 19 +++++++----
src/const.h | 3 --
src/potentials.c | 19 ++++++-----
src/potentials.h | 22 ++++--------
9 files changed, 56 insertions(+), 116 deletions(-)
diff --git a/examples/Disk-Patch/GravityOnly/README b/examples/Disk-Patch/GravityOnly/README
index 00673d274a..5bf2638fc5 100644
--- a/examples/Disk-Patch/GravityOnly/README
+++ b/examples/Disk-Patch/GravityOnly/README
@@ -1,13 +1,10 @@
Setup for a potential of a patch disk, see Creasey, Theuns &
Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+The density is given by
rho(z) = (Sigma/2b) / cosh^2(z/b)
-where Sigma is the surface density, and b the scale height
-the corresponding force is
-dphi/dz = 2 pi G Sigma tanh(z/b)
-
-which satifies d^2phi/dz^2 = 4 pi G rho
-
-
-
+where Sigma is the surface density, and b the scale height.
+The corresponding force is
+dphi/dz = 2 pi G Sigma tanh(z/b),
+which satifies d^2phi/dz^2 = 4 pi G rho.
diff --git a/examples/Disk-Patch/GravityOnly/disk-patch.yml b/examples/Disk-Patch/GravityOnly/disk-patch.yml
index 37752a7a82..78b42e7835 100644
--- a/examples/Disk-Patch/GravityOnly/disk-patch.yml
+++ b/examples/Disk-Patch/GravityOnly/disk-patch.yml
@@ -34,28 +34,10 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: Disk-Patch.hdf5 # The file to read
- h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
- shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
- shift_y: 0.
- shift_z: 0.
# External potential parameters
-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 # controls time step
-
-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
-
Disk-PatchPotential:
- surface_density: 10.
- scale_height: 100.
- z_disk: 300.
- timestep_mult: 0.03
+ surface_density: 10.
+ scale_height: 100.
+ z_disk: 300.
+ timestep_mult: 0.03
diff --git a/examples/Disk-Patch/HydroStatic/README b/examples/Disk-Patch/HydroStatic/README
index 3277115eea..5f74d73659 100644
--- a/examples/Disk-Patch/HydroStatic/README
+++ b/examples/Disk-Patch/HydroStatic/README
@@ -1,5 +1,5 @@
generate and evolve a disk-patch, where gas is in hydrostatic
-equilibrium with an imposed external gravutation force, using teh
+equilibrium with an imposed external gravutation force, using the
equations from Creasey, Theuns & Bower, 2013, MNRAS, Volume 429,
Issue 3, p.1922-1948
@@ -21,5 +21,4 @@ dynamica times
Then, use the output as new initial conditions, run using
#define EXTERNAL_POTENTIAL_DISK_PATCH
only, and see whether it stays in hydrotatic equilibrium
-Use disk-patch.yml as parameter file
-
+Use disk-patch.yml as parameter file.
diff --git a/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml b/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
index 584a8df9c2..ebf0467585 100644
--- a/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
+++ b/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
@@ -21,12 +21,12 @@ Statistics:
Snapshots:
basename: Disk-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 12. # Time difference between consecutive outputs (in internal units)
+ delta_time: 12. # 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: 0.1 # The tolerance for the targetted number of neighbours.
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
@@ -34,28 +34,11 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: Disk-Patch.hdf5 # The file to read
- h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
- shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
- shift_y: 0.
- shift_z: 0.
# External potential parameters
-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 # controls time step
-
-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
-
Disk-PatchPotential:
- surface_density: 10.
- scale_height: 100.
- z_disk: 200.
- timestep_mult: 0.03
+ surface_density: 10.
+ scale_height: 100.
+ z_disk: 200.
+ timestep_mult: 0.03
+ growth_time: 5.
diff --git a/examples/Disk-Patch/HydroStatic/disk-patch.yml b/examples/Disk-Patch/HydroStatic/disk-patch.yml
index f7df277501..55df81a08d 100644
--- a/examples/Disk-Patch/HydroStatic/disk-patch.yml
+++ b/examples/Disk-Patch/HydroStatic/disk-patch.yml
@@ -19,14 +19,14 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
- basename: Disk-Patch-dynamic # Common part of the name of output files
- time_first: 968. # Time of the first output (in internal units)
- delta_time: 24. # Time difference between consecutive outputs (in internal units)
+ basename: Disk-Patch-dynamic # Common part of the name of output files
+ time_first: 968. # Time of the first output (in internal units)
+ delta_time: 24. # 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: 0.1 # The tolerance for the targetted number of neighbours.
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
@@ -34,28 +34,10 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: Disk-Patch-dynamic.hdf5 # The file to read
- h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
- shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
- shift_y: 0.
- shift_z: 0.
# External potential parameters
-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 # controls time step
-
-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
-
Disk-PatchPotential:
- surface_density: 10.
- scale_height: 100.
- z_disk: 200.
- timestep_mult: 0.03
+ surface_density: 10.
+ scale_height: 100.
+ z_disk: 200.
+ timestep_mult: 0.03
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index edb0885d62..0a15f5e162 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -66,15 +66,22 @@ DomainDecomposition:
# Point mass external potentials
PointMass:
- position_x: 50. # location of external point mass in internal units
+ position_x: 50. # location of external point mass (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
+ mass: 1e10 # mass of external point mass (internal units)
+ timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
IsothermalPotential:
- position_x: 100. # Location of centre of isothermal potential in internal units
+ position_x: 100. # Location of centre of isothermal potential (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
+ vrot: 200. # Rotation speed of isothermal potential (internal units)
+ timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
+
+Disk-PatchPotential:
+ surface_density: 10. # Surface density of the disk (internal units)
+ scale_height: 100. # Scale height of the disk (internal units)
+ z_disk: 200. # Disk height (internal units)
+ timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
+ growth_time: 5. # (Optional) Time for the disk to grow to its final size (multiple of the dynamical time)
diff --git a/src/const.h b/src/const.h
index 6d97634360..cbaa0054fe 100644
--- a/src/const.h
+++ b/src/const.h
@@ -76,9 +76,6 @@
#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//#define EXTERNAL_POTENTIAL_DISK_PATCH
-/* Add viscuous force to gas particles to speed-up glass making for disk-patch
- * ICs */
-//#define EXTERNAL_POTENTIAL_DISK_PATCH_ICS
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/potentials.c b/src/potentials.c
index 0f656ed544..dd7aed8712 100644
--- a/src/potentials.c
+++ b/src/potentials.c
@@ -76,6 +76,8 @@ void potential_init(const struct swift_params* parameter_file,
parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
potential->disk_patch_potential.timestep_mult = parser_get_param_double(
parameter_file, "Disk-PatchPotential:timestep_mult");
+ potential->disk_patch_potential.growth_time = parser_get_opt_param_double(
+ parameter_file, "Disk-PatchPotential:growth_time", 0.);
potential->disk_patch_potential.dynamical_time =
sqrt(potential->disk_patch_potential.scale_height /
(phys_const->const_newton_G *
@@ -94,7 +96,7 @@ void potential_print(const struct external_potential* potential) {
message(
"Point mass properties are (x,y,z) = (%e, %e, %e), M = %e timestep "
- "multiplier = %e",
+ "multiplier = %e.",
potential->point_mass.x, potential->point_mass.y, potential->point_mass.z,
potential->point_mass.mass, potential->point_mass.timestep_mult);
@@ -104,24 +106,25 @@ void potential_print(const struct external_potential* potential) {
message(
"Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e "
- "timestep multiplier= %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 */
+
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+
message(
"Disk-patch potential properties are surface_density = %e disk height= "
- "%e scale height= %e timestep multiplier= %e",
+ "%e scale height = %e timestep multiplier = %e.",
potential->disk_patch_potential.surface_density,
potential->disk_patch_potential.z_disk,
potential->disk_patch_potential.scale_height,
potential->disk_patch_potential.timestep_mult);
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
- message(
- "Disk-patch potential: imposing growth of gravity over time, and adding "
- "viscous force to gravity");
-#endif
+
+ if (potential->disk_patch_potential.growth_time > 0.)
+ message("Disk will grow for %f dynamiccal times.",
+ potential->disk_patch_potential.growth_time);
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
}
diff --git a/src/potentials.h b/src/potentials.h
index 1771e5ea2f..4e33338dac 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -61,6 +61,7 @@ struct external_potential {
double scale_height;
double z_disk;
double dynamical_time;
+ double growth_time;
double timestep_mult;
} disk_patch_potential;
#endif
@@ -148,30 +149,19 @@ external_gravity_disk_patch_potential(
g->a_grav[0] += 0;
g->a_grav[1] += 0;
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
float reduction_factor = 1.;
- if (time < 5 * t_dyn) reduction_factor = time / (5 * t_dyn);
-#else
- const float reduction_factor = 1.;
-#endif
+ if (time < potential->disk_patch_potential.growth_time * t_dyn)
+ reduction_factor =
+ time / (potential->disk_patch_potential.growth_time * t_dyn);
const float z_accel =
reduction_factor * 2 * G_newton * M_PI *
potential->disk_patch_potential.surface_density *
tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
- if (dz > 0)
- g->a_grav[2] -=
- z_accel /
- G_newton; /* returned acceleraton is multiplied by G later on */
+ /* Accelerations. Note that they are multiplied by G later on */
+ if (dz > 0) g->a_grav[2] -= z_accel / G_newton;
if (dz < 0) g->a_grav[2] += z_accel / G_newton;
-
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
- /* TT: add viscous drag */
- g->a_grav[0] -= g->v_full[0] / (2 * t_dyn) / G_newton;
- g->a_grav[1] -= g->v_full[1] / (2 * t_dyn) / G_newton;
- g->a_grav[2] -= g->v_full[2] / (2 * t_dyn) / G_newton;
-#endif
}
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
--
GitLab