From 4fc9291edd0f3200a8f0891ac603d4c60aaaf339 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 12:33:20 +0000 Subject: [PATCH 01/54] Updated dt_alpha to be more consistent with other schemes. --- src/runner.c | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/runner.c b/src/runner.c index dfbb9e214..98985207f 100644 --- a/src/runner.c +++ b/src/runner.c @@ -1254,7 +1254,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) { struct xpart *restrict xparts = c->hydro.xparts; const int count = c->hydro.count; const struct engine *e = r->e; - const integertime_t ti_end = e->ti_current; + const integertime_t ti_current = e->ti_current; const int with_cosmology = (e->policy & engine_policy_cosmology); const double time_base = e->time_base; const struct cosmology *cosmo = e->cosmology; @@ -1289,9 +1289,14 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) { * This is the physical time between the start and end of the time-step * without any scale-factor powers. */ double dt_alpha; + if (with_cosmology) { const integertime_t ti_step = get_integer_timestep(p->time_bin); - dt_alpha = cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end); + const integertime_t ti_begin = + get_integer_time_begin(ti_current - 1, p->time_bin); + + dt_alpha = + cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step); } else { dt_alpha = get_timestep(p->time_bin, time_base); } @@ -1470,13 +1475,16 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { * artificial viscosity and thermal conduction terms) */ const int with_cosmology = (e->policy & engine_policy_cosmology); const double time_base = e->time_base; - const integertime_t ti_end = e->ti_current; + const integertime_t ti_current = e->ti_current; double dt_alpha; if (with_cosmology) { const integertime_t ti_step = get_integer_timestep(p->time_bin); + const integertime_t ti_begin = + get_integer_time_begin(ti_current - 1, p->time_bin); + dt_alpha = - cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end); + cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step); } else { dt_alpha = get_timestep(p->time_bin, time_base); } @@ -1614,13 +1622,17 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { * the evolution of alpha factors (i.e. those involved in the artificial * viscosity and thermal conduction terms) */ const int with_cosmology = (e->policy & engine_policy_cosmology); - const integertime_t ti_end = e->ti_current; const double time_base = e->time_base; + const integertime_t ti_current = e->ti_current; double dt_alpha; if (with_cosmology) { const integertime_t ti_step = get_integer_timestep(p->time_bin); - dt_alpha = cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end); + const integertime_t ti_begin = + get_integer_time_begin(ti_current - 1, p->time_bin); + + dt_alpha = + cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step); } else { dt_alpha = get_timestep(p->time_bin, time_base); } -- GitLab From 4d5da992b465f2a0d5178474a94ce9d83aa7a07c Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 14:54:58 +0000 Subject: [PATCH 02/54] Changed 3.f to hydro_dimension --- src/hydro/Gadget2/hydro.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index bcef24808..f05e3229a 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -562,7 +562,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->density.rot_v[2] * p->density.rot_v[2]); /* Compute the norm of div v including the Hubble flow term */ - const float div_physical_v = p->density.div_v + 3.f * cosmo->H; + const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H; const float abs_div_physical_v = fabsf(div_physical_v); /* Compute the pressure */ -- GitLab From 0d1ebebab226160950a4a9283bfe04b572b79ca8 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 16:02:58 +0000 Subject: [PATCH 03/54] Added hydro_dimension in place of 3.f --- src/hydro/Minimal/hydro.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 1d5d25f73..e2fd00695 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -561,7 +561,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->density.rot_v[2] * p->density.rot_v[2]); /* Compute the norm of div v including the Hubble flow term */ - const float div_physical_v = p->density.div_v + 3.f * cosmo->H; + const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H; const float abs_div_physical_v = fabsf(div_physical_v); /* Compute the pressure */ -- GitLab From ff5b717ce72a9909c43cd0616312ed54573423bc Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 16:25:03 +0000 Subject: [PATCH 04/54] Fixes hubble term in Pressure-Energy schemes for DivV --- src/hydro/PressureEnergy/hydro.h | 4 ++-- src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h | 12 ++++++++---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h index 0e53b5f2e..3fe3c805f 100644 --- a/src/hydro/PressureEnergy/hydro.h +++ b/src/hydro/PressureEnergy/hydro.h @@ -530,8 +530,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( /* Finish calculation of the velocity divergence, including hubble flow term */ - p->density.div_v *= - h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension; + p->density.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2; + p->density.div_v += cosmo->H * hydro_dimension; } /** diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h index 9f981569e..819abeec9 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h @@ -526,8 +526,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv; /* Finish calculation of the velocity divergence */ - p->density.div_v *= - h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension; + p->density.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2; + p->density.div_v += cosmo->H * hydro_dimension; } /** @@ -614,7 +614,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Artificial viscosity updates */ - const float inverse_tau = hydro_props->viscosity.length * soundspeed * h_inv; + /* Need to divide by a^2 here because this is a timescale and we need to + introduce the time-integration factor to balance out the scale-free + nature of the artificial viscosity (i.e. dt has a^2 built in) */ + const float inverse_tau = + (hydro_props->viscosity.length * cosmo->a2_inv) * soundspeed * h_inv; const float source_term = -1.f * min(p->density.div_v, 0.f); /* Compute da/dt */ @@ -869,4 +873,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) { __attribute__((always_inline)) INLINE static void hydro_remove_part( const struct part *p, const struct xpart *xp) {} -#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */ +#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */ \ No newline at end of file -- GitLab From b83abf424dbda881b33351cae26bedb5fe21417b Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 16:25:18 +0000 Subject: [PATCH 05/54] Fixes cosmological integration of the AV and AD in ANARCHY --- src/hydro/AnarchyPU/hydro.h | 44 +++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 37d182889..0403608fe 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -541,8 +541,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv; /* Finish calculation of the velocity divergence */ - p->viscosity.div_v *= - h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension; + p->viscosity.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2; + p->viscosity.div_v += cosmo->H * hydro_dimension; } /** @@ -694,18 +694,37 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Here we need to update the artificial viscosity */ /* Timescale for decay */ - const float tau = - p->h / (2.f * p->viscosity.v_sig * hydro_props->viscosity.length); - /* Construct time differential of div.v implicitly */ - const float div_v_dt = + + const float a2 = cosmo->a * cosmo->a; + const float a_fac_soundspeed = cosmo->a_factor_sound_speed; + + /* We integrate the whole AV scheme in physical units as cancelling the cosmology terms + is just a nightmare, especially for non-standard gas_gamma. */ + + const float h_physical = p->h * cosmo->a; + const float v_sig_physical = p->viscosity.v_sig / a_fac_soundspeed; + const float div_v_physical = + (p->viscosity.div_v - cosmo->H * hydro_dimension) * a2; + const float div_v_prev_physical = + (p->viscosity.div_v_previous_step - cosmo->H * hydro_dimension) * a2; + const float dt_alpha_physical = dt_alpha / a2; + + const float tau = h_physical / (2.f * v_sig_physical * hydro_props->viscosity.length); + + /* Construct time differential of div.v implicitly following the ANARCHY spec */ + + float div_v_dt = dt_alpha == 0.f ? 0.f - : (p->viscosity.div_v - p->viscosity.div_v_previous_step) / dt_alpha; + : (div_v_physical - div_v_prev_physical) / dt_alpha_physical; + /* Construct the source term for the AV; if shock detected this is _positive_ * as div_v_dt should be _negative_ before the shock hits */ - const float S = p->h * p->h * max(0.f, -1.f * div_v_dt); - const float v_sig_square = p->viscosity.v_sig * p->viscosity.v_sig; + const float S = h_physical * h_physical * max(0.f, -1.f * div_v_dt); + const float v_sig_square = v_sig_physical * v_sig_physical; + /* Calculate the current appropriate value of the AV based on the above */ + /* Note this is v_sig_physical squared, not comoving */ const float alpha_loc = hydro_props->viscosity.alpha_max * S / (v_sig_square + S); @@ -717,7 +736,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float alpha_dt = (alpha_loc - p->viscosity.alpha) / tau; /* Finally, we can update the actual value of the alpha */ - p->viscosity.alpha += alpha_dt * dt_alpha; + p->viscosity.alpha += alpha_dt * dt_alpha_physical; } if (p->viscosity.alpha < hydro_props->viscosity.alpha_min) { @@ -732,7 +751,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ /* alpha_diff_dt is cosmology-less */ - /* Evolution term: following Schaller+ 2015 */ + /* Evolution term: following Schaller+ 2015. This is actually cosmology-less and + so requires no conversion to physical. */ float alpha_diff_dt = hydro_props->diffusion.beta * p->h * p->diffusion.laplace_u / sqrt_u; /* Decay term: not documented in Schaller+ 2015 but was present @@ -741,7 +761,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( (p->diffusion.alpha - hydro_props->diffusion.alpha_min) / tau; float new_diffusion_alpha = p->diffusion.alpha; - new_diffusion_alpha += alpha_diff_dt * dt_alpha; + new_diffusion_alpha += alpha_diff_dt * dt_alpha_physical; /* Consistency checks to ensure min < alpha < max */ if (new_diffusion_alpha > hydro_props->diffusion.alpha_max) { -- GitLab From 5d5433c05420ff4844f3a85dc1df59e07c39b75e Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 12:33:20 +0000 Subject: [PATCH 06/54] Updated dt_alpha to be more consistent with other schemes. --- src/runner.c | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/runner.c b/src/runner.c index 719c193e3..24956dc70 100644 --- a/src/runner.c +++ b/src/runner.c @@ -1259,7 +1259,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) { struct xpart *restrict xparts = c->hydro.xparts; const int count = c->hydro.count; const struct engine *e = r->e; - const integertime_t ti_end = e->ti_current; + const integertime_t ti_current = e->ti_current; const int with_cosmology = (e->policy & engine_policy_cosmology); const double time_base = e->time_base; const struct cosmology *cosmo = e->cosmology; @@ -1294,9 +1294,14 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) { * This is the physical time between the start and end of the time-step * without any scale-factor powers. */ double dt_alpha; + if (with_cosmology) { const integertime_t ti_step = get_integer_timestep(p->time_bin); - dt_alpha = cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end); + const integertime_t ti_begin = + get_integer_time_begin(ti_current - 1, p->time_bin); + + dt_alpha = + cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step); } else { dt_alpha = get_timestep(p->time_bin, time_base); } @@ -1475,13 +1480,16 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { * artificial viscosity and thermal conduction terms) */ const int with_cosmology = (e->policy & engine_policy_cosmology); const double time_base = e->time_base; - const integertime_t ti_end = e->ti_current; + const integertime_t ti_current = e->ti_current; double dt_alpha; if (with_cosmology) { const integertime_t ti_step = get_integer_timestep(p->time_bin); + const integertime_t ti_begin = + get_integer_time_begin(ti_current - 1, p->time_bin); + dt_alpha = - cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end); + cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step); } else { dt_alpha = get_timestep(p->time_bin, time_base); } @@ -1619,13 +1627,17 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { * the evolution of alpha factors (i.e. those involved in the artificial * viscosity and thermal conduction terms) */ const int with_cosmology = (e->policy & engine_policy_cosmology); - const integertime_t ti_end = e->ti_current; const double time_base = e->time_base; + const integertime_t ti_current = e->ti_current; double dt_alpha; if (with_cosmology) { const integertime_t ti_step = get_integer_timestep(p->time_bin); - dt_alpha = cosmology_get_delta_time(cosmo, ti_end - ti_step, ti_end); + const integertime_t ti_begin = + get_integer_time_begin(ti_current - 1, p->time_bin); + + dt_alpha = + cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step); } else { dt_alpha = get_timestep(p->time_bin, time_base); } -- GitLab From 0c8488f731ba7c72110f6712f67e05f9d3f6658e Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 14:54:58 +0000 Subject: [PATCH 07/54] Changed 3.f to hydro_dimension --- src/hydro/Gadget2/hydro.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index bcef24808..f05e3229a 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -562,7 +562,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->density.rot_v[2] * p->density.rot_v[2]); /* Compute the norm of div v including the Hubble flow term */ - const float div_physical_v = p->density.div_v + 3.f * cosmo->H; + const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H; const float abs_div_physical_v = fabsf(div_physical_v); /* Compute the pressure */ -- GitLab From 884a8c0cb3dda2b6e27ddc3bc2575b5093abaede Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 16:02:58 +0000 Subject: [PATCH 08/54] Added hydro_dimension in place of 3.f --- src/hydro/Minimal/hydro.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 1d5d25f73..e2fd00695 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -561,7 +561,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->density.rot_v[2] * p->density.rot_v[2]); /* Compute the norm of div v including the Hubble flow term */ - const float div_physical_v = p->density.div_v + 3.f * cosmo->H; + const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H; const float abs_div_physical_v = fabsf(div_physical_v); /* Compute the pressure */ -- GitLab From 91e7802c4e28cdc9981055cd171f8abed974a967 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 16:25:03 +0000 Subject: [PATCH 09/54] Fixes hubble term in Pressure-Energy schemes for DivV --- src/hydro/PressureEnergy/hydro.h | 4 ++-- src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h | 12 ++++++++---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h index 0e53b5f2e..3fe3c805f 100644 --- a/src/hydro/PressureEnergy/hydro.h +++ b/src/hydro/PressureEnergy/hydro.h @@ -530,8 +530,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( /* Finish calculation of the velocity divergence, including hubble flow term */ - p->density.div_v *= - h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension; + p->density.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2; + p->density.div_v += cosmo->H * hydro_dimension; } /** diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h index 9f981569e..819abeec9 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h @@ -526,8 +526,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv; /* Finish calculation of the velocity divergence */ - p->density.div_v *= - h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension; + p->density.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2; + p->density.div_v += cosmo->H * hydro_dimension; } /** @@ -614,7 +614,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Artificial viscosity updates */ - const float inverse_tau = hydro_props->viscosity.length * soundspeed * h_inv; + /* Need to divide by a^2 here because this is a timescale and we need to + introduce the time-integration factor to balance out the scale-free + nature of the artificial viscosity (i.e. dt has a^2 built in) */ + const float inverse_tau = + (hydro_props->viscosity.length * cosmo->a2_inv) * soundspeed * h_inv; const float source_term = -1.f * min(p->density.div_v, 0.f); /* Compute da/dt */ @@ -869,4 +873,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) { __attribute__((always_inline)) INLINE static void hydro_remove_part( const struct part *p, const struct xpart *xp) {} -#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */ +#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */ \ No newline at end of file -- GitLab From 1615b7250502058352a48088e2237286f33a5aad Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 20 Mar 2019 16:25:18 +0000 Subject: [PATCH 10/54] Fixes cosmological integration of the AV and AD in ANARCHY --- src/hydro/AnarchyPU/hydro.h | 44 +++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 37d182889..0403608fe 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -541,8 +541,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv; /* Finish calculation of the velocity divergence */ - p->viscosity.div_v *= - h_inv_dim_plus_one * rho_inv * a_inv2 + cosmo->H * hydro_dimension; + p->viscosity.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2; + p->viscosity.div_v += cosmo->H * hydro_dimension; } /** @@ -694,18 +694,37 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Here we need to update the artificial viscosity */ /* Timescale for decay */ - const float tau = - p->h / (2.f * p->viscosity.v_sig * hydro_props->viscosity.length); - /* Construct time differential of div.v implicitly */ - const float div_v_dt = + + const float a2 = cosmo->a * cosmo->a; + const float a_fac_soundspeed = cosmo->a_factor_sound_speed; + + /* We integrate the whole AV scheme in physical units as cancelling the cosmology terms + is just a nightmare, especially for non-standard gas_gamma. */ + + const float h_physical = p->h * cosmo->a; + const float v_sig_physical = p->viscosity.v_sig / a_fac_soundspeed; + const float div_v_physical = + (p->viscosity.div_v - cosmo->H * hydro_dimension) * a2; + const float div_v_prev_physical = + (p->viscosity.div_v_previous_step - cosmo->H * hydro_dimension) * a2; + const float dt_alpha_physical = dt_alpha / a2; + + const float tau = h_physical / (2.f * v_sig_physical * hydro_props->viscosity.length); + + /* Construct time differential of div.v implicitly following the ANARCHY spec */ + + float div_v_dt = dt_alpha == 0.f ? 0.f - : (p->viscosity.div_v - p->viscosity.div_v_previous_step) / dt_alpha; + : (div_v_physical - div_v_prev_physical) / dt_alpha_physical; + /* Construct the source term for the AV; if shock detected this is _positive_ * as div_v_dt should be _negative_ before the shock hits */ - const float S = p->h * p->h * max(0.f, -1.f * div_v_dt); - const float v_sig_square = p->viscosity.v_sig * p->viscosity.v_sig; + const float S = h_physical * h_physical * max(0.f, -1.f * div_v_dt); + const float v_sig_square = v_sig_physical * v_sig_physical; + /* Calculate the current appropriate value of the AV based on the above */ + /* Note this is v_sig_physical squared, not comoving */ const float alpha_loc = hydro_props->viscosity.alpha_max * S / (v_sig_square + S); @@ -717,7 +736,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float alpha_dt = (alpha_loc - p->viscosity.alpha) / tau; /* Finally, we can update the actual value of the alpha */ - p->viscosity.alpha += alpha_dt * dt_alpha; + p->viscosity.alpha += alpha_dt * dt_alpha_physical; } if (p->viscosity.alpha < hydro_props->viscosity.alpha_min) { @@ -732,7 +751,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ /* alpha_diff_dt is cosmology-less */ - /* Evolution term: following Schaller+ 2015 */ + /* Evolution term: following Schaller+ 2015. This is actually cosmology-less and + so requires no conversion to physical. */ float alpha_diff_dt = hydro_props->diffusion.beta * p->h * p->diffusion.laplace_u / sqrt_u; /* Decay term: not documented in Schaller+ 2015 but was present @@ -741,7 +761,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( (p->diffusion.alpha - hydro_props->diffusion.alpha_min) / tau; float new_diffusion_alpha = p->diffusion.alpha; - new_diffusion_alpha += alpha_diff_dt * dt_alpha; + new_diffusion_alpha += alpha_diff_dt * dt_alpha_physical; /* Consistency checks to ensure min < alpha < max */ if (new_diffusion_alpha > hydro_props->diffusion.alpha_max) { -- GitLab From 17d77af75b166384cd970a5c99f76b6cedd0c76c Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Thu, 21 Mar 2019 10:01:01 +0000 Subject: [PATCH 11/54] Added printing of lambda constant in internal units --- src/cooling/const_lambda/cooling.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 1e50c162d..0d15c297f 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -352,6 +352,10 @@ static INLINE void cooling_print_backend( "* " "cm^3]", cooling->lambda_nH2_cgs); + + message( + "Lambda/n_H^2=%g [internal units]", cooling->lambda_nH2_cgs * cooling->conv_factor_energy_rate_from_cgs + ); if (cooling->cooling_tstep_mult == FLT_MAX) message("Cooling function time-step size is unlimited"); -- GitLab From 7d4fbe6015ed540f898f2630b80b46583546c6d7 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Fri, 22 Mar 2019 15:00:38 +0000 Subject: [PATCH 12/54] Added new example to test cosmology of cooling --- .../cooling_redshift_dependence_low_z.yml | 80 +++++++++++ .../cooling_redshift_dependence_no_z.yml | 72 ++++++++++ .../CoolingRedshiftDependence/getGlass.sh | 2 + .../CoolingRedshiftDependence/makeIC.py | 125 ++++++++++++++++++ .../CoolingRedshiftDependence/plotSolution.py | 111 ++++++++++++++++ .../Cooling/CoolingRedshiftDependence/run.sh | 22 +++ 6 files changed, 412 insertions(+) create mode 100644 examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml create mode 100644 examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml create mode 100755 examples/Cooling/CoolingRedshiftDependence/getGlass.sh create mode 100644 examples/Cooling/CoolingRedshiftDependence/makeIC.py create mode 100644 examples/Cooling/CoolingRedshiftDependence/plotSolution.py create mode 100755 examples/Cooling/CoolingRedshiftDependence/run.sh diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml new file mode 100644 index 000000000..c7c7d65e2 --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml @@ -0,0 +1,80 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun + UnitLength_in_cgs: 3.08567758e24 # 1 Mpc + UnitVelocity_in_cgs: 1e5 # 1 km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Cosmological parameters +Cosmology: + h: 0.6777 # Reduced Hubble constant + a_begin: 0.5 # Initial scale-factor of the simulation (z = 1.0) + a_end: 0.5061559 # Final scale factor of the simulation (~ 100 myr) + Omega_m: 0.307 # Matter density parameter + Omega_lambda: 0.693 # Dark-energy density parameter + Omega_b: 0.0455 # Baryon density parameter + +# Parameters governing the time integration +TimeIntegration: + dt_min: 1e-14 + dt_max: 5e-3 + +# Parameters governing the snapshots +Snapshots: + basename: redshift_dependence_low_z + delta_time: 1.000122373748838 + scale_factor_first: 0.5 + compression: 4 + +# Parameters governing the conserved quantities statistics +Statistics: + scale_factor_first: 0.5 + delta_time: 1.1 + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + minimal_temperature: 100 # K + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./ics_low_z.hdf5 # The file to read + periodic: 1 + +# Parameters for the EAGLE chemistry +EAGLEChemistry: # Solar abundances + init_abundance_metal: 0. + init_abundance_Hydrogen: 0.752 + init_abundance_Helium: 0.248 + init_abundance_Carbon: 0. + init_abundance_Nitrogen: 0. + init_abundance_Oxygen: 0. + init_abundance_Neon: 0. + init_abundance_Magnesium: 0. + init_abundance_Silicon: 0. + init_abundance_Iron: 0. + +# Parameters for the EAGLE cooling +EAGLECooling: + dir_name: ./coolingtables/ + H_reion_z: 11.5 + H_reion_eV_p_H: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_eV_p_H: 2.0 + +# Parameters for the EAGLE "equation of state" +EAGLEEntropyFloor: + Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in. + Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin. + Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor + Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in. + Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. + Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor + +LambdaCooling: + lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] \ No newline at end of file diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml new file mode 100644 index 000000000..ba6fe62d4 --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml @@ -0,0 +1,72 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun + UnitLength_in_cgs: 3.08567758e24 # 1 Mpc + UnitVelocity_in_cgs: 1e5 # 1 km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters governing the time integration +TimeIntegration: + time_begin: 0. + time_end: 1.02e-4 # ~ 100 myr + dt_min: 1e-14 + dt_max: 5e-8 + +# Parameters governing the snapshots +Snapshots: + basename: redshift_dependence_no_z + delta_time: 1.02e-6 + compression: 4 + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1.02e-6 + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + minimal_temperature: 100 # K + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./ics_no_z.hdf5 # The file to read + periodic: 1 + +# Parameters for the EAGLE chemistry +EAGLEChemistry: # Solar abundances + init_abundance_metal: 0. + init_abundance_Hydrogen: 0.752 + init_abundance_Helium: 0.248 + init_abundance_Carbon: 0. + init_abundance_Nitrogen: 0. + init_abundance_Oxygen: 0. + init_abundance_Neon: 0. + init_abundance_Magnesium: 0. + init_abundance_Silicon: 0. + init_abundance_Iron: 0. + +# Parameters for the EAGLE cooling +EAGLECooling: + dir_name: ./coolingtables/ + H_reion_z: 11.5 + H_reion_eV_p_H: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_eV_p_H: 2.0 + +# Parameters for the EAGLE "equation of state" +EAGLEEntropyFloor: + Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in. + Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin. + Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor + Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in. + Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. + Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor + +LambdaCooling: + lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] + cooling_tstep_mult: 1e5 # (Optional) Dimensionless pre-factor for the time-step condition. diff --git a/examples/Cooling/CoolingRedshiftDependence/getGlass.sh b/examples/Cooling/CoolingRedshiftDependence/getGlass.sh new file mode 100755 index 000000000..01b4474ac --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/getGlass.sh @@ -0,0 +1,2 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/gravity_glassCube_32.hdf5 diff --git a/examples/Cooling/CoolingRedshiftDependence/makeIC.py b/examples/Cooling/CoolingRedshiftDependence/makeIC.py new file mode 100644 index 000000000..47e3540d5 --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/makeIC.py @@ -0,0 +1,125 @@ +""" +Creates a redshift-dependent cooling box such that it always has the same +_physical_ density at the given redshift. +""" + +from swiftsimio import Writer +from swiftsimio.units import cosmo_units + +from unyt import mh, cm, s, K, Mpc, kb +import numpy as np + +import h5py + +# Physics parameters. +boxsize = 1.0 * Mpc +physical_density = 1.0 * mh / cm ** 3 +mu_hydrogen = 0.5 # Fully ionised +temperature = 1e7 * K +gamma = 5.0 / 3.0 + + +def get_coordinates(glass_filename: str) -> np.array: + """ + Gets the coordinates from the glass file. + """ + + with h5py.File(glass_filename, "r") as handle: + coordinates = handle["PartType1/Coordinates"][...] + + return coordinates + + +def generate_ics(redshift: float, filename: str, glass_filename: str) -> None: + """ + Generate initial conditions for the CoolingRedshiftDependence example. + """ + + scale_factor = 1 / (1 + redshift) + energy_factor = scale_factor ** (-3.0 * (gamma - 1)) + comoving_boxsize = boxsize * scale_factor + + glass_coordinates = get_coordinates(glass_filename) + number_of_particles = len(glass_coordinates) + + gas_particle_mass = physical_density * (comoving_boxsize ** 3) / number_of_particles + + writer = Writer(cosmo_units, comoving_boxsize) + + writer.gas.coordinates = glass_coordinates * comoving_boxsize + + writer.gas.velocities = np.zeros_like(glass_coordinates) * cm / s + + writer.gas.masses = np.ones(number_of_particles, dtype=float) * gas_particle_mass + + writer.gas.internal_energy = ( + np.ones(number_of_particles, dtype=float) + * energy_factor + * (3.0 / 2.0) + * (temperature * kb) + / (mu_hydrogen * mh) + ) + + writer.gas.generate_smoothing_lengths(boxsize=comoving_boxsize, dimension=3) + + writer.write(filename) + + return + + +if __name__ == "__main__": + """ + Sets up the initial parameters. + """ + + import argparse as ap + + parser = ap.ArgumentParser( + description=""" + Sets up the initial conditions for the cooling test. Takes two + redshifts, and produces two files: ics_high_z.hdf5 and + ics_low_z.hdf5. + """ + ) + + parser.add_argument( + "-a", + "--high", + help="The high redshift to generate initial conditions for. Default: 10.0", + default=10, + type=float, + ) + + parser.add_argument( + "-b", + "--low", + help="The low redshift to generate initial conditions for. Default: 1.0", + default=1, + type=float, + ) + + parser.add_argument( + "-n", + "--nocosmo", + help="Generate a non-cosmological box? Default: Truthy", + default=1, + type=bool, + ) + + parser.add_argument( + "-g", + "--glass", + help="Glass filename. Default: gravity_glassCube_32.hdf5", + type=str, + default="gravity_glassCube_32.hdf5", + ) + + args = parser.parse_args() + + generate_ics(args.low, filename="ics_low_z.hdf5", glass_filename=args.glass) + generate_ics(args.high, filename="ics_high_z.hdf5", glass_filename=args.glass) + + if args.nocosmo: + generate_ics(1.0, filename="ics_no_z.hdf5", glass_filename=args.glass) + + exit(0) diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py new file mode 100644 index 000000000..44427be76 --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py @@ -0,0 +1,111 @@ +""" +Plots the mean temperature. +""" + +import matplotlib.pyplot as plt + +from swiftsimio import load + +from unyt import Myr, K, mh, cm + +import numpy as np + + +def setup_axes(): + """ + Sets up the axes and returns fig, ax + """ + + fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(4, 8), sharex=True) + + ax = ax.flatten() + + ax[0].semilogy() + + ax[1].set_xlabel("Simulation time elapsed [Myr]") + ax[0].set_ylabel("Temperature of Universe [K]") + ax[1].set_ylabel("Physical Density of Universe [$n_H$ cm$^{-3}$]") + + ax[0].set_xlim(0, 50) + + fig.tight_layout() + + return fig, ax + + +def get_data(handle: float, n_snaps: int): + """ + Returns the elapsed simulation time, temperature, and density + """ + + t0 = 0.0 + + times = [] + temps = [] + densities = [] + + for snap in range(n_snaps): + data = load(f"{handle}_{snap:04d}.hdf5") + + if snap == 0: + t0 = data.metadata.t.to(Myr).value + + times.append(data.metadata.t.to(Myr).value - t0) + temps.append(np.mean(data.gas.temperature.to(K).value)) + densities.append(np.mean(data.gas.density.to(mh / cm**3).value) / (data.metadata.scale_factor**3) ) + + return times, temps, densities + + +def plot_single_data(handle: str, n_snaps: int, label: str, ax: plt.axes): + """ + Takes the a single simulation and plots it on the axes. + """ + + data = get_data(handle, n_snaps) + + ax[0].plot( + data[0], data[1], + label=label, + marker="o", + ms=2 + ) + + ax[1].plot( + data[0], data[2], + label=label, + marker="o", + ms=2 + ) + + + return + + +def make_plot(handles, names, n_snaps=100): + """ + Makes the whole plot and returns the fig, ax objects. + """ + + fig, ax = setup_axes() + + for handle, name in zip(handles, names): + plot_single_data(handle, n_snaps, name, ax) + + ax[0].legend() + + return fig, ax + + +if __name__ == "__main__": + """ + Plot everything! + """ + + handles = ["redshift_dependence_no_z", "redshift_dependence_low_z"] + names = ["No Cosmology", "Low Redshift"] + + fig, ax = make_plot(handles, names) + + fig.savefig("redshift_dependence.png", dpi=300) + diff --git a/examples/Cooling/CoolingRedshiftDependence/run.sh b/examples/Cooling/CoolingRedshiftDependence/run.sh new file mode 100755 index 000000000..25e600f65 --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/run.sh @@ -0,0 +1,22 @@ +#!/bin/bash + +# Generate the initial conditions if they are not present. +if [ ! -e gravity_glassCube_32.hdf5 ] +then + echo "Fetching initial gravity glass file for the constant cosmological box example..." + ./getGlass.sh +fi + +# Fetch the cooling tables +if [ ! -e ics_no_z.hdf5 ] +then + echo "Generating initial conditions for the uniform cosmo box example..." + python3 makeIC.py +fi + +swift_location="../../../build/anarchy-const-lambda/examples/swift" + +# Run SWIFT +$swift_location --hydro --cosmology --cooling --threads=4 cooling_redshift_dependence_low_z.yml 2>&1 | tee output_low.log +$swift_location --hydro --cooling --threads=4 cooling_redshift_dependence_no_z.yml 2>&1 | tee output_no.log + -- GitLab From d65deba5b730dc3b38f15ac7d3d913368dab7889 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Fri, 22 Mar 2019 17:30:25 +0000 Subject: [PATCH 13/54] Minor improvements to the new cooling test --- .../cooling_redshift_dependence_high_z.yml | 80 +++++++++++++++++++ .../cooling_redshift_dependence_low_z.yml | 10 +-- .../cooling_redshift_dependence_no_z.yml | 1 - .../CoolingRedshiftDependence/makeIC.py | 18 ++--- .../CoolingRedshiftDependence/plotSolution.py | 6 +- .../Cooling/CoolingRedshiftDependence/run.sh | 9 ++- .../ComovingSodShock_3D/sodShock.yml | 43 +++++++++- 7 files changed, 144 insertions(+), 23 deletions(-) create mode 100644 examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml new file mode 100644 index 000000000..e9fb6d995 --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml @@ -0,0 +1,80 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun + UnitLength_in_cgs: 3.08567758e24 # 1 Mpc + UnitVelocity_in_cgs: 1e5 # 1 km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Cosmological parameters +Cosmology: + h: 0.6777 # Reduced Hubble constant + a_begin: 0.5 # Initial scale-factor of the simulation (z = 1.0) + a_end: 0.5061559 # Final scale factor of the simulation (~ 100 myr) + Omega_m: 0.307 # Matter density parameter + Omega_lambda: 0.693 # Dark-energy density parameter + Omega_b: 0.0455 # Baryon density parameter + +# Parameters governing the time integration +TimeIntegration: + dt_min: 1e-14 + dt_max: 5e-3 + +# Parameters governing the snapshots +Snapshots: + basename: redshift_dependence_high_z + delta_time: 1.000122373748838 + scale_factor_first: 0.5 + compression: 4 + +# Parameters governing the conserved quantities statistics +Statistics: + scale_factor_first: 0.5 + delta_time: 1.1 + +# Parameters for the hydrodynamics scheme +SPH: + resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). + CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. + minimal_temperature: 100 # K + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./ics_high_z.hdf5 # The file to read + periodic: 1 + +# Parameters for the EAGLE chemistry +EAGLEChemistry: # Solar abundances + init_abundance_metal: 0. + init_abundance_Hydrogen: 0.752 + init_abundance_Helium: 0.248 + init_abundance_Carbon: 0. + init_abundance_Nitrogen: 0. + init_abundance_Oxygen: 0. + init_abundance_Neon: 0. + init_abundance_Magnesium: 0. + init_abundance_Silicon: 0. + init_abundance_Iron: 0. + +# Parameters for the EAGLE cooling +EAGLECooling: + dir_name: ./coolingtables/ + H_reion_z: 11.5 + H_reion_eV_p_H: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_eV_p_H: 2.0 + +# Parameters for the EAGLE "equation of state" +EAGLEEntropyFloor: + Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in. + Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin. + Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor + Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in. + Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. + Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor + +LambdaCooling: + lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] \ No newline at end of file diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml index c7c7d65e2..159a5e3f4 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml @@ -9,8 +9,8 @@ InternalUnitSystem: # Cosmological parameters Cosmology: h: 0.6777 # Reduced Hubble constant - a_begin: 0.5 # Initial scale-factor of the simulation (z = 1.0) - a_end: 0.5061559 # Final scale factor of the simulation (~ 100 myr) + a_begin: 0.99009 # Initial scale-factor of the simulation (z = 0.01) + a_end: 0.99700 # Final scale factor of the simulation (~ 100 myr) Omega_m: 0.307 # Matter density parameter Omega_lambda: 0.693 # Dark-energy density parameter Omega_b: 0.0455 # Baryon density parameter @@ -23,13 +23,13 @@ TimeIntegration: # Parameters governing the snapshots Snapshots: basename: redshift_dependence_low_z - delta_time: 1.000122373748838 - scale_factor_first: 0.5 + delta_time: 1.0000695206950205 + scale_factor_first: 0.99009 compression: 4 # Parameters governing the conserved quantities statistics Statistics: - scale_factor_first: 0.5 + scale_factor_first: 0.99009 delta_time: 1.1 # Parameters for the hydrodynamics scheme diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml index ba6fe62d4..8b00ca0ac 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml @@ -69,4 +69,3 @@ EAGLEEntropyFloor: LambdaCooling: lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] - cooling_tstep_mult: 1e5 # (Optional) Dimensionless pre-factor for the time-step condition. diff --git a/examples/Cooling/CoolingRedshiftDependence/makeIC.py b/examples/Cooling/CoolingRedshiftDependence/makeIC.py index 47e3540d5..489be7636 100644 --- a/examples/Cooling/CoolingRedshiftDependence/makeIC.py +++ b/examples/Cooling/CoolingRedshiftDependence/makeIC.py @@ -36,13 +36,12 @@ def generate_ics(redshift: float, filename: str, glass_filename: str) -> None: """ scale_factor = 1 / (1 + redshift) - energy_factor = scale_factor ** (-3.0 * (gamma - 1)) - comoving_boxsize = boxsize * scale_factor + comoving_boxsize = boxsize / scale_factor glass_coordinates = get_coordinates(glass_filename) number_of_particles = len(glass_coordinates) - gas_particle_mass = physical_density * (comoving_boxsize ** 3) / number_of_particles + gas_particle_mass = physical_density * (boxsize ** 3) / number_of_particles writer = Writer(cosmo_units, comoving_boxsize) @@ -52,10 +51,9 @@ def generate_ics(redshift: float, filename: str, glass_filename: str) -> None: writer.gas.masses = np.ones(number_of_particles, dtype=float) * gas_particle_mass + # Leave in physical units; handled by boxsize change. writer.gas.internal_energy = ( np.ones(number_of_particles, dtype=float) - * energy_factor - * (3.0 / 2.0) * (temperature * kb) / (mu_hydrogen * mh) ) @@ -85,16 +83,16 @@ if __name__ == "__main__": parser.add_argument( "-a", "--high", - help="The high redshift to generate initial conditions for. Default: 10.0", - default=10, + help="The high redshift to generate initial conditions for. Default: 1.0", + default=1, type=float, ) parser.add_argument( "-b", "--low", - help="The low redshift to generate initial conditions for. Default: 1.0", - default=1, + help="The low redshift to generate initial conditions for. Default: 0.01", + default=0.01, type=float, ) @@ -120,6 +118,6 @@ if __name__ == "__main__": generate_ics(args.high, filename="ics_high_z.hdf5", glass_filename=args.glass) if args.nocosmo: - generate_ics(1.0, filename="ics_no_z.hdf5", glass_filename=args.glass) + generate_ics(0.0, filename="ics_no_z.hdf5", glass_filename=args.glass) exit(0) diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py index 44427be76..d8a4ddc0f 100644 --- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py +++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py @@ -26,7 +26,7 @@ def setup_axes(): ax[0].set_ylabel("Temperature of Universe [K]") ax[1].set_ylabel("Physical Density of Universe [$n_H$ cm$^{-3}$]") - ax[0].set_xlim(0, 50) + ax[0].set_xlim(0, 100) fig.tight_layout() @@ -102,8 +102,8 @@ if __name__ == "__main__": Plot everything! """ - handles = ["redshift_dependence_no_z", "redshift_dependence_low_z"] - names = ["No Cosmology", "Low Redshift"] + handles = ["redshift_dependence_no_z", "redshift_dependence_low_z", "redshift_dependence_high_z"] + names = ["No Cosmology", "Low Redshift ($z=0.01$)", "High Redshift ($z=1$)"] fig, ax = make_plot(handles, names) diff --git a/examples/Cooling/CoolingRedshiftDependence/run.sh b/examples/Cooling/CoolingRedshiftDependence/run.sh index 25e600f65..ba399811e 100755 --- a/examples/Cooling/CoolingRedshiftDependence/run.sh +++ b/examples/Cooling/CoolingRedshiftDependence/run.sh @@ -14,9 +14,16 @@ then python3 makeIC.py fi -swift_location="../../../build/anarchy-const-lambda/examples/swift" +swift_location="../../../build/gadget-const-lambda/examples/swift" + +rm redshift_dependence_*_z_*.hdf5 # Run SWIFT +$swift_location --hydro --cosmology --cooling --threads=4 cooling_redshift_dependence_high_z.yml 2>&1 | tee output_high.log +mv timesteps_4.txt timesteps_high.txt $swift_location --hydro --cosmology --cooling --threads=4 cooling_redshift_dependence_low_z.yml 2>&1 | tee output_low.log +mv timesteps_4.txt timesteps_low.txt $swift_location --hydro --cooling --threads=4 cooling_redshift_dependence_no_z.yml 2>&1 | tee output_no.log +mv timesteps_4.txt timesteps_no.txt +python3 plotSolution.py diff --git a/examples/Cosmology/ComovingSodShock_3D/sodShock.yml b/examples/Cosmology/ComovingSodShock_3D/sodShock.yml index 2d7a5727c..f7cf51ce0 100644 --- a/examples/Cosmology/ComovingSodShock_3D/sodShock.yml +++ b/examples/Cosmology/ComovingSodShock_3D/sodShock.yml @@ -1,14 +1,14 @@ # Define the system of units to use internally. InternalUnitSystem: - UnitMass_in_cgs: 2.94e55 # Grams + UnitMass_in_cgs: 2.94e45 # Grams UnitLength_in_cgs: 3.086e18 # pc - UnitVelocity_in_cgs: 1. # km per s + UnitVelocity_in_cgs: 1e10 # km per s UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin # Parameters governing the time integration TimeIntegration: - dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). + dt_min: 1e-18 # The minimal time-step size of the simulation (in internal units). dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units). # Parameters governing the snapshots @@ -16,11 +16,13 @@ Snapshots: basename: sodShock # Common part of the name of output files time_first: 0. # Time of the first output (in internal units) delta_time: 1.06638 # Time difference between consecutive outputs (in internal units) + # scale_factor_first: 1.0 scale_factor_first: 0.001 compression: 1 # Parameters governing the conserved quantities statistics Statistics: + scale_factor_first: 1.0 delta_time: 1.02 # Time between statistics output # Parameters for the hydrodynamics scheme @@ -41,3 +43,38 @@ Cosmology: a_begin: 0.001 a_end: 0.00106638 +# Impose primoridal metallicity +EAGLEChemistry: + init_abundance_metal: 0. + init_abundance_Hydrogen: 0.752 + init_abundance_Helium: 0.248 + init_abundance_Carbon: 0.0 + init_abundance_Nitrogen: 0.0 + init_abundance_Oxygen: 0.0 + init_abundance_Neon: 0.0 + init_abundance_Magnesium: 0.0 + init_abundance_Silicon: 0.0 + init_abundance_Iron: 0.0 + +EAGLECooling: + dir_name: ./coolingtables/ + H_reion_z: 11.5 + H_reion_eV_p_H: 2.0 + He_reion_z_centre: 3.5 + He_reion_z_sigma: 0.5 + He_reion_eV_p_H: 2.0 + +# Parameters for the EAGLE "equation of state" +EAGLEEntropyFloor: + Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in. + Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin. + Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor + Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. + Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in. + Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin. + Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor + +LambdaCooling: + lambda_nH2_cgs: 1e-48 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3]) + cooling_tstep_mult: 1.0 # (Optional) Dimensionless pre-factor for the time-step condition. -- GitLab From 242dda4f172c38746f174dbdfc162a27da9186a7 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Mon, 25 Mar 2019 20:17:23 +0000 Subject: [PATCH 14/54] Changed initial density --- examples/Cooling/CoolingRedshiftDependence/makeIC.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/makeIC.py b/examples/Cooling/CoolingRedshiftDependence/makeIC.py index 489be7636..61b875668 100644 --- a/examples/Cooling/CoolingRedshiftDependence/makeIC.py +++ b/examples/Cooling/CoolingRedshiftDependence/makeIC.py @@ -13,7 +13,7 @@ import h5py # Physics parameters. boxsize = 1.0 * Mpc -physical_density = 1.0 * mh / cm ** 3 +physical_density = 0.1 * mh / cm ** 3 mu_hydrogen = 0.5 # Fully ionised temperature = 1e7 * K gamma = 5.0 / 3.0 @@ -54,6 +54,7 @@ def generate_ics(redshift: float, filename: str, glass_filename: str) -> None: # Leave in physical units; handled by boxsize change. writer.gas.internal_energy = ( np.ones(number_of_particles, dtype=float) + * 3.0 / 2.0 * (temperature * kb) / (mu_hydrogen * mh) ) -- GitLab From 6169996ded76da0c66c3953aab44affe83d28962 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Mon, 25 Mar 2019 20:17:38 +0000 Subject: [PATCH 15/54] Changed to using a data directory --- .../cooling_redshift_dependence_high_z.yml | 9 +++++++-- .../cooling_redshift_dependence_low_z.yml | 17 +++++++++++------ .../cooling_redshift_dependence_no_z.yml | 7 ++++++- 3 files changed, 24 insertions(+), 9 deletions(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml index e9fb6d995..13ed73554 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml @@ -22,7 +22,7 @@ TimeIntegration: # Parameters governing the snapshots Snapshots: - basename: redshift_dependence_high_z + basename: data/redshift_dependence_high_z delta_time: 1.000122373748838 scale_factor_first: 0.5 compression: 4 @@ -77,4 +77,9 @@ EAGLEEntropyFloor: Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor LambdaCooling: - lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] \ No newline at end of file + lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] + +ConstCooling: + cooling_rate: 1. # Cooling rate (du/dt) (internal units) + min_energy: 1. # Minimal internal energy per unit mass (internal units) + cooling_tstep_mult: 1. # Dimensionless pre-factor for the time-step condition \ No newline at end of file diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml index 159a5e3f4..5b336bc4b 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml @@ -1,8 +1,8 @@ # Define the system of units to use internally. InternalUnitSystem: - UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun - UnitLength_in_cgs: 3.08567758e24 # 1 Mpc - UnitVelocity_in_cgs: 1e5 # 1 km/s + UnitMass_in_cgs: 1.98848e33 # 10^10 M_sun + UnitLength_in_cgs: 3.08567758e21 # 1 kpc + UnitVelocity_in_cgs: 1 # 1 km/s UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin @@ -22,7 +22,7 @@ TimeIntegration: # Parameters governing the snapshots Snapshots: - basename: redshift_dependence_low_z + basename: data/redshift_dependence_low_z delta_time: 1.0000695206950205 scale_factor_first: 0.99009 compression: 4 @@ -44,7 +44,7 @@ InitialConditions: periodic: 1 # Parameters for the EAGLE chemistry -EAGLEChemistry: # Solar abundances +EAGLEChemistry: # Primordial init_abundance_metal: 0. init_abundance_Hydrogen: 0.752 init_abundance_Helium: 0.248 @@ -77,4 +77,9 @@ EAGLEEntropyFloor: Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor LambdaCooling: - lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] \ No newline at end of file + lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] + +ConstCooling: + cooling_rate: 1. # Cooling rate (du/dt) (internal units) + min_energy: 1. # Minimal internal energy per unit mass (internal units) + cooling_tstep_mult: 1. # Dimensionless pre-factor for the time-step condition diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml index 8b00ca0ac..35242f553 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml @@ -15,7 +15,7 @@ TimeIntegration: # Parameters governing the snapshots Snapshots: - basename: redshift_dependence_no_z + basename: data/redshift_dependence_no_z delta_time: 1.02e-6 compression: 4 @@ -69,3 +69,8 @@ EAGLEEntropyFloor: LambdaCooling: lambda_nH2_cgs: 2.13744785e-23 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3] + +ConstCooling: + cooling_rate: 1. # Cooling rate (du/dt) (internal units) + min_energy: 1. # Minimal internal energy per unit mass (internal units) + cooling_tstep_mult: 1. # Dimensionless pre-factor for the time-step condition \ No newline at end of file -- GitLab From 1fbf625795efe2ced4d6789e1e7b7a62a47a6fef Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Mon, 25 Mar 2019 20:18:34 +0000 Subject: [PATCH 16/54] Moved to using data directory --- examples/Cooling/CoolingRedshiftDependence/run.sh | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/run.sh b/examples/Cooling/CoolingRedshiftDependence/run.sh index ba399811e..c9bf2aea7 100755 --- a/examples/Cooling/CoolingRedshiftDependence/run.sh +++ b/examples/Cooling/CoolingRedshiftDependence/run.sh @@ -14,16 +14,18 @@ then python3 makeIC.py fi -swift_location="../../../build/gadget-const-lambda/examples/swift" +swift_location="../../swift" -rm redshift_dependence_*_z_*.hdf5 +rm data/redshift_dependence_*_z_*.hdf5 + +mkdir data # Run SWIFT -$swift_location --hydro --cosmology --cooling --threads=4 cooling_redshift_dependence_high_z.yml 2>&1 | tee output_high.log +$swift_location --hydro --cosmology --cooling --limiter --threads=4 cooling_redshift_dependence_high_z.yml 2>&1 | tee output_high.log mv timesteps_4.txt timesteps_high.txt -$swift_location --hydro --cosmology --cooling --threads=4 cooling_redshift_dependence_low_z.yml 2>&1 | tee output_low.log +$swift_location --hydro --cosmology --cooling --limiter --threads=4 cooling_redshift_dependence_low_z.yml 2>&1 | tee output_low.log mv timesteps_4.txt timesteps_low.txt -$swift_location --hydro --cooling --threads=4 cooling_redshift_dependence_no_z.yml 2>&1 | tee output_no.log +$swift_location --hydro --cooling --limiter --threads=4 cooling_redshift_dependence_no_z.yml 2>&1 | tee output_no.log mv timesteps_4.txt timesteps_no.txt python3 plotSolution.py -- GitLab From 61a81fb1def2b8b1a685f046a7f86b596b2f840a Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Mon, 25 Mar 2019 20:18:48 +0000 Subject: [PATCH 17/54] Added plot of diffusion and viscosity parameters. --- .../HydroTests/SedovBlast_3D/plotSolution.py | 48 ++++++++++++++++--- .../HydroTests/SodShock_3D/plotSolution.py | 46 +++++++++++++++--- 2 files changed, 82 insertions(+), 12 deletions(-) diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py index ad34695d3..a87e3cd24 100644 --- a/examples/HydroTests/SedovBlast_3D/plotSolution.py +++ b/examples/HydroTests/SedovBlast_3D/plotSolution.py @@ -86,6 +86,18 @@ S = sim["/PartType0/Entropy"][:] P = sim["/PartType0/Pressure"][:] rho = sim["/PartType0/Density"][:] +try: + diffusion = sim["/PartType0/Diffusion"][:] + plot_diffusion = True +except: + plot_diffusion = False + +try: + viscosity = sim["/PartType0/Viscosity"][:] + plot_viscosity = True +except: + plot_viscosity = False + # Bin te data r_bin_edge = np.arange(0., 0.5, 0.01) r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1]) @@ -105,6 +117,16 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin**2) S_sigma_bin = np.sqrt(S2_bin - S_bin**2) u_sigma_bin = np.sqrt(u2_bin - u_bin**2) +if plot_diffusion: + alpha_diff_bin,_,_ = stats.binned_statistic(x, diffusion, statistic='mean', bins=r_bin_edge) + alpha2_diff_bin,_,_ = stats.binned_statistic(x, diffusion**2, statistic='mean', bins=r_bin_edge) + alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2) + +if plot_viscosity: + alpha_visc_bin,_,_ = stats.binned_statistic(x, viscosity, statistic='mean', bins=r_bin_edge) + alpha2_visc_bin,_,_ = stats.binned_statistic(x, viscosity**2, statistic='mean', bins=r_bin_edge) + alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2) + # Now, work our the solution.... @@ -274,13 +296,27 @@ ylim(-2, 22) # Entropy profile --------------------------------- subplot(235) -plot(r, S, '.', color='r', ms=0.5, alpha=0.2) -plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2) -errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) -xlabel("${\\rm{Radius}}~r$", labelpad=0) -ylabel("${\\rm{Entropy}}~S$", labelpad=0) xlim(0, 1.3 * r_shock) -ylim(-5, 50) +xlabel("${\\rm{Radius}}~r$", labelpad=0) + + +if plot_diffusion or plot_viscosity: + if plot_diffusion: + plot(x, diffusion, ".", color='r', ms=0.5, alpha=0.2) + errorbar(r_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion") + + if plot_viscosity: + plot(x, viscosity, ".", color='g', ms=0.5, alpha=0.2) + errorbar(r_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity") + + ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0) + legend() +else: + plot(r, S, '.', color='r', ms=0.5, alpha=0.2) + plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2) + errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) + ylabel("${\\rm{Entropy}}~S$", labelpad=0) + ylim(-5, 50) # Information ------------------------------------- subplot(236, frameon=False) diff --git a/examples/HydroTests/SodShock_3D/plotSolution.py b/examples/HydroTests/SodShock_3D/plotSolution.py index 6da7193bc..7cec619ef 100644 --- a/examples/HydroTests/SodShock_3D/plotSolution.py +++ b/examples/HydroTests/SodShock_3D/plotSolution.py @@ -84,6 +84,18 @@ S = sim["/PartType0/Entropy"][:] P = sim["/PartType0/Pressure"][:] rho = sim["/PartType0/Density"][:] +try: + diffusion = sim["/PartType0/Diffusion"][:] + plot_diffusion = True +except: + plot_diffusion = False + +try: + viscosity = sim["/PartType0/Viscosity"][:] + plot_viscosity = True +except: + plot_viscosity = False + x_min = -1. x_max = 1. x += x_min @@ -108,6 +120,15 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin**2) S_sigma_bin = np.sqrt(S2_bin - S_bin**2) u_sigma_bin = np.sqrt(u2_bin - u_bin**2) +if plot_diffusion: + alpha_diff_bin,_,_ = stats.binned_statistic(x, diffusion, statistic='mean', bins=x_bin_edge) + alpha2_diff_bin,_,_ = stats.binned_statistic(x, diffusion**2, statistic='mean', bins=x_bin_edge) + alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2) + +if plot_viscosity: + alpha_visc_bin,_,_ = stats.binned_statistic(x, viscosity, statistic='mean', bins=x_bin_edge) + alpha2_visc_bin,_,_ = stats.binned_statistic(x, viscosity**2, statistic='mean', bins=x_bin_edge) + alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2) # Analytic solution c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave @@ -282,13 +303,26 @@ ylim(0.8, 2.2) # Entropy profile --------------------------------- subplot(235) -plot(x, S, '.', color='r', ms=0.5, alpha=0.2) -plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2) -errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) -xlabel("${\\rm{Position}}~x$", labelpad=0) -ylabel("${\\rm{Entropy}}~S$", labelpad=0) xlim(-0.5, 0.5) -ylim(0.8, 3.8) +xlabel("${\\rm{Position}}~x$", labelpad=0) + +if plot_diffusion or plot_viscosity: + if plot_diffusion: + plot(x, diffusion, ".", color='r', ms=0.5, alpha=0.2) + errorbar(x_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion") + + if plot_viscosity: + plot(x, viscosity, ".", color='g', ms=0.5, alpha=0.2) + errorbar(x_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity") + + ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0) + legend() +else: + plot(x, S, '.', color='r', ms=0.5, alpha=0.2) + plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2) + errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) + ylabel("${\\rm{Entropy}}~S$", labelpad=0) + ylim(0.8, 3.8) # Information ------------------------------------- subplot(236, frameon=False) -- GitLab From 0eaff002fd3ae36d1b6d46ee49de0aac7f0b49c6 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Mon, 25 Mar 2019 20:19:30 +0000 Subject: [PATCH 18/54] Added steps to plotSolution --- .../CoolingRedshiftDependence/plotSolution.py | 34 +++++++++++++++---- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py index d8a4ddc0f..d7f4fe3b7 100644 --- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py +++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py @@ -45,7 +45,7 @@ def get_data(handle: float, n_snaps: int): densities = [] for snap in range(n_snaps): - data = load(f"{handle}_{snap:04d}.hdf5") + data = load(f"data/{handle}_{snap:04d}.hdf5") if snap == 0: t0 = data.metadata.t.to(Myr).value @@ -57,7 +57,17 @@ def get_data(handle: float, n_snaps: int): return times, temps, densities -def plot_single_data(handle: str, n_snaps: int, label: str, ax: plt.axes): +def get_n_steps(timesteps_filename: str) -> int: + """ + Reads the number of steps from the timesteps file. + """ + + data = np.genfromtxt(timesteps_filename) + + return int(data[-1][0]) + + +def plot_single_data(handle: str, n_snaps: int, label: str, ax: plt.axes, n_steps:int = 0): """ Takes the a single simulation and plots it on the axes. """ @@ -73,7 +83,7 @@ def plot_single_data(handle: str, n_snaps: int, label: str, ax: plt.axes): ax[1].plot( data[0], data[2], - label=label, + label=f"Steps: {n_steps}", marker="o", ms=2 ) @@ -82,17 +92,26 @@ def plot_single_data(handle: str, n_snaps: int, label: str, ax: plt.axes): return -def make_plot(handles, names, n_snaps=100): +def make_plot(handles, names, timestep_names, n_snaps=100): """ Makes the whole plot and returns the fig, ax objects. """ fig, ax = setup_axes() - for handle, name in zip(handles, names): - plot_single_data(handle, n_snaps, name, ax) + for handle, name, timestep_name in zip(handles, names, timestep_names): + try: + n_steps = get_n_steps(timestep_name) + except: + n_steps = "Unknown" + + try: + plot_single_data(handle, n_snaps, name, ax, n_steps=n_steps) + except: + pass ax[0].legend() + ax[1].legend() return fig, ax @@ -104,8 +123,9 @@ if __name__ == "__main__": handles = ["redshift_dependence_no_z", "redshift_dependence_low_z", "redshift_dependence_high_z"] names = ["No Cosmology", "Low Redshift ($z=0.01$)", "High Redshift ($z=1$)"] + timestep_names = ["timesteps_no.txt", "timesteps_low.txt", "timesteps_high.txt"] - fig, ax = make_plot(handles, names) + fig, ax = make_plot(handles, names, timestep_names) fig.savefig("redshift_dependence.png", dpi=300) -- GitLab From fb670c26086d5622201a7f593313ccd41db929ef Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Tue, 26 Mar 2019 11:11:23 +0000 Subject: [PATCH 19/54] Used more physically meaningful variable names in anarchy-pu sph --- src/hydro/AnarchyPU/hydro.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 0403608fe..998c00996 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -218,6 +218,7 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) { /* Compute the sound speed -- see theory section for justification */ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */ + const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho); return square_rooted; @@ -233,7 +234,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_soundspeed(const struct part *restrict p, const struct cosmology *cosmo) { - return cosmo->a_factor_sound_speed * p->force.soundspeed; + return cosmo->a_factor_sound_speed * hydro_get_comoving_soundspeed(p); } /** @@ -703,13 +704,14 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float h_physical = p->h * cosmo->a; const float v_sig_physical = p->viscosity.v_sig / a_fac_soundspeed; + const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo); const float div_v_physical = (p->viscosity.div_v - cosmo->H * hydro_dimension) * a2; const float div_v_prev_physical = (p->viscosity.div_v_previous_step - cosmo->H * hydro_dimension) * a2; const float dt_alpha_physical = dt_alpha / a2; - const float tau = h_physical / (2.f * v_sig_physical * hydro_props->viscosity.length); + const float sound_crossing_time_inverse = soundspeed_physical / h_physical; /* Construct time differential of div.v implicitly following the ANARCHY spec */ @@ -733,7 +735,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->viscosity.alpha = alpha_loc; } else { /* Integrate the alpha forward in time to decay back to alpha = 0 */ - const float alpha_dt = (alpha_loc - p->viscosity.alpha) / tau; + const float alpha_dt = (alpha_loc - p->viscosity.alpha) * sound_crossing_time_inverse * hydro_props->viscosity.length; /* Finally, we can update the actual value of the alpha */ p->viscosity.alpha += alpha_dt * dt_alpha_physical; @@ -748,6 +750,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Now for the diffusive alpha */ + const float diffusion_timescale_physical_inverse = v_sig_physical / h_physical; + const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ /* alpha_diff_dt is cosmology-less */ @@ -758,7 +762,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Decay term: not documented in Schaller+ 2015 but was present * in the original EAGLE code and in Schaye+ 2015 */ alpha_diff_dt -= - (p->diffusion.alpha - hydro_props->diffusion.alpha_min) / tau; + (p->diffusion.alpha - hydro_props->diffusion.alpha_min) * diffusion_timescale_physical_inverse; float new_diffusion_alpha = p->diffusion.alpha; new_diffusion_alpha += alpha_diff_dt * dt_alpha_physical; -- GitLab From 9a1a9b600d72d6106489561881620e4ab061c8a6 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Tue, 26 Mar 2019 12:52:00 +0000 Subject: [PATCH 20/54] Updated to use "correct" alpha --- src/hydro/AnarchyPU/hydro.h | 16 ++++++++++------ src/hydro/AnarchyPU/hydro_iact.h | 12 ++++++------ 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 998c00996..cba932b89 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -735,16 +735,20 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->viscosity.alpha = alpha_loc; } else { /* Integrate the alpha forward in time to decay back to alpha = 0 */ - const float alpha_dt = (alpha_loc - p->viscosity.alpha) * sound_crossing_time_inverse * hydro_props->viscosity.length; - - /* Finally, we can update the actual value of the alpha */ - p->viscosity.alpha += alpha_dt * dt_alpha_physical; + p->viscosity.alpha = + alpha_loc + (p->viscosity.alpha - alpha_loc) * + expf(-dt_alpha_physical * sound_crossing_time_inverse * + hydro_props->viscosity.length); } if (p->viscosity.alpha < hydro_props->viscosity.alpha_min) { p->viscosity.alpha = hydro_props->viscosity.alpha_min; } + if (p->viscosity.alpha > hydro_props->viscosity.alpha_max) { + p->viscosity.alpha = hydro_props->viscosity.alpha_max; + } + /* Set our old div_v to the one for the next loop */ p->viscosity.div_v_previous_step = p->viscosity.div_v; @@ -761,8 +765,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( hydro_props->diffusion.beta * p->h * p->diffusion.laplace_u / sqrt_u; /* Decay term: not documented in Schaller+ 2015 but was present * in the original EAGLE code and in Schaye+ 2015 */ - alpha_diff_dt -= - (p->diffusion.alpha - hydro_props->diffusion.alpha_min) * diffusion_timescale_physical_inverse; + alpha_diff_dt -= (p->diffusion.alpha - hydro_props->diffusion.alpha_min) * + diffusion_timescale_physical_inverse; float new_diffusion_alpha = p->diffusion.alpha; new_diffusion_alpha += alpha_diff_dt * dt_alpha_physical; diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index c214db3b0..870f30acd 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -204,10 +204,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient( (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; - const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx); + const float dv_dx_factor = min(0, dv_dx); const float new_v_sig = - pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor; + 0.5 * (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); @@ -257,10 +257,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; - const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx); + const float dv_dx_factor = min(0, dv_dx); const float new_v_sig = - pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor; + 0.5 * (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); @@ -353,7 +353,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( const float rho_ij = rhoi + rhoj; const float alpha = pi->viscosity.alpha + pj->viscosity.alpha; const float visc = - -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; + -1.f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; /* Convolve with the kernel */ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; @@ -484,7 +484,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( const float rho_ij = rhoi + rhoj; const float alpha = pi->viscosity.alpha + pj->viscosity.alpha; const float visc = - -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; + -1.f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; /* Convolve with the kernel */ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; -- GitLab From 75ba2e13c06960b5d0eb239c0b22b27304c037eb Mon Sep 17 00:00:00 2001 From: JBorrow Date: Wed, 27 Mar 2019 16:02:35 +0000 Subject: [PATCH 21/54] Added viscosity coefficients to plotting Gresho --- .../GreshoVortex_3D/plotSolution.py | 53 ++++++++++++++++--- 1 file changed, 45 insertions(+), 8 deletions(-) diff --git a/examples/HydroTests/GreshoVortex_3D/plotSolution.py b/examples/HydroTests/GreshoVortex_3D/plotSolution.py index 8fddf148b..545440c99 100644 --- a/examples/HydroTests/GreshoVortex_3D/plotSolution.py +++ b/examples/HydroTests/GreshoVortex_3D/plotSolution.py @@ -108,6 +108,18 @@ u = sim["/PartType0/InternalEnergy"][:] S = sim["/PartType0/Entropy"][:] P = sim["/PartType0/Pressure"][:] +try: + diffusion = sim["/PartType0/Diffusion"][:] + plot_diffusion = True +except: + plot_diffusion = False + +try: + viscosity = sim["/PartType0/Viscosity"][:] + plot_viscosity = True +except: + plot_viscosity = False + # Bin te data r_bin_edge = np.arange(0., 1., 0.02) r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1]) @@ -127,6 +139,15 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin**2) S_sigma_bin = np.sqrt(S2_bin - S_bin**2) u_sigma_bin = np.sqrt(u2_bin - u_bin**2) +if plot_diffusion: + alpha_diff_bin,_,_ = stats.binned_statistic(r, diffusion, statistic='mean', bins=r_bin_edge) + alpha2_diff_bin,_,_ = stats.binned_statistic(r, diffusion**2, statistic='mean', bins=r_bin_edge) + alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2) + +if plot_viscosity: + alpha_visc_bin,_,_ = stats.binned_statistic(r, viscosity, statistic='mean', bins=r_bin_edge) + alpha2_visc_bin,_,_ = stats.binned_statistic(r, viscosity**2, statistic='mean', bins=r_bin_edge) + alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2) # Plot the interesting quantities figure() @@ -188,16 +209,32 @@ ylim(7.3, 9.1) # Radial entropy profile -------------------------------- subplot(235) - -plot(r, S, '.', color='r', ms=0.5) -plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2) -errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) -plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2) -plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2) xlabel("${\\rm{Radius}}~r$", labelpad=0) -ylabel("${\\rm{Entropy}}~S$", labelpad=0) + xlim(0, R_max) -ylim(4.9 + P0, P0 + 6.1) + +xlabel("${\\rm{Radius}}~r$", labelpad=0) + + +if plot_diffusion or plot_viscosity: + if plot_diffusion: + plot(r, diffusion, ".", color='r', ms=0.5, alpha=0.2) + errorbar(r_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion") + + if plot_viscosity: + plot(r, viscosity, ".", color='g', ms=0.5, alpha=0.2) + errorbar(r_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity") + + ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0) + legend() +else: + plot(r, S, '.', color='r', ms=0.5) + plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2) + errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) + plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2) + plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2) + ylabel("${\\rm{Entropy}}~S$", labelpad=0) + ylim(4.9 + P0, P0 + 6.1) # Image -------------------------------------------------- #subplot(234) -- GitLab From 48f4f09f9d36b3cd49d74e1edce15e1c54238c47 Mon Sep 17 00:00:00 2001 From: JBorrow Date: Wed, 27 Mar 2019 16:03:12 +0000 Subject: [PATCH 22/54] Fixed bug with Sedov plotting --- .../HydroTests/SedovBlast_3D/plotSolution.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py index a87e3cd24..b0f2e0844 100644 --- a/examples/HydroTests/SedovBlast_3D/plotSolution.py +++ b/examples/HydroTests/SedovBlast_3D/plotSolution.py @@ -24,7 +24,7 @@ # Parameters rho_0 = 1. # Background Density P_0 = 1.e-6 # Background Pressure -E_0 = 1. # Energy of the explosion +E_0 = 1.0 # Energy of the explosion gas_gamma = 5./3. # Gas polytropic index @@ -118,13 +118,13 @@ S_sigma_bin = np.sqrt(S2_bin - S_bin**2) u_sigma_bin = np.sqrt(u2_bin - u_bin**2) if plot_diffusion: - alpha_diff_bin,_,_ = stats.binned_statistic(x, diffusion, statistic='mean', bins=r_bin_edge) - alpha2_diff_bin,_,_ = stats.binned_statistic(x, diffusion**2, statistic='mean', bins=r_bin_edge) + alpha_diff_bin,_,_ = stats.binned_statistic(r, diffusion, statistic='mean', bins=r_bin_edge) + alpha2_diff_bin,_,_ = stats.binned_statistic(r, diffusion**2, statistic='mean', bins=r_bin_edge) alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2) if plot_viscosity: - alpha_visc_bin,_,_ = stats.binned_statistic(x, viscosity, statistic='mean', bins=r_bin_edge) - alpha2_visc_bin,_,_ = stats.binned_statistic(x, viscosity**2, statistic='mean', bins=r_bin_edge) + alpha_visc_bin,_,_ = stats.binned_statistic(r, viscosity, statistic='mean', bins=r_bin_edge) + alpha2_visc_bin,_,_ = stats.binned_statistic(r, viscosity**2, statistic='mean', bins=r_bin_edge) alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2) @@ -296,17 +296,16 @@ ylim(-2, 22) # Entropy profile --------------------------------- subplot(235) -xlim(0, 1.3 * r_shock) xlabel("${\\rm{Radius}}~r$", labelpad=0) if plot_diffusion or plot_viscosity: if plot_diffusion: - plot(x, diffusion, ".", color='r', ms=0.5, alpha=0.2) + plot(r, diffusion, ".", color='r', ms=0.5, alpha=0.2) errorbar(r_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion") if plot_viscosity: - plot(x, viscosity, ".", color='g', ms=0.5, alpha=0.2) + plot(r, viscosity, ".", color='g', ms=0.5, alpha=0.2) errorbar(r_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity") ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0) @@ -318,6 +317,7 @@ else: ylabel("${\\rm{Entropy}}~S$", labelpad=0) ylim(-5, 50) +xlim(0, 1.3 * r_shock) # Information ------------------------------------- subplot(236, frameon=False) -- GitLab From 6b4d59acea7b6e3e4468e875ef0dc8642aba8913 Mon Sep 17 00:00:00 2001 From: JBorrow Date: Wed, 27 Mar 2019 16:04:58 +0000 Subject: [PATCH 23/54] Minor bugfixes in ANARCHY scheme --- src/hydro/AnarchyPU/hydro.h | 8 +++++--- src/hydro/AnarchyPU/hydro_iact.h | 24 ++++++++++++++---------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index cba932b89..1fc4e9d8e 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -218,7 +218,7 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) { /* Compute the sound speed -- see theory section for justification */ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */ - + const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho); return square_rooted; @@ -702,7 +702,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* We integrate the whole AV scheme in physical units as cancelling the cosmology terms is just a nightmare, especially for non-standard gas_gamma. */ - const float h_physical = p->h * cosmo->a; + /* We use in this function that h is the radius of support */ + const float h_physical = p->h * cosmo->a * kernel_gamma; const float v_sig_physical = p->viscosity.v_sig / a_fac_soundspeed; const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo); const float div_v_physical = @@ -723,7 +724,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Construct the source term for the AV; if shock detected this is _positive_ * as div_v_dt should be _negative_ before the shock hits */ const float S = h_physical * h_physical * max(0.f, -1.f * div_v_dt); - const float v_sig_square = v_sig_physical * v_sig_physical; + /* 0.25 factor comes from our definition of v_sig */ + const float v_sig_square = 0.25 * v_sig_physical * v_sig_physical; /* Calculate the current appropriate value of the AV based on the above */ /* Note this is v_sig_physical squared, not comoving */ diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index 870f30acd..a7fc22b79 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -200,14 +200,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient( /* We need to construct the maximal signal velocity between our particle * and all of it's neighbours */ + const float r = sqrtf(r2); + const float r_inv = 1.f / r; + const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; - const float dv_dx_factor = min(0, dv_dx); + const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx); const float new_v_sig = - 0.5 * (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor; + pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor * r_inv; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); @@ -217,14 +220,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient( /* Need to get some kernel values F_ij = wi_dx */ float wi, wi_dx, wj, wj_dx; - const float r = sqrtf(r2); const float ui = r / hi; const float uj = r / hj; kernel_deval(ui, &wi, &wi_dx); kernel_deval(uj, &wj, &wj_dx); - const float delta_u_factor = (pi->u - pj->u) / r; + const float delta_u_factor = (pi->u - pj->u) * r_inv; pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho; pj->diffusion.laplace_u -= pi->mass * delta_u_factor * wj_dx / pi->rho; } @@ -253,14 +255,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( /* We need to construct the maximal signal velocity between our particle * and all of it's neighbours */ + const float r = sqrtf(r2); + const float r_inv = 1.f / r; + const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; - const float dv_dx_factor = min(0, dv_dx); + const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx); const float new_v_sig = - 0.5 * (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor; + pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor * r_inv; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); @@ -269,12 +274,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( /* Need to get some kernel values F_ij = wi_dx */ float wi, wi_dx; - const float r = sqrtf(r2); const float ui = r / hi; kernel_deval(ui, &wi, &wi_dx); - const float delta_u_factor = (pi->u - pj->u) / r; + const float delta_u_factor = (pi->u - pj->u) * r_inv; pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho; } @@ -385,7 +389,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( wj_dr * dvdr * r_inv; /* Viscosity term */ - const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble; + const float visc_du_term = 0.5 * visc_acc_term * dvdr_Hubble; /* Diffusion term */ const float v_diff = @@ -509,7 +513,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( wi_dr * dvdr * r_inv; /* Viscosity term */ - const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble; + const float visc_du_term = 0.5 * visc_acc_term * dvdr_Hubble; /* Diffusion term */ const float v_diff = -- GitLab From 469fb7b5126e7b736c167df9fb453682f3dca5e0 Mon Sep 17 00:00:00 2001 From: JBorrow Date: Wed, 27 Mar 2019 16:05:48 +0000 Subject: [PATCH 24/54] Updated ANARCHY defaults --- src/hydro_properties.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hydro_properties.c b/src/hydro_properties.c index ff578aec1..9edd72f77 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -50,7 +50,7 @@ /* This nasty #ifdef is only temporary until we separate the viscosity * and hydro components. If it is not removed by July 2019, shout at JB. */ #define hydro_props_default_viscosity_alpha_min \ - 0.01f /* values taken from Schaller+ 2015 */ + 0.1f /* values NOT the same as Schaller+ 2015 */ #define hydro_props_default_viscosity_alpha_max \ 2.0f /* values taken from Schaller+ 2015 */ #define hydro_props_default_viscosity_length \ -- GitLab From 77ba30bd19c9765baef03b66b6ca7a93cdb917eb Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 27 Mar 2019 17:47:43 +0000 Subject: [PATCH 25/54] Attempting to fix soundspeed --- src/hydro/AnarchyPU/hydro.h | 3 +-- src/hydro/AnarchyPU/hydro_iact.h | 24 ++++++++++++++++++++---- src/hydro_properties.c | 3 +++ 3 files changed, 24 insertions(+), 6 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 1fc4e9d8e..dbd78eee5 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -697,14 +697,13 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Timescale for decay */ const float a2 = cosmo->a * cosmo->a; - const float a_fac_soundspeed = cosmo->a_factor_sound_speed; /* We integrate the whole AV scheme in physical units as cancelling the cosmology terms is just a nightmare, especially for non-standard gas_gamma. */ /* We use in this function that h is the radius of support */ const float h_physical = p->h * cosmo->a * kernel_gamma; - const float v_sig_physical = p->viscosity.v_sig / a_fac_soundspeed; + const float v_sig_physical = p->viscosity.v_sig / cosmo->a_factor_sound_speed; const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo); const float div_v_physical = (p->viscosity.div_v - cosmo->H * hydro_dimension) * a2; diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index a7fc22b79..887783be4 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -203,14 +203,22 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient( const float r = sqrtf(r2); const float r_inv = 1.f / r; + /* Cosmology terms for the signal velocity */ + const float fac_mu = pow_three_gamma_minus_five_over_two(a); + const float a2_Hubble = a * a * H; + const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; - const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx); + /* Add Hubble flow */ + const float dv_dx_Hubble = dv_dx + a2_Hubble * r2; + + const float dv_dx_factor = min(0, const_viscosity_beta * fac_mu * dv_dx_Hubble); + /* Store this as a _physical_ variable because of weird combinations of a-factors */ const float new_v_sig = - pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor * r_inv; + (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor * r_inv; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); @@ -257,15 +265,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( const float r = sqrtf(r2); const float r_inv = 1.f / r; + + /* Cosmology terms for the signal velocity */ + const float fac_mu = pow_three_gamma_minus_five_over_two(a); + const float a2_Hubble = a * a * H; const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; - const float dv_dx_factor = min(0, const_viscosity_beta * dv_dx); + /* Add Hubble flow */ + const float dv_dx_Hubble = dv_dx + a2_Hubble * r2; + + const float dv_dx_factor = min(0, const_viscosity_beta * fac_mu * dv_dx_Hubble); + /* Store this as a _physical_ variable because of weird combinations of a-factors */ const float new_v_sig = - pi->force.soundspeed + pj->force.soundspeed - dv_dx_factor * r_inv; + (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor * r_inv; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); diff --git a/src/hydro_properties.c b/src/hydro_properties.c index 9edd72f77..7d2b03be7 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -49,6 +49,9 @@ #ifdef ANARCHY_PU_SPH /* This nasty #ifdef is only temporary until we separate the viscosity * and hydro components. If it is not removed by July 2019, shout at JB. */ +#undef hydro_props_default_viscosity_alpha +#define hydro_props_default_viscosity_alpha \ + 0.1f /* Use a very low initial AV paramater for hydrodynamics tests */ #define hydro_props_default_viscosity_alpha_min \ 0.1f /* values NOT the same as Schaller+ 2015 */ #define hydro_props_default_viscosity_alpha_max \ -- GitLab From e250e7988baffb0555f93b13f4f53afbdefa5f85 Mon Sep 17 00:00:00 2001 From: JBorrow Date: Thu, 28 Mar 2019 21:34:31 +0000 Subject: [PATCH 26/54] Slowly converging on cosmology fixes --- src/hydro/AnarchyPU/hydro.h | 34 +++++++++----------------------- src/hydro/AnarchyPU/hydro_iact.h | 30 +++++++++++++++++----------- 2 files changed, 27 insertions(+), 37 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index dbd78eee5..26845d928 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -694,22 +694,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Here we need to update the artificial viscosity */ - /* Timescale for decay */ - - const float a2 = cosmo->a * cosmo->a; - - /* We integrate the whole AV scheme in physical units as cancelling the cosmology terms - is just a nightmare, especially for non-standard gas_gamma. */ - /* We use in this function that h is the radius of support */ const float h_physical = p->h * cosmo->a * kernel_gamma; - const float v_sig_physical = p->viscosity.v_sig / cosmo->a_factor_sound_speed; + const float v_sig_physical = p->viscosity.v_sig * cosmo->a_factor_sound_speed; const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo); - const float div_v_physical = - (p->viscosity.div_v - cosmo->H * hydro_dimension) * a2; - const float div_v_prev_physical = - (p->viscosity.div_v_previous_step - cosmo->H * hydro_dimension) * a2; - const float dt_alpha_physical = dt_alpha / a2; const float sound_crossing_time_inverse = soundspeed_physical / h_physical; @@ -718,16 +706,17 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( float div_v_dt = dt_alpha == 0.f ? 0.f - : (div_v_physical - div_v_prev_physical) / dt_alpha_physical; + : (p->viscosity.div_v - p->viscosity.div_v_previous_step) / dt_alpha; /* Construct the source term for the AV; if shock detected this is _positive_ * as div_v_dt should be _negative_ before the shock hits */ const float S = h_physical * h_physical * max(0.f, -1.f * div_v_dt); - /* 0.25 factor comes from our definition of v_sig */ + /* 0.25 factor comes from our definition of v_sig (sum of soundspeeds rather + * than mean). */ + /* Note this is v_sig_physical squared, not comoving */ const float v_sig_square = 0.25 * v_sig_physical * v_sig_physical; /* Calculate the current appropriate value of the AV based on the above */ - /* Note this is v_sig_physical squared, not comoving */ const float alpha_loc = hydro_props->viscosity.alpha_max * S / (v_sig_square + S); @@ -735,10 +724,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Reset the value of alpha to the appropriate value */ p->viscosity.alpha = alpha_loc; } else { - /* Integrate the alpha forward in time to decay back to alpha = 0 */ + /* Integrate the alpha forward in time to decay back to alpha = alpha_loc */ p->viscosity.alpha = alpha_loc + (p->viscosity.alpha - alpha_loc) * - expf(-dt_alpha_physical * sound_crossing_time_inverse * + expf(-dt_alpha * sound_crossing_time_inverse * hydro_props->viscosity.length); } @@ -746,10 +735,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->viscosity.alpha = hydro_props->viscosity.alpha_min; } - if (p->viscosity.alpha > hydro_props->viscosity.alpha_max) { - p->viscosity.alpha = hydro_props->viscosity.alpha_max; - } - /* Set our old div_v to the one for the next loop */ p->viscosity.div_v_previous_step = p->viscosity.div_v; @@ -759,18 +744,17 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ - /* alpha_diff_dt is cosmology-less */ /* Evolution term: following Schaller+ 2015. This is actually cosmology-less and so requires no conversion to physical. */ float alpha_diff_dt = - hydro_props->diffusion.beta * p->h * p->diffusion.laplace_u / sqrt_u; + hydro_props->diffusion.beta * h_physical * cosmo->a_factor_sound_speed * p->diffusion.laplace_u / sqrt_u; /* Decay term: not documented in Schaller+ 2015 but was present * in the original EAGLE code and in Schaye+ 2015 */ alpha_diff_dt -= (p->diffusion.alpha - hydro_props->diffusion.alpha_min) * diffusion_timescale_physical_inverse; float new_diffusion_alpha = p->diffusion.alpha; - new_diffusion_alpha += alpha_diff_dt * dt_alpha_physical; + new_diffusion_alpha += alpha_diff_dt * dt_alpha / (cosmo->a * cosmo->a); /* Consistency checks to ensure min < alpha < max */ if (new_diffusion_alpha > hydro_props->diffusion.alpha_max) { diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index 887783be4..6d49b8d40 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -202,23 +202,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient( const float r = sqrtf(r2); const float r_inv = 1.f / r; + const float ci = pi->force.soundspeed; + const float cj = pi->force.soundspeed; /* Cosmology terms for the signal velocity */ const float fac_mu = pow_three_gamma_minus_five_over_two(a); const float a2_Hubble = a * a * H; - const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] + + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; /* Add Hubble flow */ - const float dv_dx_Hubble = dv_dx + a2_Hubble * r2; - const float dv_dx_factor = min(0, const_viscosity_beta * fac_mu * dv_dx_Hubble); + const float dvdr_Hubble = dvdr + a2_Hubble * r2; + /* Are the particles moving towards each others ? */ + const float omega_ij = min(dvdr_Hubble, 0.f); + const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ - /* Store this as a _physical_ variable because of weird combinations of a-factors */ - const float new_v_sig = - (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor * r_inv; + /* Signal velocity */ + const float new_v_sig = ci + cj - const_viscosity_beta * mu_ij; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); @@ -265,23 +268,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( const float r = sqrtf(r2); const float r_inv = 1.f / r; + const float ci = pi->force.soundspeed; + const float cj = pi->force.soundspeed; /* Cosmology terms for the signal velocity */ const float fac_mu = pow_three_gamma_minus_five_over_two(a); const float a2_Hubble = a * a * H; - const float dv_dx = (pi->v[0] - pj->v[0]) * dx[0] + + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + (pi->v[2] - pj->v[2]) * dx[2]; /* Add Hubble flow */ - const float dv_dx_Hubble = dv_dx + a2_Hubble * r2; - const float dv_dx_factor = min(0, const_viscosity_beta * fac_mu * dv_dx_Hubble); + const float dvdr_Hubble = dvdr + a2_Hubble * r2; + /* Are the particles moving towards each others ? */ + const float omega_ij = min(dvdr_Hubble, 0.f); + const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ - /* Store this as a _physical_ variable because of weird combinations of a-factors */ - const float new_v_sig = - (pi->force.soundspeed + pj->force.soundspeed) - dv_dx_factor * r_inv; + /* Signal velocity */ + const float new_v_sig = ci + cj - const_viscosity_beta * mu_ij; /* Update if we need to */ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig); -- GitLab From 8b6871f0d6c0e71ac6485533d87eabed0758e526 Mon Sep 17 00:00:00 2001 From: JBorrow Date: Fri, 29 Mar 2019 17:15:37 +0000 Subject: [PATCH 27/54] Fixes to the diffusion cosmology dependence. --- src/hydro/AnarchyPU/hydro.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 26845d928..d4a86085e 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -747,14 +747,14 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Evolution term: following Schaller+ 2015. This is actually cosmology-less and so requires no conversion to physical. */ float alpha_diff_dt = - hydro_props->diffusion.beta * h_physical * cosmo->a_factor_sound_speed * p->diffusion.laplace_u / sqrt_u; + hydro_props->diffusion.beta * h_physical * p->diffusion.laplace_u * cosmo->a_factor_sound_speed / (sqrt_u * cosmo->a * cosmo->a); /* Decay term: not documented in Schaller+ 2015 but was present * in the original EAGLE code and in Schaye+ 2015 */ alpha_diff_dt -= (p->diffusion.alpha - hydro_props->diffusion.alpha_min) * diffusion_timescale_physical_inverse; float new_diffusion_alpha = p->diffusion.alpha; - new_diffusion_alpha += alpha_diff_dt * dt_alpha / (cosmo->a * cosmo->a); + new_diffusion_alpha += alpha_diff_dt * dt_alpha; /* Consistency checks to ensure min < alpha < max */ if (new_diffusion_alpha > hydro_props->diffusion.alpha_max) { -- GitLab From 356bd46925ca54148a70995e9e8764747a2b1d6c Mon Sep 17 00:00:00 2001 From: JBorrow Date: Fri, 29 Mar 2019 17:16:03 +0000 Subject: [PATCH 28/54] Plotted 100x diffusion coefficient to fit both nicely on page. --- .../ComovingSodShock_3D/plotSolution.py | 47 ++++++++++++++++--- 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py b/examples/Cosmology/ComovingSodShock_3D/plotSolution.py index f05a385e8..d85f4be9a 100644 --- a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py +++ b/examples/Cosmology/ComovingSodShock_3D/plotSolution.py @@ -88,6 +88,18 @@ S = sim["/PartType0/Entropy"][:] P = sim["/PartType0/Pressure"][:] rho = sim["/PartType0/Density"][:] +try: + diffusion = sim["/PartType0/Diffusion"][:] + plot_diffusion = True +except: + plot_diffusion = False + +try: + viscosity = sim["/PartType0/Viscosity"][:] + plot_viscosity = True +except: + plot_viscosity = False + x_min = -1. x_max = 1. x += x_min @@ -112,6 +124,15 @@ P_sigma_bin = np.sqrt(P2_bin - P_bin**2) S_sigma_bin = np.sqrt(S2_bin - S_bin**2) u_sigma_bin = np.sqrt(u2_bin - u_bin**2) +if plot_diffusion: + alpha_diff_bin,_,_ = stats.binned_statistic(x, diffusion, statistic='mean', bins=x_bin_edge) + alpha2_diff_bin,_,_ = stats.binned_statistic(x, diffusion**2, statistic='mean', bins=x_bin_edge) + alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2) + +if plot_viscosity: + alpha_visc_bin,_,_ = stats.binned_statistic(x, viscosity, statistic='mean', bins=x_bin_edge) + alpha2_visc_bin,_,_ = stats.binned_statistic(x, viscosity**2, statistic='mean', bins=x_bin_edge) + alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2) # Analytic solution c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave @@ -285,13 +306,26 @@ ylim(0.8, 2.2) # Entropy profile --------------------------------- subplot(235) -plot(x, S, '.', color='r', ms=0.5, alpha=0.2) -plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2) -errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) -xlabel("${\\rm{Position}}~x$", labelpad=0) -ylabel("${\\rm{Entropy}}~S$", labelpad=0) xlim(-0.5, 0.5) -ylim(0.8, 3.8) +xlabel("${\\rm{Position}}~x$", labelpad=0) + +if plot_diffusion or plot_viscosity: + if plot_diffusion: + plot(x, diffusion * 100, ".", color='r', ms=0.5, alpha=0.2) + errorbar(x_bin, alpha_diff_bin * 100, yerr=alpha_diff_sigma_bin * 100, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion (100x)") + + if plot_viscosity: + plot(x, viscosity, ".", color='g', ms=0.5, alpha=0.2) + errorbar(x_bin, alpha_visc_bin, yerr=alpha_visc_sigma_bin, fmt=".", ms=8.0, color='y', lw=1.2, label="Viscosity") + + ylabel("${\\rm{Rate~Coefficient}}~\\alpha$", labelpad=0) + legend() +else: + plot(x, S, '.', color='r', ms=0.5, alpha=0.2) + plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2) + errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2) + ylabel("${\\rm{Entropy}}~S$", labelpad=0) + ylim(0.8, 3.8) # Information ------------------------------------- subplot(236, frameon=False) @@ -312,5 +346,4 @@ ylim(0, 1) xticks([]) yticks([]) -tight_layout() savefig("SodShock.png", dpi=200) -- GitLab From f80db4ad52a7ddcf12117756b8d2d36102234c5d Mon Sep 17 00:00:00 2001 From: JBorrow Date: Fri, 29 Mar 2019 17:16:31 +0000 Subject: [PATCH 29/54] Fixed units in comoving sodshock example --- examples/Cosmology/ComovingSodShock_3D/sodShock.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/Cosmology/ComovingSodShock_3D/sodShock.yml b/examples/Cosmology/ComovingSodShock_3D/sodShock.yml index f7cf51ce0..b5b2ff8fd 100644 --- a/examples/Cosmology/ComovingSodShock_3D/sodShock.yml +++ b/examples/Cosmology/ComovingSodShock_3D/sodShock.yml @@ -1,8 +1,8 @@ # Define the system of units to use internally. InternalUnitSystem: - UnitMass_in_cgs: 2.94e45 # Grams + UnitMass_in_cgs: 2.94e55 # Grams UnitLength_in_cgs: 3.086e18 # pc - UnitVelocity_in_cgs: 1e10 # km per s + UnitVelocity_in_cgs: 1 # km per s UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin -- GitLab From 2a16d47bd10d48381353f6ea47975cf5e9d5da7d Mon Sep 17 00:00:00 2001 From: JBorrow Date: Fri, 29 Mar 2019 17:17:18 +0000 Subject: [PATCH 30/54] Minor change to sodshock plotting script --- examples/HydroTests/SodShock_3D/plotSolution.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/HydroTests/SodShock_3D/plotSolution.py b/examples/HydroTests/SodShock_3D/plotSolution.py index 7cec619ef..69b2fe488 100644 --- a/examples/HydroTests/SodShock_3D/plotSolution.py +++ b/examples/HydroTests/SodShock_3D/plotSolution.py @@ -308,8 +308,8 @@ xlabel("${\\rm{Position}}~x$", labelpad=0) if plot_diffusion or plot_viscosity: if plot_diffusion: - plot(x, diffusion, ".", color='r', ms=0.5, alpha=0.2) - errorbar(x_bin, alpha_diff_bin, yerr=alpha_diff_sigma_bin, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion") + plot(x, diffusion * 100, ".", color='r', ms=0.5, alpha=0.2) + errorbar(x_bin, alpha_diff_bin * 100, yerr=alpha_diff_sigma_bin * 100, fmt=".", ms=8.0, color='b', lw=1.2, label="Diffusion (100x)") if plot_viscosity: plot(x, viscosity, ".", color='g', ms=0.5, alpha=0.2) -- GitLab From ac4f739121f3dfed0432aaf9bf9366a1dc892be1 Mon Sep 17 00:00:00 2001 From: JBorrow Date: Fri, 29 Mar 2019 17:17:37 +0000 Subject: [PATCH 31/54] Made UniformBox_3D work with python3 --- examples/HydroTests/UniformBox_3D/makeICbig.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/HydroTests/UniformBox_3D/makeICbig.py b/examples/HydroTests/UniformBox_3D/makeICbig.py index bd5cf627f..8841b6016 100644 --- a/examples/HydroTests/UniformBox_3D/makeICbig.py +++ b/examples/HydroTests/UniformBox_3D/makeICbig.py @@ -45,7 +45,7 @@ internalEnergy = P / ((gamma - 1.)*rho) n_iterations = numPart / N remainder = numPart % N -print "Writing", numPart, "in", n_iterations, "iterations of size", N, "and a remainder of", remainder +print("Writing", numPart, "in", n_iterations, "iterations of size", N, "and a remainder of", remainder) #--------------------------------------------------- @@ -55,8 +55,8 @@ file = h5py.File(fileName, 'w') # Header grp = file.create_group("/Header") grp.attrs["BoxSize"] = boxSize -grp.attrs["NumPart_Total"] = [numPart % (long(1)<<32), 0, 0, 0, 0, 0] -grp.attrs["NumPart_Total_HighWord"] = [numPart / (long(1)<<32), 0, 0, 0, 0, 0] +grp.attrs["NumPart_Total"] = [numPart % (int(1)<<32), 0, 0, 0, 0, 0] +grp.attrs["NumPart_Total_HighWord"] = [numPart / (int(1)<<32), 0, 0, 0, 0, 0] grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0] grp.attrs["Time"] = 0.0 grp.attrs["NumFilesPerSnapshot"] = 1 @@ -120,7 +120,7 @@ for n in range(n_iterations): coords = zeros((1,3)) offset += N - print "Done", offset,"/", numPart, "(%.1f %%)"%(100*(float)(offset)/numPart) + print("Done", offset,"/", numPart, "(%.1f %%)"%(100*(float)(offset)/numPart)) # And now, the remainder v = zeros((remainder, 3)) @@ -150,7 +150,7 @@ coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L) coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L) ods_x[offset:offset+remainder,:] = coords -print "Done", offset+remainder,"/", numPart +print("Done", offset+remainder,"/", numPart) -- GitLab From e7f02f30472eb946f03ec761f92a6fc23726f7b3 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Fri, 29 Mar 2019 17:20:08 +0000 Subject: [PATCH 32/54] Code formatting --- src/cooling/const_lambda/cooling.h | 7 +++---- src/hydro/AnarchyPU/hydro.h | 17 ++++++++++------- src/hydro/AnarchyPU/hydro_iact.h | 14 +++++++------- 3 files changed, 20 insertions(+), 18 deletions(-) diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 0d15c297f..f5816b53e 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -352,10 +352,9 @@ static INLINE void cooling_print_backend( "* " "cm^3]", cooling->lambda_nH2_cgs); - - message( - "Lambda/n_H^2=%g [internal units]", cooling->lambda_nH2_cgs * cooling->conv_factor_energy_rate_from_cgs - ); + + message("Lambda/n_H^2=%g [internal units]", + cooling->lambda_nH2_cgs * cooling->conv_factor_energy_rate_from_cgs); if (cooling->cooling_tstep_mult == FLT_MAX) message("Cooling function time-step size is unlimited"); diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index d4a86085e..683917143 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -218,7 +218,7 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) { /* Compute the sound speed -- see theory section for justification */ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */ - + const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho); return square_rooted; @@ -701,7 +701,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float sound_crossing_time_inverse = soundspeed_physical / h_physical; - /* Construct time differential of div.v implicitly following the ANARCHY spec */ + /* Construct time differential of div.v implicitly following the ANARCHY spec + */ float div_v_dt = dt_alpha == 0.f @@ -740,14 +741,16 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Now for the diffusive alpha */ - const float diffusion_timescale_physical_inverse = v_sig_physical / h_physical; + const float diffusion_timescale_physical_inverse = + v_sig_physical / h_physical; const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ - /* Evolution term: following Schaller+ 2015. This is actually cosmology-less and - so requires no conversion to physical. */ - float alpha_diff_dt = - hydro_props->diffusion.beta * h_physical * p->diffusion.laplace_u * cosmo->a_factor_sound_speed / (sqrt_u * cosmo->a * cosmo->a); + /* Evolution term: following Schaller+ 2015. This is actually cosmology-less + and so requires no conversion to physical. */ + float alpha_diff_dt = hydro_props->diffusion.beta * h_physical * + p->diffusion.laplace_u * cosmo->a_factor_sound_speed / + (sqrt_u * cosmo->a * cosmo->a); /* Decay term: not documented in Schaller+ 2015 but was present * in the original EAGLE code and in Schaye+ 2015 */ alpha_diff_dt -= (p->diffusion.alpha - hydro_props->diffusion.alpha_min) * diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index 6d49b8d40..d04cd10f3 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -205,13 +205,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient( const float ci = pi->force.soundspeed; const float cj = pi->force.soundspeed; - /* Cosmology terms for the signal velocity */ + /* Cosmology terms for the signal velocity */ const float fac_mu = pow_three_gamma_minus_five_over_two(a); const float a2_Hubble = a * a * H; const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + - (pi->v[1] - pj->v[1]) * dx[1] + - (pi->v[2] - pj->v[2]) * dx[2]; + (pi->v[1] - pj->v[1]) * dx[1] + + (pi->v[2] - pj->v[2]) * dx[2]; /* Add Hubble flow */ @@ -270,14 +270,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( const float r_inv = 1.f / r; const float ci = pi->force.soundspeed; const float cj = pi->force.soundspeed; - - /* Cosmology terms for the signal velocity */ + + /* Cosmology terms for the signal velocity */ const float fac_mu = pow_three_gamma_minus_five_over_two(a); const float a2_Hubble = a * a * H; const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + - (pi->v[1] - pj->v[1]) * dx[1] + - (pi->v[2] - pj->v[2]) * dx[2]; + (pi->v[1] - pj->v[1]) * dx[1] + + (pi->v[2] - pj->v[2]) * dx[2]; /* Add Hubble flow */ -- GitLab From 2510cdb20592bd62a6a49324b84cfb41cc6b1eb4 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Fri, 29 Mar 2019 17:21:11 +0000 Subject: [PATCH 33/54] Added comments explaining cosmology --- src/hydro/AnarchyPU/hydro.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 683917143..7f6428c66 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -746,8 +746,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ - /* Evolution term: following Schaller+ 2015. This is actually cosmology-less - and so requires no conversion to physical. */ + /* Evolution term: following Schaller+ 2015. This is made up of several + cosmology factors: physical h, sound speed from u / sqrt(u), and the + 1 / a^2 coming from the laplace operator. */ float alpha_diff_dt = hydro_props->diffusion.beta * h_physical * p->diffusion.laplace_u * cosmo->a_factor_sound_speed / (sqrt_u * cosmo->a * cosmo->a); -- GitLab From 1835710d6d6fc1e4ba95ce561c0e946281c898fb Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Sun, 31 Mar 2019 09:57:45 +0100 Subject: [PATCH 34/54] Changed variable name from h -> kernel support --- src/hydro/AnarchyPU/hydro.h | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 7f6428c66..0ac52165e 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -695,11 +695,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Here we need to update the artificial viscosity */ /* We use in this function that h is the radius of support */ - const float h_physical = p->h * cosmo->a * kernel_gamma; + const float kernel_support_physical = p->h * cosmo->a * kernel_gamma; const float v_sig_physical = p->viscosity.v_sig * cosmo->a_factor_sound_speed; const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo); - const float sound_crossing_time_inverse = soundspeed_physical / h_physical; + const float sound_crossing_time_inverse = + soundspeed_physical / kernel_support_physical; /* Construct time differential of div.v implicitly following the ANARCHY spec */ @@ -711,7 +712,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Construct the source term for the AV; if shock detected this is _positive_ * as div_v_dt should be _negative_ before the shock hits */ - const float S = h_physical * h_physical * max(0.f, -1.f * div_v_dt); + const float S = kernel_support_physical * kernel_support_physical * + max(0.f, -1.f * div_v_dt); /* 0.25 factor comes from our definition of v_sig (sum of soundspeeds rather * than mean). */ /* Note this is v_sig_physical squared, not comoving */ @@ -742,14 +744,14 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Now for the diffusive alpha */ const float diffusion_timescale_physical_inverse = - v_sig_physical / h_physical; + v_sig_physical / kernel_support_physical; const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ /* Evolution term: following Schaller+ 2015. This is made up of several cosmology factors: physical h, sound speed from u / sqrt(u), and the 1 / a^2 coming from the laplace operator. */ - float alpha_diff_dt = hydro_props->diffusion.beta * h_physical * + float alpha_diff_dt = hydro_props->diffusion.beta * kernel_support_physical * p->diffusion.laplace_u * cosmo->a_factor_sound_speed / (sqrt_u * cosmo->a * cosmo->a); /* Decay term: not documented in Schaller+ 2015 but was present -- GitLab From 911a3fedaf1dd60ae26b1b3feaf1e60bc0bcf689 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Sun, 31 Mar 2019 10:00:59 +0100 Subject: [PATCH 35/54] Fixed incorrect soundspeed in gradient loop for anarchy --- src/hydro/AnarchyPU/hydro_iact.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index d04cd10f3..cc2128b03 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -203,7 +203,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient( const float r = sqrtf(r2); const float r_inv = 1.f / r; const float ci = pi->force.soundspeed; - const float cj = pi->force.soundspeed; + const float cj = pj->force.soundspeed; /* Cosmology terms for the signal velocity */ const float fac_mu = pow_three_gamma_minus_five_over_two(a); @@ -269,7 +269,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient( const float r = sqrtf(r2); const float r_inv = 1.f / r; const float ci = pi->force.soundspeed; - const float cj = pi->force.soundspeed; + const float cj = pj->force.soundspeed; /* Cosmology terms for the signal velocity */ const float fac_mu = pow_three_gamma_minus_five_over_two(a); @@ -369,7 +369,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ /* Compute sound speeds and signal velocity */ - const float v_sig = 0.5 * (pi->viscosity.v_sig + pj->viscosity.v_sig); + const float v_sig = 0.5f * (pi->viscosity.v_sig + pj->viscosity.v_sig); /* Balsara term */ const float balsara_i = pi->force.balsara; @@ -411,12 +411,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( wj_dr * dvdr * r_inv; /* Viscosity term */ - const float visc_du_term = 0.5 * visc_acc_term * dvdr_Hubble; + const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble; /* Diffusion term */ const float v_diff = max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f); - const float alpha_diff = 0.5 * (pi->diffusion.alpha + pj->diffusion.alpha); + const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha); /* wi_dx + wj_dx / 2 is F_ij */ const float diff_du_term = alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij; @@ -500,7 +500,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ /* Compute sound speeds and signal velocity */ - const float v_sig = 0.5 * (pi->viscosity.v_sig + pj->viscosity.v_sig); + const float v_sig = 0.5f * (pi->viscosity.v_sig + pj->viscosity.v_sig); /* Balsara term */ const float balsara_i = pi->force.balsara; @@ -535,12 +535,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( wi_dr * dvdr * r_inv; /* Viscosity term */ - const float visc_du_term = 0.5 * visc_acc_term * dvdr_Hubble; + const float visc_du_term = f * visc_acc_term * dvdr_Hubble; /* Diffusion term */ const float v_diff = max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f); - const float alpha_diff = 0.5 * (pi->diffusion.alpha + pj->diffusion.alpha); + const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha); /* wi_dx + wj_dx / 2 is F_ij */ const float diff_du_term = alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij; -- GitLab From 2399e7823b4a5c79573711c843150338c359d24b Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Sun, 31 Mar 2019 10:03:09 +0100 Subject: [PATCH 36/54] Fixed accidental deletion of a constant in anarchy --- src/hydro/AnarchyPU/hydro_iact.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index cc2128b03..199617a8f 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -535,7 +535,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( wi_dr * dvdr * r_inv; /* Viscosity term */ - const float visc_du_term = f * visc_acc_term * dvdr_Hubble; + const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble; /* Diffusion term */ const float v_diff = -- GitLab From 91d71b1bde4d17081ae9a6b884ac6fe30921db57 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Tue, 2 Apr 2019 11:19:13 +0100 Subject: [PATCH 37/54] Attempted a minor fix to the viscosity --- src/hydro/AnarchyPU/hydro_iact.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index 199617a8f..d091ebebc 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -379,7 +379,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( const float rho_ij = rhoi + rhoj; const float alpha = pi->viscosity.alpha + pj->viscosity.alpha; const float visc = - -1.f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; + -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; /* Convolve with the kernel */ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; @@ -510,7 +510,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( const float rho_ij = rhoi + rhoj; const float alpha = pi->viscosity.alpha + pj->viscosity.alpha; const float visc = - -1.f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; + -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; /* Convolve with the kernel */ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; -- GitLab From b7cf1576678ea17f1425c2abec7c3f37b3320201 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Tue, 2 Apr 2019 13:46:30 +0100 Subject: [PATCH 38/54] Updated ANARCHY defaults --- src/hydro_properties.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hydro_properties.c b/src/hydro_properties.c index 7d2b03be7..eb506839c 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -53,11 +53,11 @@ #define hydro_props_default_viscosity_alpha \ 0.1f /* Use a very low initial AV paramater for hydrodynamics tests */ #define hydro_props_default_viscosity_alpha_min \ - 0.1f /* values NOT the same as Schaller+ 2015 */ + 0.0f /* values NOT the same as Schaller+ 2015 */ #define hydro_props_default_viscosity_alpha_max \ 2.0f /* values taken from Schaller+ 2015 */ #define hydro_props_default_viscosity_length \ - 0.01f /* values taken from Schaller+ 2015 */ + 0.25f /* values taken from Schaller+ 2015 */ #else #define hydro_props_default_viscosity_alpha_min \ 0.1f /* values taken from (price,2004), not used in legacy gadget mode */ -- GitLab From 4efc2605473e2db2b6e6187bf062db2fb098b4ab Mon Sep 17 00:00:00 2001 From: JBorrow Date: Wed, 3 Apr 2019 12:41:46 +0100 Subject: [PATCH 39/54] Updated min-dt values --- .../InteractingBlastWaves_1D/interactingBlastWaves.yml | 2 +- examples/HydroTests/Noh_1D/noh.yml | 4 ++-- examples/HydroTests/Noh_2D/noh.yml | 2 +- examples/HydroTests/SedovBlast_2D/sedov.yml | 2 +- examples/HydroTests/SedovBlast_3D/sedov.yml | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml b/examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml index c4960dfa2..aaad57dc0 100644 --- a/examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml +++ b/examples/HydroTests/InteractingBlastWaves_1D/interactingBlastWaves.yml @@ -10,7 +10,7 @@ InternalUnitSystem: TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). time_end: 3.8e-2 # The end time of the simulation (in internal units). - dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). + dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). dt_max: 1e-5 # The maximal time-step size of the simulation (in internal units). # Parameters governing the snapshots diff --git a/examples/HydroTests/Noh_1D/noh.yml b/examples/HydroTests/Noh_1D/noh.yml index 58e13ddda..3a7e328df 100644 --- a/examples/HydroTests/Noh_1D/noh.yml +++ b/examples/HydroTests/Noh_1D/noh.yml @@ -10,7 +10,7 @@ InternalUnitSystem: TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). time_end: 0.6 # The end time of the simulation (in internal units). - dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). + dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units). # Parameters governing the snapshots @@ -33,4 +33,4 @@ InitialConditions: file_name: ./noh.hdf5 # The file to read periodic: 1 - \ No newline at end of file + diff --git a/examples/HydroTests/Noh_2D/noh.yml b/examples/HydroTests/Noh_2D/noh.yml index eaf991631..76ac8eb91 100644 --- a/examples/HydroTests/Noh_2D/noh.yml +++ b/examples/HydroTests/Noh_2D/noh.yml @@ -10,7 +10,7 @@ InternalUnitSystem: TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). time_end: 0.6 # The end time of the simulation (in internal units). - dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). + dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units). # Parameters governing the snapshots diff --git a/examples/HydroTests/SedovBlast_2D/sedov.yml b/examples/HydroTests/SedovBlast_2D/sedov.yml index b4252581d..208e4ac72 100644 --- a/examples/HydroTests/SedovBlast_2D/sedov.yml +++ b/examples/HydroTests/SedovBlast_2D/sedov.yml @@ -10,7 +10,7 @@ InternalUnitSystem: TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). time_end: 5e-2 # The end time of the simulation (in internal units). - dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). + dt_min: 1e-9 # 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). # Parameters governing the snapshots diff --git a/examples/HydroTests/SedovBlast_3D/sedov.yml b/examples/HydroTests/SedovBlast_3D/sedov.yml index 19e8c7253..00419ba94 100644 --- a/examples/HydroTests/SedovBlast_3D/sedov.yml +++ b/examples/HydroTests/SedovBlast_3D/sedov.yml @@ -10,7 +10,7 @@ InternalUnitSystem: TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). time_end: 5e-2 # The end time of the simulation (in internal units). - dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). + dt_min: 1e-9 # 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). # Parameters governing the snapshots -- GitLab From d69708edbee30a60c4ad081b1c564538571b7ad5 Mon Sep 17 00:00:00 2001 From: JBorrow Date: Wed, 3 Apr 2019 12:43:14 +0100 Subject: [PATCH 40/54] Minor changes to documentation --- doc/RTD/source/GettingStarted/running_example.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/RTD/source/GettingStarted/running_example.rst b/doc/RTD/source/GettingStarted/running_example.rst index 9dfbdd8c8..a5614fb2a 100644 --- a/doc/RTD/source/GettingStarted/running_example.rst +++ b/doc/RTD/source/GettingStarted/running_example.rst @@ -11,7 +11,7 @@ as ``wget`` for grabbing the glass). .. code-block:: bash - cd examples/SodShock_3D + cd examples/HydroTests/SodShock_3D ./getGlass.sh python makeIC.py ../swift --hydro --threads=4 sodShock.yml -- GitLab From 48f58459b7630d87474a2861d1535c0334f0c27b Mon Sep 17 00:00:00 2001 From: JBorrow Date: Wed, 3 Apr 2019 12:53:25 +0100 Subject: [PATCH 41/54] Added anarchy sph to docs --- doc/RTD/source/HydroSchemes/anarchy_sph.rst | 33 +++++++++++++++++++++ doc/RTD/source/HydroSchemes/index.rst | 1 + 2 files changed, 34 insertions(+) create mode 100644 doc/RTD/source/HydroSchemes/anarchy_sph.rst diff --git a/doc/RTD/source/HydroSchemes/anarchy_sph.rst b/doc/RTD/source/HydroSchemes/anarchy_sph.rst new file mode 100644 index 000000000..8d0928003 --- /dev/null +++ b/doc/RTD/source/HydroSchemes/anarchy_sph.rst @@ -0,0 +1,33 @@ +.. ANARCHY-SPH + Josh Borrow 5th April 2018 + +ANARCHY-PU SPH +============== + +This scheme is similar to the one used in the EAGLE code. This scheme +includes: + ++ Durier & Dalla Vecchia (2012) time-step limiter ++ Pressure-Energy SPH ++ Thermal diffusion following Price (2012) ++ A simplified version of the 'Inviscid SPH' artificial viscosity + (Cullen & Denhen 2010). + +More information will be made available in a forthcoming publication. + +The scheme as-implemented in SWIFT is slightly different to the one +implemented in the original EAGLE code: + ++ Pressure-Energy SPH is used instead of Pressure-Entropy SPH ++ Artificial viscosity coefficients have changed -- from minimal + value of 0.1 to 0.0, and from length of 0.1 to 0.25. This + is based on performance of hydrodynamics tests in SWIFT and may + be to do with our choice of smoothing length definition. ++ Recommended kernel changed from Wendland-C2 (with 100 Ngb) to + Quintic Spline (with ~82 Ngb). + + +.. code-block:: bash + + ./configure --with-hydro=anarchy-pu --with-kernel=quintic-spline --disable-hand-vec + diff --git a/doc/RTD/source/HydroSchemes/index.rst b/doc/RTD/source/HydroSchemes/index.rst index 462bb7378..7331331bd 100644 --- a/doc/RTD/source/HydroSchemes/index.rst +++ b/doc/RTD/source/HydroSchemes/index.rst @@ -17,6 +17,7 @@ schemes available in SWIFT, as well as how to implement your own. minimal_sph planetary hopkins_sph + anarchy_sph gizmo adding_your_own -- GitLab From fa1a141da984014ae24750036c55aeec900b045e Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 3 Apr 2019 14:25:41 +0100 Subject: [PATCH 42/54] Added build folder --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 798b187d4..e0e0be678 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ config.h.in config.sub ltmain.sh libtool +build src/version_string.h swift*.tar.gz -- GitLab From 617afd4fe30fff9afffd64c4d54f597d672334f5 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 3 Apr 2019 16:38:09 +0100 Subject: [PATCH 43/54] Increase the maximal time-step size of the non-cosmological cooling example. --- .../cooling_redshift_dependence_no_z.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml index 35242f553..8bce64984 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml @@ -11,7 +11,7 @@ TimeIntegration: time_begin: 0. time_end: 1.02e-4 # ~ 100 myr dt_min: 1e-14 - dt_max: 5e-8 + dt_max: 5e-5 # Parameters governing the snapshots Snapshots: -- GitLab From 25ad838f154681d2dfba68f50656b6496791f2b2 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 3 Apr 2019 16:46:41 +0100 Subject: [PATCH 44/54] Working version of the cosmological constant lambda cooling model. --- src/cooling/const_lambda/cooling.h | 77 ++++++++++++++++-------------- 1 file changed, 41 insertions(+), 36 deletions(-) diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index f5816b53e..8e53c044e 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -123,10 +123,10 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( if (dt == 0.) return; /* Current energy */ - const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo); + // const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo); - /* Current du_dt in physical coordinates (internal units) */ - const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo); + /* Y' */ + const float hydro_du_dt_com = hydro_get_comoving_internal_energy_dt(p); /* Calculate cooling du_dt (in cgs units) */ const double cooling_du_dt_cgs = @@ -136,45 +136,49 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( float cooling_du_dt = cooling_du_dt_cgs * cooling->conv_factor_energy_rate_from_cgs; - /* Add cosmological term */ - cooling_du_dt *= cosmo->a * cosmo->a; + /* Add cosmological term to get Y_cooling' */ + cooling_du_dt *= cosmo->a * cosmo->a / cosmo->a_factor_internal_energy; - float total_du_dt = hydro_du_dt + cooling_du_dt; + /* Y_total' */ + float total_du_dt = hydro_du_dt_com + cooling_du_dt; /* We now need to check that we are not going to go below any of the limits */ - /* Limit imposed by the entropy floor */ - const float A_floor = entropy_floor(p, cosmo, floor_props); - const float rho = hydro_get_physical_density(p, cosmo); - const float u_floor = gas_internal_energy_from_entropy(rho, A_floor); - - /* Absolute minimum */ - const float u_minimal = hydro_props->minimal_internal_energy; - - /* Largest of both limits */ - const float u_limit = max(u_minimal, u_floor); - - /* First, check whether we may end up below the minimal energy after - * this step 1/2 kick + another 1/2 kick that could potentially be for - * a time-step twice as big. We hence check for 1.5 delta_t. */ - if (u_old + total_du_dt * 1.5 * dt_therm < u_limit) { - total_du_dt = (u_limit - u_old) / (1.5f * dt_therm); - } - - /* Second, check whether the energy used in the prediction could get negative. - * We need to check for the 1/2 dt kick followed by a full time-step drift - * that could potentially be for a time-step twice as big. We hence check - * for 2.5 delta_t but this time against 0 energy not the minimum */ - if (u_old + total_du_dt * 2.5 * dt_therm < 0.) { - total_du_dt = -u_old / ((2.5f + 0.0001f) * dt_therm); - } + /* /\* Limit imposed by the entropy floor *\/ */ + /* const float A_floor = entropy_floor(p, cosmo, floor_props); */ + /* const float rho = hydro_get_physical_density(p, cosmo); */ + /* const float u_floor = gas_internal_energy_from_entropy(rho, A_floor); */ + + /* /\* Absolute minimum *\/ */ + /* const float u_minimal = hydro_props->minimal_internal_energy; */ + + /* /\* Largest of both limits *\/ */ + /* const float u_limit = max(u_minimal, u_floor); */ + + /* /\* First, check whether we may end up below the minimal energy after */ + /* * this step 1/2 kick + another 1/2 kick that could potentially be for */ + /* * a time-step twice as big. We hence check for 1.5 delta_t. *\/ */ + /* if (u_old + total_du_dt * 1.5 * dt_therm < u_limit) { */ + /* total_du_dt = (u_limit - u_old) / (1.5f * dt_therm); */ + /* } */ + + /* /\* Second, check whether the energy used in the prediction could get + * negative. */ + /* * We need to check for the 1/2 dt kick followed by a full time-step drift + */ + /* * that could potentially be for a time-step twice as big. We hence check + */ + /* * for 2.5 delta_t but this time against 0 energy not the minimum *\/ */ + /* if (u_old + total_du_dt * 2.5 * dt_therm < 0.) { */ + /* total_du_dt = -u_old / ((2.5f + 0.0001f) * dt_therm); */ + /* } */ /* Update the internal energy time derivative */ - hydro_set_physical_internal_energy_dt(p, cosmo, total_du_dt); + hydro_set_comoving_internal_energy_dt(p, total_du_dt); /* Store the radiated energy (assuming dt will not change) */ - xp->cooling_data.radiated_energy += - -hydro_get_mass(p) * (total_du_dt - hydro_du_dt) * dt_therm; + /* xp->cooling_data.radiated_energy += */ + /* -hydro_get_mass(p) * (total_du_dt - hydro_du_dt) * dt_therm; */ } /** @@ -208,8 +212,9 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( cooling_rate_cgs(cosmo, hydro_props, cooling, p); /* Convert to internal units */ - const float cooling_du_dt = - cooling_du_dt_cgs * cooling->conv_factor_energy_rate_from_cgs; + const float cooling_du_dt = cooling_du_dt_cgs * + cooling->conv_factor_energy_rate_from_cgs / + cosmo->a_factor_internal_energy; /* If we are close to (or below) the limit, we ignore the condition */ if (u < 1.01f * hydro_props->minimal_internal_energy || cooling_du_dt == 0.f) -- GitLab From 429a7cc90561303517cba2cfb17c7943d0e10484 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 3 Apr 2019 17:03:01 +0100 Subject: [PATCH 45/54] Update the const-lambda cooling scheme to obey the minimal temperature floor. --- src/cooling/const_lambda/cooling.h | 69 +++++++++++++++--------------- 1 file changed, 35 insertions(+), 34 deletions(-) diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 8e53c044e..51ac66074 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -122,8 +122,8 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( /* Nothing to do here? */ if (dt == 0.) return; - /* Current energy */ - // const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo); + /* Current energy (in internal units) */ + const float u_old_com = hydro_get_comoving_internal_energy(p, xp); /* Y' */ const float hydro_du_dt_com = hydro_get_comoving_internal_energy_dt(p); @@ -133,52 +133,53 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( cooling_rate_cgs(cosmo, hydro_props, cooling, p); /* Convert to internal units */ - float cooling_du_dt = + const float cooling_du_dt_physical = cooling_du_dt_cgs * cooling->conv_factor_energy_rate_from_cgs; /* Add cosmological term to get Y_cooling' */ - cooling_du_dt *= cosmo->a * cosmo->a / cosmo->a_factor_internal_energy; + const float cooling_du_dt = cooling_du_dt_physical * cosmo->a * cosmo->a / + cosmo->a_factor_internal_energy; /* Y_total' */ float total_du_dt = hydro_du_dt_com + cooling_du_dt; /* We now need to check that we are not going to go below any of the limits */ - /* /\* Limit imposed by the entropy floor *\/ */ - /* const float A_floor = entropy_floor(p, cosmo, floor_props); */ - /* const float rho = hydro_get_physical_density(p, cosmo); */ - /* const float u_floor = gas_internal_energy_from_entropy(rho, A_floor); */ - - /* /\* Absolute minimum *\/ */ - /* const float u_minimal = hydro_props->minimal_internal_energy; */ - - /* /\* Largest of both limits *\/ */ - /* const float u_limit = max(u_minimal, u_floor); */ - - /* /\* First, check whether we may end up below the minimal energy after */ - /* * this step 1/2 kick + another 1/2 kick that could potentially be for */ - /* * a time-step twice as big. We hence check for 1.5 delta_t. *\/ */ - /* if (u_old + total_du_dt * 1.5 * dt_therm < u_limit) { */ - /* total_du_dt = (u_limit - u_old) / (1.5f * dt_therm); */ - /* } */ - - /* /\* Second, check whether the energy used in the prediction could get - * negative. */ - /* * We need to check for the 1/2 dt kick followed by a full time-step drift - */ - /* * that could potentially be for a time-step twice as big. We hence check - */ - /* * for 2.5 delta_t but this time against 0 energy not the minimum *\/ */ - /* if (u_old + total_du_dt * 2.5 * dt_therm < 0.) { */ - /* total_du_dt = -u_old / ((2.5f + 0.0001f) * dt_therm); */ - /* } */ + /* Limit imposed by the entropy floor (comoving) + * (Recall entropy is the same in physical and comoving frames) */ + const float A_floor_com = entropy_floor(p, cosmo, floor_props); + const float rho_com = hydro_get_comoving_density(p); + const float u_floor_com = + gas_internal_energy_from_entropy(rho_com, A_floor_com); + + /* Absolute minimum */ + const float u_minimal_com = + hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; + + /* Largest of both limits */ + const float u_limit_com = max(u_minimal_com, u_floor_com); + + /* First, check whether we may end up below the minimal energy after + * this step 1/2 kick + another 1/2 kick that could potentially be for + * a time-step twice as big. We hence check for 1.5 delta_t. */ + if (u_old_com + total_du_dt * 1.5 * dt_therm < u_limit_com) { + total_du_dt = (u_limit_com - u_old_com) / (1.5f * dt_therm); + } + + /* Second, check whether the energy used in the prediction could get negative. + * We need to check for the 1/2 dt kick followed by a full time-step drift + * that could potentially be for a time-step twice as big. We hence check + * for 2.5 delta_t but this time against 0 energy not the minimum */ + if (u_old_com + total_du_dt * 2.5 * dt_therm < 0.) { + total_du_dt = -u_old_com / ((2.5f + 0.0001f) * dt_therm); + } /* Update the internal energy time derivative */ hydro_set_comoving_internal_energy_dt(p, total_du_dt); /* Store the radiated energy (assuming dt will not change) */ - /* xp->cooling_data.radiated_energy += */ - /* -hydro_get_mass(p) * (total_du_dt - hydro_du_dt) * dt_therm; */ + xp->cooling_data.radiated_energy += + -hydro_get_mass(p) * cooling_du_dt_physical * dt; } /** -- GitLab From d0f76d64c27324bcf5068210b78ac39bacda9548 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 3 Apr 2019 17:23:07 +0100 Subject: [PATCH 46/54] Add the radiated energy to the fields output by the const-lambda cooling. --- src/cooling/const_lambda/cooling_io.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h index 9437f0f94..2e2ba799a 100644 --- a/src/cooling/const_lambda/cooling_io.h +++ b/src/cooling/const_lambda/cooling_io.h @@ -76,7 +76,10 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles( UNIT_CONV_TEMPERATURE, parts, xparts, convert_part_T); - return 1; + list[1] = io_make_output_field("RadiatedEnergy", FLOAT, 1, UNIT_CONV_ENERGY, + xparts, cooling_data.radiated_energy); + + return 2; } #endif /* SWIFT_COOLING_CONST_LAMBDA_IO_H */ -- GitLab From 389bfb4372fef915909547142002d87fd67467ce Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 3 Apr 2019 17:56:42 +0100 Subject: [PATCH 47/54] Added panel to show energy radiated --- .../CoolingRedshiftDependence/plotSolution.py | 136 ++++++++++++++---- 1 file changed, 111 insertions(+), 25 deletions(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py index d7f4fe3b7..7b750a240 100644 --- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py +++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py @@ -6,29 +6,69 @@ import matplotlib.pyplot as plt from swiftsimio import load -from unyt import Myr, K, mh, cm +from unyt import Myr, K, mh, cm, erg import numpy as np +def get_data_dump(metadata): + """ + Gets a big data dump from the SWIFT metadata + """ + + try: + viscosity = metadata.viscosity_info + except: + viscosity = "No info" + + try: + diffusion = metadata.diffusion_info + except: + diffusion = "No info" + + output = ( + "$\\bf{SWIFT}$\n" + + metadata.code_info + + "\n\n" + + "$\\bf{Compiler}$\n" + + metadata.compiler_info + + "\n\n" + + "$\\bf{Hydrodynamics}$\n" + + metadata.hydro_info + + "\n\n" + + "$\\bf{Viscosity}$\n" + + viscosity + + "\n\n" + + "$\\bf{Diffusion}$\n" + + diffusion + ) + + return output + + def setup_axes(): """ Sets up the axes and returns fig, ax """ - fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(4, 8), sharex=True) + fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8, 8)) ax = ax.flatten() ax[0].semilogy() + ax[2].semilogy() ax[1].set_xlabel("Simulation time elapsed [Myr]") + ax[2].set_xlabel("Simulation time elapsed [Myr]") + ax[0].set_ylabel("Temperature of Universe [K]") ax[1].set_ylabel("Physical Density of Universe [$n_H$ cm$^{-3}$]") + ax[2].set_ylabel("Physical Energy [erg]") - ax[0].set_xlim(0, 100) + for a in ax[:-1]: + a.set_xlim(0, 100) - fig.tight_layout() + fig.tight_layout(pad=0.5) return fig, ax @@ -43,6 +83,8 @@ def get_data(handle: float, n_snaps: int): times = [] temps = [] densities = [] + energies = [] + radiated_energies = [] for snap in range(n_snaps): data = load(f"data/{handle}_{snap:04d}.hdf5") @@ -52,9 +94,18 @@ def get_data(handle: float, n_snaps: int): times.append(data.metadata.t.to(Myr).value - t0) temps.append(np.mean(data.gas.temperature.to(K).value)) - densities.append(np.mean(data.gas.density.to(mh / cm**3).value) / (data.metadata.scale_factor**3) ) + densities.append( + np.mean(data.gas.density.to(mh / cm ** 3).value) + / (data.metadata.scale_factor ** 3) + ) + + energies.append( + np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value) + * data.metadata.scale_factor ** (2) + ) + radiated_energies.append(np.mean(data.gas.radiated_energy.to(erg).value)) - return times, temps, densities + return times, temps, densities, energies, radiated_energies def get_n_steps(timesteps_filename: str) -> int: @@ -67,28 +118,49 @@ def get_n_steps(timesteps_filename: str) -> int: return int(data[-1][0]) -def plot_single_data(handle: str, n_snaps: int, label: str, ax: plt.axes, n_steps:int = 0): +def plot_single_data( + handle: str, n_snaps: int, label: str, ax: plt.axes, n_steps: int = 0, run: int = 0 +): """ Takes the a single simulation and plots it on the axes. """ data = get_data(handle, n_snaps) - ax[0].plot( - data[0], data[1], - label=label, - marker="o", - ms=2 - ) + ax[0].plot(data[0], data[1], label=label, marker="o", ms=2, c=f"C{run}") ax[1].plot( - data[0], data[2], - label=f"Steps: {n_steps}", - marker="o", - ms=2 + data[0], data[2], label=f"Steps: {n_steps}", marker="o", ms=2, c=f"C{run}" + ) + + if run == 0: + label_energy = "Particle Energy" + label_radiated = "Radiated Energy" + label_sum = "Sum" + else: + label_energy = "" + label_radiated = "" + label_sum = "" + + ax[2].plot( + data[0], data[3], label=label_energy, ls="dotted", C=f"C{run}" + ) + + ax[2].plot( + data[0], + data[4], + label=label_radiated, + ls="dashed", + C=f"C{run}", + ) + + ax[2].plot( + data[0], + [x + y for x, y in zip(*data[3:5])], + label=label_sum, + C=f"C{run}", ) - return @@ -99,21 +171,32 @@ def make_plot(handles, names, timestep_names, n_snaps=100): fig, ax = setup_axes() + run = 0 + for handle, name, timestep_name in zip(handles, names, timestep_names): try: n_steps = get_n_steps(timestep_name) except: n_steps = "Unknown" - try: - plot_single_data(handle, n_snaps, name, ax, n_steps=n_steps) - except: - pass + plot_single_data(handle, n_snaps, name, ax, n_steps=n_steps, run=run) + run += 1 ax[0].legend() ax[1].legend() + ax[2].legend() - return fig, ax + info_axis = ax[-1] + + info = get_data_dump(load(f"data/{handles[0]}_0000.hdf5").metadata) + + info_axis.text( + 0.5, 0.45, info, ha="center", va="center", fontsize=7, transform=info_axis.transAxes + ) + + info_axis.axis("off") + + return fig, ax if __name__ == "__main__": @@ -121,11 +204,14 @@ if __name__ == "__main__": Plot everything! """ - handles = ["redshift_dependence_no_z", "redshift_dependence_low_z", "redshift_dependence_high_z"] + handles = [ + "redshift_dependence_no_z", + "redshift_dependence_low_z", + "redshift_dependence_high_z", + ] names = ["No Cosmology", "Low Redshift ($z=0.01$)", "High Redshift ($z=1$)"] timestep_names = ["timesteps_no.txt", "timesteps_low.txt", "timesteps_high.txt"] fig, ax = make_plot(handles, names, timestep_names) fig.savefig("redshift_dependence.png", dpi=300) - -- GitLab From fb614b04aab9f51c15f686392772e525d2226a45 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Wed, 3 Apr 2019 18:24:23 +0100 Subject: [PATCH 48/54] Made having the energies radiated optional --- .../CoolingRedshiftDependence/plotSolution.py | 38 +++++++++---------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py index 7b750a240..2be73a165 100644 --- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py +++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py @@ -99,11 +99,14 @@ def get_data(handle: float, n_snaps: int): / (data.metadata.scale_factor ** 3) ) - energies.append( - np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value) - * data.metadata.scale_factor ** (2) - ) - radiated_energies.append(np.mean(data.gas.radiated_energy.to(erg).value)) + try: + energies.append( + np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value) + * data.metadata.scale_factor ** (2) + ) + radiated_energies.append(np.mean(data.gas.radiated_energy.to(erg).value)) + except AttributeError: + continue return times, temps, densities, energies, radiated_energies @@ -142,23 +145,12 @@ def plot_single_data( label_radiated = "" label_sum = "" - ax[2].plot( - data[0], data[3], label=label_energy, ls="dotted", C=f"C{run}" - ) + ax[2].plot(data[0], data[3], label=label_energy, ls="dotted", C=f"C{run}") - ax[2].plot( - data[0], - data[4], - label=label_radiated, - ls="dashed", - C=f"C{run}", - ) + ax[2].plot(data[0], data[4], label=label_radiated, ls="dashed", C=f"C{run}") ax[2].plot( - data[0], - [x + y for x, y in zip(*data[3:5])], - label=label_sum, - C=f"C{run}", + data[0], [x + y for x, y in zip(*data[3:5])], label=label_sum, C=f"C{run}" ) return @@ -191,7 +183,13 @@ def make_plot(handles, names, timestep_names, n_snaps=100): info = get_data_dump(load(f"data/{handles[0]}_0000.hdf5").metadata) info_axis.text( - 0.5, 0.45, info, ha="center", va="center", fontsize=7, transform=info_axis.transAxes + 0.5, + 0.45, + info, + ha="center", + va="center", + fontsize=7, + transform=info_axis.transAxes, ) info_axis.axis("off") -- GitLab From 89f3d1c61f171df6ef23bffd301799dddb9ff986 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Thu, 4 Apr 2019 09:53:19 +0100 Subject: [PATCH 49/54] Use the drifted internal energy for the temperature calculation in the i/o. --- src/cooling/const_lambda/cooling.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 51ac66074..3518fe21f 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -276,7 +276,7 @@ INLINE static float cooling_get_temperature( const double mu_ionised = hydro_props->mu_ionised; /* Particle temperature */ - const double u = hydro_get_physical_internal_energy(p, xp, cosmo); + const double u = hydro_get_drifted_physical_internal_energy(p, cosmo); /* Temperature over mean molecular weight */ const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B; -- GitLab From 233ab52cac06b6444fb18502a93313da0f9acf51 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Thu, 4 Apr 2019 16:38:45 +0100 Subject: [PATCH 50/54] Fix the script for non completed runs. --- .../CoolingRedshiftDependence/plotSolution.py | 70 +++++++++++-------- 1 file changed, 39 insertions(+), 31 deletions(-) diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py index 2be73a165..9cde36cb0 100644 --- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py +++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py @@ -87,25 +87,30 @@ def get_data(handle: float, n_snaps: int): radiated_energies = [] for snap in range(n_snaps): - data = load(f"data/{handle}_{snap:04d}.hdf5") - - if snap == 0: - t0 = data.metadata.t.to(Myr).value + try: + data = load(f"data/{handle}_{snap:04d}.hdf5") - times.append(data.metadata.t.to(Myr).value - t0) - temps.append(np.mean(data.gas.temperature.to(K).value)) - densities.append( - np.mean(data.gas.density.to(mh / cm ** 3).value) - / (data.metadata.scale_factor ** 3) - ) + if snap == 0: + t0 = data.metadata.t.to(Myr).value - try: - energies.append( - np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value) - * data.metadata.scale_factor ** (2) + times.append(data.metadata.t.to(Myr).value - t0) + temps.append(np.mean(data.gas.temperature.to(K).value)) + densities.append( + np.mean(data.gas.density.to(mh / cm ** 3).value) + / (data.metadata.scale_factor ** 3) ) - radiated_energies.append(np.mean(data.gas.radiated_energy.to(erg).value)) - except AttributeError: + + try: + energies.append( + np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value) + * data.metadata.scale_factor ** (2) + ) + radiated_energies.append( + np.mean(data.gas.radiated_energy.to(erg).value) + ) + except AttributeError: + continue + except OSError: continue return times, temps, densities, energies, radiated_energies @@ -147,11 +152,11 @@ def plot_single_data( ax[2].plot(data[0], data[3], label=label_energy, ls="dotted", C=f"C{run}") - ax[2].plot(data[0], data[4], label=label_radiated, ls="dashed", C=f"C{run}") + # ax[2].plot(data[0], data[4], label=label_radiated, ls="dashed", C=f"C{run}") - ax[2].plot( - data[0], [x + y for x, y in zip(*data[3:5])], label=label_sum, C=f"C{run}" - ) + # ax[2].plot( + # data[0], [x + y for x, y in zip(*data[3:5])], label=label_sum, C=f"C{run}" + # ) return @@ -180,17 +185,20 @@ def make_plot(handles, names, timestep_names, n_snaps=100): info_axis = ax[-1] - info = get_data_dump(load(f"data/{handles[0]}_0000.hdf5").metadata) - - info_axis.text( - 0.5, - 0.45, - info, - ha="center", - va="center", - fontsize=7, - transform=info_axis.transAxes, - ) + try: + info = get_data_dump(load(f"data/{handles[0]}_0000.hdf5").metadata) + + info_axis.text( + 0.5, + 0.45, + info, + ha="center", + va="center", + fontsize=7, + transform=info_axis.transAxes, + ) + except OSError: + pass info_axis.axis("off") -- GitLab From 6ed712bb9e3c5e071f52e7daa98263b7a7424928 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Sat, 6 Apr 2019 09:59:42 +0100 Subject: [PATCH 51/54] Added data gitignore to the CoolingRedshiftDependence example --- examples/Cooling/CoolingRedshiftDependence/.gitignore | 1 + 1 file changed, 1 insertion(+) create mode 100644 examples/Cooling/CoolingRedshiftDependence/.gitignore diff --git a/examples/Cooling/CoolingRedshiftDependence/.gitignore b/examples/Cooling/CoolingRedshiftDependence/.gitignore new file mode 100644 index 000000000..60baa9cb8 --- /dev/null +++ b/examples/Cooling/CoolingRedshiftDependence/.gitignore @@ -0,0 +1 @@ +data/* -- GitLab From 3dc745d0f6fba89f0589021b95d0ca14cca5e2db Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Sat, 6 Apr 2019 12:26:56 +0100 Subject: [PATCH 52/54] Clarified comment for the comoving nature of du/dt in the const-lambda cooling model. --- src/cooling/const_lambda/cooling.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 3518fe21f..fddc9a620 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -125,7 +125,8 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( /* Current energy (in internal units) */ const float u_old_com = hydro_get_comoving_internal_energy(p, xp); - /* Y' */ + /* Y' = RHS of the comoving equation for du/dt that will be integrated + forward in time using dt_therm */ const float hydro_du_dt_com = hydro_get_comoving_internal_energy_dt(p); /* Calculate cooling du_dt (in cgs units) */ -- GitLab From 8f7de5b9216870d1d08e477395674d8947dbc518 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Mon, 8 Apr 2019 14:22:40 -0400 Subject: [PATCH 53/54] Clarified comment in Morris&Monaghan --- src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h index 819abeec9..92d0fa84e 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h @@ -614,9 +614,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Artificial viscosity updates */ - /* Need to divide by a^2 here because this is a timescale and we need to - introduce the time-integration factor to balance out the scale-free - nature of the artificial viscosity (i.e. dt has a^2 built in) */ + /* We divide by a^2 here to make this transform under cosmology the + * same as the velocity (which in SWIFT has an extra 1/a^2 factor. + * See the cosmology theory documents for more information. */ const float inverse_tau = (hydro_props->viscosity.length * cosmo->a2_inv) * soundspeed * h_inv; const float source_term = -1.f * min(p->density.div_v, 0.f); @@ -873,4 +873,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) { __attribute__((always_inline)) INLINE static void hydro_remove_part( const struct part *p, const struct xpart *xp) {} -#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */ \ No newline at end of file +#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H */ -- GitLab From 6a5a87d9650d26fba82744bf69bd8b54c317ce2a Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Mon, 8 Apr 2019 14:26:51 -0400 Subject: [PATCH 54/54] Added TODO comment in the morris & monaghan for cosmology calculations --- src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h index 92d0fa84e..be6a789eb 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h @@ -614,6 +614,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Artificial viscosity updates */ + /* TODO: Actually work out why this cosmology factor is correct + * and update the SPH / cosmology theory documents. */ + /* We divide by a^2 here to make this transform under cosmology the * same as the velocity (which in SWIFT has an extra 1/a^2 factor. * See the cosmology theory documents for more information. */ -- GitLab