Merge branch 'debug' into IO

MintonGroup · Sep 16, 2021 · 21b859a · 21b859a
2 parents aaa2f67 + 58583ab
commit 21b859a
Show file tree

Hide file tree

Showing 11 changed files with 109 additions and 58 deletions.
diff --git a/Makefile b/Makefile
@@ -59,25 +59,30 @@ SWIFTEST_MODULES =   swiftest_globals.f90 \
 
 include Makefile.Defines
 
+MKL_ROOT = /apps/spack/bell/apps/intel-parallel-studio/cluster.2019.5-intel-19.0.5-4brgqlf/mkl/lib
+IMKL = -I$(MKLROOT)/include
+LMKL = -L$(MKLROOT)/lib/intel64 -qopt-matmul 
+
 MODULES         = $(SWIFTEST_MODULES) $(USER_MODULES) 
 
-.PHONY : all mod lib libdir drivers bin clean force 
+.PHONY : all mod lib libdir fast drivers bin clean force 
 
 % : %.f90 force
-	$(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include  $< -o $@ \
-	  -L$(SWIFTEST_HOME)/lib -lswiftest -L$(NETCDF_FORTRAN_HOME)/lib -lnetcdf -lnetcdff
+	$(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) $< -o $@ \
+	  -L$(SWIFTEST_HOME)/lib -lswiftest -L$(NETCDF_FORTRAN_HOME)/lib -lnetcdf -lnetcdff $(LMKL)
 	$(INSTALL_PROGRAM) $@ $(SWIFTEST_HOME)/bin
 	rm -f $@
 
 all:
 	cd $(SWIFTEST_HOME); \
 	  make mod; \
 	  make lib; \
+	  make fast; \
 	  make drivers; \
 
 mod:
 	cd $(SWIFTEST_HOME)/src/modules/; \
-	  $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -c $(MODULES); \
+	  $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c $(MODULES); \
 	  $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o; \
 	  $(INSTALL_DATA) *.mod *.smod $(SWIFTEST_HOME)/include; \
 	  rm -f *.o *.mod  *.smod
@@ -93,11 +98,6 @@ lib:
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make libdir
-	cd $(SWIFTEST_HOME)/src/fraggle; \
-	  rm -f Makefile.Defines Makefile; \
-	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
-	  ln -s $(SWIFTEST_HOME)/Makefile .; \
-	  make libdir
 	cd $(SWIFTEST_HOME)/src/gr; \
 	  rm -f Makefile.Defines Makefile; \
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
@@ -143,11 +143,6 @@ lib:
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make libdir
-	cd $(SWIFTEST_HOME)/src/util; \
-	  rm -f Makefile.Defines Makefile; \
-	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
-	  ln -s $(SWIFTEST_HOME)/Makefile .; \
-	  make libdir
 	cd $(SWIFTEST_HOME)/src/whm; \
 	  rm -f Makefile.Defines Makefile; \
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
@@ -174,8 +169,39 @@ lib:
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make libdir
 
+fast:
+	cd $(SWIFTEST_HOME)/src/fraggle; \
+	  rm -f Makefile.Defines Makefile; \
+	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
+	  ln -s $(SWIFTEST_HOME)/Makefile .; \
+	  make fastdir
+
+	cd $(SWIFTEST_HOME)/src/util; \
+	  rm -f Makefile.Defines Makefile; \
+	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
+	  ln -s $(SWIFTEST_HOME)/Makefile .; \
+	  make fastdir
+
+	cd $(SWIFTEST_HOME)/src/rmvs; \
+		$(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c rmvs_encounter_check.f90; \
+		$(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \
+		$(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \
+		rm -f *.o *.smod	
+
+	cd $(SWIFTEST_HOME)/src/symba; \
+		$(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c symba_encounter_check.f90; \
+		$(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \
+		$(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \
+		rm -f *.o *.smod
+
 libdir:
-	$(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -c *.f90; \
+	$(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c *.f90; \
+	$(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \
+	$(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \
+	rm -f *.o *.smod
+
+fastdir:
+	$(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c *.f90; \
 	$(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \
 	$(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \
 	rm -f *.o *.smod

diff --git a/Makefile.Defines b/Makefile.Defines
@@ -54,10 +54,10 @@ VTUNE_FLAGS = -g -O2 -qopt-report=5 -simd -shared-intel -qopenmp -debug inline-d
 IDEBUG = -O0 -init=snan,arrays  -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
 STRICTREAL = -fp-model strict -prec-div -prec-sqrt -assume protect-parens 
 SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except 
-PAR = -qopenmp #-parallel  #Something goes wrong in SyMBA at the moment with auto-paralellization enabled
+PAR = -qopenmp -parallel  
 HEAPARR = -heap-arrays 4194304 
 OPTREPORT = -qopt-report=5 
-IPRODUCTION = -init=snan,arrays  -no-wrap-margin -O3 $(PAR) $(SIMDVEC) $(STRICTREAL) #$(HEAPARR)
+IPRODUCTION = -no-wrap-margin -O3 $(PAR) $(SIMDVEC) #$(HEAPARR)
 
 #gfortran flags
 GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none
@@ -66,8 +66,10 @@ GMEM =  -fsanitize-address-use-after-scope -fstack-check  -fsanitize=bounds-stri
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
 GPRODUCTION = -O2 -ffree-line-length-none $(GPAR)
 
+
 #FFLAGS  = $(IDEBUG) $(SIMDVEC) $(PAR)
-FFLAGS = $(IPRODUCTION) $(OPTREPORT) #$(ADVIXE_FLAGS)
+FFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT) 
+FFASTFLAGS = $(IPRODUCTION) -fp-model fast $(OPTREPORT)
 FORTRAN = ifort
 #AR     = xiar
 

diff --git a/examples/symba_mars_disk/aescattermovie.py b/examples/symba_mars_disk/aescattermovie.py
@@ -79,8 +79,7 @@ def setup_plot(self):
     def data_stream(self, frame=0):
         while True:
             d = self.ds.isel(time=frame)
-
-            d = d.where(d['radius'] < self.Rcb, drop=True)
+            d = d.where(np.invert(np.isnan(d['a'])), drop=True)
             d['radmarker'] = (d['radius'] / self.Rcb) * radscale
             radius = d['radmarker'].values
 

diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in
@@ -1,17 +1,17 @@
 !Parameter file for the SyMBA-RINGMOONS test
 T0             0.0 
-TSTOP          1.0e12 
+TSTOP          60000.0
 DT             600.0
 CB_IN          cb.in
 PL_IN          mars.in
 TP_IN          tp.in
 IN_TYPE        ASCII
-ISTEP_OUT      1
-ISTEP_DUMP     1
+ISTEP_OUT      100
+ISTEP_DUMP     100
 BIN_OUT        bin.nc 
 PARTICLE_OUT   particle.dat
-OUT_TYPE       NETCDF_DOUBLE
-OUT_FORM       XVEL
+OUT_TYPE       REAL8
+OUT_FORM       XV
 OUT_STAT       REPLACE
 CHK_CLOSE      yes
 CHK_RMIN       3389500.0

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/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py
@@ -694,8 +694,8 @@ def swiftest2xr(param):
         ds = xr.combine_by_coords([cbds, plds, tpds])
 
     elif ((param['OUT_TYPE'] == 'NETCDF_DOUBLE') or (param['OUT_TYPE'] == 'NETCDF_FLOAT')):
-        print('\nCreating Dataset')
-        ds = xr.open_dataset(param['BIN_OUT'])
+        print('\nCreating Dataset from NetCDF file')
+        ds = xr.open_dataset(param['BIN_OUT'], mask_and_scale=False)
         ds = clean_string_values(param, ds)
     else:
         print(f"Error encountered. OUT_TYPE {param['OUT_TYPE']} not recognized.")

diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90
@@ -545,14 +545,18 @@ function radial_objective_function(v_r_mag_input) result(fval)
             integer(I4B)                          :: i
             real(DP), dimension(:,:), allocatable :: v_shift
             real(DP), dimension(frag%nbody)       :: kearr
-            real(DP)                              :: keo, ke_radial
+            real(DP)                              :: keo, ke_radial, rotmag2, vmag2
 
             associate(nfrag => frag%nbody)
                allocate(v_shift, mold=frag%vb)
                v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, frag%v_r_unit, frag%v_t_mag, frag%v_t_unit, frag%mass, frag%vbcom) 
-               do concurrent(i = 1:nfrag)
-                  kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * dot_product(frag%rot(:, i), frag%rot(:, i)) + dot_product(v_shift(:, i), v_shift(:, i)))
+               !$omp do simd 
+               do i = 1,nfrag
+                  rotmag2 = frag%rot(1,i)**2 + frag%rot(2,i)**2 + frag%rot(3,i)**2
+                  vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2
+                  kearr(i) = frag%mass(i) * (frag%Ip(3, i) * frag%radius(i)**2 * rotmag2 + vmag2) 
                end do
+               !$omp end do simd
                keo = 2 * frag%ke_budget - sum(kearr(:))
                ke_radial = frag%ke_budget - frag%ke_orbit - frag%ke_spin
                ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for

diff --git a/src/kick/kick.f90 b/src/kick/kick.f90
@@ -66,7 +66,7 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius,
       ahi(:,:) = 0.0_DP
       ahj(:,:) = 0.0_DP
 
-      !$omp parallel do simd default(private) schedule(static)&
+      !$omp parallel do default(private) schedule(static)&
       !$omp shared(nplpl, k_plpl, x, Gmass, radius) &
       !$omp lastprivate(rji2, rlim2, xr, yr, zr) &
       !$omp reduction(+:ahi) &
@@ -81,8 +81,7 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius,
          rlim2 = (radius(i) + radius(j))**2
          if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
       end do
-      !$omp end parallel do simd
-
+      !$omp end parallel do 
 
       do concurrent(i = 1:npl)
          acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i)
@@ -113,7 +112,8 @@ module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
       integer(I4B) :: i, j
 
       !$omp parallel do default(private) schedule(static)&
-      !$omp shared(npl, ntp, lmask, xtp, xpl, acc)
+      !$omp shared(npl, ntp, lmask, xtp, xpl) &
+      !$omp reduction(-:acc)
       do i = 1, ntp
          if (lmask(i)) then
             do j = 1, npl

diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90
@@ -38,17 +38,22 @@ module subroutine netcdf_initialize_output(self, param)
       !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
       !!
       !! Initialize a NetCDF file system and defines all variables.
+      use, intrinsic :: ieee_arithmetic
       implicit none
       ! Arguments
       class(netcdf_parameters),   intent(inout) :: self    !! Parameters used to identify a particular NetCDF dataset
       class(swiftest_parameters), intent(in)    :: param           !! Current run configuration parameters 
       ! Internals
       logical :: fileExists
-      integer(I4B) :: old_mode
+      integer(I4B) :: old_mode, nvar, varid, vartype
+      real(DP) :: dfill
+      real(SP) :: sfill
+
+      dfill = ieee_value(dfill, IEEE_QUIET_NAN)
+      sfill = ieee_value(sfill, IEEE_QUIET_NAN)
 
       !! Create the new output file, deleting any previously existing output file of the same name
       call check( nf90_create(param%outfile, NF90_NETCDF4, self%ncid) )
-      call check( nf90_set_fill(self%ncid, nf90_nofill, old_mode) )
 
       ! Define the NetCDF dimensions with particle name as the record dimension
       call check( nf90_def_dim(self%ncid, ID_DIMNAME, NF90_UNLIMITED, self%id_dimid) )     ! 'x' dimension
@@ -138,6 +143,22 @@ module subroutine netcdf_initialize_output(self, param)
       call check( nf90_def_var(self%ncid, DISCARD_VHZ_VARNAME, self%out_type, self%id_dimid,  self%discard_vhz_varid) )
       call check( nf90_def_var(self%ncid, DISCARD_BODY_ID_VARNAME, NF90_INT, self%id_dimid, self%discard_body_id_varid) )
 
+      ! Set fill mode to NaN for all variables
+      call check( nf90_inquire(self%ncid, nVariables=nvar) )
+      do varid = 1, nvar
+         call check( nf90_inquire_variable(self%ncid, varid, xtype=vartype) )
+         select case(vartype)
+         case(NF90_INT)
+            call check( nf90_def_var_fill(self%ncid, varid, 0, NF90_FILL_INT) )
+         case(NF90_FLOAT)
+            call check( nf90_def_var_fill(self%ncid, varid, 0, sfill) )
+         case(NF90_DOUBLE)
+            call check( nf90_def_var_fill(self%ncid, varid, 0, dfill) )
+         case(NF90_CHAR)
+            call check( nf90_def_var_fill(self%ncid, varid, 0, 0) )
+         end select
+      end do
+
       return
    end subroutine netcdf_initialize_output
 
@@ -150,11 +171,8 @@ module subroutine netcdf_open(self, param)
       ! Arguments
       class(netcdf_parameters),   intent(inout) :: self   !! Parameters used to identify a particular NetCDF dataset
       class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters
-      ! Internals
-      integer(I4B) :: old_mode
 
       call check( nf90_open(param%outfile, nf90_write, self%ncid) )
-      call check( nf90_set_fill(self%ncid, nf90_nofill, old_mode) )
 
       call check( nf90_inq_varid(self%ncid, TIME_DIMNAME, self%time_varid))
       call check( nf90_inq_varid(self%ncid, ID_DIMNAME, self%id_varid))
@@ -249,14 +267,15 @@ module subroutine netcdf_write_frame_base(self, iu, param)
       class(netcdf_parameters),   intent(inout) :: iu     !! Parameters used to identify a particular NetCDF dataset
       class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters
       ! Internals
-      integer(I4B)                              :: i, j, tslot, strlen, idslot
+      integer(I4B)                              :: i, j, tslot, strlen, idslot, old_mode
       integer(I4B), dimension(:), allocatable   :: ind
       character(len=:), allocatable             :: charstring
 
       call self%write_particle_info(iu)
 
       tslot = int(param%ioutput, kind=I4B) + 1
 
+      call check( nf90_set_fill(iu%ncid, nf90_nofill, old_mode) )
       select type(self)
          class is (swiftest_body)
          associate(n => self%nbody)
@@ -332,6 +351,7 @@ module subroutine netcdf_write_frame_base(self, iu, param)
          end if
 
       end select
+      call check( nf90_set_fill(iu%ncid, old_mode, old_mode) )
 
       return
    end subroutine netcdf_write_frame_base
@@ -346,13 +366,14 @@ module subroutine netcdf_write_particle_info_base(self, iu)
       class(swiftest_base),       intent(in)    :: self   !! Swiftest particle object
       class(netcdf_parameters),   intent(inout) :: iu     !! Parameters used to identify a particular NetCDF dataset
       ! Internals
-      integer(I4B)                              :: i, j, tslot, strlen, idslot
+      integer(I4B)                              :: i, j, tslot, strlen, idslot, old_mode
       integer(I4B), dimension(:), allocatable   :: ind
       character(len=:), allocatable             :: charstring
       character(len=NAMELEN)                    :: emptystr, lenstr
       character(len=:), allocatable :: fmtlabel
 
       ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables
+      call check( nf90_set_fill(iu%ncid, nf90_nofill, old_mode) )
       write(lenstr, *) NAMELEN
       fmtlabel = "(A" // trim(adjustl(lenstr)) // ")"
       write(emptystr, fmtlabel) " "
@@ -448,6 +469,7 @@ module subroutine netcdf_write_particle_info_base(self, iu)
          call check( nf90_put_var(iu%ncid, iu%discard_vhz_varid, self%info%discard_vh(3), start=[idslot]) )
       end select
 
+      call check( nf90_set_fill(iu%ncid, old_mode, old_mode) )
       return
    end subroutine netcdf_write_particle_info_base
 
@@ -464,12 +486,11 @@ module subroutine netcdf_write_hdr_system(self, iu, param)
       class(netcdf_parameters),     intent(inout) :: iu    !! Parameters used to for writing a NetCDF dataset to file
       class(swiftest_parameters),   intent(in)    :: param !! Current run configuration parameters
       ! Internals
-      integer(I4B) :: tslot, old_mode
+      integer(I4B) :: tslot
 
       tslot = int(param%ioutput, kind=I4B) + 1
 
       call check( nf90_open(param%outfile, nf90_write, iu%ncid) )
-      call check( nf90_set_fill(iu%ncid, nf90_nofill, old_mode) )
 
       call check( nf90_put_var(iu%ncid, iu%time_varid, param%t, start=[tslot]) )
       call check( nf90_put_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]) )

diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90
@@ -284,9 +284,14 @@ module pure subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhil
       integer(I4B), intent(in)  :: irec
       logical,      intent(out) :: lencounter, lvdotr
       ! Internals
-      real(DP)     :: r2crit
+      real(DP)     :: r2crit, rshell_irec
+      integer(I4B) :: i
 
-      r2crit = (rhill1 + rhill2)*RHSCALE*(RSHELL**(irec))
+      rshell_irec = 1._DP
+      do i = 1, irec
+         rshell_irec = rshell_irec * RSHELL
+      end do
+      r2crit = (rhill1 + rhill2) * RHSCALE * rshell_irec
       r2crit = r2crit**2
       call rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit, lencounter, lvdotr)
 

diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90
@@ -139,19 +139,13 @@ subroutine util_get_energy_potential(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mas
          pecb(i) = -GMcb * mass(i) / norm2(xb(:,i)) 
       end do
 
-      ! Do the potential energy between pairs of massive bodies
-      do concurrent (k = 1:nplpl)
+      !$omp parallel do default(private) schedule(static)&
+      !$omp shared(nplpl, k_plpl, xb, mass, Gmass, pepl, lstatpl, lmask)
+      do k = 1, nplpl
          i = k_plpl(1,k)
          j = k_plpl(2,k)
          lstatpl(k) = (lmask(i) .and. lmask(j))
-      end do
-
-      !$omp parallel do default(private) schedule(static)&
-      !$omp shared(nplpl, k_plpl, xb, mass, Gmass, pepl, lstatpl)
-      do k = 1, nplpl
          if (lstatpl(k)) then
-            i = k_plpl(1,k)
-            j = k_plpl(2,k)
             pepl(k) = -(Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j))
          else
             pepl(k) = 0.0_DP