From 9eee64132f10511b3ca911ffbf039366b72cf42a Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 22 Feb 2024 16:00:59 -0500 Subject: [PATCH] Updated the swiftest_discard_system method to save all discard bodies to the collisions.nc file --- src/collision/collision_io.f90 | 111 +++++++++++++++------------- src/collision/collision_resolve.f90 | 6 +- src/helio/helio_util.f90 | 1 + src/rmvs/rmvs_discard.f90 | 1 + src/swiftest/swiftest_discard.f90 | 78 +++++++++++++------ src/swiftest/swiftest_module.f90 | 78 +++++++++++-------- src/swiftest/swiftest_util.f90 | 52 ++++++++++++- src/symba/symba_discard.f90 | 19 +++-- src/symba/symba_util.f90 | 38 +--------- src/whm/whm_util.f90 | 1 + 10 files changed, 227 insertions(+), 158 deletions(-) diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index bcb7ba058..78d4b8d0f 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -467,62 +467,69 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) if (allocated(tp)) deallocate(tp) select case(stage) case(1) - if (.not. allocated(before%pl)) cycle - allocate(pl, source=before%pl) + if (allocated(before%pl)) allocate(pl, source=before%pl) if (allocated(before%tp)) allocate(tp, source=before%tp) case(2) - if (.not. allocated(after%pl)) cycle - allocate(pl, source=after%pl) + if (allocated(after%pl)) allocate(pl, source=after%pl) if (allocated(after%tp)) allocate(tp, source=after%tp) end select - npl = pl%nbody - - ! This ensures that there first idslot will have the first body in it, not id 0 which is the default for a new - ! idvals array - if (.not.allocated(nc%idvals)) allocate(nc%idvals, source=pl%id) - do i = 1, npl - call nc%find_idslot(pl%id(i), idslot) - call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), & - "collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: pl") - charstring = trim(adjustl(pl%info(i)%name)) - call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], & - count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: pl") - charstring = trim(adjustl(pl%info(i)%particle_type)) - call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], & - count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: pl") - call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], & - count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: pl") - call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], & - count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: pl") - call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), & - "collision_io_netcdf_write_frame_snapshot nf90_put_var Gmass_varid: pl") - call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), & - "collision_io_netcdf_write_frame_snapshot nf90_put_var radius_varid: pl") - if (param%lrotation) then - call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], & - count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid: pl") - call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i)*RAD2DEG, & - start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), & - "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid: pl") - end if - end do - - ntp = tp%nbody - do i = 1, ntp - call nc%find_idslot(tp%id(i), idslot) - call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[ idslot ]), & - "collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: tp" ) - charstring = trim(adjustl(tp%info(i)%name)) - call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], & - count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: tp" ) - charstring = trim(adjustl(tp%info(i)%particle_type)) - call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], & - count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: tp" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1, idslot, stage, eslot], & - count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: tp" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1, idslot, stage, eslot], & - count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: tp" ) - end do + + if (.not. (allocated(pl) .or. allocated(tp))) cycle + + if (allocated(pl)) then + npl = pl%nbody + ! This ensures that there first idslot will have the first body in it, not id 0 which is the default for a new + ! idvals array + if (.not.allocated(nc%idvals)) allocate(nc%idvals, source=pl%id) + do i = 1, npl + call nc%find_idslot(pl%id(i), idslot) + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: pl") + charstring = trim(adjustl(pl%info(i)%name)) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], & + count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: pl") + charstring = trim(adjustl(pl%info(i)%particle_type)) + call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], & + count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: pl") + call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], & + count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: pl") + call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], & + count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: pl") + call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var Gmass_varid: pl") + call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var radius_varid: pl") + if (param%lrotation) then + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], & + count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid: pl") + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i)*RAD2DEG, & + start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid: pl") + end if + end do + end if + + if (allocated(tp)) then + ntp = tp%nbody + ! This ensures that there first idslot will have the first body in it, not id 0 which is the default for a new + ! idvals array + if (.not.allocated(nc%idvals)) allocate(nc%idvals, source=tp%id) + do i = 1, ntp + call nc%find_idslot(tp%id(i), idslot) + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[ idslot ]), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var id_varid: tp" ) + charstring = trim(adjustl(tp%info(i)%name)) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], & + count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var name_varid: tp" ) + charstring = trim(adjustl(tp%info(i)%particle_type)) + call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], & + count=[NAMELEN, 1, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var particle_type_varid: tp" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1, idslot, stage, eslot], & + count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rh_varid: tp" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1, idslot, stage, eslot], & + count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid: tp" ) + end do + end if end do end select end select diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index cd1661dc3..bed6f5b62 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -532,7 +532,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) allocate(plsub, mold=pl) call pl%spill(plsub, lmask, ldestructive=.false.) - call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)]) + ! call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)]) ! Save the before/after snapshots select type(before => collider%before) @@ -644,7 +644,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) ! Destroy the collision list now that the collisions are resolved call plpl_collision%setup(0_I8B) - if ((nbody_system%pl_adds%nbody == 0) .and. (nbody_system%pl_discards%nbody == 0)) exit + if ((nbody_system%pl_adds%nbody == 0) .and. (.not.any(pl%ldiscard(:)))) exit if (allocated(idnew)) deallocate(idnew) nnew = nbody_system%pl_adds%nbody allocate(idnew, source=nbody_system%pl_adds%id) @@ -743,7 +743,7 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) ! Restructure the massive bodies based on the outcome of the collision call tp%rearray(nbody_system, param) - ! Discard the collider + ! Check for discards call nbody_system%tp%discard(nbody_system, param) associate(idx1 => pltp_collision%index1, idx2 => pltp_collision%index2) diff --git a/src/helio/helio_util.f90 b/src/helio/helio_util.f90 index 52b588148..c528d7669 100644 --- a/src/helio/helio_util.f90 +++ b/src/helio/helio_util.f90 @@ -28,6 +28,7 @@ module subroutine helio_util_setup_initialize_system(self, system_history, param call self%tp%h2b(self%cb) ! Make sure that the discard list gets allocated initially + call self%pl_discards%setup(0, param) call self%tp_discards%setup(0, param) call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 0b8b7d9ba..9d3e2ac1f 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -44,6 +44,7 @@ module subroutine rmvs_discard_tp(self, nbody_system, param) call swiftest_io_log_one_message(COLLISION_LOG_OUT,message) tp%ldiscard(i) = .true. tp%lmask(i) = .false. + pl%ldiscard(iplperP) = .true. call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_rh=tp%rh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=pl%id(iplperP)) end if diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index b96112e62..e9a19e374 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/src/swiftest/swiftest_discard.f90 @@ -21,33 +21,80 @@ module subroutine swiftest_discard_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals logical :: lpl_discards, ltp_discards, lpl_check, ltp_check + logical, dimension(:), allocatable :: ldiscard + integer(I4B) :: i, nstart, nend, nsub + character(len=STRMAX) :: idstr + class(swiftest_pl), allocatable :: plsub + class(swiftest_tp), allocatable :: tpsub lpl_check = allocated(self%pl_discards) ltp_check = allocated(self%tp_discards) - associate(nbody_system => self,tp => self%tp,pl => self%pl,tp_discards => self%tp_discards,pl_discards => self%pl_discards) + associate(nbody_system => self,tp => self%tp,pl => self%pl,tp_discards => self%tp_discards,pl_discards => self%pl_discards, & + npl => self%pl%nbody, ntp => self%tp%nbody, t => self%t, collision_history => self%collision_history, & + collider => self%collider) lpl_discards = .false. ltp_discards = .false. if (lpl_check .and. pl%nbody > 0) then pl%ldiscard = pl%status(:) /= ACTIVE call pl%discard(nbody_system, param) - lpl_discards = (pl_discards%nbody > 0) + lpl_discards = any(pl%ldiscard(1:npl)) end if if (ltp_check .and. tp%nbody > 0) then tp%ldiscard = tp%status(:) /= ACTIVE call tp%discard(nbody_system, param) - ltp_discards = (tp_discards%nbody > 0) + ltp_discards = any(tp%ldiscard(1:ntp)) + lpl_discards = any(pl%ldiscard(1:npl)) end if if (ltp_discards.or.lpl_discards) then - if (lpl_discards) then - if (param%lenergy) call self%conservation_report(param, lterminal=.false.) - call pl_discards%setup(0,param) - end if if (ltp_discards) then + allocate(ldiscard, source=tp%ldiscard(:)) + allocate(tpsub, mold=tp) + call tp%spill(tpsub, ldiscard, ldestructive=.true.) + nsub = tpsub%nbody + nstart = tp_discards%nbody + 1 + nend = tp_discards%nbody + nsub + call tp_discards%append(tpsub, lsource_mask=[(.true., i = 1, nsub)]) + deallocate(ldiscard) + select type(before => collider%before) + class is (swiftest_nbody_system) + if (allocated(before%tp)) deallocate(before%tp) + allocate(before%tp, source=tp_discards) + end select call tp_discards%setup(0,param) end if + + if (lpl_discards) then ! In the base integrators, massive bodies are not true discards. The discard is + ! simply used to trigger a snapshot. + if (param%lenergy) call self%conservation_report(param, lterminal=.false.) + allocate(ldiscard, source=pl%ldiscard(:)) + allocate(plsub, mold=pl) + call pl%spill(plsub, ldiscard, ldestructive=.false.) + nsub = plsub%nbody + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + nsub + call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nsub)]) + deallocate(ldiscard) + pl%ldiscard(1:npl) = .false. + ! Save the before snapshots + select type(before => collider%before) + class is (swiftest_nbody_system) + if (allocated(before%pl)) deallocate(before%pl) + allocate(before%pl, source=pl_discards) + end select + call pl_discards%setup(0,param) + end if + ! Advance the collision id number and save it + collider%maxid_collision = max(collider%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id)) + collider%maxid_collision = collider%maxid_collision + 1 + collider%collision_id = collider%maxid_collision + collider%impactors%regime = COLLRESOLVE_REGIME_MERGE + write(idstr,*) collider%collision_id + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr))) + + call collision_history%take_snapshot(param,nbody_system, t, "particle") end if end associate @@ -86,15 +133,11 @@ module subroutine swiftest_discard_tp(self, nbody_system, param) class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameter - ! Internals - logical, dimension(self%nbody) :: ldiscard - integer(I4B) :: i, nstart, nend, nsub - class(swiftest_tp), allocatable :: tpsub if (self%nbody == 0) return associate(tp => self, ntp => self%nbody, cb => nbody_system%cb, pl => nbody_system%pl, npl => nbody_system%pl%nbody, & - tp_discards => nbody_system%tp_discards) + tp_discards => nbody_system%tp_discards, pl_discards => nbody_system%pl_discards) if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. & (param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then @@ -113,15 +156,6 @@ module subroutine swiftest_discard_tp(self, nbody_system, param) if (param%lclose) then call swiftest_discard_pl_tp(tp, nbody_system, param) end if - if (any(tp%ldiscard(1:ntp))) then - ldiscard(1:ntp) = tp%ldiscard(1:ntp) - allocate(tpsub, mold=tp) - call tp%spill(tpsub, ldiscard, ldestructive=.false.) - nsub = tpsub%nbody - nstart = tp_discards%nbody + 1 - nend = tp_discards%nbody + nsub - call tp_discards%append(tpsub, lsource_mask=[(.true., i = 1, nsub)]) - end if end associate return @@ -242,7 +276,7 @@ subroutine swiftest_discard_peri_tp(tp, nbody_system, param) call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) tp%ldiscard(i) = .true. call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), & - discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) + discard_vh=tp%vh(:,i), discard_body_id=cb%id) end if end if end if diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 6350a42dd..e513623b7 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -10,27 +10,31 @@ module swiftest !! author: David A. Minton !! - !! This module serves to combine all of the Swiftest project modules under a single umbrella so that they can be accessed from individual submodule implementations - !! with a simple "use swiftest" line. + !! This module serves to combine all of the Swiftest project modules under a single umbrella so that they can be accessed from + !! individual submodule implementations with a simple "use swiftest" line. !! - !! The project structure is divided into a heirarchy of modules. The lowest level of the heirarchy are the modules called in the "use" statements below. Next the - !! "swiftest" !! modules (this one), and finally each individual integrator (and potential future integrators) sit at the top. This structure is a consequence of two - !! competing constraints: - !! 1) The desire that much of the basic functionality of the code is modular, such that new functionality can be easily added without altering too much of the basic code. + !! The project structure is divided into a heirarchy of modules. The lowest level of the heirarchy are the modules called in the + !! "use" statements below. Next the "swiftest" modules (this one), and finally each individual integrator (and potential future + !! integrators) sit at the top. This structure is a consequence of two competing constraints: + !! 1) The desire that much of the basic functionality of the code is modular, such that new functionality can be easily added + !! without altering too much of the basic code. !! 2) Adhering to Modern Fortran's typing rules. !! - !! A set of "base" types is defined in the base module. These define classes of objects, (i.e. central body, massive body, and test particles) and other major types - !! used throughout the project. However, none of the derived data types are defined with concrete type-bound procedures attached (only abstract procedures). - !! However, the *interfaces* of type-bound procedures are defined using the base types as arguments. Because of the typing rules of Modern Fortran's type-bound procedure overrides, any non-pass arguments - !! (i.e. arguments not named self) must be identical in all extended types. Because some of the basic functionality in the project is split across multiple modules, - !! we cannot define type-bound procedures in base class objects until the all interfaces are defined. In order to avoid these dependency issues and not end up with a - !! massive base class with every possibly type-bound procedure interface in the project (thus reducing the modularity of the project), the type-bound procedures are added - !! to the base types here. + !! A set of "base" types is defined in the base module. These define classes of objects, (i.e. central body, massive body, and + !! test particles) and other major types used throughout the project. However, none of the derived data types are defined with + !! concrete type-bound procedures attached (only abstract procedures). However, the *interfaces* of type-bound procedures are + !! defined using the base types as arguments. Because of the typing rules of Modern Fortran's type-bound procedure overrides, any + !! non-pass arguments(i.e. arguments not named self) must be identical in all extended types. Because some of the basic + !! functionality in the project is split across multiple modules, we cannot define type-bound procedures in base class objects + !! until the all interfaces are defined. In order to avoid these dependency issues and not end up with a massive base class with + !! every possibly type-bound procedure interface in the project (thus reducing the modularity of the project), the type-bound + !! procedures are added to the base types here. !! - !! Structuring this code this way adds somewhat to the verbosity of the code. The main thing that has to happen is that for any procedures where one wishes to make use of an - !! type-bound procedures defined for arguments at the swiftest-type level or higher, but that are passsed to base-level procedures, must have their arguments wrapped in - !! a select type(...); class is(...) construct in order to "reveal" the procedures. This is done throughout the project at the beginning of many procedures (along with - !! copious amounts of associate(...) statements, in order to help with code readibility) + !! Structuring this code this way adds somewhat to the verbosity of the code. The main thing that has to happen is that for any + !! procedures where one wishes to make use of an type-bound procedures defined for arguments at the swiftest-type level or + !! higher, but that are passsed to base-level procedures, must have their arguments wrapped in a select type(...); class is(...) + !! construct in order to "reveal" the procedures. This is done throughout the project at the beginning of many procedures (along + !! with copious amounts of associate(...) statements, in order to help with code readibility) !! !! Adapted from David E. Kaufmann's Swifter routine: module_swifter.f90 use globals @@ -54,10 +58,13 @@ module swiftest type, extends(netcdf_parameters) :: swiftest_netcdf_parameters contains - procedure :: initialize => swiftest_io_netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object - procedure :: get_valid_masks => swiftest_io_netcdf_get_valid_masks !! Gets logical masks indicating which bodies are valid pl and tp type at the current time - procedure :: open => swiftest_io_netcdf_open !! Opens a NetCDF file and does the variable inquiries to activate variable ids - procedure :: flush => swiftest_io_netcdf_flush !! Flushes a NetCDF file by closing it then opening it again + procedure :: initialize => swiftest_io_netcdf_initialize_output !! Initialize a set of parameters used to identify a + !! NetCDF output object + procedure :: get_valid_masks => swiftest_io_netcdf_get_valid_masks !! Gets logical masks indicating which bodies are valid + !! pl and tp type at the current time + procedure :: open => swiftest_io_netcdf_open !! Opens a NetCDF file and does the variable inquiries to + !! activate variable ids + procedure :: flush => swiftest_io_netcdf_flush !! Flushes a NetCDF file by closing it then opening it again #ifdef COARRAY procedure :: coclone => swiftest_coarray_coclone_nc #endif @@ -68,15 +75,19 @@ module swiftest class(swiftest_netcdf_parameters), allocatable :: nc !! NetCDF object attached to this storage object contains procedure :: dump => swiftest_io_dump_storage !! Dumps storage object contents to file - procedure :: dealloc => swiftest_util_dealloc_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 - procedure :: get_index_values => swiftest_util_get_vals_storage !! Gets the unique values of the indices of a storage object (i.e. body id or time value) - procedure :: make_index_map => swiftest_util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - procedure :: take_snapshot => swiftest_util_snapshot_system !! Takes a snapshot of the nbody_system for later file storage + procedure :: dealloc => swiftest_util_dealloc_storage !! Resets a storage object by deallocating all items and + !! resetting the frame counter to 0 + procedure :: get_index_values => swiftest_util_get_vals_storage !! Gets the unique values of the indices of a storage object + !! (i.e. body id or time value) + procedure :: make_index_map => swiftest_util_index_map_storage !! Maps body id values to storage index values so we don't + !! have to use unlimited dimensions for id + procedure :: take_snapshot => swiftest_util_snapshot_system !! Takes a snapshot of the nbody_system for later file storage final :: swiftest_final_storage end type swiftest_storage - ! The following extended types or their children should be used, where possible, as the base of any types defined in additional modules, such as new integrators. + ! The following extended types or their children should be used, where possible, as the base of any types defined in additional + ! modules, such as new integrators. type, extends(base_parameters) :: swiftest_parameters contains procedure :: dump => swiftest_io_dump_param @@ -95,7 +106,8 @@ module swiftest contains procedure :: dealloc => swiftest_util_dealloc_kin !! Deallocates all allocatable arrays #ifdef COARRAY - procedure :: coclone => swiftest_coarray_coclone_kin !! Clones the image 1 body object to all other images in the coarray structure. + procedure :: coclone => swiftest_coarray_coclone_kin !! Clones the image 1 body object to all other images in the coarray + !! structure. #endif final :: swiftest_final_kin !! Finalizes the Swiftest kinship object - deallocates all allocatables end type swiftest_kinship @@ -103,8 +115,10 @@ module swiftest type, extends(base_particle_info) :: swiftest_particle_info character(len=NAMELEN) :: name !! Non-unique name - character(len=NAMELEN) :: particle_type !! String containing a description of the particle type (e.g. Central Body, Massive Body, Test Particle) - character(len=NAMELEN) :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.) + character(len=NAMELEN) :: particle_type !! String containing a description of the particle type (e.g. Central Body, + !! Massive Body, Test Particle) + character(len=NAMELEN) :: origin_type !! String containing a description of the origin of the particle (e.g. Initial + !! Conditions, Supercatastrophic, Disruption, etc.) real(DP) :: origin_time !! The time of the particle's formation integer(I4B) :: collision_id !! The ID of the collision that formed the particle real(DP), dimension(NDIM) :: origin_rh !! The heliocentric distance vector at the time of the particle's formation @@ -115,8 +129,10 @@ module swiftest 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 :: copy => swiftest_util_copy_particle_info !! Copies one set of information object components into another, component-by-component - procedure :: set_value => swiftest_util_set_particle_info !! Sets one or more values of the particle information metadata object + procedure :: copy => swiftest_util_copy_particle_info !! Copies one set of information object components into another, + !! component-by-component + procedure :: set_value => swiftest_util_set_particle_info !! Sets one or more values of the particle information metadata + !! object end type swiftest_particle_info diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 67cac063a..bcc4a5398 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1730,7 +1730,12 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) end if ! Reindex the new list of bodies - call pl%sort("mass", ascending=.false.) + select type(pl) + class is (helio_pl) + call pl%sort("mass", ascending=.false.) + class is (whm_pl) + call pl%sort("ir3h", ascending=.false.) + end select call pl%flatten(param) call pl%set_rhill(cb) @@ -1886,7 +1891,7 @@ module subroutine swiftest_util_rearray_tp(self, nbody_system, param) return end if - ! Reset all of thes tatus flags for the remaining bodies + ! Reset all of the status flags for the remaining bodies tp%status(1:ntp) = ACTIVE do i = 1, ntp call tp%info(i)%set_value(status="ACTIVE") @@ -2439,6 +2444,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) allocate(helio_cb :: nbody_system%cb) allocate(helio_pl :: nbody_system%pl) allocate(helio_tp :: nbody_system%tp) + allocate(helio_pl :: nbody_system%pl_discards) allocate(helio_tp :: nbody_system%tp_discards) end select param%collision_model = "MERGE" @@ -2453,6 +2459,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) allocate(whm_cb :: nbody_system%cb) allocate(whm_pl :: nbody_system%pl) allocate(whm_tp :: nbody_system%tp) + allocate(whm_pl :: nbody_system%pl_discards) allocate(whm_tp :: nbody_system%tp_discards) end select param%collision_model = "MERGE" @@ -2463,6 +2470,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) allocate(rmvs_cb :: nbody_system%cb) allocate(rmvs_pl :: nbody_system%pl) allocate(rmvs_tp :: nbody_system%tp) + allocate(rmvs_pl :: nbody_system%pl_discards) allocate(rmvs_tp :: nbody_system%tp_discards) end select param%collision_model = "MERGE" @@ -2543,8 +2551,13 @@ module subroutine swiftest_util_setup_initialize_system(self, system_history, pa class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object class(swiftest_storage), allocatable, intent(inout) :: system_history !! Stores the system history between output dumps class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals + type(encounter_storage) :: encounter_history + type(collision_storage) :: collision_history + + call encounter_history%setup(4096) + call collision_history%setup(4096) + if (allocated(system_history)) then call system_history%dealloc() deallocate(system_history) @@ -2577,6 +2590,39 @@ module subroutine swiftest_util_setup_initialize_system(self, system_history, pa call nbody_system%initialize_output_file(nc, param) call nc%close() + allocate(collision_basic :: nbody_system%collider) + call nbody_system%collider%setup(nbody_system) + + if (param%lenc_save_trajectory .or. param%lenc_save_closest) then + allocate(encounter_netcdf_parameters :: encounter_history%nc) + select type(nc => encounter_history%nc) + class is (encounter_netcdf_parameters) + nc%file_name = ENCOUNTER_OUTFILE + if (.not.param%lrestart) then + call nc%initialize(param) + call nc%close() + end if + end select + allocate(nbody_system%encounter_history, source=encounter_history) + end if + + allocate(collision_netcdf_parameters :: collision_history%nc) + select type(nc => collision_history%nc) + class is (collision_netcdf_parameters) + nc%file_name = COLLISION_OUTFILE + if (param%lrestart) then + call nc%open(param) ! This will find the nc%max_idslot variable + else + call nc%initialize(param) + end if + call nc%close() + nbody_system%collider%maxid_collision = nc%max_idslot + end select + + allocate(nbody_system%collision_history, source=collision_history) + + nbody_system%collider%max_rot = MAX_ROT_SI * param%TU2S + end associate return diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 1d86fb637..7c95da30a 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -233,17 +233,16 @@ subroutine symba_discard_nonplpl(pl, nbody_system, param) call symba_discard_cb_pl(pl, nbody_system, param) end if if (param%qmin >= 0.0_DP) call symba_discard_peri_pl(pl, nbody_system, param) - if (any(pl%ldiscard(1:npl))) then - ldiscard(1:npl) = pl%ldiscard(1:npl) + ! if (any(pl%ldiscard(1:npl))) then + ! ldiscard(1:npl) = pl%ldiscard(1:npl) - allocate(plsub, mold=pl) - call pl%spill(plsub, ldiscard, ldestructive=.false.) - nsub = plsub%nbody - nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + nsub - call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nsub)]) - - end if + ! allocate(plsub, mold=pl) + ! call pl%spill(plsub, ldiscard, ldestructive=.false.) + ! nsub = plsub%nbody + ! nstart = pl_discards%nbody + 1 + ! nend = pl_discards%nbody + nsub + ! call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nsub)]) + ! end if end associate return diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 8e4b45e37..38b13978a 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -277,12 +277,7 @@ module subroutine symba_util_setup_initialize_system(self, system_history, param class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody_system object class(swiftest_storage), allocatable, intent(inout) :: system_history !! Stores the system history between output dumps class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - type(encounter_storage) :: encounter_history - type(collision_storage) :: collision_history - call encounter_history%setup(4096) - call collision_history%setup(4096) ! Call parent method associate(nbody_system => self) call helio_util_setup_initialize_system(nbody_system, system_history, param) @@ -291,32 +286,7 @@ module subroutine symba_util_setup_initialize_system(self, system_history, param call nbody_system%plpl_collision%setup(0_I8B) call nbody_system%pltp_collision%setup(0_I8B) - if (param%lenc_save_trajectory .or. param%lenc_save_closest) then - allocate(encounter_netcdf_parameters :: encounter_history%nc) - select type(nc => encounter_history%nc) - class is (encounter_netcdf_parameters) - nc%file_name = ENCOUNTER_OUTFILE - if (.not.param%lrestart) then - call nc%initialize(param) - call nc%close() - end if - end select - allocate(nbody_system%encounter_history, source=encounter_history) - end if - - allocate(collision_netcdf_parameters :: collision_history%nc) - select type(nc => collision_history%nc) - class is (collision_netcdf_parameters) - nc%file_name = COLLISION_OUTFILE - if (param%lrestart) then - call nc%open(param) ! This will find the nc%max_idslot variable - else - call nc%initialize(param) - end if - call nc%close() - end select - allocate(nbody_system%collision_history, source=collision_history) - + if (allocated(nbody_system%collider)) deallocate(nbody_system%collider) select case(param%collision_model) case("MERGE") allocate(collision_basic :: nbody_system%collider) @@ -327,12 +297,6 @@ module subroutine symba_util_setup_initialize_system(self, system_history, param end select call nbody_system%collider%setup(nbody_system) - nbody_system%collider%max_rot = MAX_ROT_SI * param%TU2S - select type(nc => collision_history%nc) - class is (collision_netcdf_parameters) - nbody_system%collider%maxid_collision = nc%max_idslot - end select - end associate return diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index cb461fb03..5b1218550 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -216,6 +216,7 @@ module subroutine whm_util_setup_initialize_system(self, system_history, param) call self%pl%flatten(param) ! Make sure that the discard list gets allocated initially + call self%pl_discards%setup(0, param) call self%tp_discards%setup(0, param) call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb)