From 71cc5878b2793e3a6d1ede8c35848bbe226c8cee Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 30 Aug 2021 12:21:34 -0400 Subject: [PATCH] Added discard info and rearranged code to properly record discard information. Still hunting for mass loss bug --- python/swiftest/swiftest/io.py | 8 +- src/discard/discard.f90 | 5 + src/io/io.f90 | 56 ++++---- src/modules/swiftest_classes.f90 | 239 +++++++++++++++++-------------- src/modules/swiftest_globals.f90 | 119 ++++++++------- src/modules/symba_classes.f90 | 2 +- src/netcdf/netcdf.f90 | 82 ++++++++--- src/rmvs/rmvs_discard.f90 | 1 + src/setup/setup.f90 | 51 +++---- src/symba/symba_collision.f90 | 86 ++++++----- src/symba/symba_discard.f90 | 4 + src/symba/symba_io.f90 | 14 +- src/symba/symba_util.f90 | 1 + src/util/util_copy.f90 | 19 ++- src/util/util_set.f90 | 56 ++++++++ 15 files changed, 460 insertions(+), 283 deletions(-) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 0397335d6..a82c5f043 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -718,8 +718,12 @@ def clean_string_values(param, ds): ------- ds : xarray dataset with the strings cleaned up """ - ds['name'] = ds['name'].str.decode(encoding='utf-8') - ds['particle_type'] = ds['particle_type'].str.decode(encoding='utf-8') + if 'name' in ds: + ds['name'] = ds['name'].str.decode(encoding='utf-8') + if 'particle_type' in ds: + ds['particle_type'] = ds['particle_type'].str.decode(encoding='utf-8') + if 'status' in ds: + ds['status'] = ds['status'].str.decode(encoding='utf-8') if 'origin_type' in ds: ds['origin_type'] = ds['origin_type'].str.decode(encoding='utf-8') return ds diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index cf38c2bbd..e0dc8813e 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -128,6 +128,7 @@ subroutine discard_cb_tp(tp, system, param) write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. + call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then tp%status(i) = DISCARDED_RMIN write(idstr, *) tp%id(i) @@ -135,6 +136,7 @@ subroutine discard_cb_tp(tp, system, param) write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. + call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) vb2 = dot_product(tp%vb(:, i), tp%vb(:, i)) @@ -146,6 +148,7 @@ subroutine discard_cb_tp(tp, system, param) write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. + call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i)) end if end if end if @@ -195,6 +198,7 @@ subroutine discard_peri_tp(tp, system, param) write(timestr, *) param%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " perihelion distance too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. + call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) end if end if end if @@ -243,6 +247,7 @@ subroutine discard_pl_tp(tp, system, param) // " too close to massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstrj)) & // " at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. + call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) exit end if end do diff --git a/src/io/io.f90 b/src/io/io.f90 index a7cea0d7f..a7e1d5d9f 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -938,7 +938,7 @@ module subroutine io_read_in_body(self, param) ! Internals integer(I4B), parameter :: LUN = 7 !! Unit number of input file integer(I4B) :: iu = LUN - integer(I4B) :: nbody + integer(I4B) :: i, nbody logical :: is_ascii, is_pl character(len=:), allocatable :: infile real(DP) :: t @@ -974,6 +974,9 @@ module subroutine io_read_in_body(self, param) ierr = self%read_frame(iu, param) self%status(:) = ACTIVE self%lmask(:) = .true. + do i = 1, nbody + call self%info(i)%set_value(status="ACTIVE") + end do end if close(iu, err = 667, iomsg = errmsg) @@ -1547,7 +1550,7 @@ module subroutine io_write_discard(self, param) implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B), parameter :: LUN = 40 integer(I4B) :: i @@ -1561,9 +1564,17 @@ module subroutine io_write_discard(self, param) class(swiftest_body), allocatable :: pltemp character(len=STRMAX) :: errmsg, out_stat - if (param%discard_out == "") return - associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody) + + ! Record the discarded body metadata information to file + if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + call param%nciu%open(param) + call tp_discards%write_particle_info(param%nciu) + call param%nciu%close(param) + end if + + if (param%discard_out == "") return + if (nsp == 0) return if (lfirst) then out_stat = param%out_stat @@ -1589,22 +1600,22 @@ module subroutine io_write_discard(self, param) write(LUN, VECFMT, err = 667, iomsg = errmsg) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i) end do if (param%lbig_discard) then - if (param%lgr) then - allocate(pltemp, source = pl) - call pltemp%pv2v(param) - allocate(vh, source = pltemp%vh) - deallocate(pltemp) - else - allocate(vh, source = pl%vh) - end if + if (param%lgr) then + allocate(pltemp, source = pl) + call pltemp%pv2v(param) + allocate(vh, source = pltemp%vh) + deallocate(pltemp) + else + allocate(vh, source = pl%vh) + end if - write(LUN, NPLFMT) npl - do i = 1, npl - write(LUN, PLNAMEFMT, err = 667, iomsg = errmsg) pl%id(i), pl%Gmass(i), pl%radius(i) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i) - write(LUN, VECFMT, err = 667, iomsg = errmsg) vh(1, i), vh(2, i), vh(3, i) - end do - deallocate(vh) + write(LUN, NPLFMT) npl + do i = 1, npl + write(LUN, PLNAMEFMT, err = 667, iomsg = errmsg) pl%id(i), pl%Gmass(i), pl%radius(i) + write(LUN, VECFMT, err = 667, iomsg = errmsg) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i) + write(LUN, VECFMT, err = 667, iomsg = errmsg) vh(1, i), vh(2, i), vh(3, i) + end do + deallocate(vh) end if close(LUN) end associate @@ -1878,13 +1889,6 @@ module subroutine io_write_frame_system(self, param) end if end select - select type(param) - class is (symba_parameters) - param%nciu%ltrack_origin = param%lfragmentation - class default - param%nciu%ltrack_origin = .false. - end select - select case(param%out_stat) case('APPEND') call param%nciu%open(param) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index ec1a3e9df..faffcfb62 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -9,65 +9,74 @@ module swiftest_classes type :: netcdf_parameters - integer(I4B) :: out_type !! NetCDF output type (will be assigned either NF90_DOUBLE or NF90_FLOAT, depending on the user parameter) - integer(I4B) :: ncid !! NetCDF ID for the output file - integer(I4B) :: dimids(3) !! Dimensions of the NetCDF file - integer(I4B) :: time_dimid !! NetCDF ID for the time dimension - integer(I4B) :: id_dimid !! NetCDF ID for the particle id dimension - integer(I4B) :: str_dimid !! NetCDF ID for the character string dimension - integer(I4B) :: time_varid !! NetCDF ID for the time variable - integer(I4B) :: id_varid !! NetCDF ID for the particle name variable - integer(I4B) :: name_varid !! NetCDF ID for the namevariable - integer(I4B) :: ptype_varid !! NetCDF ID for the particle type variable - integer(I4B) :: npl_varid !! NetCDF ID for the number of active massive bodies variable - integer(I4B) :: ntp_varid !! NetCDF ID for the number of active test particles variable - integer(I4B) :: a_varid !! NetCDF ID for the semimajor axis variable - integer(I4B) :: e_varid !! NetCDF ID for the eccentricity variable - integer(I4B) :: inc_varid !! NetCDF ID for the inclination variable - integer(I4B) :: capom_varid !! NetCDF ID for the long. asc. node variable - integer(I4B) :: omega_varid !! NetCDF ID for the arg. periapsis variable - integer(I4B) :: capm_varid !! NetCDF ID for the mean anomaly variable - integer(I4B) :: xhx_varid !! NetCDF ID for the heliocentric position x variable - integer(I4B) :: xhy_varid !! NetCDF ID for the heliocentric position y variable - integer(I4B) :: xhz_varid !! NetCDF ID for the heliocentric position z variable - integer(I4B) :: vhx_varid !! NetCDF ID for the heliocentric velocity x variable - integer(I4B) :: vhy_varid !! NetCDF ID for the heliocentric velocity y variable - integer(I4B) :: vhz_varid !! NetCDF ID for the heliocentric velocity z variable - integer(I4B) :: Gmass_varid !! NetCDF ID for the mass variable - integer(I4B) :: rhill_varid !! NetCDF ID for the hill radius variable - integer(I4B) :: radius_varid !! NetCDF ID for the radius variable - integer(I4B) :: Ip1_varid !! NetCDF ID for the axis 1 principal moment of inertia variable - integer(I4B) :: Ip2_varid !! NetCDF ID for the axis 2 principal moment of inertia variable - integer(I4B) :: Ip3_varid !! NetCDF ID for the axis 3 principal moment of inertia variable - integer(I4B) :: rotx_varid !! NetCDF ID for the rotation x variable - integer(I4B) :: roty_varid !! NetCDF ID for the rotation y variable - integer(I4B) :: rotz_varid !! NetCDF ID for the rotation z variable - integer(I4B) :: k2_varid !! NetCDF ID for the Love number variable - integer(I4B) :: Q_varid !! NetCDF ID for the energy dissipation variable - integer(I4B) :: KE_orb_varid !! NetCDF ID for the system orbital kinetic energy variable - integer(I4B) :: KE_spin_varid !! NetCDF ID for the system spin kinetic energy variable - integer(I4B) :: PE_varid !! NetCDF ID for the system potential energy variable - integer(I4B) :: L_orbx_varid !! NetCDF ID for the system orbital angular momentum x variable - integer(I4B) :: L_orby_varid !! NetCDF ID for the system orbital angular momentum y variable - integer(I4B) :: L_orbz_varid !! NetCDF ID for the system orbital angular momentum z variable - integer(I4B) :: L_spinx_varid !! NetCDF ID for the system spin angular momentum x variable - integer(I4B) :: L_spiny_varid !! NetCDF ID for the system spin angular momentum y variable - integer(I4B) :: L_spinz_varid !! NetCDF ID for the system spin angular momentum z variable - integer(I4B) :: L_escapex_varid !! NetCDF ID for the escaped angular momentum x variable - integer(I4B) :: L_escapey_varid !! NetCDF ID for the escaped angular momentum x variable - integer(I4B) :: L_escapez_varid !! NetCDF ID for the escaped angular momentum x variable - integer(I4B) :: Ecollisions_varid !! NetCDF ID for the energy lost in collisions variable - integer(I4B) :: Euntracked_varid !! NetCDF ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) - integer(I4B) :: GMescape_varid !! NetCDF ID for the G*Mass of bodies that escape the system - logical :: ltrack_origin = .false. !! Indicate whether to track particle origin (SyMBA with Fragmentation) - integer(I4B) :: origin_type_varid !! NetCDF ID for the origin type - integer(I4B) :: origin_time_varid !! NetCDF ID for the origin type - integer(I4B) :: origin_xhx_varid !! NetCDF ID for the origin xh x component - integer(I4B) :: origin_xhy_varid !! NetCDF ID for the origin xh y component - integer(I4B) :: origin_xhz_varid !! NetCDF ID for the origin xh z component - integer(I4B) :: origin_vhx_varid !! NetCDF ID for the origin xh x component - integer(I4B) :: origin_vhy_varid !! NetCDF ID for the origin xh y component - integer(I4B) :: origin_vhz_varid !! NetCDF ID for the origin xh z component + integer(I4B) :: out_type !! NetCDF output type (will be assigned either NF90_DOUBLE or NF90_FLOAT, depending on the user parameter) + integer(I4B) :: ncid !! NetCDF ID for the output file + integer(I4B) :: dimids(3) !! Dimensions of the NetCDF file + integer(I4B) :: time_dimid !! NetCDF ID for the time dimension + integer(I4B) :: id_dimid !! NetCDF ID for the particle id dimension + integer(I4B) :: str_dimid !! NetCDF ID for the character string dimension + integer(I4B) :: time_varid !! NetCDF ID for the time variable + integer(I4B) :: id_varid !! NetCDF ID for the particle name variable + integer(I4B) :: name_varid !! NetCDF ID for the namevariable + integer(I4B) :: ptype_varid !! NetCDF ID for the particle type variable + integer(I4B) :: npl_varid !! NetCDF ID for the number of active massive bodies variable + integer(I4B) :: ntp_varid !! NetCDF ID for the number of active test particles variable + integer(I4B) :: a_varid !! NetCDF ID for the semimajor axis variable + integer(I4B) :: e_varid !! NetCDF ID for the eccentricity variable + integer(I4B) :: inc_varid !! NetCDF ID for the inclination variable + integer(I4B) :: capom_varid !! NetCDF ID for the long. asc. node variable + integer(I4B) :: omega_varid !! NetCDF ID for the arg. periapsis variable + integer(I4B) :: capm_varid !! NetCDF ID for the mean anomaly variable + integer(I4B) :: xhx_varid !! NetCDF ID for the heliocentric position x variable + integer(I4B) :: xhy_varid !! NetCDF ID for the heliocentric position y variable + integer(I4B) :: xhz_varid !! NetCDF ID for the heliocentric position z variable + integer(I4B) :: vhx_varid !! NetCDF ID for the heliocentric velocity x variable + integer(I4B) :: vhy_varid !! NetCDF ID for the heliocentric velocity y variable + integer(I4B) :: vhz_varid !! NetCDF ID for the heliocentric velocity z variable + integer(I4B) :: Gmass_varid !! NetCDF ID for the mass variable + integer(I4B) :: rhill_varid !! NetCDF ID for the hill radius variable + integer(I4B) :: radius_varid !! NetCDF ID for the radius variable + integer(I4B) :: Ip1_varid !! NetCDF ID for the axis 1 principal moment of inertia variable + integer(I4B) :: Ip2_varid !! NetCDF ID for the axis 2 principal moment of inertia variable + integer(I4B) :: Ip3_varid !! NetCDF ID for the axis 3 principal moment of inertia variable + integer(I4B) :: rotx_varid !! NetCDF ID for the rotation x variable + integer(I4B) :: roty_varid !! NetCDF ID for the rotation y variable + integer(I4B) :: rotz_varid !! NetCDF ID for the rotation z variable + integer(I4B) :: k2_varid !! NetCDF ID for the Love number variable + integer(I4B) :: Q_varid !! NetCDF ID for the energy dissipation variable + integer(I4B) :: KE_orb_varid !! NetCDF ID for the system orbital kinetic energy variable + integer(I4B) :: KE_spin_varid !! NetCDF ID for the system spin kinetic energy variable + integer(I4B) :: PE_varid !! NetCDF ID for the system potential energy variable + integer(I4B) :: L_orbx_varid !! NetCDF ID for the system orbital angular momentum x variable + integer(I4B) :: L_orby_varid !! NetCDF ID for the system orbital angular momentum y variable + integer(I4B) :: L_orbz_varid !! NetCDF ID for the system orbital angular momentum z variable + integer(I4B) :: L_spinx_varid !! NetCDF ID for the system spin angular momentum x variable + integer(I4B) :: L_spiny_varid !! NetCDF ID for the system spin angular momentum y variable + integer(I4B) :: L_spinz_varid !! NetCDF ID for the system spin angular momentum z variable + integer(I4B) :: L_escapex_varid !! NetCDF ID for the escaped angular momentum x variable + integer(I4B) :: L_escapey_varid !! NetCDF ID for the escaped angular momentum x variable + integer(I4B) :: L_escapez_varid !! NetCDF ID for the escaped angular momentum x variable + integer(I4B) :: Ecollisions_varid !! NetCDF ID for the energy lost in collisions variable + integer(I4B) :: Euntracked_varid !! NetCDF ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + integer(I4B) :: GMescape_varid !! NetCDF ID for the G*Mass of bodies that escape the system + integer(I4B) :: status_varid !! NetCDF ID for the status variable + integer(I4B) :: origin_type_varid !! NetCDF ID for the origin type + integer(I4B) :: origin_time_varid !! NetCDF ID for the origin time + integer(I4B) :: origin_xhx_varid !! NetCDF ID for the origin xh x component + integer(I4B) :: origin_xhy_varid !! NetCDF ID for the origin xh y component + integer(I4B) :: origin_xhz_varid !! NetCDF ID for the origin xh z component + integer(I4B) :: origin_vhx_varid !! NetCDF ID for the origin xh x component + integer(I4B) :: origin_vhy_varid !! NetCDF ID for the origin xh y component + integer(I4B) :: origin_vhz_varid !! NetCDF ID for the origin xh z component + integer(I4B) :: discard_time_varid !! NetCDF ID for the time of discard variable + integer(I4B) :: discard_xhx_varid !! NetCDF ID for the heliocentric position of the body at the time of discard x variable + integer(I4B) :: discard_xhy_varid !! NetCDF ID for the heliocentric position of the body at the time of discard y variable + integer(I4B) :: discard_xhz_varid !! NetCDF ID for the heliocentric position of the body at the time of discard z variable + integer(I4B) :: discard_vhx_varid !! NetCDF ID for the heliocentric velocity of the body at the time of discard x variable + integer(I4B) :: discard_vhy_varid !! NetCDF ID for the heliocentric velocity of the body at the time of discard y variable + integer(I4B) :: discard_vhz_varid !! NetCDF ID for the heliocentric velocity of the body at the time of discard z variable + integer(I4B) :: discard_body_id_varid !! NetCDF ID for the id of the other body involved in the discard + contains procedure :: close => netcdf_close !! Closes an open NetCDF file procedure :: initialize => netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object @@ -172,9 +181,10 @@ module swiftest_classes real(DP), dimension(NDIM) :: discard_vh !! The heliocentric velocity vector at the time of the particle's discard integer(I4B) :: discard_body_id !! The id of the other body involved in the discard (0 if no other body involved) contains - procedure :: dump => io_dump_particle_info !! Dumps contents of particle information to file - procedure :: read_in => io_read_in_particle_info !! Read in a particle information object from an open file - procedure :: copy => util_copy_particle_info !! Copies one set of information object components into another, component-by-component + procedure :: dump => io_dump_particle_info !! Dumps contents of particle information to file + procedure :: read_in => io_read_in_particle_info !! Read in a particle information object from an open file + procedure :: copy => util_copy_particle_info !! Copies one set of information object components into another, component-by-component + procedure :: set_value => util_set_particle_info !! Sets one or more values of the particle information metadata object end type swiftest_particle_info !******************************************************************************************************************************** @@ -234,29 +244,29 @@ module swiftest_classes !> An abstract class for a generic collection of Swiftest bodies type, abstract, extends(swiftest_base) :: swiftest_body !! Superclass that defines the generic elements of a Swiftest particle - logical :: lfirst = .true. !! Run the current step as a first - integer(I4B) :: nbody = 0 !! Number of bodies - type(swiftest_particle_info), dimension(:), allocatable :: info !! Particle metadata information - integer(I4B), dimension(:), allocatable :: id !! External identifier (unique) - integer(I4B), dimension(:), allocatable :: status !! An integrator-specific status indicator - logical, dimension(:), allocatable :: ldiscard !! Body should be discarded - logical, dimension(:), allocatable :: lmask !! Logical mask used to select a subset of bodies when performing certain operations (drift, kick, accel, etc.) - real(DP), dimension(:), allocatable :: mu !! G * (Mcb + [m]) - real(DP), dimension(:,:), allocatable :: xh !! Swiftestcentric position - real(DP), dimension(:,:), allocatable :: vh !! Swiftestcentric velocity - real(DP), dimension(:,:), allocatable :: xb !! Barycentric position - real(DP), dimension(:,:), allocatable :: vb !! Barycentric velocity - real(DP), dimension(:,:), allocatable :: ah !! Total heliocentric acceleration - real(DP), dimension(:,:), allocatable :: aobl !! Barycentric accelerations of bodies due to central body oblatenes - real(DP), dimension(:,:), allocatable :: atide !! Tanngential component of acceleration of bodies due to tides - real(DP), dimension(:,:), allocatable :: agr !! Acceleration due to post-Newtonian correction - real(DP), dimension(:), allocatable :: ir3h !! Inverse heliocentric radius term (1/rh**3) - real(DP), dimension(:), allocatable :: a !! Semimajor axis (pericentric distance for a parabolic orbit) - real(DP), dimension(:), allocatable :: e !! Eccentricity - real(DP), dimension(:), allocatable :: inc !! Inclination - real(DP), dimension(:), allocatable :: capom !! Longitude of ascending node - real(DP), dimension(:), allocatable :: omega !! Argument of pericenter - real(DP), dimension(:), allocatable :: capm !! Mean anomaly + logical :: lfirst = .true. !! Run the current step as a first + integer(I4B) :: nbody = 0 !! Number of bodies + type(swiftest_particle_info), dimension(:), allocatable :: info !! Particle metadata information + integer(I4B), dimension(:), allocatable :: id !! External identifier (unique) + integer(I4B), dimension(:), allocatable :: status !! An integrator-specific status indicator + logical, dimension(:), allocatable :: ldiscard !! Body should be discarded + logical, dimension(:), allocatable :: lmask !! Logical mask used to select a subset of bodies when performing certain operations (drift, kick, accel, etc.) + real(DP), dimension(:), allocatable :: mu !! G * (Mcb + [m]) + real(DP), dimension(:,:), allocatable :: xh !! Swiftestcentric position + real(DP), dimension(:,:), allocatable :: vh !! Swiftestcentric velocity + real(DP), dimension(:,:), allocatable :: xb !! Barycentric position + real(DP), dimension(:,:), allocatable :: vb !! Barycentric velocity + real(DP), dimension(:,:), allocatable :: ah !! Total heliocentric acceleration + real(DP), dimension(:,:), allocatable :: aobl !! Barycentric accelerations of bodies due to central body oblatenes + real(DP), dimension(:,:), allocatable :: atide !! Tanngential component of acceleration of bodies due to tides + real(DP), dimension(:,:), allocatable :: agr !! Acceleration due to post-Newtonian correction + real(DP), dimension(:), allocatable :: ir3h !! Inverse heliocentric radius term (1/rh**3) + real(DP), dimension(:), allocatable :: a !! Semimajor axis (pericentric distance for a parabolic orbit) + real(DP), dimension(:), allocatable :: e !! Eccentricity + real(DP), dimension(:), allocatable :: inc !! Inclination + real(DP), dimension(:), allocatable :: capom !! Longitude of ascending node + real(DP), dimension(:), allocatable :: omega !! Argument of pericenter + real(DP), dimension(:), allocatable :: capm !! Mean anomaly !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the !! component list, such as setup_body and util_spill contains @@ -341,10 +351,10 @@ module swiftest_classes !> An abstract class for a generic collection of Swiftest test particles type, abstract, extends(swiftest_body) :: swiftest_tp !! Superclass that defines the generic elements of a Swiftest test particle - integer(I4B), dimension(:), allocatable :: isperi !! Perihelion passage flag - real(DP), dimension(:), allocatable :: peri !! Perihelion distance - real(DP), dimension(:), allocatable :: atp !! Semimajor axis following perihelion passage - integer(I4B), dimension(:,:), allocatable :: k_pltp !! Index array used to convert flattened the body-body comparison upper triangular matrix + integer(I4B), dimension(:), allocatable :: isperi !! Perihelion passage flag + real(DP), dimension(:), allocatable :: peri !! Perihelion distance + real(DP), dimension(:), allocatable :: atp !! Semimajor axis following perihelion passage + integer(I4B), dimension(:,:), allocatable :: k_pltp !! Index array used to convert flattened the body-body comparison upper triangular matrix integer(I8B) :: npltp !! Number of pl-tp comparisons in the flattened upper triangular matrix !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the !! component list, such as setup_tp and util_spill_tp @@ -376,22 +386,22 @@ module swiftest_classes !> An abstract class for a basic Swiftest nbody system type, abstract :: swiftest_nbody_system !! This superclass contains a minimial system of a set of test particles (tp), massive bodies (pl), and a central body (cb) - class(swiftest_cb), allocatable :: cb !! Central body data structure - class(swiftest_pl), allocatable :: pl !! Massive body data structure - class(swiftest_tp), allocatable :: tp !! Test particle data structure - class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure - class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure - real(DP) :: GMtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion - real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy - real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy - real(DP) :: pe = 0.0_DP !! System potential energy - real(DP) :: te = 0.0_DP !! System total energy - 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 - 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 + class(swiftest_cb), allocatable :: cb !! Central body data structure + class(swiftest_pl), allocatable :: pl !! Massive body data structure + class(swiftest_tp), allocatable :: tp !! Test particle data structure + class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure + class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure + real(DP) :: GMtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion + real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy + real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy + real(DP) :: pe = 0.0_DP !! System potential energy + real(DP) :: te = 0.0_DP !! System total energy + 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 + 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 contains !> Each integrator will have its own version of the step procedure(abstract_step_system), deferred :: step @@ -783,7 +793,7 @@ end subroutine io_read_particle_info_system module subroutine io_write_discard(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine io_write_discard module subroutine io_toupper(string) @@ -1122,7 +1132,6 @@ module subroutine util_append_arr_logical(arr, source, nold, nsrc, lsource_mask) integer(I4B), intent(in) :: nold, nsrc !! Extend of the old array and the source array, respectively logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_logical - end interface interface @@ -1432,6 +1441,22 @@ module subroutine util_set_mu_tp(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine util_set_mu_tp + module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) + implicit none + class(swiftest_particle_info), intent(inout) :: self + character(len=*), intent(in), optional :: name !! Non-unique name + character(len=*), intent(in), optional :: particle_type !! String containing a description of the particle type (e.g. Central Body, Massive Body, Test Particle) + character(len=*), intent(in), optional :: status !! Particle status description: Active, Merged, Fragmented, etc. + character(len=*), intent(in), optional :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.) + real(DP), intent(in), optional :: origin_time !! The time of the particle's formation + real(DP), dimension(:), intent(in), optional :: origin_xh !! The heliocentric distance vector at the time of the particle's formation + real(DP), dimension(:), intent(in), optional :: origin_vh !! The heliocentric velocity vector at the time of the particle's formation + real(DP), intent(in), optional :: discard_time !! The time of the particle's discard + real(DP), dimension(:), intent(in), optional :: discard_xh !! The heliocentric distance vector at the time of the particle's discard + real(DP), dimension(:), intent(in), optional :: discard_vh !! The heliocentric velocity vector at the time of the particle's discard + integer(I4B), intent(in), optional :: discard_body_id !! The id of the other body involved in the discard (0 if no other body involved) + end subroutine util_set_particle_info + module subroutine util_set_rhill(self,cb) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index a9dd2d762..72d204a49 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -58,6 +58,7 @@ module swiftest_globals character(*), parameter :: CB_TYPE_NAME = "Central Body" character(*), parameter :: PL_TYPE_NAME = "Massive Body" character(*), parameter :: TP_TYPE_NAME = "Test Particle" + character(*), parameter :: PL_TINY_TYPE_NAME = "Semi-Interacting Massive Body" ! OpenMP Parameters integer(I4B) :: nthreads = 1 !! Number of OpenMP threads @@ -135,60 +136,68 @@ module swiftest_globals real(DP), parameter :: einsteinC = 299792458.0_DP !! Speed of light in SI units !> NetCDF variable names and constants - character(*), parameter :: NETCDF_OUTFILE = 'bin.nc' !! Default output file name - character(*), parameter :: TIME_DIMNAME = "time" !! NetCDF name of the time dimension - character(*), parameter :: ID_DIMNAME = "id" !! NetCDF name of the particle id dimension - character(*), parameter :: STR_DIMNAME = "str" !! NetCDF name of the particle id dimension - character(*), parameter :: PTYPE_VARNAME = "particle_type" ! NetCDF name of the particle type variable - character(*), parameter :: NAME_VARNAME = "name" !! NetCDF name of the particle name variable - character(*), parameter :: NPL_VARNAME = "npl" !! NetCDF name of the number of active massive bodies variable - character(*), parameter :: NTP_VARNAME = "ntp" !! NetCDF name of the number of active test particles variable - character(*), parameter :: A_VARNAME = "a" !! NetCDF name of the semimajor axis variable - character(*), parameter :: E_VARNAME = "e" !! NetCDF name of the eccentricity variable - character(*), parameter :: INC_VARNAME = "inc" !! NetCDF name of the inclination variable - character(*), parameter :: CAPOM_VARNAME = "capom" !! NetCDF name of the long. asc. node variable - character(*), parameter :: OMEGA_VARNAME = "omega" !! NetCDF name of the arg. periapsis variable - character(*), parameter :: CAPM_VARNAME = "capm" !! NetCDF name of the mean anomaly variable - character(*), parameter :: XHX_VARNAME = "xhx" !! NetCDF name of the heliocentric position x variable - character(*), parameter :: XHY_VARNAME = "xhy" !! NetCDF name of the heliocentric position y variable - character(*), parameter :: XHZ_VARNAME = "xhz" !! NetCDF name of the heliocentric position z variable - character(*), parameter :: VHX_VARNAME = "vhx" !! NetCDF name of the heliocentric velocity x variable - character(*), parameter :: VHY_VARNAME = "vhy" !! NetCDF name of the heliocentric velocity y variable - character(*), parameter :: VHZ_VARNAME = "vhz" !! NetCDF name of the heliocentric velocity z variable - character(*), parameter :: GMASS_VARNAME = "Gmass" !! NetCDF name of the mass variable - character(*), parameter :: RHILL_VARNAME = "rhill" !! NetCDF name of the hill radius variable - character(*), parameter :: RADIUS_VARNAME = "radius" !! NetCDF name of the radius variable - character(*), parameter :: IP1_VARNAME = "Ip1" !! NetCDF name of the axis 1 principal moment of inertial variable - character(*), parameter :: IP2_VARNAME = "Ip2" !! NetCDF name of the axis 2 principal moment of inertial variable - character(*), parameter :: IP3_VARNAME = "Ip3" !! NetCDF name of the axis 3 principal moment of inertial variable - character(*), parameter :: ROTX_VARNAME = "rotx" !! NetCDF name of the rotation x variable - character(*), parameter :: ROTY_VARNAME = "roty" !! NetCDF name of the rotation y variable - character(*), parameter :: ROTZ_VARNAME = "rotz" !! NetCDF name of the rotation z variable - character(*), parameter :: K2_VARNAME = "k2" !! NetCDF name of the Love number variable - character(*), parameter :: Q_VARNAME = "Q" !! NetCDF name of the energy dissipation variable - character(*), parameter :: KE_ORB_VARNAME = "KE_orb" !! NetCDF name of the system orbital kinetic energy variable - character(*), parameter :: KE_SPIN_VARNAME = "KE_spin" !! NetCDF name of the system spin kinetic energy variable - character(*), parameter :: PE_VARNAME = "PE" !! NetCDF name of the system potential energy variable - character(*), parameter :: L_ORBX_VARNAME = "L_orbx" !! NetCDF name of the orbital angular momentum x variable - character(*), parameter :: L_ORBY_VARNAME = "L_orby" !! NetCDF name of the orbital angular momentum y variable - character(*), parameter :: L_ORBZ_VARNAME = "L_orbz" !! NetCDF name of the orbital angular momentum z variable - character(*), parameter :: L_SPINX_VARNAME = "L_spinx" !! NetCDF name of the spin angular momentum x variable - character(*), parameter :: L_SPINY_VARNAME = "L_spiny" !! NetCDF name of the spin angular momentum y variable - character(*), parameter :: L_SPINZ_VARNAME = "L_spinz" !! NetCDF name of the spin angular momentum z variable - character(*), parameter :: L_ESCAPEX_VARNAME = "L_escapex" !! NetCDF name of the escaped angular momentum x variable - character(*), parameter :: L_ESCAPEY_VARNAME = "L_escapey" !! NetCDF name of the escaped angular momentum y variable - character(*), parameter :: L_ESCAPEZ_VARNAME = "L_escapez" !! NetCDF name of the escaped angular momentum z variable - character(*), parameter :: ECOLLISIONS_VARNAME = "Ecollisions" !! NetCDF name of the escaped angular momentum y variable - character(*), parameter :: EUNTRACKED_VARNAME = "Euntracked" !! NetCDF name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) - character(*), parameter :: GMESCAPE_VARNAME = "GMescape" !! NetCDF name of the G*Mass of bodies that escape the system - character(*), parameter :: ORIGIN_TYPE_VARNAME = "origin_type" - character(*), parameter :: ORIGIN_TIME_VARNAME = "origin_time" - character(*), parameter :: ORIGIN_XHX_VARNAME = "origin_xhx" - character(*), parameter :: ORIGIN_XHY_VARNAME = "origin_xhy" - character(*), parameter :: ORIGIN_XHZ_VARNAME = "origin_xhz" - character(*), parameter :: ORIGIN_VHX_VARNAME = "origin_vhx" - character(*), parameter :: ORIGIN_VHY_VARNAME = "origin_vhy" - character(*), parameter :: ORIGIN_VHZ_VARNAME = "origin_vhz" - character(*), parameter :: PL_TINY_TYPE_NAME = "Semi-Interacting Massive Body" + character(*), parameter :: NETCDF_OUTFILE = 'bin.nc' !! Default output file name + character(*), parameter :: TIME_DIMNAME = "time" !! NetCDF name of the time dimension + character(*), parameter :: ID_DIMNAME = "id" !! NetCDF name of the particle id dimension + character(*), parameter :: STR_DIMNAME = "str" !! NetCDF name of the particle id dimension + character(*), parameter :: PTYPE_VARNAME = "particle_type" !! NetCDF name of the particle type variable + character(*), parameter :: NAME_VARNAME = "name" !! NetCDF name of the particle name variable + character(*), parameter :: NPL_VARNAME = "npl" !! NetCDF name of the number of active massive bodies variable + character(*), parameter :: NTP_VARNAME = "ntp" !! NetCDF name of the number of active test particles variable + character(*), parameter :: A_VARNAME = "a" !! NetCDF name of the semimajor axis variable + character(*), parameter :: E_VARNAME = "e" !! NetCDF name of the eccentricity variable + character(*), parameter :: INC_VARNAME = "inc" !! NetCDF name of the inclination variable + character(*), parameter :: CAPOM_VARNAME = "capom" !! NetCDF name of the long. asc. node variable + character(*), parameter :: OMEGA_VARNAME = "omega" !! NetCDF name of the arg. periapsis variable + character(*), parameter :: CAPM_VARNAME = "capm" !! NetCDF name of the mean anomaly variable + character(*), parameter :: XHX_VARNAME = "xhx" !! NetCDF name of the heliocentric position x variable + character(*), parameter :: XHY_VARNAME = "xhy" !! NetCDF name of the heliocentric position y variable + character(*), parameter :: XHZ_VARNAME = "xhz" !! NetCDF name of the heliocentric position z variable + character(*), parameter :: VHX_VARNAME = "vhx" !! NetCDF name of the heliocentric velocity x variable + character(*), parameter :: VHY_VARNAME = "vhy" !! NetCDF name of the heliocentric velocity y variable + character(*), parameter :: VHZ_VARNAME = "vhz" !! NetCDF name of the heliocentric velocity z variable + character(*), parameter :: GMASS_VARNAME = "Gmass" !! NetCDF name of the mass variable + character(*), parameter :: RHILL_VARNAME = "rhill" !! NetCDF name of the hill radius variable + character(*), parameter :: RADIUS_VARNAME = "radius" !! NetCDF name of the radius variable + character(*), parameter :: IP1_VARNAME = "Ip1" !! NetCDF name of the axis 1 principal moment of inertial variable + character(*), parameter :: IP2_VARNAME = "Ip2" !! NetCDF name of the axis 2 principal moment of inertial variable + character(*), parameter :: IP3_VARNAME = "Ip3" !! NetCDF name of the axis 3 principal moment of inertial variable + character(*), parameter :: ROTX_VARNAME = "rotx" !! NetCDF name of the rotation x variable + character(*), parameter :: ROTY_VARNAME = "roty" !! NetCDF name of the rotation y variable + character(*), parameter :: ROTZ_VARNAME = "rotz" !! NetCDF name of the rotation z variable + character(*), parameter :: K2_VARNAME = "k2" !! NetCDF name of the Love number variable + character(*), parameter :: Q_VARNAME = "Q" !! NetCDF name of the energy dissipation variable + character(*), parameter :: KE_ORB_VARNAME = "KE_orb" !! NetCDF name of the system orbital kinetic energy variable + character(*), parameter :: KE_SPIN_VARNAME = "KE_spin" !! NetCDF name of the system spin kinetic energy variable + character(*), parameter :: PE_VARNAME = "PE" !! NetCDF name of the system potential energy variable + character(*), parameter :: L_ORBX_VARNAME = "L_orbx" !! NetCDF name of the orbital angular momentum x variable + character(*), parameter :: L_ORBY_VARNAME = "L_orby" !! NetCDF name of the orbital angular momentum y variable + character(*), parameter :: L_ORBZ_VARNAME = "L_orbz" !! NetCDF name of the orbital angular momentum z variable + character(*), parameter :: L_SPINX_VARNAME = "L_spinx" !! NetCDF name of the spin angular momentum x variable + character(*), parameter :: L_SPINY_VARNAME = "L_spiny" !! NetCDF name of the spin angular momentum y variable + character(*), parameter :: L_SPINZ_VARNAME = "L_spinz" !! NetCDF name of the spin angular momentum z variable + character(*), parameter :: L_ESCAPEX_VARNAME = "L_escapex" !! NetCDF name of the escaped angular momentum x variable + character(*), parameter :: L_ESCAPEY_VARNAME = "L_escapey" !! NetCDF name of the escaped angular momentum y variable + character(*), parameter :: L_ESCAPEZ_VARNAME = "L_escapez" !! NetCDF name of the escaped angular momentum z variable + character(*), parameter :: ECOLLISIONS_VARNAME = "Ecollisions" !! NetCDF name of the escaped angular momentum y variable + character(*), parameter :: EUNTRACKED_VARNAME = "Euntracked" !! NetCDF name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + character(*), parameter :: GMESCAPE_VARNAME = "GMescape" !! NetCDF name of the G*Mass of bodies that escape the system + character(*), parameter :: STATUS_VARNAME = "status" !! NetCDF name of the current status of the body variable (includes discard type) + character(*), parameter :: ORIGIN_TYPE_VARNAME = "origin_type" !! NetCDF name of the origin type variable (Initial Conditions, Disruption, etc.) + character(*), parameter :: ORIGIN_TIME_VARNAME = "origin_time" !! NetCDF name of the time of origin variable + character(*), parameter :: ORIGIN_XHX_VARNAME = "origin_xhx" !! NetCDF name of the heliocentric position of the body at the time of origin x variable + character(*), parameter :: ORIGIN_XHY_VARNAME = "origin_xhy" !! NetCDF name of the heliocentric position of the body at the time of origin y variable + character(*), parameter :: ORIGIN_XHZ_VARNAME = "origin_xhz" !! NetCDF name of the heliocentric position of the body at the time of origin z variable + character(*), parameter :: ORIGIN_VHX_VARNAME = "origin_vhx" !! NetCDF name of the heliocentric velocity of the body at the time of origin x variable + character(*), parameter :: ORIGIN_VHY_VARNAME = "origin_vhy" !! NetCDF name of the heliocentric velocity of the body at the time of origin y variable + character(*), parameter :: ORIGIN_VHZ_VARNAME = "origin_vhz" !! NetCDF name of the heliocentric velocity of the body at the time of origin z variable + character(*), parameter :: DISCARD_TIME_VARNAME = "discard_time" !! NetCDF name of the time of discard variable + character(*), parameter :: DISCARD_XHX_VARNAME = "discard_xhx" !! NetCDF name of the heliocentric position of the body at the time of discard x variable + character(*), parameter :: DISCARD_XHY_VARNAME = "discard_xhy" !! NetCDF name of the heliocentric position of the body at the time of discard y variable + character(*), parameter :: DISCARD_XHZ_VARNAME = "discard_xhz" !! NetCDF name of the heliocentric position of the body at the time of discard z variable + character(*), parameter :: DISCARD_VHX_VARNAME = "discard_vhx" !! NetCDF name of the heliocentric velocity of the body at the time of discard x variable + character(*), parameter :: DISCARD_VHY_VARNAME = "discard_vhy" !! NetCDF name of the heliocentric velocity of the body at the time of discard y variable + character(*), parameter :: DISCARD_VHZ_VARNAME = "discard_vhz" !! NetCDF name of the heliocentric velocity of the body at the time of discard z variable + character(*), parameter :: DISCARD_BODY_ID_VARNAME = "discard_body_id" !! NetCDF name of the id of the other body involved in the discard end module swiftest_globals diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 60dbd4feb..ca51def4c 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -369,7 +369,7 @@ module subroutine symba_io_write_discard(self, param) use swiftest_classes, only : swiftest_parameters implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine symba_io_write_discard module subroutine symba_kick_getacch_int_pl(self) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index e37703782..6672ee7d0 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -120,16 +120,23 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_var(self%ncid, GMESCAPE_VARNAME, self%out_type, self%time_dimid, self%GMescape_varid) ) end if - if (self%ltrack_origin) then - call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%origin_type_varid) ) - call check( nf90_def_var(self%ncid, ORIGIN_TIME_VARNAME, self%out_type, self%id_dimid, self%origin_time_varid) ) - call check( nf90_def_var(self%ncid, ORIGIN_XHX_VARNAME, self%out_type, self%id_dimid, self%origin_xhx_varid) ) - call check( nf90_def_var(self%ncid, ORIGIN_XHY_VARNAME, self%out_type, self%id_dimid, self%origin_xhy_varid) ) - call check( nf90_def_var(self%ncid, ORIGIN_XHZ_VARNAME, self%out_type, self%id_dimid, self%origin_xhz_varid) ) - call check( nf90_def_var(self%ncid, ORIGIN_VHX_VARNAME, self%out_type, self%id_dimid, self%origin_vhx_varid) ) - call check( nf90_def_var(self%ncid, ORIGIN_VHY_VARNAME, self%out_type, self%id_dimid, self%origin_vhy_varid) ) - call check( nf90_def_var(self%ncid, ORIGIN_VHZ_VARNAME, self%out_type, self%id_dimid, self%origin_vhz_varid) ) - end if + call check( nf90_def_var(self%ncid, STATUS_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%status_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%origin_type_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_TIME_VARNAME, self%out_type, self%id_dimid, self%origin_time_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_XHX_VARNAME, self%out_type, self%id_dimid, self%origin_xhx_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_XHY_VARNAME, self%out_type, self%id_dimid, self%origin_xhy_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_XHZ_VARNAME, self%out_type, self%id_dimid, self%origin_xhz_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_VHX_VARNAME, self%out_type, self%id_dimid, self%origin_vhx_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_VHY_VARNAME, self%out_type, self%id_dimid, self%origin_vhy_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_VHZ_VARNAME, self%out_type, self%id_dimid, self%origin_vhz_varid) ) + call check( nf90_def_var(self%ncid, DISCARD_TIME_VARNAME, self%out_type, self%id_dimid, self%discard_time_varid) ) + call check( nf90_def_var(self%ncid, DISCARD_XHX_VARNAME, self%out_type, self%id_dimid, self%discard_xhx_varid) ) + call check( nf90_def_var(self%ncid, DISCARD_XHY_VARNAME, self%out_type, self%id_dimid, self%discard_xhy_varid) ) + call check( nf90_def_var(self%ncid, DISCARD_XHZ_VARNAME, self%out_type, self%id_dimid, self%discard_xhz_varid) ) + call check( nf90_def_var(self%ncid, DISCARD_VHX_VARNAME, self%out_type, self%id_dimid, self%discard_vhx_varid) ) + call check( nf90_def_var(self%ncid, DISCARD_VHY_VARNAME, self%out_type, self%id_dimid, self%discard_vhy_varid) ) + 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) ) return end subroutine netcdf_initialize_output @@ -208,16 +215,24 @@ module subroutine netcdf_open(self, param) call check( nf90_inq_varid(self%ncid, GMESCAPE_VARNAME, self%GMescape_varid) ) end if - if (self%ltrack_origin) then - call check( nf90_inq_varid(self%ncid, ORIGIN_TYPE_VARNAME, self%origin_type_varid)) - call check( nf90_inq_varid(self%ncid, ORIGIN_TIME_VARNAME, self%origin_time_varid)) - call check( nf90_inq_varid(self%ncid, ORIGIN_XHX_VARNAME, self%origin_xhx_varid)) - call check( nf90_inq_varid(self%ncid, ORIGIN_XHY_VARNAME, self%origin_xhy_varid)) - call check( nf90_inq_varid(self%ncid, ORIGIN_XHZ_VARNAME, self%origin_xhz_varid)) - call check( nf90_inq_varid(self%ncid, ORIGIN_VHX_VARNAME, self%origin_vhx_varid)) - call check( nf90_inq_varid(self%ncid, ORIGIN_VHY_VARNAME, self%origin_vhy_varid)) - call check( nf90_inq_varid(self%ncid, ORIGIN_VHZ_VARNAME, self%origin_vhz_varid)) - end if + call check( nf90_inq_varid(self%ncid, STATUS_VARNAME, self%status_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_TYPE_VARNAME, self%origin_type_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_TIME_VARNAME, self%origin_time_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_XHX_VARNAME, self%origin_xhx_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_XHY_VARNAME, self%origin_xhy_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_XHZ_VARNAME, self%origin_xhz_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_VHX_VARNAME, self%origin_vhx_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_VHY_VARNAME, self%origin_vhy_varid)) + call check( nf90_inq_varid(self%ncid, ORIGIN_VHZ_VARNAME, self%origin_vhz_varid)) + + call check( nf90_inq_varid(self%ncid, DISCARD_TIME_VARNAME, self%discard_time_varid)) + call check( nf90_inq_varid(self%ncid, DISCARD_XHX_VARNAME, self%discard_xhx_varid)) + call check( nf90_inq_varid(self%ncid, DISCARD_XHY_VARNAME, self%discard_xhy_varid)) + call check( nf90_inq_varid(self%ncid, DISCARD_XHZ_VARNAME, self%discard_xhz_varid)) + call check( nf90_inq_varid(self%ncid, DISCARD_VHX_VARNAME, self%discard_vhx_varid)) + call check( nf90_inq_varid(self%ncid, DISCARD_VHY_VARNAME, self%discard_vhy_varid)) + call check( nf90_inq_varid(self%ncid, DISCARD_VHZ_VARNAME, self%discard_vhz_varid)) + call check( nf90_inq_varid(self%ncid, DISCARD_BODY_ID_VARNAME, self%discard_body_id_varid)) return end subroutine netcdf_open @@ -344,7 +359,6 @@ module subroutine netcdf_write_particle_info_base(self, iu) integer(I4B), dimension(:), allocatable :: ind character(len=:), allocatable :: charstring - select type(self) class is (swiftest_body) associate(n => self%nbody) @@ -364,9 +378,14 @@ module subroutine netcdf_write_particle_info_base(self, iu) strlen = len(charstring) call check( nf90_put_var(iu%ncid, iu%ptype_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + charstring = trim(adjustl(self%info(j)%status)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%status_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + charstring = trim(adjustl(self%info(j)%origin_type)) strlen = len(charstring) call check( nf90_put_var(iu%ncid, iu%origin_type_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + call check( nf90_put_var(iu%ncid, iu%origin_time_varid, self%info(j)%origin_time, start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_xhx_varid, self%info(j)%origin_xh(1), start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_xhy_varid, self%info(j)%origin_xh(2), start=[idslot]) ) @@ -374,6 +393,14 @@ module subroutine netcdf_write_particle_info_base(self, iu) call check( nf90_put_var(iu%ncid, iu%origin_vhx_varid, self%info(j)%origin_vh(1), start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_vhy_varid, self%info(j)%origin_vh(2), start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_vhz_varid, self%info(j)%origin_vh(3), start=[idslot]) ) + + call check( nf90_put_var(iu%ncid, iu%discard_time_varid, self%info(j)%discard_time, start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_xhx_varid, self%info(j)%discard_xh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_xhy_varid, self%info(j)%discard_xh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_xhz_varid, self%info(j)%discard_xh(3), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_vhx_varid, self%info(j)%discard_vh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_vhy_varid, self%info(j)%discard_vh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_vhz_varid, self%info(j)%discard_vh(3), start=[idslot]) ) end do end associate @@ -389,9 +416,14 @@ module subroutine netcdf_write_particle_info_base(self, iu) strlen = len(charstring) call check( nf90_put_var(iu%ncid, iu%ptype_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + charstring = trim(adjustl(self%info%status)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%status_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + charstring = trim(adjustl(self%info%origin_type)) strlen = len(charstring) call check( nf90_put_var(iu%ncid, iu%origin_type_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + call check( nf90_put_var(iu%ncid, iu%origin_time_varid, self%info%origin_time, start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_xhx_varid, self%info%origin_xh(1), start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_xhy_varid, self%info%origin_xh(2), start=[idslot]) ) @@ -399,6 +431,14 @@ module subroutine netcdf_write_particle_info_base(self, iu) call check( nf90_put_var(iu%ncid, iu%origin_vhx_varid, self%info%origin_vh(1), start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_vhy_varid, self%info%origin_vh(2), start=[idslot]) ) call check( nf90_put_var(iu%ncid, iu%origin_vhz_varid, self%info%origin_vh(3), start=[idslot]) ) + + call check( nf90_put_var(iu%ncid, iu%discard_time_varid, self%info%discard_time, start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_xhx_varid, self%info%discard_xh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_xhy_varid, self%info%discard_xh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_xhz_varid, self%info%discard_xh(3), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_vhx_varid, self%info%discard_vh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_vhy_varid, self%info%discard_vh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%discard_vhz_varid, self%info%discard_vh(3), start=[idslot]) ) end select return diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 724c107be..b179178dc 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -34,6 +34,7 @@ module subroutine rmvs_discard_tp(self, system, param) // ") is too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. + call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(iplperP)) end if end if end associate diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 4477fd788..82e549ce5 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -137,33 +137,17 @@ module subroutine setup_initialize_particle_info_system(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - logical :: ltrack_origin integer(I4B) :: i associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) - if (param%nciu%ltrack_origin) then - cb%info%origin_type = "Initial conditions" - cb%info%origin_time = param%t0 - cb%info%origin_xh(:) = 0.0_DP - cb%info%origin_vh(:) = 0.0_DP - do i = 1, self%pl%nbody - pl%info(i)%origin_type = "Initial conditions" - pl%info(i)%origin_time = param%t0 - pl%info(i)%origin_xh(:) = self%pl%xh(:,i) - pl%info(i)%origin_vh(:) = self%pl%vh(:,i) - end do - do i = 1, self%tp%nbody - tp%info(i)%origin_type = "Initial conditions" - tp%info(i)%origin_time = param%t0 - tp%info(i)%origin_xh(:) = self%tp%xh(:,i) - tp%info(i)%origin_vh(:) = self%tp%vh(:,i) - end do - end if - - cb%info%particle_type = CB_TYPE_NAME - if (npl > 0) pl%info(1:npl)%particle_type = PL_TYPE_NAME - if (ntp > 0) tp%info(1:ntp)%particle_type = TP_TYPE_NAME + call cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP]) + do i = 1, self%pl%nbody + call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%pl%xh(:,i), origin_vh=self%pl%vh(:,i)) + end do + do i = 1, self%tp%nbody + call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%tp%xh(:,i), origin_vh=self%tp%vh(:,i)) + end do end associate @@ -219,6 +203,7 @@ module subroutine setup_body(self, n, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter ! Internals integer(I4B) :: i + self%nbody = n if (n <= 0) return self%lfirst = .true. @@ -250,14 +235,22 @@ module subroutine setup_body(self, n, param) allocate(self%lmask(n)) self%id(:) = 0 - self%info(:)%name = "UNNAMED" - self%info(:)%particle_type = "UKNOWN" - self%info(:)%origin_type = "UNKNOWN" - self%info(:)%origin_time = -1.0_DP do i = 1, n - self%info(i)%origin_xh(:) = 0.0_DP - self%info(i)%origin_vh(:) = 0.0_DP + call self%info(i)%set_value(& + name = "UNNAMED", & + particle_type = "UNKNOWN", & + status = "INACTIVE", & + origin_type = "UNKNOWN", & + origin_time = -huge(1.0_DP), & + origin_xh = [0.0_DP, 0.0_DP, 0.0_DP], & + origin_vh = [0.0_DP, 0.0_DP, 0.0_DP], & + discard_time = -huge(1.0_DP), & + discard_xh = [0.0_DP, 0.0_DP, 0.0_DP], & + discard_vh = [0.0_DP, 0.0_DP, 0.0_DP], & + discard_body_id = -1 & + ) end do + self%status(:) = INACTIVE self%lmask(:) = .false. self%ldiscard(:) = .false. diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 2bf7387d6..922cadc09 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -315,6 +315,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad end do status = MERGED + call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) end select @@ -875,7 +876,7 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, real(DP), dimension(:,:), intent(in) :: xb_frag, vb_frag, rot_frag !! Fragment barycentric position, barycentric velocity, and rotation vectors integer(I4B), intent(in) :: status !! Status flag to assign to adds ! Internals - integer(I4B) :: i, ibiggest, nstart, nend, nfamily, nfrag + integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, nfamily, nfrag logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' @@ -885,26 +886,16 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, class is (symba_pl) select type(pl_discards => system%pl_discards) class is (symba_merger) - associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb) + associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody) ! Add the family bodies to the subtraction list nfamily = size(family(:)) nfrag = size(m_frag(:)) - lmask(:) = .false. - lmask(family(:)) = .true. - pl%status(family(:)) = MERGED - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nfamily - call pl_discards%append(pl, lmask) - pl%ldiscard(family(:)) = .true. - pl%lcollision(family(:)) = .true. - - ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = nfamily - + ! Setup new bodies allocate(plnew, mold=pl) call plnew%setup(nfrag, param) - ibiggest = family(maxloc(pl%Gmass(family(:)), dim=1)) + ibiggest = family(maxloc(pl%Gmass(family(:)), dim=1)) + ismallest = family(minloc(pl%Gmass(family(:)), dim=1)) ! Copy over identification, information, and physical properties of the new bodies from the fragment list plnew%id(1:nfrag) = id_frag(1:nfrag) @@ -923,40 +914,54 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, select case(status) case(DISRUPTION) - plnew%info(1:nfrag)%origin_type = "Disruption" plnew%status(1:nfrag) = NEW_PARTICLE - plnew%info(1:nfrag)%origin_time = param%t do i = 1, nfrag write(newname, FRAGFMT) id_frag(i) - plnew%info(i)%name = newname - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i)) + end do + do i = 1, nfamily + if (family(i) == ibiggest) then + iother = ismallest + else + iother = ibiggest + end if + call pl%info(family(i))%set_value(status="Disruption", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(SUPERCATASTROPHIC) - plnew%info(1:nfrag)%origin_type = "Supercatastrophic" plnew%status(1:nfrag) = NEW_PARTICLE - plnew%info(1:nfrag)%origin_time = param%t do i = 1, nfrag write(newname, FRAGFMT) id_frag(i) - plnew%info(i)%name = newname - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i)) + end do + do i = 1, nfamily + if (family(i) == ibiggest) then + iother = ismallest + else + iother = ibiggest + end if + call pl%info(family(i))%set_value(status="Supercatastrophic", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(HIT_AND_RUN_DISRUPT) call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE - plnew%status(2:nfrag) = NEW_PARTICLE - plnew%info(2:nfrag)%origin_type = "Hit and run fragment" - plnew%info(2:nfrag)%origin_time = param%t do i = 2, nfrag write(newname, FRAGFMT) id_frag(i) - plnew%info(i)%name = newname - plnew%info(i)%origin_xh(:) = plnew%xh(:,i) - plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + call plnew%info(i)%set_value(origin_type="Hit and run fragmention", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i)) end do + do i = 1, nfamily + if (family(i) == ibiggest) cycle + iother = ibiggest + call pl%info(family(i))%set_value(status="Hit and run fragmention", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + end do case(MERGED) call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE + do i = 1, nfamily + if (family(i) == ibiggest) cycle + + iother = ibiggest + call pl%info(family(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + end do end select if (param%lrotation) then @@ -978,13 +983,28 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, plnew%levelg(1:nfrag) = pl%levelg(ibiggest) plnew%levelm(1:nfrag) = pl%levelm(ibiggest) - ! Append the new merged body to the list and record how many we made + ! Append the new merged body to the list nstart = pl_adds%nbody + 1 nend = pl_adds%nbody + plnew%nbody call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) + ! Record how many bodies were added in this event pl_adds%ncomp(nstart:nend) = plnew%nbody - call plnew%setup(0, param) + + ! Add the discarded bodies to the discard list + pl%status(family(:)) = MERGED + lmask(:) = .false. + lmask(family(:)) = .true. + pl%ldiscard(family(:)) = .true. + pl%lcollision(family(:)) = .true. + + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nfamily + call pl_discards%append(pl, lmask) + + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = nfamily + deallocate(plnew) end associate end select diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index f812dc664..877f7cc95 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -37,6 +37,7 @@ subroutine symba_discard_cb_pl(pl, system, param) write(idstr, *) pl%id(i) write(timestr, *) param%t write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. pl%lcollision(i) = .false. @@ -44,6 +45,7 @@ subroutine symba_discard_cb_pl(pl, system, param) write(idstr, *) pl%id(i) write(timestr, *) param%t write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) vb2 = dot_product(pl%vb(:,i), pl%vb(:,i)) @@ -55,6 +57,7 @@ subroutine symba_discard_cb_pl(pl, system, param) write(idstr, *) pl%id(i) write(timestr, *) param%t write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) end if end if end if @@ -275,6 +278,7 @@ subroutine symba_discard_peri_pl(pl, system, param) write(timestr, *) param%t write(idstr, *) pl%id(i) write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ") perihelion distance too small at t = " // trim(adjustl(timestr)) + call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=system%cb%id) end if end if end if diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 8b64950e9..71717aabe 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -182,7 +182,7 @@ end subroutine symba_io_param_writer module subroutine symba_io_write_discard(self, param) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B), parameter :: LUN = 40 integer(I4B) :: iadd, isub, j, nsub, nadd @@ -196,13 +196,21 @@ module subroutine symba_io_write_discard(self, param) class(swiftest_body), allocatable :: pltemp character(STRMAX) :: errmsg, out_stat - if (param%discard_out == "") return - associate(pl => self%pl, npl => self%pl%nbody, pl_adds => self%pl_adds) + if (self%tp_discards%nbody > 0) call io_write_discard(self, param) select type(pl_discards => self%pl_discards) class is (symba_merger) if (pl_discards%nbody == 0) return + + ! Record the discarded body metadata information to file + if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + call param%nciu%open(param) + call pl_discards%write_particle_info(param%nciu) + call param%nciu%close(param) + end if + + if (param%discard_out == "") return if (lfirst) then out_stat = param%out_stat else diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index c804aebb7..ce29e22aa 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -454,6 +454,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! Reset all of the status flags for this body where(pl%status(1:npl) /= INACTIVE) pl%status(1:npl) = ACTIVE + pl%info(1:npl)%status = "ACTIVE" pl%ldiscard(1:npl) = .false. pl%lcollision(1:npl) = .false. pl%lmask(1:npl) = .true. diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index 0d29cc92f..60d93e45e 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -39,12 +39,19 @@ module subroutine util_copy_particle_info(self, source) class(swiftest_particle_info), intent(inout) :: self class(swiftest_particle_info), intent(in) :: source - self%name = source%name - self%particle_type = source%particle_type - self%origin_type = source%origin_type - self%origin_time = source%origin_time - self%origin_xh(:) = source%origin_xh(:) - self%origin_vh(:) = source%origin_vh(:) + call self%set_value(& + name = source%name, & + particle_type = source%particle_type, & + status = source%status, & + origin_type = source%origin_type, & + origin_time = source%origin_time, & + origin_xh = source%origin_xh(:), & + origin_vh = source%origin_vh(:), & + discard_time = source%discard_time, & + discard_xh = source%discard_xh(:), & + discard_vh = source%discard_vh(:), & + discard_body_id = source%discard_body_id & + ) return end subroutine util_copy_particle_info diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index 4e686eec8..48c0006c2 100644 --- a/src/util/util_set.f90 +++ b/src/util/util_set.f90 @@ -98,6 +98,62 @@ module subroutine util_set_mu_tp(self, cb) return end subroutine util_set_mu_tp + module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) + !! author: David A. Minton + !! + !! Sets one or more values of the particle information metadata object + implicit none + ! Arguments + class(swiftest_particle_info), intent(inout) :: self + character(len=*), intent(in), optional :: name !! Non-unique name + character(len=*), intent(in), optional :: particle_type !! String containing a description of the particle type (e.g. Central Body, Massive Body, Test Particle) + character(len=*), intent(in), optional :: status !! Particle status description: ACTIVE, MERGED, FRAGMENTED, etc. + character(len=*), intent(in), optional :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.) + real(DP), intent(in), optional :: origin_time !! The time of the particle's formation + real(DP), dimension(:), intent(in), optional :: origin_xh !! The heliocentric distance vector at the time of the particle's formation + real(DP), dimension(:), intent(in), optional :: origin_vh !! The heliocentric velocity vector at the time of the particle's formation + real(DP), intent(in), optional :: discard_time !! The time of the particle's discard + real(DP), dimension(:), intent(in), optional :: discard_xh !! The heliocentric distance vector at the time of the particle's discard + real(DP), dimension(:), intent(in), optional :: discard_vh !! The heliocentric velocity vector at the time of the particle's discard + integer(I4B), intent(in), optional :: discard_body_id !! The id of the other body involved in the discard (0 if no other body involved) + + if (present(name)) then + self%name = name + end if + if (present(particle_type)) then + self%particle_type = particle_type + end if + if (present(origin_type)) then + self%origin_type = origin_type + end if + if (present(origin_time)) then + self%origin_time = origin_time + end if + if (present(origin_xh)) then + self%origin_xh(:) = origin_xh(:) + end if + if (present(origin_vh)) then + self%origin_vh(:) = origin_vh(:) + end if + if (present(discard_time)) then + self%discard_time = discard_time + end if + if (present(status)) then + self%status = status + end if + if (present(discard_xh)) then + self%discard_xh(:) = discard_xh(:) + end if + if (present(discard_vh)) then + self%discard_vh(:) = discard_vh(:) + end if + if (present(discard_body_id)) then + self%discard_body_id = discard_body_id + end if + + return + end subroutine util_set_particle_info + module subroutine util_set_rhill(self,cb) !! author: David A. Minton