From 600c1b0f67f7805385d794e28a3f075feb634908 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 11 Jan 2023 20:30:19 -0500 Subject: [PATCH 1/4] Fixed another string variable NetCDF bug --- src/encounter/encounter_io.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 1333c05aa..50180113c 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -232,7 +232,7 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param) ! Internals integer(I4B) :: i, idslot, old_mode, npl, ntp - character(len=:), allocatable :: charstring + character(len=STRMAX) :: charstring select type(param) class is (swiftest_parameters) From acb580e8c642124784eab2b361acf86f3c7d28b7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 11 Jan 2023 20:30:40 -0500 Subject: [PATCH 2/4] Fixed potential issues with underflow and NaN values for degenerate orbits --- src/swiftest/swiftest_gr.f90 | 4 ++++ src/swiftest/swiftest_orbel.f90 | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/src/swiftest/swiftest_gr.f90 b/src/swiftest/swiftest_gr.f90 index 3274f218e..1c37d3da4 100644 --- a/src/swiftest/swiftest_gr.f90 +++ b/src/swiftest/swiftest_gr.f90 @@ -200,6 +200,10 @@ pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv) pv(1:NDIM) = vh(1:NDIM) ! Initial guess rterm = 3 * mu / norm2(rh(:)) v2 = dot_product(vh(:), vh(:)) + if (v2 < TINY(1.0_DP)) then + pv(:) = 0.0_DP + return + end if do n = 1, MAXITER pv2 = dot_product(pv(:), pv(:)) G = 1.0_DP - inv_c2 * (0.5_DP * pv2 + rterm) diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90 index a18ef2a19..bf72c91b2 100644 --- a/src/swiftest/swiftest_orbel.f90 +++ b/src/swiftest/swiftest_orbel.f90 @@ -939,6 +939,11 @@ pure module subroutine swiftest_orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, in capom = 0.0_DP omega = 0.0_DP capm = 0.0_DP + varpi = 0.0_DP + lam = 0.0_DP + f = 0.0_DP + cape = 0.0_DP + capf = 0.0_DP x = [px, py, pz] v = [vx, vy, vz] r = .mag. x(:) From 86e7812285dcc3c095722be91d89d6c32716b307 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 11 Jan 2023 20:40:50 -0500 Subject: [PATCH 3/4] Fixed possible units issue in the computation of Qloss --- src/collision/collision_regime.f90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index 6e5199ca3..4de99a844 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -65,7 +65,7 @@ subroutine collision_regime_LS12(impactors, nbody_system, param) integer(I4B) :: jtarg, jproj real(DP), dimension(2) :: radius_si, mass_si, density_si real(DP) :: min_mfrag_si, Mcb_si - real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si, runit + real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si real(DP) :: mlr, mslr, mtot, dentot integer(I4B), parameter :: NMASS_DIST = 3 !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder) @@ -149,7 +149,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, real(DP), parameter :: CRUFU = 2.0_DP - 3 * MU_BAR ! central potential variable from Rufu and Aharonson (2019) real(DP), parameter :: SUPERCAT_QRATIO = 1.8_DP ! See Section 4.1 of LS12 ! Internals - real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta + real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta, ke, pe, Qmerge real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat real(DP) :: Mint, Mtot, Mtmp @@ -213,13 +213,16 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3) bcrit = rad1 / (rad1 + rad2) Qloss = 0.0_DP - U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 + U_binding = (3 * GC * Mtot**2) / (5 * Rp) ! LS12 eq. 27 + ke = 0.5_DP * Mtot * Vimp**2 + pe = - GC * m1 * m2 / norm2(rh2 - rh1) + Qmerge = ke + pe + U_binding if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP - Qloss = 0.0_DP + Qloss = Qmerge call swiftest_io_log_one_message(COLLISION_LOG_OUT, & "Fragments would have mass below the minimum. Converting this collision into a merger.") else @@ -227,18 +230,18 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP - Qloss = 0.0_DP + Qloss = Qmerge else if (Vimp < Verosion) then if (b < bcrit) then regime = COLLRESOLVE_REGIME_MERGE !partial accretion regime" Mlr = Mtot Mslr = 0.0_DP - Qloss = 0.0_DP + Qloss = Qmerge else if ((b > bcrit) .and. (Vimp < Vhr)) then regime = COLLRESOLVE_REGIME_MERGE ! graze and merge Mlr = Mtot Mslr = 0.0_DP - Qloss = 0.0_DP + Qloss = Qmerge else Mlr = m1 Mslr = max(calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star), min_mfrag) @@ -250,7 +253,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, regime = COLLRESOLVE_REGIME_MERGE !cratering regime" Mlr = Mtot Mslr = 0.0_DP - Qloss = 0.0_DP + Qloss = Qmerge else Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag) ! LS12 eq (37) regime = COLLRESOLVE_REGIME_DISRUPTION !disruption From 337b4bd0286b3f1689c55afce8971396804bf30f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 11 Jan 2023 22:54:18 -0500 Subject: [PATCH 4/4] Improved the convergence of Fraggle on the energy constraint. Fixed issue where, because the data file is not written to every time there's a collision, particles can come and go before ever being recorded, leaving gaps in the name dimension of the NetCDF file. --- python/swiftest/swiftest/io.py | 1 + src/collision/collision_resolve.f90 | 1 - src/fraggle/fraggle_generate.f90 | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index e76c772aa..d7367774a 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -820,6 +820,7 @@ def process_netcdf_input(ds, param): elif param['OUT_TYPE'] == "NETCDF_FLOAT": ds = fix_types(ds,ftype=np.float32) + ds = ds.where(ds['id']>=0, drop=True) return ds def swiftest2xr(param, verbose=True): diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 07bdc6fb3..52948c4d2 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -325,7 +325,6 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) integer(I4B) :: merge_text_length character(len=NAMELEN) :: merge_text character(len=NAMELEN) :: newname, origin_type - character(len=STRMAX) :: message real(DP) :: volume select type(nbody_system) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index e7b0aabba..611a34eb3 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -420,7 +420,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP), parameter :: vmax_initial_factor = 5.0_DP ! For the initial "guess" of fragment velocities, this is the maximum velocity relative to escape velocity that the fragments will have integer(I4B), parameter :: MAXLOOP = 200 integer(I4B), parameter :: MAXTRY = 20 - real(DP), parameter :: SUCCESS_METRIC = 1.0_DP + real(DP), parameter :: SUCCESS_METRIC = 1.0e-3_DP class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message