Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Restarts from NetCDF files fixed
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Oct 20, 2021
1 parent f5807f2 commit a90c125
Show file tree
Hide file tree
Showing 6 changed files with 133 additions and 68 deletions.
68 changes: 45 additions & 23 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,29 +46,29 @@ module subroutine io_conservation_report(self, param, lterminal)
Lorbit_now(:) = system%Lorbit(:)
Lspin_now(:) = system%Lspin(:)
Eorbit_now = ke_orbit_now + ke_spin_now + pe_now
Ltot_now(:) = system%Ltot(:) + param%Lescape(:)
GMtot_now = system%GMtot + param%GMescape
Ltot_now(:) = system%Ltot(:) + system%Lescape(:)
GMtot_now = system%GMtot + system%GMescape

if (param%lfirstenergy) then
param%Eorbit_orig = Eorbit_now
param%GMtot_orig = GMtot_now
param%Lorbit_orig(:) = Lorbit_now(:)
param%Lspin_orig(:) = Lspin_now(:)
param%Ltot_orig(:) = Ltot_now(:)
system%Eorbit_orig = Eorbit_now
system%GMtot_orig = GMtot_now
system%Lorbit_orig(:) = Lorbit_now(:)
system%Lspin_orig(:) = Lspin_now(:)
system%Ltot_orig(:) = Ltot_now(:)
param%lfirstenergy = .false.
end if

if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then
write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, param%Ecollisions, Ltot_now, GMtot_now
write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, system%Ecollisions, Ltot_now, GMtot_now
close(EGYIU, err = 667, iomsg = errmsg)
end if

if (.not.param%lfirstenergy) then
Lerror = norm2(Ltot_now(:) - param%Ltot_orig(:)) / norm2(param%Ltot_orig(:))
Eorbit_error = (Eorbit_now - param%Eorbit_orig) / abs(param%Eorbit_orig)
Ecoll_error = param%Ecollisions / abs(param%Eorbit_orig)
Etotal_error = (Eorbit_now - param%Ecollisions - param%Eorbit_orig - param%Euntracked) / abs(param%Eorbit_orig)
Merror = (GMtot_now - param%GMtot_orig) / param%GMtot_orig
Lerror = norm2(Ltot_now(:) - system%Ltot_orig(:)) / norm2(system%Ltot_orig(:))
Eorbit_error = (Eorbit_now - system%Eorbit_orig) / abs(system%Eorbit_orig)
Ecoll_error = system%Ecollisions / abs(system%Eorbit_orig)
Etotal_error = (Eorbit_now - system%Ecollisions - system%Eorbit_orig - system%Euntracked) / abs(system%Eorbit_orig)
Merror = (GMtot_now - system%GMtot_orig) / system%GMtot_orig
if (lterminal) write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror
if (abs(Merror) > 100 * epsilon(Merror)) then
write(*,*) "Severe error! Mass not conserved! Halting!"
Expand Down Expand Up @@ -281,6 +281,17 @@ module subroutine io_dump_system(self, param)
dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx)))
dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx)))
dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx)))

dump_param%Eorbit_orig = self%Eorbit_orig
dump_param%GMtot_orig = self%GMtot_orig
dump_param%Ltot_orig(:) = self%Ltot_orig(:)
dump_param%Lorbit_orig(:) = self%Lorbit_orig(:)
dump_param%Lspin_orig(:) = self%Lspin_orig(:)
dump_param%GMescape = self%GMescape
dump_param%Ecollisions = self%Ecollisions
dump_param%Euntracked = self%Euntracked
dump_param%Lescape(:) = self%Lescape

else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
dump_param%in_type = NETCDF_DOUBLE_TYPE
dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx)))
Expand Down Expand Up @@ -995,15 +1006,17 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)

if (param%lenergy) then
call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit)
call io_param_writer_one("EORBIT_ORIG", param%Eorbit_orig, unit)
call io_param_writer_one("GMTOT_ORIG", param%GMtot_orig, unit)
call io_param_writer_one("LTOT_ORIG", param%Ltot_orig(:), unit)
call io_param_writer_one("LORBIT_ORIG", param%Lorbit_orig(:), unit)
call io_param_writer_one("LSPIN_ORIG", param%Lspin_orig(:), unit)
call io_param_writer_one("LESCAPE", param%Lescape(:), unit)
call io_param_writer_one("GMESCAPE",param%GMescape, unit)
call io_param_writer_one("ECOLLISIONS",param%Ecollisions, unit)
call io_param_writer_one("EUNTRACKED",param%Euntracked, unit)
if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then
call io_param_writer_one("EORBIT_ORIG", param%Eorbit_orig, unit)
call io_param_writer_one("GMTOT_ORIG", param%GMtot_orig, unit)
call io_param_writer_one("LTOT_ORIG", param%Ltot_orig(:), unit)
call io_param_writer_one("LORBIT_ORIG", param%Lorbit_orig(:), unit)
call io_param_writer_one("LSPIN_ORIG", param%Lspin_orig(:), unit)
call io_param_writer_one("LESCAPE", param%Lescape(:), unit)
call io_param_writer_one("GMESCAPE",param%GMescape, unit)
call io_param_writer_one("ECOLLISIONS",param%Ecollisions, unit)
call io_param_writer_one("EUNTRACKED",param%Euntracked, unit)
end if
end if
call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit)
call io_param_writer_one("MAXID",param%maxid, unit)
Expand Down Expand Up @@ -1372,13 +1385,22 @@ module subroutine io_read_in_system(self, param)
allocate(tmp_param, source=param)
tmp_param%outfile = param%in_netcdf
tmp_param%out_form = param%in_form
ierr = self%read_frame(tmp_param%nciu, tmp_param)
ierr = self%read_frame(tmp_param%nciu, tmp_param)
deallocate(tmp_param)
if (ierr /=0) call util_exit(FAILURE)
else
call self%cb%read_in(param)
call self%pl%read_in(param)
call self%tp%read_in(param)
! Copy over param file variable inputs
self%Eorbit_orig = param%Eorbit_orig
self%GMtot_orig = param%GMtot_orig
self%Ltot_orig(:) = param%Ltot_orig(:)
self%Lorbit_orig(:) = param%Lorbit_orig(:)
self%Lspin_orig(:) = param%Lspin_orig(:)
self%Lescape(:) = param%Lescape(:)
self%Ecollisions = param%Ecollisions
self%Euntracked = param%Euntracked
end if

return
Expand Down
2 changes: 1 addition & 1 deletion src/main/swiftest_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,8 @@ program swiftest_driver
old_t_final = nbody_system%get_old_t_final_netcdf(param)
else
old_t_final = t0
if (istep_out > 0) call nbody_system%write_frame(param)
if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum
if (istep_out > 0) call nbody_system%write_frame(param)
end if

!> Define the maximum number of threads
Expand Down
16 changes: 12 additions & 4 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,12 +152,12 @@ module swiftest_classes

! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values)
real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy
real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass
real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass
real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector
real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum
real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector
real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping)
real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions
real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies
logical :: lfirstenergy = .true. !! This is the first time computing energe
Expand Down Expand Up @@ -424,6 +424,15 @@ module swiftest_classes
real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! System orbital angular momentum vector
real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! System spin angular momentum vector
real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector
real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy
real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass
real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector
real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum
real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector
real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping)
real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions
real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies
logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated
!! separately from massive bodies. Massive body variables are saved at half steps, and passed to
!! the test particles
Expand All @@ -442,7 +451,6 @@ module swiftest_classes
procedure :: read_frame_netcdf => netcdf_read_frame_system !! Read in a frame of input data from file
procedure :: write_frame_netcdf => netcdf_write_frame_system !! Write a frame of input data from file
procedure :: write_hdr_bin => io_write_hdr_system !! Write a header for an output frame in Fortran binary format
!procedure :: read_hdr_bin => io_read_hdr !! Read a header for an output frame in Fortran binary format
procedure :: read_hdr_netcdf => netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format
procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format
procedure :: read_in => io_read_in_system !! Reads the initial conditions for an nbody system
Expand Down Expand Up @@ -961,7 +969,7 @@ end subroutine netcdf_flush

module function netcdf_get_old_t_final_system(self, param) result(old_t_final)
implicit none
class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP) :: old_t_final !! Final time from last run
end function netcdf_get_old_t_final_system
Expand Down
93 changes: 64 additions & 29 deletions src/netcdf/netcdf.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,23 +58,58 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final)
!!
implicit none
! Arguments
class(swiftest_nbody_system), intent(in) :: self
class(swiftest_nbody_system), intent(inout) :: self
class(swiftest_parameters), intent(inout) :: param
! Result
real(DP) :: old_t_final
! Internals
integer(I4B) :: itmax
real(DP), dimension(:), allocatable :: tvals

integer(I4B) :: itmax, idmax
real(DP), dimension(:), allocatable :: vals
real(DP), dimension(1) :: val
real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, Ltmp

call param%nciu%open(param)
call check( nf90_inquire_dimension(param%nciu%ncid, param%nciu%time_dimid, len=itmax) )
allocate(tvals(itmax))
call check( nf90_get_var(param%nciu%ncid, param%nciu%time_varid, tvals(:)) )
call check( nf90_inquire_dimension(param%nciu%ncid, param%nciu%id_dimid, len=idmax) )
allocate(vals(idmax))
call check( nf90_get_var(param%nciu%ncid, param%nciu%time_varid, val, start=[itmax], count=[1]) )

old_t_final = val(1)

if (param%lenergy) then
call check( nf90_get_var(param%nciu%ncid, param%nciu%KE_orb_varid, val, start=[1], count=[1]) )
KE_orb_orig = val(1)

call check( nf90_get_var(param%nciu%ncid, param%nciu%KE_spin_varid, val, start=[1], count=[1]) )
KE_spin_orig = val(1)

call check( nf90_get_var(param%nciu%ncid, param%nciu%PE_varid, val, start=[1], count=[1]) )
PE_orig = val(1)

old_t_final = tvals(itmax)
self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig

call check( nf90_get_var(param%nciu%ncid, param%nciu%L_orbx_varid, val, start=[1], count=[1]) )
self%Lorbit_orig(1) = val(1)
call check( nf90_get_var(param%nciu%ncid, param%nciu%L_orby_varid, val, start=[1], count=[1]) )
self%Lorbit_orig(2) = val(1)
call check( nf90_get_var(param%nciu%ncid, param%nciu%L_orbz_varid, val, start=[1], count=[1]) )
self%Lorbit_orig(3) = val(1)

call check( nf90_get_var(param%nciu%ncid, param%nciu%L_spinx_varid, val, start=[1], count=[1]) )
self%Lspin_orig(1) = val(1)
call check( nf90_get_var(param%nciu%ncid, param%nciu%L_spiny_varid, val, start=[1], count=[1]) )
self%Lspin_orig(2) = val(1)
call check( nf90_get_var(param%nciu%ncid, param%nciu%L_spinz_varid, val, start=[1], count=[1]) )
self%Lspin_orig(3) = val(1)

self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:)

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))

end if

deallocate(tvals)
deallocate(vals)

return
end function netcdf_get_old_t_final_system
Expand Down Expand Up @@ -607,21 +642,21 @@ module subroutine netcdf_read_hdr_system(self, iu, param)
call check( nf90_get_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]) )

if (param%lenergy) then
call check( nf90_get_var(iu%ncid, iu%KE_orb_varid, self%ke_orbit, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%KE_spin_varid, self%ke_spin, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%PE_varid, self%pe, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_orbx_varid, self%Lorbit(1), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_orby_varid, self%Lorbit(2), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_orbz_varid, self%Lorbit(3), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_spinx_varid, self%Lspin(1), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_spiny_varid, self%Lspin(2), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_spinz_varid, self%Lspin(3), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_escapex_varid, param%Lescape(1), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_escapey_varid, param%Lescape(2), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_escapez_varid, param%Lescape(3), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%Ecollisions_varid, param%Ecollisions, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%Euntracked_varid, param%Euntracked, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%GMescape_varid, param%GMescape, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%KE_orb_varid, self%ke_orbit, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%KE_spin_varid, self%ke_spin, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%PE_varid, self%pe, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_orbx_varid, self%Lorbit(1), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_orby_varid, self%Lorbit(2), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_orbz_varid, self%Lorbit(3), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_spinx_varid, self%Lspin(1), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_spiny_varid, self%Lspin(2), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_spinz_varid, self%Lspin(3), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_escapex_varid, self%Lescape(1), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_escapey_varid, self%Lescape(2), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%L_escapez_varid, self%Lescape(3), start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%Ecollisions_varid, self%Ecollisions, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%Euntracked_varid, self%Euntracked, start=[tslot]) )
call check( nf90_get_var(iu%ncid, iu%GMescape_varid, self%GMescape, start=[tslot]) )
end if

return
Expand Down Expand Up @@ -1077,12 +1112,12 @@ module subroutine netcdf_write_hdr_system(self, iu, param)
call check( nf90_put_var(iu%ncid, iu%L_spinx_varid, self%Lspin(1), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_spiny_varid, self%Lspin(2), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_spinz_varid, self%Lspin(3), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_escapex_varid, param%Lescape(1), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_escapey_varid, param%Lescape(2), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_escapez_varid, param%Lescape(3), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%Ecollisions_varid, param%Ecollisions, start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%Euntracked_varid, param%Euntracked, start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%GMescape_varid, param%GMescape, start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_escapex_varid, self%Lescape(1), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_escapey_varid, self%Lescape(2), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%L_escapez_varid, self%Lescape(3), start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%Ecollisions_varid, self%Ecollisions, start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%Euntracked_varid, self%Euntracked, start=[tslot]) )
call check( nf90_put_var(iu%ncid, iu%GMescape_varid, self%GMescape, start=[tslot]) )
end if

return
Expand Down
Loading

0 comments on commit a90c125

Please sign in to comment.