Merge branch 'debug' into OOPTides

MintonGroup · Aug 20, 2021 · 486531c · 486531c
2 parents 6ce85dc + 728cacf
commit 486531c
Show file tree

Hide file tree

Showing 11 changed files with 205 additions and 52 deletions.
diff --git a/Makefile.Defines b/Makefile.Defines
@@ -61,17 +61,17 @@ OPTREPORT = -qopt-report=5
 #gfortran flags
 GDEBUG = -g -Og -fbacktrace -fbounds-check 
 GPRODUCTION = -O3
-GPAR = -fopenmp -ftree-parallelize-loops=4
+GPAR = -fopenmp #-ftree-parallelize-loops=4
 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak 
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
 
 #FFLAGS  = $(IDEBUG) $(HEAPARR) $(SIMDVEC) $(PAR)
-FFLAGS  = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(PAR) $(SIMDVEC) $(HEAPARR) 
-FORTRAN = ifort
+#FFLAGS  = -init=snan,arrays -no-wrap-margin -O0 -g -traceback $(STRICTREAL) $(PAR) $(SIMDVEC) $(HEAPARR) 
+#FORTRAN = ifort
 #AR     = xiar
 
-#FORTRAN = gfortran
-#FFLAGS =  -ffree-line-length-none $(GPRODUCTION) $(GPAR) #$(GDEBUG) $(GMEM)
+FORTRAN = gfortran
+FFLAGS =  -ffree-line-length-none $(GPAR) $(GPRODUCTION) #$(GMEM)
 AR     = ar
 
 # DO NOT include in CFLAGS the "-c" option to compile object only

diff --git a/examples/in_form_el/cb.in b/examples/in_form_el/cb.in
@@ -0,0 +1,5 @@
+0
+0.0002959122081920778
+0.004650467260962157
+4.7535806948127355e-12
+-2.2473967953572827e-18
diff --git a/examples/in_form_el/param.in b/examples/in_form_el/param.in
@@ -0,0 +1,35 @@
+! VERSION        Swiftest parameter input from file param.in
+T0               0.0
+TSTOP            365250000000.0
+DT               5.0
+ISTEP_OUT        7305000
+ISTEP_DUMP       7305000
+OUT_FORM         EL
+OUT_TYPE         REAL8
+OUT_STAT         REPLACE
+IN_TYPE          ASCII
+PL_IN            pl.in
+TP_IN            tp.in
+CB_IN            cb.in
+BIN_OUT          bin.dat
+CHK_QMIN         0.004650467260962157
+CHK_RMIN         0.004650467260962157
+CHK_RMAX         1000.0
+CHK_EJECT        1000.0
+CHK_QMIN_COORD   HELIO
+CHK_QMIN_RANGE   0.004650467260962157 1000.0
+MU2KG            1.988409870698051e+30
+TU2S             86400.0
+DU2M             149597870700.0
+IN_FORM          EL
+ENC_OUT          enc.dat
+EXTRA_FORCE      NO
+DISCARD_OUT      discard.out
+BIG_DISCARD      NO
+CHK_CLOSE        YES
+RHILL_PRESENT    YES
+FRAGMENTATION    NO
+ROTATION         NO
+TIDES            NO
+ENERGY           NO
+GR               YES
diff --git a/examples/in_form_el/pl.in b/examples/in_form_el/pl.in
@@ -0,0 +1,33 @@
+8
+1 4.9125474498983625e-11 0.0014751256791533402
+1.6306381826061646e-05
+0.38710383731080894 0.24071190040066506 8.640070276952633
+30.683987164281632 57.949871930463694 151.08868778597147
+2 7.243452483873647e-10 0.0067591154147628296
+4.0453784346544176e-05
+0.7233315774753523 0.02498272794314777 2.4772002857938187
+125.76075677093282 6.742571718202682 350.2972657114627
+3 8.997011382166019e-10 0.010044745408888887
+4.25875607065041e-05
+0.9999904186618421 0.05232751342856539 2.0722848768764717
+113.84624589363827 312.12453079628943 148.9860363695549
+4 9.549535102761465e-11 0.007246758104365299
+2.2657408050928896e-05
+1.5236293140305446 0.005382832120155795 4.468216946001687
+325.6098790099478 152.40711385323957 110.04936074450407
+5 2.825345908631355e-07 0.3552714410085414
+0.0004673261703049093
+5.203095539549931 0.03063821583296288 1.3661365597246853
+113.82182400586191 303.52091573107697 96.03226696213453
+6 8.459715183006416e-08 0.43765112940564116
+0.00038925687730393614
+9.579329287028964 0.08182512713593486 2.54624964218005
+100.32672664745415 306.8523855342038 64.14978871000847
+7 1.2920249163736674e-08 0.4695173932519177
+0.00016953449859497232
+19.185526866180975 0.02725895554029442 0.6468816262831063
+75.84341829565291 203.1828812930829 337.91763726094723
+8 1.5243589003230834e-08 0.7812638143548571
+0.00016458790412449367
+30.033406636357004 0.014808837377897958 1.0172902231424474
+112.31047122167676 61.6982073585126 219.66063286932612
diff --git a/examples/in_form_el/tp.in b/examples/in_form_el/tp.in
@@ -0,0 +1,76 @@
+25
+1030101
+1.9298376430729873 0.012639757841142886 24.991702977377773
+341.74785371597466 248.1210122966632 133.62225322652563
+1030102
+1.9576788467130908 0.027970945719088838 25.55968189226752
+168.9126943096874 174.5978360268346 251.14179333855816
+1030103
+1.967626111006107 0.09226955476091864 1.5279774530000956
+343.6507149897438 15.943105536008861 203.03043240705532
+1030104
+1.9481905063928624 0.0011928078627614026 13.712706682327513
+65.70102045298484 272.5349689212673 248.3878548054673
+1030105
+1.9770133412534745 0.07760502570359305 21.113956523745326
+352.7420976555242 265.31463697173825 90.20915048805514
+1030106
+1.9218154288668319 0.024539743705313734 6.3122997131518
+96.26526838131817 245.40058043017729 203.1458881543752
+1030107
+1.8804043534140047 0.051960914508786316 5.641080857843637
+135.37074231372475 310.57579433118855 92.01889766237802
+1030108
+1.905574620628057 0.0094526649542939 27.795518144975087
+257.055868633706 262.18164024081847 174.33454380356292
+1030109
+1.9610371451706512 0.0642394856575368 11.238178270217231
+155.06613544539204 179.7982091248266 227.65351567825778
+1030110
+1.8442099887042422 0.06774466132861966 7.7755827872476635
+53.89722157548661 231.27772356054047 224.74925586092095
+1030111
+1.918614204096997 0.09311188233380624 22.656238093717935
+8.91787503082778 121.09392635435141 61.064615587866086
+1030112
+1.846083421031052 0.058817512099736215 23.17733139808657
+204.2210767808533 325.0096873835248 80.8104389106319
+1030113
+1.877211316979287 0.09545803603787674 22.75125883423326
+338.55088656488704 260.97440743659024 358.27993618149003
+1030114
+1.847699999560908 0.04422349734135026 24.400107433270254
+175.41066398268853 271.9213024922804 40.10860995923193
+1030115
+1.9693855972875431 0.01087613268562867 18.3688586644925
+288.07225817085316 351.2861239674453 138.17614312177193
+1030116
+1.8499142747340205 0.05785154439430367 18.680424125046468
+170.34777680935684 121.2501486818594 116.63207933400314
+1030117
+1.8484833703421926 0.027361686628703055 12.35340287511839
+77.72361263316543 170.96254103643864 153.80655599509745
+1030118
+1.8844066519945433 0.012253504795668403 28.11812027401767
+5.816089316051807 234.40531004184038 284.0432510837761
+1030119
+1.8321110668963587 0.05520820027904241 4.246848772229706
+195.98962278787414 196.26877791689589 296.01468603344676
+1030120
+1.9298252957677076 0.04051603959563857 11.38835718089586
+14.020062127093711 350.7131191785245 61.91415736777295
+1030121
+1.8296584872076656 0.02265081514507832 19.242848763922552
+170.21169721403766 255.88594273947075 131.4644380844818
+1030122
+1.8699958208259342 0.08248881864547655 6.775274873263883
+131.47697608382563 36.70815817760613 166.48187957273242
+1030123
+1.8427387543178881 0.07305783679736073 5.034586865380743
+311.1450141807796 214.40263993828856 332.6408661254218
+1030124
+1.8567511528637752 0.014548178017731906 8.946207161705699
+27.813385444122652 111.5675183218652 201.4816583454089
+1030125
+1.8846946839402372 0.06619567127431686 13.271555732235106
+271.1533808977903 174.59924852630596 257.0822687250529
diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py
@@ -174,17 +174,17 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds):
             p6.append(pldata[key].vectors()['vz'][0] * VCONV)
             p7.append(pldata[key].elements()['a'][0] * DCONV)
             p8.append(pldata[key].elements()['e'][0])
-            p9.append(pldata[key].elements()['incl'][0] * np.pi / 180.0)
-            p10.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0)
-            p11.append(pldata[key].elements()['w'][0] * np.pi / 180.0)
-            p12.append(pldata[key].elements()['M'][0] * np.pi / 180.0)
+            p9.append(pldata[key].elements()['incl'][0])
+            p10.append(pldata[key].elements()['Omega'][0])
+            p11.append(pldata[key].elements()['w'][0])
+            p12.append(pldata[key].elements()['M'][0])
         elif param['OUT_FORM'] == 'EL':
             p1.append(pldata[key].elements()['a'][0] * DCONV)
             p2.append(pldata[key].elements()['e'][0])
-            p3.append(pldata[key].elements()['incl'][0] * np.pi / 180.0)
-            p4.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0)
-            p5.append(pldata[key].elements()['w'][0] * np.pi / 180.0)
-            p6.append(pldata[key].elements()['M'][0] * np.pi / 180.0)
+            p3.append(pldata[key].elements()['incl'][0])
+            p4.append(pldata[key].elements()['Omega'][0])
+            p5.append(pldata[key].elements()['w'][0])
+            p6.append(pldata[key].elements()['M'][0])
             p7.append(pldata[key].vectors()['x'][0] * DCONV)
             p8.append(pldata[key].vectors()['y'][0] * DCONV)
             p9.append(pldata[key].vectors()['z'][0] * DCONV)

diff --git a/python/swiftest/swiftest/tool.py b/python/swiftest/swiftest/tool.py
@@ -8,10 +8,10 @@
 """
 
 def wrap_angle(angle):
-    while np.any(angle >= 2 * np.pi):
-        angle[angle >= 2 * np.pi] -= 2 * np.pi
+    while np.any(angle >= 360.0 ):
+        angle[angle >= 360.0] -= 360.0
     while np.any(angle < 0.0):
-        angle[angle < 0.0] += 2 * np.pi
+        angle[angle < 0.0] += 360.0
     return angle
 
 def follow_swift(ds, ifol=None, nskp=None):
@@ -50,15 +50,9 @@ def follow_swift(ds, ifol=None, nskp=None):
         intxt = input('Input the print frequency\n')
         nskp = int(intxt)
 
-    dr = 180.0 / np.pi
     fol['obar'] = fol['capom'] + fol['omega']
     fol['obar'] = fol['obar'].fillna(0)
     fol['obar'] = wrap_angle(fol['obar'])
-    fol['obar'] = fol['obar'] * dr
-    fol['inc'] = fol['inc'] * dr
-    fol['capom'] = fol['capom'] * dr
-    fol['omega'] = fol['omega'] * dr
-    fol['capm'] = fol['capm'] * dr
     fol['peri'] = fol['a'] * (1.0 - fol['e'])
     fol['apo']  = fol['a'] * (1.0 + fol['e'])
 

diff --git a/src/io/io.f90 b/src/io/io.f90
@@ -1101,6 +1101,13 @@ module function io_read_frame_body(self, iu, param) result(ierr)
 
             end do
          end select
+
+         if (param%in_form == EL) then
+            self%inc(1:n)   = self%inc(1:n) * DEG2RAD
+            self%capom(1:n) = self%capom(1:n) * DEG2RAD
+            self%omega(1:n) = self%omega(1:n) * DEG2RAD
+            self%capm(1:n)  = self%capm(1:n) * DEG2RAD
+         end if
       end associate
 
       ierr = 0
@@ -1205,11 +1212,6 @@ module function io_read_frame_system(self, iu, param) result(ierr)
          goto 667
       end if
 
-      if (param%in_form == EL) then
-         call self%pl%el2xv(self%cb)
-         call self%tp%el2xv(self%cb)
-      end if
-
       return
 
       667 continue
@@ -1482,10 +1484,10 @@ module subroutine io_write_frame_body(self, iu, param)
          case (EL) 
             write(iu, err = 667, iomsg = errmsg) self%a(1:n)
             write(iu, err = 667, iomsg = errmsg) self%e(1:n)
-            write(iu, err = 667, iomsg = errmsg) self%inc(1:n)
-            write(iu, err = 667, iomsg = errmsg) self%capom(1:n)
-            write(iu, err = 667, iomsg = errmsg) self%omega(1:n)
-            write(iu, err = 667, iomsg = errmsg) self%capm(1:n)
+            write(iu, err = 667, iomsg = errmsg) self%inc(1:n) * RAD2DEG
+            write(iu, err = 667, iomsg = errmsg) self%capom(1:n) * RAD2DEG
+            write(iu, err = 667, iomsg = errmsg) self%omega(1:n) * RAD2DEG
+            write(iu, err = 667, iomsg = errmsg) self%capm(1:n) * RAD2DEG
          case (XV)
             write(iu, err = 667, iomsg = errmsg) self%xh(1, 1:n)
             write(iu, err = 667, iomsg = errmsg) self%xh(2, 1:n)

diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90
@@ -22,7 +22,8 @@ module swiftest_globals
    real(DP), parameter :: PI3BY2 = 4.712388980384689857693965074919254326296_DP !! Definition of /(3 \pi / 2\)
    real(DP), parameter :: TWOPI  = 6.283185307179586476925286766559005768394_DP !! Definition of 2 \pi
    real(DP), parameter :: THIRD  = 0.333333333333333333333333333333333333333_DP !! Definition of 1 / 3
-   real(DP), parameter :: DEGRAD = 180.0_DP/PI !! Definition of conversion factor from degrees to radians
+   real(DP), parameter :: DEG2RAD = PI / 180.0_DP !! Definition of conversion factor from degrees to radians
+   real(DP), parameter :: RAD2DEG = 180.0_DP / PI !! Definition of conversion factor from degrees to radians
 
    integer(I4B), parameter :: LOWERCASE_BEGIN  = iachar('a') !! ASCII character set parameter for lower to upper conversion - start of lowercase
    integer(I4B), parameter :: LOWERCASE_END    = iachar('z') !! ASCII character set parameter for lower to upper conversion - end of lowercase

diff --git a/src/setup/setup.f90 b/src/setup/setup.f90
@@ -141,6 +141,10 @@ module subroutine setup_initialize_system(self, param)
       call self%set_msys()
       call self%pl%set_mu(self%cb) 
       call self%tp%set_mu(self%cb) 
+      if (param%in_form == EL) then
+         call self%pl%el2xv(self%cb)
+         call self%tp%el2xv(self%cb)
+      end if
       call self%pl%eucl_index()
       if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb)
       self%pl%lfirst = param%lfirstkick

diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90
@@ -417,22 +417,24 @@ module subroutine symba_util_rearray_pl(self, system, param)
       class(symba_parameters),   intent(in)    :: param  !! Current run configuration parameters
       ! Internals
       class(symba_pl), allocatable :: tmp !! The discarded body list.
-      integer(I4B) :: i, j, k
+      integer(I4B) :: i, j, k, npl
       logical, dimension(:), allocatable :: lmask
       class(symba_plplenc), allocatable :: plplenc_old
       logical :: lencounter
 
-      associate(pl => self, pl_adds => system%pl_adds)
+      associate(pl => self, pl_adds => system%pl_adds, nadd => system%pl_adds%nbody)
 
+         npl = pl%nbody
          ! Deallocate any temporary variables
          if (allocated(pl%xbeg)) deallocate(pl%xbeg)
          if (allocated(pl%xend)) deallocate(pl%xend)
 
          ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere
-         allocate(lmask, source=pl%ldiscard(:))
-         lmask(:) = lmask(:) .or. pl%status(:) == INACTIVE
+         allocate(lmask(npl))
+         lmask(1:npl) = pl%ldiscard(1:npl) .or. pl%status(1:npl) == INACTIVE
          allocate(tmp, mold=self)
          call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.)
+         npl = pl%nbody
          call tmp%setup(0,param)
          deallocate(tmp)
          deallocate(lmask)
@@ -442,39 +444,40 @@ module subroutine symba_util_rearray_pl(self, system, param)
          call plplenc_old%copy(system%plplenc_list)
 
          ! Add in any new bodies
-         if (pl_adds%nbody > 0) then
+         if (nadd > 0) then
             ! Append the adds to the main pl object
-            call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)])
+            call pl%append(pl_adds, lsource_mask=[(.true., i=1, nadd)])
+            npl = pl%nbody
 
-            allocate(lmask(pl%nbody)) 
-            lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE
+            allocate(lmask(npl)) 
+            lmask(1:npl) = pl%status(1:npl) == NEW_PARTICLE
 
-            call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask))
+            call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, npl)], lmask))
 
             deallocate(lmask)
          end if
 
          ! Reset all of the status flags for this body
-         where(pl%status(1:pl%nbody) /= INACTIVE) 
-            pl%status(1:pl%nbody) = ACTIVE
-            pl%ldiscard(1:pl%nbody) = .false.
-            pl%lcollision(1:pl%nbody) = .false.
-            pl%lmtiny(1:pl%nbody) = pl%Gmass(1:pl%nbody) > param%GMTINY
-            pl%lmask(1:pl%nbody) = .true.
+         where(pl%status(1:npl) /= INACTIVE) 
+            pl%status(1:npl) = ACTIVE
+            pl%ldiscard(1:npl) = .false.
+            pl%lcollision(1:npl) = .false.
+            pl%lmtiny(1:npl) = pl%Gmass(1:npl) > param%GMTINY
+            pl%lmask(1:npl) = .true.
          elsewhere
-            pl%ldiscard(1:pl%nbody) = .true.
-            pl%lmask(1:pl%nbody) = .false.
+            pl%ldiscard(1:npl) = .true.
+            pl%lmask(1:npl) = .false.
          end where
 
-         pl%nplm = count(pl%lmtiny(1:pl%nbody) .and. pl%lmask(1:pl%nbody))
+         pl%nplm = count(pl%lmtiny(1:npl) .and. pl%lmask(1:npl))
 
          ! Reindex the new list of bodies 
          call pl%sort("mass", ascending=.false.)
          call pl%eucl_index()
 
          ! Reset the kinship trackers
-         pl%kin(1:pl%nbody)%nchild = 0
-         pl%kin(1:pl%nbody)%parent = [(i, i=1, pl%nbody)]
+         pl%kin(1:npl)%nchild = 0
+         pl%kin(1:npl)%parent = [(i, i=1, npl)]
 
          ! Re-build the encounter list
          lencounter = pl%encounter_check(system, param%dt, 0)