Merge branch 'debug'

MintonGroup · Aug 20, 2021 · f7a2ff0 · f7a2ff0
2 parents a2e0f4e + a621047
commit f7a2ff0
Show file tree

Hide file tree

Showing 29 changed files with 699 additions and 1,620 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/Makefile.Defines b/Makefile.Defines
@@ -54,24 +54,26 @@ VTUNE_FLAGS = -g -O2 -vec -simd -shared-intel -qopenmp -debug inline-debug-info
 IDEBUG = -O0 -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -init=snan,arrays 
 STRICTREAL = -fp-model strict -fp-model no-except -prec-div -prec-sqrt -assume protect-parens 
 SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -prec-div -prec-sqrt -assume protect-parens
-PAR = -qopenmp #-parallel  Something goes wrong in SyMBA at the moment with auto-paralellization enabled
+PAR = -qopenmp #-parallel  #Something goes wrong in SyMBA at the moment with auto-paralellization enabled
 HEAPARR = -heap-arrays 1048576 
 OPTREPORT = -qopt-report=5 
+IPRODUCTION = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(PAR) $(SIMDVEC) $(HEAPARR) 
 
 #gfortran flags
-GDEBUG = -g -Og -fbacktrace -fbounds-check 
-GPRODUCTION = -O3
+GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none
 GPAR = -fopenmp #-ftree-parallelize-loops=4
 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak 
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
+GPRODUCTION = -O3 -ffree-line-length-none $(GPAR)
 
 #FFLAGS  = $(IDEBUG) $(HEAPARR) $(SIMDVEC) $(PAR)
-#FFLAGS  = -init=snan,arrays -no-wrap-margin -O0 -g -traceback $(STRICTREAL) $(PAR) $(SIMDVEC) $(HEAPARR) 
+#FFLAGS  = $(IPRODUCTION)
 #FORTRAN = ifort
 #AR     = xiar
 
 FORTRAN = gfortran
-FFLAGS =  -ffree-line-length-none $(GPAR) $(GPRODUCTION) #$(GMEM)
+#FFLAGS =  $(GDEBUG) #$(GMEM) #$(GPAR)
+FFLAGS =  $(GPRODUCTION)
 AR     = ar
 
 # DO NOT include in CFLAGS the "-c" option to compile object only

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_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,8 @@ BIG_DISCARD    no
 RHILL_PRESENT  yes
 GMTINY         1000.0
 ENERGY         yes
-FRAGMENTATION  yes 
+ENERGY_OUT     energy.dat
+FRAGMENTATION  yes
 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
@@ -100,7 +100,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
 
       call reset_fragments()
       call define_coordinate_system()
-      call construct_temporary_system()
 
       ! Calculate the initial energy of the system without the collisional family
       call calculate_system_energy(linclude_fragments=.false.)
@@ -150,25 +149,25 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       end do
       call restore_scale_factors()
 
-      write(*,        "(' -------------------------------------------------------------------------------------')")
-      write(*,        "('  Final diagnostic')")
-      write(*,        "(' -------------------------------------------------------------------------------------')")
-      call calculate_system_energy(linclude_fragments=.true.)
-      if (lfailure) then
-         write(*,*) "symba_frag_pos failed after: ",try," tries"
-         do ii = 1, nfrag
-            vb_frag(:, ii) = vcom(:)
-         end do
-      else
-         write(*,*) "symba_frag_pos succeeded after: ",try," tries"
-         write(*,        "(' dL_tot should be very small' )")
-         write(*,fmtlabel) ' dL_tot      |', dLmag / Lmag_before
-         write(*,        "(' dE_tot should be negative and equal to Qloss' )")
-         write(*,fmtlabel) ' dE_tot      |', dEtot / abs(Etot_before)
-         write(*,fmtlabel) ' Qloss       |', -Qloss / abs(Etot_before)
-         write(*,fmtlabel) ' dE - Qloss  |', (Etot_after - Etot_before + Qloss) / abs(Etot_before)
-      end if
-      write(*,        "(' -------------------------------------------------------------------------------------')")
+      ! write(*,        "(' -------------------------------------------------------------------------------------')")
+      ! write(*,        "('  Final diagnostic')")
+      ! write(*,        "(' -------------------------------------------------------------------------------------')")
+      ! call calculate_system_energy(linclude_fragments=.true.)
+      ! if (lfailure) then
+      !    write(*,*) "symba_frag_pos failed after: ",try," tries"
+      !    do ii = 1, nfrag
+      !       vb_frag(:, ii) = vcom(:)
+      !    end do
+      ! else
+      !    write(*,*) "symba_frag_pos succeeded after: ",try," tries"
+      !    write(*,        "(' dL_tot should be very small' )")
+      !    write(*,fmtlabel) ' dL_tot      |', dLmag / Lmag_before
+      !    write(*,        "(' dE_tot should be negative and equal to Qloss' )")
+      !    write(*,fmtlabel) ' dE_tot      |', dEtot / abs(Etot_before)
+      !    write(*,fmtlabel) ' Qloss       |', -Qloss / abs(Etot_before)
+      !    write(*,fmtlabel) ' dE - Qloss  |', (Etot_after - Etot_before + Qloss) / abs(Etot_before)
+      ! end if
+      ! write(*,        "(' -------------------------------------------------------------------------------------')")
 
       call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)  ! Save the current halting modes so we can turn them off temporarily
 
@@ -355,7 +354,7 @@ subroutine construct_temporary_system()
 
             associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb)
                if (size(lexclude) /= npl + nfrag) then 
-                  allocate(lexclude_tmp(npl_new))
+                  allocate(lexclude_tmp(npl + nfrag))
                   lexclude_tmp(1:npl) = lexclude(1:npl)
                   call move_alloc(lexclude_tmp, lexclude)
                end if
@@ -364,8 +363,9 @@ subroutine construct_temporary_system()
                elsewhere
                   lexclude(1:npl) = .false. 
                end where
-               lexclude(npl+1:npl_new) = .true.
+               lexclude(npl+1:(npl + nfrag)) = .true.
                if (allocated(tmpparam)) deallocate(tmpparam) 
+               if (allocated(tmpsys)) deallocate(tmpsys)
                allocate(tmpparam, source=param)
                call setup_construct_system(tmpsys, param)
                call tmpsys%tp%setup(0, param)
@@ -388,8 +388,11 @@ subroutine add_fragments_to_tmpsys()
             implicit none
             ! Internals
             integer(I4B) :: i
+            class(swiftest_pl), allocatable :: pl_discards
+            logical, dimension(:), allocatable :: lexclude_tmp
 
             associate(pl => system%pl, npl => system%pl%nbody)
+               npl_new = npl + nfrag
 
                tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag)
                tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * tmpparam%GU
@@ -407,11 +410,20 @@ subroutine add_fragments_to_tmpsys()
                ! Disable the collisional family for subsequent energy calculations and coordinate shifts
                lexclude(family(:)) = .true.
                lexclude(npl+1:npl_new) = .false.
-               where(lexclude(:)) 
-                  tmpsys%pl%status(:) = INACTIVE
+               where(lexclude(1:npl_new)) 
+                  tmpsys%pl%status(1:npl_new) = INACTIVE
                elsewhere
-                  tmpsys%pl%status(:) = ACTIVE
+                  tmpsys%pl%status(1:npl_new) = ACTIVE
                end where
+               allocate(pl_discards, mold=tmpsys%pl)
+               call tmpsys%pl%spill(pl_discards, lspill_list=lexclude(1:npl_new), ldestructive=.true.)
+               npl_new = count(.not.lexclude(:))
+
+               if (size(lexclude) /= npl_new) then 
+                  allocate(lexclude_tmp(npl_new))
+                  call move_alloc(lexclude_tmp, lexclude)
+               end if
+               lexclude(1:npl_new) = .false.
 
             end associate
 
@@ -432,49 +444,30 @@ 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
             ! is in the main program. This will temporarily expand the planet list in a temporary system object called tmpsys to feed it into symba_energy
             associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb)
 
-               where (lexclude(1:npl_new))
-                  tmpsys%pl%status(1:npl_new) = INACTIVE
-               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 construct_temporary_system()
+               if (linclude_fragments) call add_fragments_to_tmpsys()
+
+               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/io/io.f90 b/src/io/io.f90
@@ -905,6 +905,7 @@ module subroutine io_read_in_cb(self, param)
             read(iu, *, err = 667, iomsg = errmsg) self%Ip(1), self%Ip(2), self%Ip(3)
             read(iu, *, err = 667, iomsg = errmsg) self%rot(1), self%rot(2), self%rot(3)
          end if
+         ierr = 0
       else
          open(unit = iu, file = param%incbfile, status = 'old', form = 'UNFORMATTED', err = 667, iomsg = errmsg)
          ierr = self%read_frame(iu, param)

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,13 @@ 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)
+         use swiftest_classes, only : swiftest_parameters
+         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,13 @@ 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)
+         use swiftest_classes, only : swiftest_parameters
+         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