Skip to content
Snippets Groups Projects
Commit fe50afb3 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Disk-patch now without drag force. Disk growth time-scale is now a parameter.

parent 9c308b79
No related branches found
No related tags found
1 merge request!217Disk patch
Setup for a potential of a patch disk, see Creasey, Theuns & Setup for a potential of a patch disk, see Creasey, Theuns &
Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948 Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
The density is given by
rho(z) = (Sigma/2b) / cosh^2(z/b) rho(z) = (Sigma/2b) / cosh^2(z/b)
where Sigma is the surface density, and b the scale height 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
The corresponding force is
dphi/dz = 2 pi G Sigma tanh(z/b),
which satifies d^2phi/dz^2 = 4 pi G rho.
...@@ -34,28 +34,10 @@ SPH: ...@@ -34,28 +34,10 @@ SPH:
# Parameters related to the initial conditions # Parameters related to the initial conditions
InitialConditions: InitialConditions:
file_name: Disk-Patch.hdf5 # The file to read 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 # 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: Disk-PatchPotential:
surface_density: 10. surface_density: 10.
scale_height: 100. scale_height: 100.
z_disk: 300. z_disk: 300.
timestep_mult: 0.03 timestep_mult: 0.03
generate and evolve a disk-patch, where gas is in hydrostatic 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, equations from Creasey, Theuns & Bower, 2013, MNRAS, Volume 429,
Issue 3, p.1922-1948 Issue 3, p.1922-1948
...@@ -21,5 +21,4 @@ dynamica times ...@@ -21,5 +21,4 @@ dynamica times
Then, use the output as new initial conditions, run using Then, use the output as new initial conditions, run using
#define EXTERNAL_POTENTIAL_DISK_PATCH #define EXTERNAL_POTENTIAL_DISK_PATCH
only, and see whether it stays in hydrotatic equilibrium only, and see whether it stays in hydrotatic equilibrium
Use disk-patch.yml as parameter file Use disk-patch.yml as parameter file.
...@@ -21,12 +21,12 @@ Statistics: ...@@ -21,12 +21,12 @@ Statistics:
Snapshots: Snapshots:
basename: Disk-Patch # Common part of the name of output files basename: Disk-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units) 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 # Parameters for the hydrodynamics scheme
SPH: SPH:
resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel). 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. 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_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). max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
...@@ -34,28 +34,11 @@ SPH: ...@@ -34,28 +34,11 @@ SPH:
# Parameters related to the initial conditions # Parameters related to the initial conditions
InitialConditions: InitialConditions:
file_name: Disk-Patch.hdf5 # The file to read 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 # 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: Disk-PatchPotential:
surface_density: 10. surface_density: 10.
scale_height: 100. scale_height: 100.
z_disk: 200. z_disk: 200.
timestep_mult: 0.03 timestep_mult: 0.03
growth_time: 5.
...@@ -19,14 +19,14 @@ Statistics: ...@@ -19,14 +19,14 @@ Statistics:
# Parameters governing the snapshots # Parameters governing the snapshots
Snapshots: Snapshots:
basename: Disk-Patch-dynamic # Common part of the name of output files basename: Disk-Patch-dynamic # Common part of the name of output files
time_first: 968. # Time of the first output (in internal units) time_first: 968. # Time of the first output (in internal units)
delta_time: 24. # Time difference between consecutive outputs (in internal units) delta_time: 24. # Time difference between consecutive outputs (in internal units)
# Parameters for the hydrodynamics scheme # Parameters for the hydrodynamics scheme
SPH: SPH:
resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel). 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. 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_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). max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
...@@ -34,28 +34,10 @@ SPH: ...@@ -34,28 +34,10 @@ SPH:
# Parameters related to the initial conditions # Parameters related to the initial conditions
InitialConditions: InitialConditions:
file_name: Disk-Patch-dynamic.hdf5 # The file to read 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 # 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: Disk-PatchPotential:
surface_density: 10. surface_density: 10.
scale_height: 100. scale_height: 100.
z_disk: 200. z_disk: 200.
timestep_mult: 0.03 timestep_mult: 0.03
...@@ -66,15 +66,22 @@ DomainDecomposition: ...@@ -66,15 +66,22 @@ DomainDecomposition:
# Point mass external potentials # Point mass external potentials
PointMass: 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_y: 50.
position_z: 50. position_z: 50.
mass: 1e10 # mass of external point mass in internal units mass: 1e10 # mass of external point mass (internal units)
timestep_mult: 0.03 # Pre-factor for the time-step condition timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
IsothermalPotential: 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_y: 100.
position_z: 100. position_z: 100.
vrot: 200. # Rotation speed of isothermal potential in internal units vrot: 200. # Rotation speed of isothermal potential (internal units)
timestep_mult: 0.03 # Pre-factor for the time-step condition 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)
...@@ -76,9 +76,6 @@ ...@@ -76,9 +76,6 @@
#define EXTERNAL_POTENTIAL_POINTMASS #define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//#define EXTERNAL_POTENTIAL_DISK_PATCH //#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 ? */ /* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS //#define SWIFT_DEBUG_CHECKS
......
...@@ -76,6 +76,8 @@ void potential_init(const struct swift_params* parameter_file, ...@@ -76,6 +76,8 @@ void potential_init(const struct swift_params* parameter_file,
parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk"); parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
potential->disk_patch_potential.timestep_mult = parser_get_param_double( potential->disk_patch_potential.timestep_mult = parser_get_param_double(
parameter_file, "Disk-PatchPotential:timestep_mult"); 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 = potential->disk_patch_potential.dynamical_time =
sqrt(potential->disk_patch_potential.scale_height / sqrt(potential->disk_patch_potential.scale_height /
(phys_const->const_newton_G * (phys_const->const_newton_G *
...@@ -94,7 +96,7 @@ void potential_print(const struct external_potential* potential) { ...@@ -94,7 +96,7 @@ void potential_print(const struct external_potential* potential) {
message( message(
"Point mass properties are (x,y,z) = (%e, %e, %e), M = %e timestep " "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.x, potential->point_mass.y, potential->point_mass.z,
potential->point_mass.mass, potential->point_mass.timestep_mult); potential->point_mass.mass, potential->point_mass.timestep_mult);
...@@ -104,24 +106,25 @@ void potential_print(const struct external_potential* potential) { ...@@ -104,24 +106,25 @@ void potential_print(const struct external_potential* potential) {
message( message(
"Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e " "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.x, potential->isothermal_potential.y,
potential->isothermal_potential.z, potential->isothermal_potential.vrot, potential->isothermal_potential.z, potential->isothermal_potential.vrot,
potential->isothermal_potential.timestep_mult); potential->isothermal_potential.timestep_mult);
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */ #endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH #ifdef EXTERNAL_POTENTIAL_DISK_PATCH
message( message(
"Disk-patch potential properties are surface_density = %e disk height= " "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.surface_density,
potential->disk_patch_potential.z_disk, potential->disk_patch_potential.z_disk,
potential->disk_patch_potential.scale_height, potential->disk_patch_potential.scale_height,
potential->disk_patch_potential.timestep_mult); potential->disk_patch_potential.timestep_mult);
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
message( if (potential->disk_patch_potential.growth_time > 0.)
"Disk-patch potential: imposing growth of gravity over time, and adding " message("Disk will grow for %f dynamiccal times.",
"viscous force to gravity"); potential->disk_patch_potential.growth_time);
#endif
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */ #endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
} }
...@@ -61,6 +61,7 @@ struct external_potential { ...@@ -61,6 +61,7 @@ struct external_potential {
double scale_height; double scale_height;
double z_disk; double z_disk;
double dynamical_time; double dynamical_time;
double growth_time;
double timestep_mult; double timestep_mult;
} disk_patch_potential; } disk_patch_potential;
#endif #endif
...@@ -148,30 +149,19 @@ external_gravity_disk_patch_potential( ...@@ -148,30 +149,19 @@ external_gravity_disk_patch_potential(
g->a_grav[0] += 0; g->a_grav[0] += 0;
g->a_grav[1] += 0; g->a_grav[1] += 0;
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
float reduction_factor = 1.; float reduction_factor = 1.;
if (time < 5 * t_dyn) reduction_factor = time / (5 * t_dyn); if (time < potential->disk_patch_potential.growth_time * t_dyn)
#else reduction_factor =
const float reduction_factor = 1.; time / (potential->disk_patch_potential.growth_time * t_dyn);
#endif
const float z_accel = const float z_accel =
reduction_factor * 2 * G_newton * M_PI * reduction_factor * 2 * G_newton * M_PI *
potential->disk_patch_potential.surface_density * potential->disk_patch_potential.surface_density *
tanh(fabs(dz) / potential->disk_patch_potential.scale_height); tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
if (dz > 0) /* Accelerations. Note that they are multiplied by G later on */
g->a_grav[2] -= if (dz > 0) g->a_grav[2] -= z_accel / G_newton;
z_accel /
G_newton; /* returned acceleraton is 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 */ #endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment