Restructured code so that the eucl method falls under util and is now…

… accessed using %index. Now include sorting in that so that pl arrays are always sorted prior to indexing, rather than remembering to do it each time.
MintonGroup · Aug 20, 2021 · 5e80def · 5e80def
1 parent 728cacf
commit 5e80def
Show file tree

Hide file tree

Showing 31 changed files with 1,627 additions and 1,585 deletions.
diff --git a/Makefile b/Makefile
@@ -91,11 +91,6 @@ lib:
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make libdir
-	cd $(SWIFTEST_HOME)/src/eucl; \
-	  rm -f Makefile.Defines Makefile; \
-	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
-	  ln -s $(SWIFTEST_HOME)/Makefile .; \
-	  make libdir
 	cd $(SWIFTEST_HOME)/src/fragmentation; \
 	  rm -f Makefile.Defines Makefile; \
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
@@ -193,7 +188,6 @@ clean:
 	cd $(SWIFTEST_HOME)/src/modules;  rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/discard; rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/drift;   rm -f Makefile.Defines Makefile *.gc*
-	cd $(SWIFTEST_HOME)/src/eucl;   rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/fragmentation;   rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/gr;          rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/helio;    rm -f Makefile.Defines Makefile *.gc*

diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb
@@ -81,8 +81,8 @@
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x2b8ff43aafd0>,\n",
-       " <matplotlib.lines.Line2D at 0x2b8ff485d650>]"
+       "[<matplotlib.lines.Line2D at 0x2b105fdb3b90>,\n",
+       " <matplotlib.lines.Line2D at 0x2b105fd97b90>]"
       ]
      },
      "execution_count": 6,

diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swifter.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swifter.in
@@ -1,9 +1,9 @@
 ! VERSION        Swifter parameter file converted from Swiftest
 T0               0.0
-TSTOP            36525.0
+TSTOP            9500.0
 DT               1.0
-ISTEP_OUT        10
-ISTEP_DUMP       10
+ISTEP_OUT        100
+ISTEP_DUMP       100
 OUT_FORM         XV
 OUT_TYPE         REAL8
 OUT_STAT         UNKNOWN

diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swiftest.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swiftest.in
@@ -1,9 +1,9 @@
 ! VERSION        Swiftest parameter input
 T0               0.0
-TSTOP            36525.0
+TSTOP            9500.0
 DT               1.0
-ISTEP_OUT        10
-ISTEP_DUMP       10
+ISTEP_OUT        100
+ISTEP_DUMP       100
 OUT_FORM         XV
 OUT_TYPE         REAL8
 OUT_STAT         UNKNOWN

diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb
diff --git a/examples/symba_chambers_2013/.gitignore b/examples/symba_chambers_2013/.gitignore
@@ -0,0 +1 @@
+savestate/*
diff --git a/examples/symba_chambers_2013/savestate2/param.in b/examples/symba_chambers_2013/savestate2/param.in
@@ -0,0 +1,32 @@
+!
+! Parameter file for Chambers 2013 in units of Msun, AU, year
+!
+T0             0.0e0
+TSTOP          1e8        ! simulation length in years
+DT             0.016              ! stepsize in years
+CB_IN          sun_MsunAUYR.in
+PL_IN          pl_chambers_2013.in
+TP_IN          tp.in
+IN_TYPE        ASCII
+ISTEP_OUT      625              ! output cadence 
+ISTEP_DUMP     625             ! system dump cadence
+BIN_OUT        bin.dat
+PARTICLE_OUT   particle.dat
+OUT_TYPE       REAL8              ! double precision real output
+OUT_FORM       EL                 ! osculating element output
+OUT_STAT       REPLACE
+CHK_CLOSE      yes                ! check for planetary close encounters
+CHK_RMAX       100000.0            ! discard outside of
+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
+MU2KG          1.98847e30         ! (M_sun-> kg)
+DU2M           1.495979e11        ! distance unit to meters (AU --> m)
+TU2S           3.15569259747e7    ! time unit to seconds (years --> seconds)
+GMTINY         1e-10		  
+ENERGY         yes
+ENERGY_OUT     energy.dat
+ROTATION       yes
+FRAGMENTATION  yes
+DISCARD_OUT    discard.out
+SEED           8 12261555 871132 92734722 21132443 36344777 4334443 219291656 3848566
diff --git a/examples/symba_chambers_2013/savestate2/pl_chambers_2013.in b/examples/symba_chambers_2013/savestate2/pl_chambers_2013.in
diff --git a/examples/symba_chambers_2013/savestate2/sun_MsunAUYR.in b/examples/symba_chambers_2013/savestate2/sun_MsunAUYR.in
@@ -0,0 +1,7 @@
+0                                     ! id
+39.476926408897626                    ! G*Mass
+0.004650467260962157                  ! Radius
+0.0 !4.7535806948127355e-12                ! J2
+0.0 !-2.2473967953572827e-18               ! J4
+0.0 0.0 0.07                          ! Principle axes moments of inertia
+11.2093063  -38.75937204  82.25088158 ! Rotation vector (rad/TU)
diff --git a/examples/symba_chambers_2013/savestate2/tp.in b/examples/symba_chambers_2013/savestate2/tp.in
@@ -0,0 +1 @@
+0
diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb
diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/swiftest_symba_vs_swifter_symba.ipynb b/examples/symba_swifter_comparison/8pl_16tp_encounters/swiftest_symba_vs_swifter_symba.ipynb
diff --git a/examples/symba_swifter_comparison/mars_disk/param.swifter.in b/examples/symba_swifter_comparison/mars_disk/param.swifter.in
@@ -1,5 +1,5 @@
 T0             0.0 
-TSTOP          3000.0 
+TSTOP          4200.0 
 DT             600.0
 PL_IN          pl.swifter.in
 TP_IN          tp.in

diff --git a/examples/symba_swifter_comparison/mars_disk/param.swiftest.in b/examples/symba_swifter_comparison/mars_disk/param.swiftest.in
@@ -1,6 +1,6 @@
 !Parameter file for the SyMBA-RINGMOONS test
 T0             0.0 
-TSTOP          3000.0 
+TSTOP          4200.0 
 DT             600.0
 PL_IN          pl.swiftest.in
 TP_IN          tp.in
@@ -24,7 +24,7 @@ BIG_DISCARD    no
 RHILL_PRESENT  yes
 GMTINY         1000.0
 ENERGY         yes
-FRAGMENTATION  yes 
+FRAGMENTATION  no  
 ROTATION       yes
 MU2KG          1.0
 DU2M           1.0

diff --git a/examples/symba_swifter_comparison/mars_disk/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/mars_disk/swiftest_vs_swifter.ipynb
diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90
@@ -432,10 +432,6 @@ subroutine calculate_system_energy(linclude_fragments)
             logical, dimension(:), allocatable :: lexclude_tmp
             logical :: lk_plpl
 
-            ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the
-            ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by
-            ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this
-            ! subroutine should be called relatively infrequently, it shouldn't matter too much.
 
             ! Build the internal planet list out of the non-excluded bodies and optionally with fragments appended. This
             ! will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it
@@ -447,34 +443,23 @@ subroutine calculate_system_energy(linclude_fragments)
                elsewhere
                   tmpsys%pl%status(1:npl_new) = ACTIVE
                end where
-
-               select type(plwksp => tmpsys%pl)
-               class is (symba_pl)
-                  select type(param)
-                  class is (symba_parameters)
-                     if (linclude_fragments) call add_fragments_to_tmpsys() ! Adds or updates the fragment properties to their current values
-                     plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale)
-                  end select
-               end select
-
+
+
+               ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the
+               ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by
+               ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this
+               ! subroutine should be called relatively infrequently, it shouldn't matter too much.
                lk_plpl = allocated(pl%k_plpl)
                if (lk_plpl) deallocate(pl%k_plpl)
 
-               call tmpsys%pl%eucl_index()
+               call tmpsys%pl%index(param)
 
                call tmpsys%get_energy_and_momentum(param) 
 
                ! Restore the big array
                deallocate(tmpsys%pl%k_plpl) 
-               select type(pl)
-               class is (symba_pl)
-                  select type(param)
-                  class is (symba_parameters)
-                     pl%nplm = count(pl%Gmass > param%Gmtiny)
-                  end select
-               end select
-
-               if (lk_plpl) call pl%eucl_index()
+
+               if (lk_plpl) call pl%index(param)
 
                ! Calculate the current fragment energy and momentum balances
                if (linclude_fragments) then

diff --git a/src/helio/helio_util.f90 b/src/helio/helio_util.f90
@@ -0,0 +1,21 @@
+submodule(helio_classes) s_helio_eucl
+   use swiftest
+contains
+
+   module subroutine helio_util_index_eucl_plpl(self, param)
+      !! author: David A. Minton
+      !!
+      !! Wrapper for the indexing method for WHM massive bodies. Sorts the massive bodies by heliocentric distance and then flattens the pl-pl upper triangular matrix
+      implicit none
+      ! Arguments
+      class(helio_pl),            intent(inout) :: self  !! Helio massive body object
+      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+
+      call self%sort("mass", ascending=.false.)
+      call util_index_eucl_plpl(self, param)
+
+      return
+   end subroutine helio_util_index_eucl_plpl
+
+end submodule s_helio_eucl
+
diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90
@@ -38,6 +38,7 @@ module helio_classes
    contains
       procedure :: drift       => helio_drift_pl           !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure :: lindrift    => helio_drift_linear_pl    !! Method for linear drift of massive bodies due to barycentric momentum of Sun
+      procedure :: index       => helio_util_index_eucl_plpl   !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
       procedure :: accel_gr    => helio_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction
       procedure :: gr_pos_kick => helio_gr_p4_pl           !! Position kick due to p**4 term in the post-Newtonian correction
       procedure :: accel       => helio_kick_getacch_pl    !! Compute heliocentric accelerations of massive bodies
@@ -212,6 +213,12 @@ module subroutine helio_step_tp(self, system, param, t, dt)
          real(DP),                     intent(in)    :: t      !! Current simulation time
          real(DP),                     intent(in)    :: dt     !! Stepsizee
       end subroutine helio_step_tp
+
+      module subroutine helio_util_index_eucl_plpl(self, param)
+         implicit none
+         class(helio_pl),            intent(inout) :: self  !! Helio massive body object
+         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+      end subroutine helio_util_index_eucl_plpl
    end interface
 
 end module helio_classes
diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
@@ -216,7 +216,7 @@ module swiftest_classes
       ! Massive body-specific concrete methods 
       ! These are concrete because they are the same implemenation for all integrators
       procedure :: discard      => discard_pl             !! Placeholder method for discarding massive bodies 
-      procedure :: eucl_index   => eucl_dist_index_plpl   !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
+      procedure :: index        => util_index_eucl_plpl   !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
       procedure :: accel_int    => kick_getacch_int_pl    !! Compute direct cross (third) term heliocentric accelerations of massive bodies
       procedure :: accel_obl    => obl_acc_pl             !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
       procedure :: setup        => setup_pl               !! A base constructor that sets the number of bodies and allocates and initializes all arrays  
@@ -452,9 +452,10 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag
          integer(I4B), intent(out)      :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
       end subroutine drift_one
 
-      module subroutine eucl_dist_index_plpl(self)
+      module subroutine util_index_eucl_plpl(self, param)
          implicit none
-         class(swiftest_pl), intent(inout) :: self  !! Swiftest massive body object
+         class(swiftest_pl),         intent(inout) :: self  !! Swiftest massive body object
+         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
       end subroutine
 
       module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, &

diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90
@@ -84,6 +84,7 @@ module symba_classes
       type(symba_particle_info), dimension(:), allocatable :: info
    contains
       procedure :: make_family     => symba_collision_make_family_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family
+      procedure :: index           => symba_util_index_eucl_plpl     !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
       procedure :: discard         => symba_discard_pl               !! Process massive body discards
       procedure :: drift           => symba_drift_pl                 !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
       procedure :: encounter_check => symba_encounter_check_pl       !! Checks if massive bodies are going through close encounters with each other
@@ -352,6 +353,12 @@ module function symba_collision_casesupercatastrophic(system, param, family, x,
          integer(I4B)                                   :: status           !! Status flag assigned to this outcome
       end function symba_collision_casesupercatastrophic
 
+      module subroutine symba_util_index_eucl_plpl(self, param)
+         implicit none
+         class(symba_pl),            intent(inout) :: self  !! SyMBA massive body object
+         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+      end subroutine symba_util_index_eucl_plpl
+
       module subroutine symba_io_write_discard(self, param)
          use swiftest_classes, only : swiftest_parameters
          implicit none

diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90
@@ -34,6 +34,7 @@ module whm_classes
       procedure :: j2h         => whm_coord_j2h_pl           !! Convert position and velcoity vectors from Jacobi to helliocentric coordinates 
       procedure :: vh2vj       => whm_coord_vh2vj_pl         !! Convert velocity vectors from heliocentric to Jacobi coordinates 
       procedure :: drift       => whm_drift_pl               !! Loop through massive bodies and call Danby drift routine to jacobi coordinates
+      procedure :: index       => whm_util_index_eucl_plpl   !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
       procedure :: accel_gr    => whm_gr_kick_getacch_pl     !! Acceleration term arising from the post-Newtonian correction
       procedure :: gr_pos_kick => whm_gr_p4_pl               !! Position kick due to p**4 term in the post-Newtonian correction
       procedure :: accel       => whm_kick_getacch_pl        !! Compute heliocentric accelerations of massive bodies
@@ -108,6 +109,13 @@ module subroutine whm_drift_pl(self, system, param, dt)
          real(DP),                     intent(in)    :: dt     !! Stepsize
       end subroutine whm_drift_pl
 
+      module subroutine whm_util_index_eucl_plpl(self, param)
+         use swiftest_classes, only : swiftest_parameters
+         implicit none
+         class(whm_pl),            intent(inout) :: self  !! WHM massive body object
+         class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters
+      end subroutine whm_util_index_eucl_plpl
+
       !> Get heliocentric accelration of massive bodies
       module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg)
          use swiftest_classes, only : swiftest_cb, swiftest_parameters

diff --git a/src/python/orbel.f90 → src/python_bindings/orbel.f90 b/src/python/orbel.f90 → src/python_bindings/orbel.f90
@@ -4,7 +4,7 @@ module orbel
     public :: xv2el
 contains
     pure elemental subroutine xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
-        use module_interfaces, only : orbel_xv2el
+        use swiftest_classes, only : orbel_xv2el
         implicit none
         ! Arguments
         real*8, intent(in) :: mu, px, py, pz, vx, vy, vz
@@ -26,7 +26,7 @@ pure elemental subroutine xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, om
         real*8, dimension(3) :: x, v
         x = [px, py, pz]
         v = [vx, vy, vz]
-        call orbel_xv2el(x, v, mu, a, e, inc, capom, omega, capm)
+        call orbel_xv2el(mu, x(:), v(:), a, e, inc, capom, omega, capm)
         return
     end subroutine xv2el
 end module orbel
diff --git a/src/setup/setup.f90 b/src/setup/setup.f90
@@ -145,7 +145,7 @@ module subroutine setup_initialize_system(self, param)
          call self%pl%el2xv(self%cb)
          call self%tp%el2xv(self%cb)
       end if
-      call self%pl%eucl_index()
+      call self%pl%index(param)
       if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb)
       self%pl%lfirst = param%lfirstkick
       self%tp%lfirst = param%lfirstkick

diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90
@@ -901,7 +901,6 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
                !Copy over or set integration parameters for new bodies
                plnew%lcollision(1:nfrag) = .false.
                plnew%ldiscard(1:nfrag) = .false.
-               plnew%lmtiny(1:nfrag) = plnew%Gmass(1:nfrag) > param%GMTINY
                plnew%levelg(1:nfrag) = pl%levelg(ibiggest)
                plnew%levelm(1:nfrag) = pl%levelm(ibiggest)
 

diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90
@@ -77,19 +77,13 @@ module subroutine symba_setup_initialize_system(self, param)
          call system%pltpenc_list%setup(0)
          call system%plplenc_list%setup(0)
          call system%plplcollision_list%setup(0)
-         select type(pl => system%pl)
-         class is (symba_pl)
-            call pl%sort("mass", ascending=.false.)
-            select type(param)
-            class is (symba_parameters)
-               pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
-               pl%nplm = count(pl%lmtiny(:))
-               if (param%lrestart) then
-                  call symba_io_read_particle(system, param)
-               else
-                  call symba_setup_initialize_particle_info(system, param) 
-               end if
-            end select
+         select type(param)
+         class is (symba_parameters)
+            if (param%lrestart) then
+               call symba_io_read_particle(system, param)
+            else
+               call symba_setup_initialize_particle_info(system, param) 
+            end if
          end select
       end associate