From 17438ebd38b28a1518401c9d71bd2b277de5a76a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 Nov 2021 12:18:19 -0500 Subject: [PATCH 1/6] Fixed typo in compiler flags --- Makefile.Defines | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index 849f205b1..d7244a886 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -51,7 +51,7 @@ COLLRESOLVE_HOME = $(ROOT_DIR)/collresolve/ ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vecabi=cmdtarget -simd -shared-intel -debug inline-debug-info -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -xhost -traceback -parallel-source-info=2 #Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune -IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margi-parallel-source-info=2 +IDEBUG = -O2 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -parallel-source-info=2 STRICTREAL = -fp-model=precise -prec-div -prec-sqrt -assume protect-parens SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma PAR = -qopenmp -parallel @@ -69,8 +69,8 @@ GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul $(PAR) -FSTRICTFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) -FFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) +FSTRICTFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) $(HEAPARR) +FFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) $(HEAPARR) # FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) -g -traceback # FFLAGS = $(IPRODUCTION) -fp-model=fast -g -traceback From e44c0d9236a5096f9805869025c9855a7c8cfe86 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 Nov 2021 15:23:19 -0500 Subject: [PATCH 2/6] Correctly read in initial values of central body mass, radius, and angular momentum when restarting from NetCDF. --- src/netcdf/netcdf.f90 | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index b99188855..fd2406720 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -66,6 +66,7 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) integer(I4B) :: itmax, idmax real(DP), dimension(:), allocatable :: vals real(DP), dimension(1) :: val + real(DP), dimension(NDIM) :: rot0 real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, Ltmp call param%nciu%open(param) @@ -74,7 +75,8 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) allocate(vals(idmax)) call check( nf90_get_var(param%nciu%ncid, param%nciu%time_varid, val, start=[1], count=[1]) ) - old_t_final = val(1) + !old_t_final = val(1) + old_t_final = param%t0 ! For NetCDF it is safe to overwrite the final t value on a restart if (param%lenergy) then call check( nf90_get_var(param%nciu%ncid, param%nciu%KE_orb_varid, val, start=[1], count=[1]) ) @@ -107,6 +109,24 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) call check( nf90_get_var(param%nciu%ncid, param%nciu%Gmass_varid, vals, start=[1,1], count=[idmax,1]) ) self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + select type(cb => self%cb) + class is (symba_cb) + cb%GM0 = vals(1) + cb%dGM = cb%Gmass - cb%GM0 + + call check( nf90_get_var(param%nciu%ncid, param%nciu%radius_varid, val, start=[1,1], count=[1,1]) ) + cb%R0 = val(1) + + call check( nf90_get_var(param%nciu%ncid, param%nciu%rotx_varid, val, start=[1,1], count=[1,1]) ) + rot0(1) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%roty_varid, val, start=[1,1], count=[1,1]) ) + rot0(2) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%rotz_varid, val, start=[1,1], count=[1,1]) ) + rot0(3) = val(1) + + cb%L0(:) = cb%Ip(3) * cb%GM0 * cb%R0**2 * rot0(:) + end select + end if deallocate(vals) From 001247ba2f8d25a72a25404fc1eacb8abc5d76d2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 Nov 2021 15:23:47 -0500 Subject: [PATCH 3/6] Misc. formatting fixes --- src/io/io.f90 | 1 - src/kick/kick.f90 | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index ebdd68531..0266af0bc 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1351,7 +1351,6 @@ subroutine io_read_in_cb(self, param) close(iu, err = 667, iomsg = errmsg) if (ierr == 0) then - if (self%j2rp2 /= 0.0_DP) param%loblatecb = .true. if (param%rmin < 0.0) param%rmin = self%radius diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 891936f20..c3e37d927 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -189,7 +189,7 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius end do !$omp end parallel do else - !$omp parallel do default(private) schedule(static)& + !$omp parallel do default(private) schedule(static)& !$omp shared(npl, nplm, x, Gmass, radius) & !$omp lastprivate(rji2, xr, yr, zr) & !$omp reduction(+:ahi) & From 4260a05c99ea75c490be16b019f4c2d34da3c217 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 Nov 2021 15:24:53 -0500 Subject: [PATCH 4/6] Fixed debug flags for non-parallel debugging --- Makefile.Defines | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index d7244a886..90a2ceb58 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -51,10 +51,10 @@ COLLRESOLVE_HOME = $(ROOT_DIR)/collresolve/ ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vecabi=cmdtarget -simd -shared-intel -debug inline-debug-info -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -xhost -traceback -parallel-source-info=2 #Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune -IDEBUG = -O2 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -parallel-source-info=2 +IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin STRICTREAL = -fp-model=precise -prec-div -prec-sqrt -assume protect-parens SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma -PAR = -qopenmp -parallel +PAR = -qopenmp -parallel -parallel-source-info=2 HEAPARR = -heap-arrays 4194304 OPTREPORT = -qopt-report=5 IPRODUCTION = -no-wrap-margin -O3 -qopt-prefetch=0 -qopt-matmul -sox $(PAR) $(SIMDVEC) #$(HEAPARR) @@ -67,12 +67,12 @@ GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporari GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include -LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul $(PAR) +LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul #$(PAR) -FSTRICTFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) $(HEAPARR) -FFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) $(HEAPARR) -# FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) -g -traceback -# FFLAGS = $(IPRODUCTION) -fp-model=fast -g -traceback +FSTRICTFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) $(HEAPARR) +FFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) $(HEAPARR) +#FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) -g -traceback +#FFLAGS = $(IPRODUCTION) -fp-model=fast -g -traceback FORTRAN = ifort AR = xiar From 1a8d18f5121fdcbe2b62741e9ac76ca6d6080764 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 Nov 2021 15:46:56 -0500 Subject: [PATCH 5/6] Fixedd calculation of initial angular momentum and change of angular momentum of central body for a restarted run --- src/io/io.f90 | 5 ++++- src/netcdf/netcdf.f90 | 29 +++++++++++++++++++++-------- 2 files changed, 25 insertions(+), 9 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 0266af0bc..b84569cca 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1360,7 +1360,10 @@ subroutine io_read_in_cb(self, param) cb%GM0 = cb%Gmass cb%dGM = 0.0_DP cb%R0 = cb%radius - cb%L0(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) + if (param%lrotation) then + cb%L0(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) + cb%dL(:) = 0.0_DP + end if end select end if return diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index fd2406720..d28fcc788 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -66,7 +66,7 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) integer(I4B) :: itmax, idmax real(DP), dimension(:), allocatable :: vals real(DP), dimension(1) :: val - real(DP), dimension(NDIM) :: rot0 + real(DP), dimension(NDIM) :: rot0, Ip0, Lnow real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, Ltmp call param%nciu%open(param) @@ -117,14 +117,27 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) call check( nf90_get_var(param%nciu%ncid, param%nciu%radius_varid, val, start=[1,1], count=[1,1]) ) cb%R0 = val(1) - call check( nf90_get_var(param%nciu%ncid, param%nciu%rotx_varid, val, start=[1,1], count=[1,1]) ) - rot0(1) = val(1) - call check( nf90_get_var(param%nciu%ncid, param%nciu%roty_varid, val, start=[1,1], count=[1,1]) ) - rot0(2) = val(1) - call check( nf90_get_var(param%nciu%ncid, param%nciu%rotz_varid, val, start=[1,1], count=[1,1]) ) - rot0(3) = val(1) + if (param%lrotation) then - cb%L0(:) = cb%Ip(3) * cb%GM0 * cb%R0**2 * rot0(:) + call check( nf90_get_var(param%nciu%ncid, param%nciu%rotx_varid, val, start=[1,1], count=[1,1]) ) + rot0(1) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%roty_varid, val, start=[1,1], count=[1,1]) ) + rot0(2) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%rotz_varid, val, start=[1,1], count=[1,1]) ) + rot0(3) = val(1) + + call check( nf90_get_var(param%nciu%ncid, param%nciu%Ip1_varid, val, start=[1,1], count=[1,1]) ) + Ip0(1) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%Ip2_varid, val, start=[1,1], count=[1,1]) ) + Ip0(2) = val(1) + call check( nf90_get_var(param%nciu%ncid, param%nciu%Ip3_varid, val, start=[1,1], count=[1,1]) ) + Ip0(3) = val(1) + + cb%L0(:) = Ip0(3) * cb%GM0 * cb%R0**2 * rot0(:) + + Lnow(:) = cb%Ip(3) * cb%Gmass * cb%radius**2 * cb%rot(:) + cb%dL(:) = Lnow(:) - cb%L0(:) + end if end select end if From 966dc068ef214a42bf0d2279c02035105358ee6e Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 Nov 2021 16:33:09 -0500 Subject: [PATCH 6/6] Improved handling of mass conservation error halting. Failed output info is now saved in a time slot after the present one in order to prevent overwriting of the bin.nc with bad data --- src/io/io.f90 | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index b84569cca..7ec71d9ae 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -72,11 +72,20 @@ module subroutine io_conservation_report(self, param, lterminal) if (lterminal) write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror if (abs(Merror) > 100 * epsilon(Merror)) then write(*,*) "Severe error! Mass not conserved! Halting!" - call pl%xv2el(cb) - call self%write_hdr(param%nciu, param) - call cb%write_frame(param%nciu, param) - call pl%write_frame(param%nciu, param) - call param%nciu%close() + if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then + write(*,*) "Merror = ", Merror + write(*,*) "GMtot_now : ",GMtot_now + write(*,*) "GMtot_orig: ",system%GMtot_orig + write(*,*) "Difference: ",GMtot_now - system%GMtot_orig + else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + ! Save the frame of data to the bin file in the slot just after the present one for diagnostics + param%ioutput = param%ioutput + 1_I8B + call pl%xv2el(cb) + call self%write_hdr(param%nciu, param) + call cb%write_frame(param%nciu, param) + call pl%write_frame(param%nciu, param) + call param%nciu%close() + end if call util_exit(FAILURE) end if end if