diff --git a/examples/Disk-Patch/GravityOnly/README b/examples/Disk-Patch/GravityOnly/README index 00673d274a50f9ab3043adb35d54629e1a993f0f..5bf2638fc5ed48ebe248773223dde888af0c3bc8 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 37752a7a82e82c5d86b2603e3a30128c4d2664c6..78b42e78356f83e80eee8e7f5f91ad7dcf90c37f 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 3277115eea12885c3136d205f341086ca5d98385..5f74d736591b0344696376b18b2c990c1d267396 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 584a8df9c297dbb1337a2cb343f5eb5d9fcdbeba..ebf04675852a7663119ed1ecfd33a05da6b7bb15 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 f7df277501825a677c9b6c36730ee0a6a1d8a0ff..55df81a08d16c6a4f39aa5e9e9205101dedaa3a9 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 edb0885d621975850a09e4298b8e035ebb45a3cd..0a15f5e1628f70306cf021c4bd05bef699d78214 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 6d97634360a1fd5fcbbc2ef329bef926971d3768..cbaa0054fe3bc9f829c59b19b71197264fd0c5af 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 0f656ed544f59e875f149b4a552671c29c2f4c7c..dd7aed8712c01921a01462bf9de713433c81f5c1 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 1771e5ea2fdf017503915f03c7241046e5d06a95..4e33338dac2e8ce2e1c5f9be9ba411fcec31bfbd 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 */