Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Refactored MTINY into GMTINY to be consistent in distinguishing G*mas…
Browse files Browse the repository at this point in the history
…s and mass terms
  • Loading branch information
daminton committed Aug 6, 2021
1 parent ac34918 commit 90faadf
Show file tree
Hide file tree
Showing 30 changed files with 72 additions and 80 deletions.
10 changes: 5 additions & 5 deletions docs/src/fragmentation.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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()
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion docs/src/io.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions docs/src/swiftest_classes.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 4 additions & 4 deletions docs/src/symba_classes.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion docs/src/symba_collision.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
12 changes: 6 additions & 6 deletions docs/src/symba_io.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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"
Expand Down
2 changes: 1 addition & 1 deletion docs/src/symba_setup.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion docs/src/symba_util.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 1 addition & 3 deletions examples/symba_energy_momentum/param.disruption_headon.in
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 1 addition & 3 deletions examples/symba_energy_momentum/param.escape.in
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 1 addition & 3 deletions examples/symba_energy_momentum/param.sun.in
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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')



Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -30,4 +30,4 @@ MU2KG 1.988409870698051e+30
DU2M 149597870700.0
TU2S 31557600.0
RHILL_PRESENT yes
MTINY 1e-12
GMTINY 1e-12
Original file line number Diff line number Diff line change
Expand Up @@ -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')



Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -29,4 +29,4 @@ MU2KG 1.988409870698051e+30
DU2M 149597870700.0
TU2S 31557600.0
RHILL_PRESENT yes
MTINY 1e-12
GMTINY 1e-12
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -32,4 +32,4 @@ ROTATION NO
TIDES NO
ENERGY NO
GR NO
MTINY 1e-12
GMTINY 1e-12
18 changes: 9 additions & 9 deletions python/swiftest/swiftest/io.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 90faadf

Please sign in to comment.