diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 1e7cd5cd6..813bd8dfb 100755 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -54,13 +54,13 @@ capom_pl = default_rng().uniform(0.0, 360.0, npl) omega_pl = default_rng().uniform(0.0, 360.0, npl) capm_pl = default_rng().uniform(0.0, 360.0, npl) -GM_pl = (np.array([6e23, 8e23, 1e24, 3e24, 5e24]) / sim.param['MU2KG']) * sim.GU -R_pl = np.full(npl, (3 * (GM_pl / sim.GU) / (4 * np.pi * density_pl)) ** (1.0 / 3.0)) -Rh_pl = a_pl * ((GM_pl) / (3 * sim.GU)) ** (1.0 / 3.0) +M_pl = np.array([6e23, 8e23, 1e24, 3e24, 5e24]) * sim.KG2MU +R_pl = np.full(npl, (3 * M_pl/ (4 * np.pi * density_pl)) ** (1.0 / 3.0)) Ip_pl = np.full((npl,3),0.4,) rot_pl = np.zeros((npl,3)) +mtiny = 1.01 * np.max(M_pl) -sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, Gmass=GM_pl, radius=R_pl, rhill=Rh_pl, Ip=Ip_pl, rot=rot_pl) +sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, mass=M_pl, radius=R_pl, Ip=Ip_pl, rot=rot_pl) # Add 10 user-defined test particles. ntp = 10 @@ -74,7 +74,7 @@ capm_tp = default_rng().uniform(0.0, 360.0, ntp) sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp) -sim.set_parameter(tstart=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=0) +sim.set_parameter(tstart=0.0, tstop=1.0e3, dt=0.01, istep_out=100, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny) # Display the run configuration parameters. sim.get_parameter() sim.save() diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 921ea8308..38b1735d7 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -123,8 +123,6 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) real(DP), dimension(NDIM) :: L_total_now, L_orbit_now, L_spin_now real(DP) :: ke_orbit_now, ke_spin_now, pe_now, E_orbit_now, be_now, be_cb_now, be_cb_orig, te_now real(DP) :: GMtot_now - character(len=STRMAX) :: errmsg - integer(I4B), parameter :: EGYIU = 72 character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5, "; DE_orbit/|E0| = ", ES12.5, "; DE_total/|E0| = ", ES12.5, "; DM/M0 = ", ES12.5)' associate(nbody_system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit) @@ -204,10 +202,6 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) end associate return - - 667 continue - write(*,*) "Error writing energy and momentum tracking file: " // trim(adjustl(errmsg)) - call base_util_exit(FAILURE) end subroutine swiftest_io_conservation_report @@ -226,7 +220,6 @@ module subroutine swiftest_io_display_run_information(self, param, integration_t real(DP) :: tfrac !! Fraction of total simulation time completed type(progress_bar), save :: pbar !! Object used to print out a progress bar character(len=64) :: pbarmessage - character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)' #ifdef COARRAY character(*), parameter :: co_statusfmt = '("Image: ",I4, "; Time = ", ES12.5, "; fraction done = ", F6.3, ' // & '"; Number of active pl, tp = ", I6, ", ", I6)' @@ -465,7 +458,6 @@ module subroutine swiftest_io_get_args(integrator, param_file_name, display_styl character(len=STRMAX), dimension(:), allocatable :: arg integer(I4B), dimension(:), allocatable :: ierr integer :: i,narg - character(len=*),parameter :: linefmt = '(A)' narg = command_argument_count() if (narg > 0) then @@ -653,6 +645,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param) ! Internals integer(I4B) :: itmax, idmax, tslot real(DP), dimension(:), allocatable :: vals + logical, dimension(:), allocatable :: plmask, tpmask real(DP), dimension(1) :: rtemp real(DP), dimension(NDIM) :: rot0, Ip0, L real(DP) :: mass0 @@ -689,7 +682,8 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param) self%L_total_orig(:) = self%L_orbit_orig(:) + self%L_spin_orig(:) call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,tslot], count=[idmax,1]), "netcdf_io_get_t0_values_system Gmass_varid" ) - self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + call nc%get_valid_masks(plmask,tpmask) + self%GMtot_orig = vals(1) + sum(vals(2:idmax), plmask(:)) cb%GM0 = vals(1) cb%dGM = cb%Gmass - cb%GM0 @@ -1030,7 +1024,7 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) end subroutine swiftest_io_netcdf_open - module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask) + module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask, plmmask, Gmtiny) !! author: David A. Minton !! !! Given an open NetCDF, returns logical masks indicating which bodies in the body arrays are active pl and tp type at the current time. @@ -1039,15 +1033,17 @@ module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask) use, intrinsic :: ieee_arithmetic implicit none ! Arguments - class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - logical, dimension(:), allocatable, intent(out) :: plmask !! Logical mask indicating which bodies are massive bodies - logical, dimension(:), allocatable, intent(out) :: tpmask !! Logical mask indicating which bodies are test particles + class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + logical, dimension(:), allocatable, intent(out) :: plmask !! Logical mask indicating which bodies are massive bodies + logical, dimension(:), allocatable, intent(out) :: tpmask !! Logical mask indicating which bodies are test particles + logical, dimension(:), allocatable, intent(out), optional :: plmmask !! Logical mask indicating which bodies are fully interacting massive bodies + real(DP), intent(in), optional :: Gmtiny !! The cutoff G*mass between semi-interacting and fully interacting massive bodies ! Internals real(DP), dimension(:), allocatable :: Gmass, a real(DP), dimension(:,:), allocatable :: rh integer(I4B), dimension(:), allocatable :: body_status logical, dimension(:), allocatable :: lvalid - integer(I4B) :: idmax, status, i + integer(I4B) :: idmax, status logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily @@ -1094,6 +1090,14 @@ module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask) ! Select only active bodies plmask(:) = plmask(:) .and. lvalid(:) tpmask(:) = tpmask(:) .and. lvalid(:) + + if (present(plmmask) .and. present(Gmtiny)) then + allocate(plmmask, source=plmask) + where(plmask(:)) + plmmask = Gmass(:) > Gmtiny + endwhere + end if + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) end associate @@ -1316,10 +1320,6 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier ierr = 0 return - - 667 continue - write(*,*) "Error reading nbody_system frame in netcdf_io_read_frame_system" - end function swiftest_io_netcdf_read_frame_system @@ -1335,21 +1335,15 @@ module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for reading a NetCDF dataset to file class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: status, idmax + integer(I4B) :: status logical, dimension(:), allocatable :: plmask, tpmask, plmmask - real(DP), dimension(:), allocatable :: Gmtemp associate(tslot => nc%tslot) call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar time_varid" ) - call nc%get_valid_masks(plmask, tpmask) - if (param%lmtiny_pl) then - idmax = size(plmask) - allocate(plmmask(idmax)) - allocate(Gmtemp(idmax)) - call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, Gmtemp, start=[1,tslot], count=[idmax,1]), "netcdf_io_read_hdr_system nf90_getvar Gmass_varid" ) - where(Gmtemp(:) /= Gmtemp(:)) Gmtemp(:) = 0.0_DP - plmmask(:) = plmask(:) .and. Gmtemp(:) > param%GMTINY + call nc%get_valid_masks(plmask, tpmask, plmmask, param%GMTINY) + else + call nc%get_valid_masks(plmask, tpmask) end if status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid) @@ -1624,7 +1618,7 @@ module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param) class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i, j, idslot, old_mode + integer(I4B) :: i, j, idslot, old_mode, tmp integer(I4B), dimension(:), allocatable :: ind real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs real(DP) :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf @@ -1729,7 +1723,7 @@ module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param) end select #endif - call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "netcdf_io_write_frame_body nf90_set_fill old_mode" ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp), "netcdf_io_write_frame_body nf90_set_fill old_mode" ) return end subroutine swiftest_io_netcdf_write_frame_body @@ -1745,7 +1739,7 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to for writing a NetCDF dataset to file class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: idslot, old_mode + integer(I4B) :: idslot, old_mode, tmp associate(tslot => nc%tslot) call self%write_info(nc, param) @@ -1766,7 +1760,7 @@ module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param) call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:) * RAD2DEG, start=[1, idslot, tslot], count=[NDIM,1,1]), "swiftest_io_netcdf_write_frame_cby nf90_put_var cb rot_varid" ) end if - call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode), "swiftest_io_netcdf_write_frame_cb nf90_set_fill old_mode" ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp), "swiftest_io_netcdf_write_frame_cb nf90_set_fill old_mode" ) end associate return @@ -1848,7 +1842,7 @@ module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i, j, idslot, old_mode + integer(I4B) :: i, j, idslot, old_mode, tmp integer(I4B), dimension(:), allocatable :: ind character(len=NAMELEN) :: charstring character(len=NAMELEN), dimension(self%nbody) :: origin_type @@ -1892,7 +1886,7 @@ module subroutine swiftest_io_netcdf_write_info_body(self, nc, param) end associate end select - call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) ) return end subroutine swiftest_io_netcdf_write_info_body @@ -1906,7 +1900,7 @@ module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param) class(swiftest_netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: idslot, old_mode + integer(I4B) :: idslot, old_mode, tmp character(len=NAMELEN) :: charstring ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables @@ -1935,7 +1929,7 @@ module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param) call netcdf_io_check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info%discard_rh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb discard_rh_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info%discard_vh(:), start=[1, idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var cb discard_vh_varid" ) end if - call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) + call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) ) return end subroutine swiftest_io_netcdf_write_info_cb @@ -1989,7 +1983,6 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i character(*),parameter :: linefmt = '(A)' !! Format code for simple text string integer(I4B) :: nseeds, nseeds_from_file logical :: seed_set = .false. !! Is the random seed set in the input file? - character(len=:), allocatable :: integrator real(DP) :: tratio, y #ifdef COARRAY type(swiftest_parameters), codimension[*], save :: coparam @@ -2517,9 +2510,6 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals - character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values - character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values - character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values integer(I4B) :: nseeds associate(param => self) @@ -2588,7 +2578,6 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i iomsg = "UDIO not implemented" end associate - 667 continue return end subroutine swiftest_io_param_writer @@ -2935,7 +2924,7 @@ module subroutine swiftest_io_read_in_system(self, nc, param) if (ierr /=0) call base_util_exit(FAILURE) end if - param%loblatecb = ((self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP)) + param%loblatecb = ((abs(self%cb%j2rp2) > 0.0_DP) .or. (abs(self%cb%j4rp4) > 0.0_DP)) if (.not.param%loblatecb) then if (allocated(self%pl%aobl)) deallocate(self%pl%aobl) if (allocated(self%tp%aobl)) deallocate(self%tp%aobl) @@ -3083,7 +3072,6 @@ module subroutine swiftest_io_read_in_param(self, param_file_name) call self%reader(LUN, iotype= "none", v_list=[""], iostat = ierr, iomsg = errmsg) if (ierr == 0) return - 667 continue write(self%display_unit,*) "Error reading parameter file: " // trim(adjustl(errmsg)) call base_util_exit(FAILURE) end subroutine swiftest_io_read_in_param @@ -3177,7 +3165,7 @@ module subroutine swiftest_io_initialize_output_file_system(self, nc, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals logical, save :: lfirst = .true. !! Flag to determine if this is the first call of this method - character(len=STRMAX) :: errmsg + character(len=2*STRMAX) :: errmsg logical :: fileExists associate (pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) @@ -3191,13 +3179,13 @@ module subroutine swiftest_io_initialize_output_file_system(self, nc, param) select case(param%out_stat) case('APPEND') if (.not.fileExists) then - errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN" + errmsg = trim(adjustl(param%outfile)) // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN" goto 667 end if call nc%open(param) case('NEW') if (fileExists) then - errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" + errmsg = trim(adjustl(param%outfile)) // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" goto 667 end if call nc%initialize(param)