diff --git a/docs/src/fragmentation.f90 b/docs/src/fragmentation.f90 index 460060183..c90f64cb4 100644 --- a/docs/src/fragmentation.f90 +++ b/docs/src/fragmentation.f90 @@ -369,7 +369,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - plwksp%nplm = count(plwksp%Gmass > param%mtiny / mscale) + plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale) end select end select call tmpsys%pl%eucl_index() @@ -381,7 +381,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - nplm = count(pl%mass > param%mtiny) + nplm = count(pl%mass > param%Gmtiny) end select end select if (lk_plpl) call pl%eucl_index() @@ -836,7 +836,7 @@ end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -857,7 +857,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v ! Arguments integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision ! Constants @@ -931,7 +931,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v Qloss = 0.0_DP U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 - if ((m1 < mtiny).or.(m2 < mtiny)) then + if ((m1 < Gmtiny).or.(m2 < Gmtiny)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP diff --git a/docs/src/io.f90 b/docs/src/io.f90 index de1226951..1e0e8d626 100644 --- a/docs/src/io.f90 +++ b/docs/src/io.f90 @@ -493,7 +493,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *) param%Ecollisions case("EUNTRACKED") read(param_value, *) param%Euntracked - case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 diff --git a/docs/src/swiftest_classes.f90 b/docs/src/swiftest_classes.f90 index 2455e77f2..051add1aa 100644 --- a/docs/src/swiftest_classes.f90 +++ b/docs/src/swiftest_classes.f90 @@ -469,11 +469,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP), intent(inout) :: Qloss end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) implicit none integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision end subroutine fragmentation_regime diff --git a/docs/src/symba_classes.f90 b/docs/src/symba_classes.f90 index 4cf69a3e1..6ecb27751 100644 --- a/docs/src/symba_classes.f90 +++ b/docs/src/symba_classes.f90 @@ -19,7 +19,7 @@ module symba_classes type, extends(swiftest_parameters) :: symba_parameters character(STRMAX) :: particle_file = PARTICLE_OUTFILE !! Name of output particle information file - real(DP) :: MTINY = -1.0_DP !! Smallest mass that is fully gravitating + real(DP) :: GMTINY = -1.0_DP !! Smallest mass that is fully gravitating integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. contains @@ -73,9 +73,9 @@ module symba_classes type, extends(helio_pl) :: symba_pl logical, dimension(:), allocatable :: lcollision !! flag indicating whether body has merged with another this time step logical, dimension(:), allocatable :: lencounter !! flag indicating whether body is part of an encounter this time step - logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the MTINY cutoff value - integer(I4B) :: nplm !! number of bodies above the MTINY limit - integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above MTINY) comparisons in the flattened upper triangular matrix + logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the GMTINY cutoff value + integer(I4B) :: nplm !! number of bodies above the GMTINY limit + integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above GMTINY) comparisons in the flattened upper triangular matrix integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with other planets this time step integer(I4B), dimension(:), allocatable :: ntpenc !! number of encounters with test particles this time step integer(I4B), dimension(:), allocatable :: levelg !! level at which this body should be moved diff --git a/docs/src/symba_collision.f90 b/docs/src/symba_collision.f90 index 952d59709..caffdf296 100644 --- a/docs/src/symba_collision.f90 +++ b/docs/src/symba_collision.f90 @@ -449,7 +449,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children Mcb_si = cb%mass * param%MU2KG - mtiny_si = (param%MTINY / param%GU) * param%MU2KG + mtiny_si = (param%GMTINY / param%GU) * param%MU2KG mass_res(:) = 0.0_DP diff --git a/docs/src/symba_io.f90 b/docs/src/symba_io.f90 index 8e6420f61..713429365 100644 --- a/docs/src/symba_io.f90 +++ b/docs/src/symba_io.f90 @@ -71,8 +71,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms case ("FRAGMENTATION") call io_toupper(param_value) if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. - case ("MTINY") - read(param_value, *) param%mtiny + case ("GMTINY") + read(param_value, *) param%Gmtiny case("SEED") read(param_value, *) nseeds_from_file ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different @@ -111,12 +111,12 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) end if - if (self%mtiny < 0.0_DP) then - write(iomsg,*) "MTINY invalid or not set: ", self%mtiny + if (self%Gmtiny < 0.0_DP) then + write(iomsg,*) "GMTINY invalid or not set: ", self%Gmtiny iostat = -1 return else - write(*,*) "MTINY = ", self%mtiny + write(*,*) "GMTINY = ", self%Gmtiny end if if (.not.self%lclose) then @@ -167,7 +167,7 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms ! Special handling is required for writing the random number seed array as its size is not known until runtime ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "MTINY"; write(param_value, Rfmt) param%mtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%Gmtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%lfragmentation) then write(param_name, Afmt) "SEED" diff --git a/docs/src/symba_setup.f90 b/docs/src/symba_setup.f90 index 021873a70..ab8b5543e 100644 --- a/docs/src/symba_setup.f90 +++ b/docs/src/symba_setup.f90 @@ -25,7 +25,7 @@ module subroutine symba_setup_initialize_system(self, param) call pl%sort("mass", ascending=.false.) select type(param) class is (symba_parameters) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) end select end select diff --git a/docs/src/symba_util.f90 b/docs/src/symba_util.f90 index 0517cff9a..911f85304 100644 --- a/docs/src/symba_util.f90 +++ b/docs/src/symba_util.f90 @@ -392,7 +392,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then call pl%sort("mass", ascending=.false.) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) call pl%eucl_index() end if diff --git a/examples/symba_energy_momentum/param.disruption_headon.in b/examples/symba_energy_momentum/param.disruption_headon.in index de3c83bea..6dbe1f788 100644 --- a/examples/symba_energy_momentum/param.disruption_headon.in +++ b/examples/symba_energy_momentum/param.disruption_headon.in @@ -17,13 +17,11 @@ CHK_RMIN 0.005 CHK_RMAX 1e6 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -!CHK_QMIN_COORD HELIO ! commented out here -!CHK_QMIN_RANGE 1.0 1000.0 ! commented out here ENC_OUT enc.disruption_headon.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.disruption_off_axis.in b/examples/symba_energy_momentum/param.disruption_off_axis.in index bb0915050..39303284e 100644 --- a/examples/symba_energy_momentum/param.disruption_off_axis.in +++ b/examples/symba_energy_momentum/param.disruption_off_axis.in @@ -22,7 +22,7 @@ DISCARD_OUT discard.disruption_off_axis.out EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.escape.in b/examples/symba_energy_momentum/param.escape.in index 2b84eb719..99d572b75 100644 --- a/examples/symba_energy_momentum/param.escape.in +++ b/examples/symba_energy_momentum/param.escape.in @@ -19,13 +19,11 @@ CHK_RMIN 0.00465047 ! check for close solar encounters in AU CHK_RMAX 10000.0 ! discard outside of CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -!CHK_QMIN_COORD HELIO ! commented out here -!CHK_QMIN_RANGE 1.0 1000.0 ! commented out here ENC_OUT enc.escape.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT no ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index 65365b120..5e26e4cd3 100644 --- a/examples/symba_energy_momentum/param.sun.in +++ b/examples/symba_energy_momentum/param.sun.in @@ -19,13 +19,11 @@ CHK_RMIN 0.005 CHK_RMAX 1e2 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -!CHK_QMIN_COORD HELIO ! commented out here -!CHK_QMIN_RANGE 1.0 1000.0 ! commented out here ENC_OUT enc.escape.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT no ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_headon.in b/examples/symba_energy_momentum/param.supercatastrophic_headon.in index 19b15de7b..3ba223ad9 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_headon.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_headon.in @@ -17,13 +17,11 @@ CHK_RMIN 0.005 CHK_RMAX 1e6 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -!CHK_QMIN_COORD HELIO ! commented out here -!CHK_QMIN_RANGE 1.0 1000.0 ! commented out here ENC_OUT enc.supercatastrophic_headon.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in index 9cd214534..d033c76f1 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_off_axis.in @@ -23,7 +23,7 @@ ENC_OUT enc.supercatastrophic_off_axis.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -MTINY 1.0e-16 +GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 TU2S 3.1556925e7 diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py index 245f5fae0..6547c2802 100755 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py @@ -180,7 +180,7 @@ print(f'DU2M {DU2M}') print(f'TU2S {TU2S}') print(f'RHILL_PRESENT yes') -print(f'MTINY 1e-12') +print(f'GMTINY 1e-12') diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in index 3050dea4a..c69ee07f9 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in @@ -30,4 +30,4 @@ MU2KG 1.988409870698051e+30 DU2M 149597870700.0 TU2S 31557600.0 RHILL_PRESENT yes -MTINY 1e-12 +GMTINY 1e-12 diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py index 6b9664542..ac55abb2b 100755 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py @@ -178,7 +178,7 @@ print(f'DU2M {DU2M}') print(f'TU2S {TU2S}') print(f'RHILL_PRESENT yes') -print(f'MTINY 1e-12') +print(f'GMTINY 1e-12') diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in index 32cfb89bd..9fb0bf743 100644 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in @@ -29,4 +29,4 @@ MU2KG 1.988409870698051e+30 DU2M 149597870700.0 TU2S 31557600.0 RHILL_PRESENT yes -MTINY 1e-12 +GMTINY 1e-12 diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py b/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py index 707c80b51..546dcbb88 100755 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/init_cond.py @@ -39,7 +39,7 @@ sim.param['GR'] = 'NO' sim.param['CHK_CLOSE'] = 'YES' sim.param['RHILL_PRESENT'] = 'YES' -sim.param['MTINY'] = 1.0e-12 +sim.param['GMTINY'] = 1.0e-12 sim.param['MU2KG'] = swiftest.MSun sim.param['TU2S'] = swiftest.JD2S diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in index c72e9f4b4..f26f33592 100644 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in @@ -32,4 +32,4 @@ ROTATION NO TIDES NO ENERGY NO GR NO -MTINY 1e-12 +GMTINY 1e-12 diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index c0cd7d61f..81840ba05 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -76,8 +76,8 @@ def read_swiftest_param(param_file_name, param): param['TIDES'] = param['TIDES'].upper() param['ENERGY'] = param['ENERGY'].upper() param['GR'] = param['GR'].upper() - if 'MTINY' in param: - param['MTINY'] = real2float(param['MTINY']) + if 'GMTINY' in param: + param['GMTINY'] = real2float(param['GMTINY']) except IOError: print(f"{param_file_name} not found.") return param @@ -900,7 +900,7 @@ def swift2swifter(swift_param, plname="", tpname="", conversion_questions={}): intxt = conversion_questions.get('BIG_DISCARD', None) if not intxt: - intxt = input("BIG_DISCARD: include data for all bodies > MTINY for each discard record? (y/N)> ") + intxt = input("BIG_DISCARD: include data for all bodies > GMTINY for each discard record? (y/N)> ") if intxt.upper() == 'Y': swifter_param['BIG_DISCARD'] = 'YES' else: @@ -1215,11 +1215,11 @@ def swifter2swiftest(swifter_param, plname="", tpname="", cbname="", conversion_ print(f"Cannot write to file {swiftest_param['CB_IN']}") return swifter_param - MTINY = conversion_questions.get('MTINY', None) - if not MTINY: - MTINY = input(f"Value of MTINY if this is a SyMBA simulation (enter nothing if this is not a SyMBA parameter file)> ") - if MTINY != '' and real2float(MTINY.strip()) > 0: - swiftest_param['MTINY'] = real2float(MTINY.strip()) + GMTINY = conversion_questions.get('GMTINY', None) + if not GMTINY: + GMTINY = input(f"Value of GMTINY if this is a SyMBA simulation (enter nothing if this is not a SyMBA parameter file)> ") + if GMTINY != '' and real2float(GMTINY.strip()) > 0: + swiftest_param['GMTINY'] = real2float(GMTINY.strip()) # Remove the unneeded parameters if 'C' in swiftest_param: @@ -1265,7 +1265,7 @@ def swift2swiftest(swift_param, plname="", tpname="", cbname="", conversion_ques def swiftest2swifter_param(swiftest_param, J2=0.0, J4=0.0): swifter_param = swiftest_param CBIN = swifter_param.pop("CB_IN", None) - MTINY = swifter_param.pop("MTINY", None) + GMTINY = swifter_param.pop("GMTINY", None) DISCARD_OUT = swifter_param.pop("DISCARD_OUT", None) MU2KG = swifter_param.pop("MU2KG", 1.0) DU2M = swifter_param.pop("DU2M", 1.0) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 460060183..c0cfd3c97 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -165,10 +165,9 @@ subroutine set_scale_factors() vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot ! Set scale factors - !! Because of the implied G, mass is actually G*mass with units of distance**3 / time**2 Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) rscale = sum(radius(:)) - mscale = sqrt(Escale * rscale) + mscale = sqrt(Escale * rscale / param%GU) vscale = sqrt(Escale / mscale) tscale = rscale / vscale Lscale = mscale * rscale * vscale @@ -346,6 +345,7 @@ subroutine calculate_system_energy(linclude_fragments) if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) + tmpsys%pl%Gmass(npl+1:npl_new) = param%GU * m_frag(1:nfrag) tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag) @@ -369,7 +369,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - plwksp%nplm = count(plwksp%Gmass > param%mtiny / mscale) + plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale) end select end select call tmpsys%pl%eucl_index() @@ -381,7 +381,7 @@ subroutine calculate_system_energy(linclude_fragments) class is (symba_pl) select type(param) class is (symba_parameters) - nplm = count(pl%mass > param%mtiny) + nplm = count(pl%Gmass > param%Gmtiny) end select end select if (lk_plpl) call pl%eucl_index() @@ -836,7 +836,7 @@ end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -857,7 +857,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v ! Arguments integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision ! Constants @@ -931,7 +931,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v Qloss = 0.0_DP U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 - if ((m1 < mtiny).or.(m2 < mtiny)) then + if ((m1 < Gmtiny).or.(m2 < Gmtiny)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime Mlr = Mtot Mslr = 0.0_DP diff --git a/src/io/io.f90 b/src/io/io.f90 index 05313275c..ebfe2b1ef 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -493,7 +493,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) read(param_value, *) param%Ecollisions case("EUNTRACKED") read(param_value, *) param%Euntracked - case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 25c18295a..7dc70a37f 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -470,11 +470,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, real(DP), intent(inout) :: Qloss end subroutine fragmentation_initialize - module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss) implicit none integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision end subroutine fragmentation_regime diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index ed495ecfa..dfe1c0326 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -19,7 +19,7 @@ module symba_classes type, extends(swiftest_parameters) :: symba_parameters character(STRMAX) :: particle_file = PARTICLE_OUTFILE !! Name of output particle information file - real(DP) :: MTINY = -1.0_DP !! Smallest mass that is fully gravitating + real(DP) :: GMTINY = -1.0_DP !! Smallest mass that is fully gravitating integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. contains @@ -73,9 +73,9 @@ module symba_classes type, extends(helio_pl) :: symba_pl logical, dimension(:), allocatable :: lcollision !! flag indicating whether body has merged with another this time step logical, dimension(:), allocatable :: lencounter !! flag indicating whether body is part of an encounter this time step - logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the MTINY cutoff value - integer(I4B) :: nplm !! number of bodies above the MTINY limit - integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above MTINY) comparisons in the flattened upper triangular matrix + logical, dimension(:), allocatable :: lmtiny !! flag indicating whether this body is below the GMTINY cutoff value + integer(I4B) :: nplm !! number of bodies above the GMTINY limit + integer(I8B) :: nplplm !! Number of body (all massive)-body (only those above GMTINY) comparisons in the flattened upper triangular matrix integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with other planets this time step integer(I4B), dimension(:), allocatable :: ntpenc !! number of encounters with test particles this time step integer(I4B), dimension(:), allocatable :: levelg !! level at which this body should be moved diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 952d59709..1910411b9 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -22,7 +22,7 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec logical, dimension(:), allocatable :: lcollision, lmask real(DP), dimension(NDIM) :: xr, vr integer(I4B) :: k - real(DP) :: rlim, mtot + real(DP) :: rlim, Gmtot logical :: isplpl if (self%nenc == 0) return @@ -55,8 +55,8 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec xr(:) = pl%xh(:, ind1(k)) - pl%xh(:, ind2(k)) vr(:) = pl%vb(:, ind1(k)) - pl%vb(:, ind2(k)) rlim = pl%radius(ind1(k)) + pl%radius(ind2(k)) - mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, self%lvdotr(k)) + Gmtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k)) end do else do concurrent(k = 1:nenc, lmask(k)) @@ -359,7 +359,7 @@ module subroutine symba_collision_make_family_pl(self, idx) p2 = pl%kin(idx(2))%parent if (p1 == p2) return ! This is a collision between to children of a shared parent. We will ignore it. - if (pl%Gmass(p1) > pl%Gmass(p2)) then + if (pl%mass(p1) > pl%mass(p2)) then index_parent = p1 index_child = p2 else @@ -449,7 +449,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children Mcb_si = cb%mass * param%MU2KG - mtiny_si = (param%MTINY / param%GU) * param%MU2KG + mtiny_si = (param%GMTINY / param%GU) * param%MU2KG mass_res(:) = 0.0_DP @@ -463,8 +463,8 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) mass_res(1) = min(max(mlr, 0.0_DP), mtot) mass_res(2) = min(max(mslr, 0.0_DP), mtot) mass_res(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - mass_res(:) = (mass_res(:) / param%MU2KG) * param%GU - Qloss = Qloss * (param%GU / param%MU2KG) * (param%TU2S / param%DU2M)**2 + mass_res(:) = (mass_res(:) / param%MU2KG) + Qloss = Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG select case (regime) case (COLLRESOLVE_REGIME_DISRUPTION) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 4401fb5db..ba972db8b 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -71,8 +71,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms case ("FRAGMENTATION") call io_toupper(param_value) if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. - case ("MTINY") - read(param_value, *) param%mtiny + case ("GMTINY") + read(param_value, *) param%Gmtiny case("SEED") read(param_value, *) nseeds_from_file ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different @@ -111,12 +111,12 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) end if - if (self%mtiny < 0.0_DP) then - write(iomsg,*) "MTINY invalid or not set: ", self%mtiny + if (self%Gmtiny < 0.0_DP) then + write(iomsg,*) "GMTINY invalid or not set: ", self%Gmtiny iostat = -1 return else - write(*,*) "MTINY = ", self%mtiny + write(*,*) "GMTINY = ", self%Gmtiny end if if (.not.self%lclose) then @@ -167,7 +167,7 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms ! Special handling is required for writing the random number seed array as its size is not known until runtime ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) - write(param_name, Afmt) "MTINY"; write(param_value, Rfmt) param%mtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%Gmtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%lfragmentation) then write(param_name, Afmt) "SEED" diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 021873a70..ab8b5543e 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -25,7 +25,7 @@ module subroutine symba_setup_initialize_system(self, param) call pl%sort("mass", ascending=.false.) select type(param) class is (symba_parameters) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) end select end select diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 2ebc11ebb..2ab088d02 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -394,7 +394,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then call pl%sort("mass", ascending=.false.) - pl%lmtiny(:) = pl%Gmass(:) > param%MTINY + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) call pl%eucl_index() end if diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index b004f84b7..fd331bcc3 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -83,7 +83,7 @@ module subroutine util_get_energy_momentum_system(self, param) ! Do the potential energy between pairs of massive bodies do k = 1, pl%nplpl associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k)) - pepl(k) = -pl%mass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) + pepl(k) = -param%GU * pl%mass(ik) * pl%mass(jk) / norm2(pl%xh(:, jk) - pl%xh(:, ik)) lstatpl(k) = (lstatus(ik) .and. lstatus(jk)) end associate end do