Merge branch 'debug' into IntelAdvisor

MintonGroup · Nov 15, 2021 · 677060f · 677060f
2 parents e4ddc20 + 966dc06
commit 677060f
Show file tree

Hide file tree

Showing 3 changed files with 53 additions and 9 deletions.
diff --git 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
@@ -1351,7 +1360,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
@@ -1361,7 +1369,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/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) &

diff --git 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, Ip0, Lnow
       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,37 @@ 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) 
+
+            if (param%lrotation) then
+
+               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
 
       deallocate(vals)