diff --git a/examples/IsothermalPotential/GravityOnly/externalGravity.yml b/examples/IsothermalPotential/GravityOnly/externalGravity.yml index 53c8a0cf20a0fd9666b5f9fb43baeff6045ce978..150b06dd0f3fbcc0f3ef60be24dbb0554bd63e4d 100644 --- a/examples/IsothermalPotential/GravityOnly/externalGravity.yml +++ b/examples/IsothermalPotential/GravityOnly/externalGravity.yml @@ -63,11 +63,12 @@ PointMass: 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 external point mass in internal units + 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 diff --git a/src/potentials.c b/src/potentials.c index d55d3d669f1a146a0932c46b551ee06f37ed3921..282ac8de9252744848f0fdc713c30e756bfba126 100644 --- a/src/potentials.c +++ b/src/potentials.c @@ -46,6 +46,8 @@ void potential_init(const struct swift_params* parameter_file, parser_get_param_double(parameter_file, "PointMass:position_z"); potential->point_mass.mass = parser_get_param_double(parameter_file, "PointMass:mass"); + potential->point_mass.timestep_mult = + parser_get_param_double(parameter_file, "PointMass:timestep_mult"); #endif /* EXTERNAL_POTENTIAL_POINTMASS */ #ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL @@ -57,6 +59,8 @@ void potential_init(const struct swift_params* parameter_file, parser_get_param_double(parameter_file, "IsothermalPotential:position_z"); potential -> isothermal_potential.vrot = parser_get_param_double(parameter_file, "IsothermalPotential:vrot"); + potential -> isothermal_potential.timestep_mult = + parser_get_param_double(parameter_file, "IsothermalPotential:timestep_mult"); #endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */ } @@ -68,13 +72,15 @@ void potential_init(const struct swift_params* parameter_file, void potential_print(const struct external_potential* potential) { #ifdef EXTERNAL_POTENTIAL_POINTMASS - message("Point mass properties are (x,y,z) = (%e, %e, %e), M = %e", + message("Point mass properties are (x,y,z) = (%e, %e, %e), M = %e timestep multiplier = %e", potential->point_mass.x, potential->point_mass.y, - potential->point_mass.z, potential->point_mass.mass); + potential->point_mass.z, potential->point_mass.mass, + potential->point_mass.timestep_mult); #endif /* EXTERNAL_POTENTIAL_POINTMASS */ #ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL - message("Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e", + message("Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e timestep multiplier= %e", 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); #endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */ } diff --git a/src/potentials.h b/src/potentials.h index c2ee511af3e390a4942c176867f64b8c7bd43fa5..8dbb60bf264009c75b97ed190219d038778f31fb 100644 --- a/src/potentials.h +++ b/src/potentials.h @@ -42,12 +42,14 @@ struct external_potential { struct { double x, y, z; double mass; + double timestep_mult; } point_mass; #endif #ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL struct { double x, y, z; double vrot; + double timestep_mult; } isothermal_potential; #endif }; @@ -83,7 +85,8 @@ __attribute__((always_inline)) // error(" dx= %e dvx= %e rinv2= %e a2= %e dota_2= %e dt= %e", dx, g->v_full[1], rinv2, a_2, dota_2, 0.03f * sqrtf(a_2 / dota_2)); - return EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR * sqrtf(a_2 / dota_2); + //return EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR * sqrtf(a_2 / dota_2); + return potential->isothermal_potential.timestep_mult * sqrtf(a_2 / dota_2); } /** * @brief Computes the gravitational acceleration of a particle due to a point @@ -142,7 +145,8 @@ external_gravity_pointmass_timestep(const struct external_potential* potential, 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 EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR * sqrtf(a_2 / dota_2); + /* return EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR * sqrtf(a_2 / dota_2); */ + return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2); } /**