From a6a96e989da46d6641d1050f54ee72d5150533cf Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 22 Feb 2024 08:54:15 -0500 Subject: [PATCH 1/6] Refactored to enforce line limit and started started to switch to tracking all discards through the collision module --- src/collision/collision_resolve.f90 | 62 +++++++++++++++++------------ src/collision/collision_util.f90 | 4 +- src/swiftest/swiftest_discard.f90 | 44 +++++++++++++------- 3 files changed, 68 insertions(+), 42 deletions(-) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 0d355f3c4..cd1661dc3 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -11,19 +11,19 @@ use swiftest contains - module subroutine collision_resolve_consolidate_impactors(self, nbody_system, param, idx_parent, lflag) !! author: David A. Minton !! - !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all impactors%id members, - !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. + !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all impactors%id + !! members, and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the + !! collisional outcome. implicit none ! Arguments - class(collision_impactors), intent(out) :: self !! Collision impactors object - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current run configuration parameters with Swiftest additions - integer(I4B), dimension(:), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical, intent(out) :: lflag !! Logical flag indicating whether a impactors%id was successfully created or not + class(collision_impactors),intent(out) :: self !! Collision impactors object + class(base_nbody_system),intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(in) :: param !! Current run configuration parameters with Swiftest additions + integer(I4B), dimension(:), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + logical, intent(out) :: lflag !! Logical flag indicating whether a impactors%id was successfully created or not ! Internals type collidx_array integer(I4B), dimension(:), allocatable :: id @@ -45,7 +45,8 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa ! If all of these bodies share a parent, but this is still a unique collision, move the last child ! out of the parent's position and make it the secondary body if (idx_parent(1) == idx_parent(2)) then - if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare) + if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring + ! of the kinship relationships, though it should be rare) lflag = .false. call pl%reset_kinship([idx_parent(1)]) return @@ -233,15 +234,16 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param) end associate end do - ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't - ! contain a parent body on the unique parent list. + ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique + ! pairs that don't contain a parent body on the unique parent list. ncollisions = nunique_parent + count(lplpl_unique_parent) collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))] ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them lplpl_collision(:) = .false. lplpl_collision(collision_idx(:)) = .true. - call self%spill(nbody_system%plpl_collision, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. + ! Extract any encounters that are not collisions from the list. + call self%spill(nbody_system%plpl_collision, lplpl_collision, ldestructive=.true.) end associate end select end select @@ -281,7 +283,8 @@ module subroutine collision_resolve_extract_pltp(self, nbody_system, param) ! Create a mask that contains only the pl-tp encounters that did not result in a collision, and then discard them lpltp_collision(:) = .false. lpltp_collision(collision_idx(:)) = .true. - call self%spill(nbody_system%pltp_collision, lpltp_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. + ! Extract any encounters that are not collisions from the list. + call self%spill(nbody_system%pltp_collision, lpltp_collision, ldestructive=.true.) end associate end select end select @@ -383,8 +386,9 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, & - collider => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments) + associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, & + pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, collider => nbody_system%collider, & + impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments) ! Add the impactors%id bodies to the subtraction list nimpactors = impactors%ncoll @@ -646,10 +650,12 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) allocate(idnew, source=nbody_system%pl_adds%id) mnew = sum(nbody_system%pl_adds%mass(:)) - ! Rearrange the arrays: Remove discarded bodies, add any new bodies, re-sort, and recompute all indices and encounter lists + ! Rearrange the arrays: Remove discarded bodies, add any new bodies, re-sort, and recompute all indices and + ! encounter lists call pl%rearray(nbody_system, param) - ! Destroy the add/discard list so that we don't append the same body multiple times if another collision is detected + ! Destroy the add/discard list so that we don't append the same body multiple times if another collision + ! is detected call nbody_system%pl_discards%setup(0, param) call nbody_system%pl_adds%setup(0, param) @@ -673,13 +679,17 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) end if - ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plpl_collision + ! Check whether or not any of the particles that were just added are themselves in a collision state. This will + ! generate a new plpl_collision call self%collision_check(nbody_system, param, t, dt, irec, lplpl_collision) if (.not.lplpl_collision) exit if (loop == MAXCASCADE) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT,"A runaway collisional cascade has been detected in collision_resolve_plpl.") - call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Consider reducing the step size or changing the parameters in the collisional model to reduce the number of fragments.") + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"A runaway collisional cascade has been detected in " // & + "collision_resolve_plpl.") + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Consider reducing the step size or changing the " // & + "parameters in the collisional model to reduce the " // & + "number of fragments.") call base_util_exit(FAILURE,unit=param%display_unit) end if end associate @@ -740,12 +750,12 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) ncollisions = pltp_collision%nenc write(timestr,*) t call swiftest_io_log_one_message(COLLISION_LOG_OUT, "") - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // & - "***********************************************************") - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision between test particle and massive body detected at time t = " // & - trim(adjustl(timestr))) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // & - "***********************************************************") + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"***********************************************************" // & + "***********************************************************") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision between test particle and massive body detected " // & + "at time t = " // trim(adjustl(timestr))) + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"***********************************************************" // & + "***********************************************************") do k = 1_I8B, ncollisions ! Advance the collision id number and save it diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index fab184281..49c3e1814 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -285,6 +285,7 @@ module subroutine collision_util_dealloc_snapshot(self) return end subroutine collision_util_dealloc_snapshot + module subroutine collision_util_dealloc_impactors(self) !! author: David A. Minton !! @@ -419,7 +420,8 @@ end subroutine collision_util_dealloc_basic module subroutine collision_util_reset_fragments(self) !! author: David A. Minton !! - !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) + !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, + !! radius, or other values that get set prior to the call to fraggle_generate) implicit none ! Arguments class(collision_fragments), intent(inout) :: self diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index 311bf21f1..b96112e62 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/src/swiftest/swiftest_discard.f90 @@ -25,7 +25,7 @@ module subroutine swiftest_discard_system(self, param) 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) lpl_discards = .false. ltp_discards = .false. if (lpl_check .and. pl%nbody > 0) then @@ -59,8 +59,8 @@ end subroutine swiftest_discard_system module subroutine swiftest_discard_pl(self, nbody_system, param) !! author: David A. Minton !! - !! Placeholder method for discarding massive bodies. This method does nothing except to ensure that the discard flag is set to false. - !! This method is intended to be overridden by more advanced integrators. + !! Placeholder method for discarding massive bodies. This method does nothing except to ensure that the discard flag is set + !! to false. This method is intended to be overridden by more advanced integrators. implicit none ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object @@ -87,14 +87,14 @@ module subroutine swiftest_discard_tp(self, nbody_system, param) 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(:), allocatable :: ldiscard - integer(I4B) :: npl, ntp + logical, dimension(self%nbody) :: ldiscard + integer(I4B) :: i, nstart, nend, nsub + class(swiftest_tp), allocatable :: tpsub if (self%nbody == 0) return - associate(tp => self, cb => nbody_system%cb, pl => nbody_system%pl) - ntp = tp%nbody - npl = pl%nbody + 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) 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 @@ -102,12 +102,25 @@ module subroutine swiftest_discard_tp(self, nbody_system, param) call tp%h2b(cb) end if - if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) call swiftest_discard_cb_tp(tp, nbody_system, param) - if (param%qmin >= 0.0_DP) call swiftest_discard_peri_tp(tp, nbody_system, param) - if (param%lclose) call swiftest_discard_pl_tp(tp, nbody_system, param) + if ((param%rmin >= 0.0_DP) .or. & + (param%rmax >= 0.0_DP) .or. & + (param%rmaxu >= 0.0_DP)) then + call swiftest_discard_cb_tp(tp, nbody_system, param) + end if + if (param%qmin >= 0.0_DP) then + call swiftest_discard_peri_tp(tp, nbody_system, param) + end if + if (param%lclose) then + call swiftest_discard_pl_tp(tp, nbody_system, param) + end if if (any(tp%ldiscard(1:ntp))) then - allocate(ldiscard, source=tp%ldiscard) - call tp%spill(nbody_system%tp_discards, ldiscard(1:ntp), ldestructive=.true.) + 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 @@ -275,8 +288,9 @@ subroutine swiftest_discard_pl_tp(tp, nbody_system, param) write(idstrj, *) pl%id(j) write(timestr, *) nbody_system%t write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & - // " at t = " // trim(adjustl(timestr)) + // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" & + // trim(adjustl(idstrj)) // ")" & + // " at t = " // trim(adjustl(timestr)) call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) tp%ldiscard(i) = .true. call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=nbody_system%t, discard_rh=tp%rh(:,i), & From b050701d821c12fd9f90ea4b694cf93de0446d82 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 22 Feb 2024 12:46:37 -0500 Subject: [PATCH 2/6] Fixed typo --- src/collision/collision_io.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index bae71daf5..3ad6e332f 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -420,7 +420,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) end if end do - ntp = pl%nbody + 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" ) From 7800b5fd73af2641f07b18f203738a973a5cb8a3 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 22 Feb 2024 13:00:55 -0500 Subject: [PATCH 3/6] Refactored to enforce line limit --- src/collision/collision_io.f90 | 317 +++++++++++++++++++++------------ 1 file changed, 208 insertions(+), 109 deletions(-) diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 3ad6e332f..bcb7ba058 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -125,7 +125,8 @@ end subroutine collision_io_netcdf_dump module subroutine collision_io_netcdf_initialize_output(self, param) !! author: David A. Minton !! - !! Initialize a NetCDF fragment history file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables. + !! Initialize a NetCDF fragment history file system. This is a simplified version of the main simulation output NetCDF file, + !! but with fewer variables. use, intrinsic :: ieee_arithmetic use netcdf implicit none @@ -163,102 +164,143 @@ module subroutine collision_io_netcdf_initialize_output(self, param) close(unit=LUN, status="delete") end if - call netcdf_io_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "collision_io_netcdf_initialize_output nf90_create" ) + call netcdf_io_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), & + "collision_io_netcdf_initialize_output nf90_create" ) nc%lfile_is_open = .true. ! Dimensions - call netcdf_io_check( nf90_def_dim(nc%id, nc%collision_id_varname, NF90_UNLIMITED, nc%collision_id_dimid), "collision_io_netcdf_initialize_output nf90_def_dim collision_id_dimid" ) ! Dimension to store individual collision events - call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), "collision_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers - call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - call netcdf_io_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_netcdf_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" + call netcdf_io_check( nf90_def_dim(nc%id, nc%collision_id_varname, NF90_UNLIMITED, nc%collision_id_dimid), & + "collision_io_netcdf_initialize_output nf90_def_dim collision_id_dimid" ) ! Dimension to store collision events + call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), & + "collision_io_netcdf_initialize_output nf90_def_dim space_dimid") ! 3D space dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), & + "collision_io_netcdf_initialize_output nf90_def_dim name_dimid") ! Dimension to store particle id numbers + call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), & + "collision_io_netcdf_initialize_output nf90_def_dim str_dimid") ! Dimension for string variables (character arrays) + call netcdf_io_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), & + "collision_io_netcdf_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables + ! (aka "before" vs. "after") ! Dimension coordinates - call netcdf_io_check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, [nc%collision_id_dimid], nc%collision_id_varid), "collision_io_netcdf_initialize_output nf90_def_var collision_id_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_netcdf_initialize_output nf90_def_var space_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_netcdf_initialize_output nf90_def_var name_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_netcdf_initialize_output nf90_def_var stage_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, & + [nc%collision_id_dimid], nc%collision_id_varid), & + "collision_io_netcdf_initialize_output nf90_def_var collision_id_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), & + "collision_io_netcdf_initialize_output nf90_def_var space_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, & + [nc%str_dimid, nc%name_dimid], nc%name_varid), & + "collision_io_netcdf_initialize_output nf90_def_var name_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, & + [nc%str_dimid, nc%stage_dimid], nc%stage_varid), & + "collision_io_netcdf_initialize_output nf90_def_var stage_varid") ! Variables - call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "collision_io_netcdf_initialize_output nf90_def_var id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), & + "collision_io_netcdf_initialize_output nf90_def_var id_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & - nc%collision_id_dimid, nc%time_varid), "collision_io_netcdf_initialize_output nf90_def_var time_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & - [nc%str_dimid, nc%collision_id_dimid], nc%regime_varid), "collision_io_netcdf_initialize_output nf90_def_var regime_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & - [ nc%collision_id_dimid], nc%Qloss_varid), "collision_io_netcdf_initialize_output nf90_def_var Qloss_varid") + nc%collision_id_dimid, nc%time_varid), & + "collision_io_netcdf_initialize_output nf90_def_var time_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & + [nc%str_dimid, nc%collision_id_dimid], nc%regime_varid), & + "collision_io_netcdf_initialize_output nf90_def_var regime_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & + [nc%collision_id_dimid], nc%Qloss_varid), & + "collision_io_netcdf_initialize_output nf90_def_var Qloss_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & - [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%ptype_varid), "collision_io_netcdf_initialize_output nf90_def_var ptype_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & + [nc%str_dimid, nc%name_dimid,nc%stage_dimid, nc%collision_id_dimid], nc%ptype_varid), & + "collision_io_netcdf_initialize_output nf90_def_var ptype_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rh_varid), "collision_io_netcdf_initialize_output nf90_def_var rh_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, & + [nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rh_varid), & + "collision_io_netcdf_initialize_output nf90_def_var rh_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%vh_varid), "collision_io_netcdf_initialize_output nf90_def_var vh_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, & + [nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%vh_varid), & + "collision_io_netcdf_initialize_output nf90_def_var vh_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Gmass_varid), "collision_io_netcdf_initialize_output nf90_def_var Gmass_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, & + [nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Gmass_varid), & + "collision_io_netcdf_initialize_output nf90_def_var Gmass_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%radius_varid), "collision_io_netcdf_initialize_output nf90_def_var radius_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, & + [nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%radius_varid), & + "collision_io_netcdf_initialize_output nf90_def_var radius_varid") if (param%lrotation) then call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Ip_varid), "collision_io_netcdf_initialize_output nf90_def_var Ip_varid") + [nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Ip_varid), & + "collision_io_netcdf_initialize_output nf90_def_var Ip_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rot_varid), "collision_io_netcdf_initialize_output nf90_def_var rot_varid") + [nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rot_varid), & + "collision_io_netcdf_initialize_output nf90_def_var rot_varid") end if if (param%lenergy) then - call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& - [ nc%stage_dimid, nc%collision_id_dimid], nc%KE_orb_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_orb_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type, & + [nc%stage_dimid, nc%collision_id_dimid], nc%KE_orb_varid), & + "collision_io_netcdf_initialize_output nf90_def_var KE_orb_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& - [ nc%stage_dimid, nc%collision_id_dimid], nc%KE_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_spin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, & + [nc%stage_dimid, nc%collision_id_dimid], nc%KE_spin_varid), & + "collision_io_netcdf_initialize_output nf90_def_var KE_spin_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& - [ nc%stage_dimid, nc%collision_id_dimid], nc%PE_varid), "collision_io_netcdf_initialize_output nf90_def_var PE_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, & + [nc%stage_dimid, nc%collision_id_dimid], nc%PE_varid), & + "collision_io_netcdf_initialize_output nf90_def_var PE_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type,& - [ nc%stage_dimid, nc%collision_id_dimid], nc%BE_varid), "collision_io_netcdf_initialize_output nf90_def_var BE_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%te_varname, nc%out_type,& - [ nc%stage_dimid, nc%collision_id_dimid], nc%TE_varid), "collision_io_netcdf_initialize_output nf90_def_var TE_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type, & + [nc%stage_dimid, nc%collision_id_dimid], nc%BE_varid), & + "collision_io_netcdf_initialize_output nf90_def_var BE_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%te_varname, nc%out_type, & + [nc%stage_dimid, nc%collision_id_dimid], nc%TE_varid), & + "collision_io_netcdf_initialize_output nf90_def_var TE_varid") if (param%lrotation) then - call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, & - [ nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_orbit_varid), "collision_io_netcdf_initialize_output nf90_def_var L_orbit_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, & + [nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_orbit_varid), & + "collision_io_netcdf_initialize_output nf90_def_var L_orbit_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& - [ nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var L_spin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_spin_varname,nc%out_type, & + [nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_spin_varid), & + "collision_io_netcdf_initialize_output nf90_def_var L_spin_varid") end if end if - call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_netcdf_initialize_output nf90_inquire nVariables" ) + call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), & + "collision_io_netcdf_initialize_output nf90_inquire nVariables") do varid = 1, nvar - call netcdf_io_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "collision_io_netcdf_initialize_output nf90_inquire_variable" ) + call netcdf_io_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), & + "collision_io_netcdf_initialize_output nf90_inquire_variable") select case(vartype) case(NF90_INT) - call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_INT" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), & + "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_INT") case(NF90_FLOAT) - call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), & + "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_FLOAT") case(NF90_DOUBLE) - call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), & + "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE") case(NF90_CHAR) - call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) + call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), & + "collision_io_netcdf_initialize_output nf90_def_var_fill NF90_CHAR") end select end do ! Take the file out of define mode - call netcdf_io_check( nf90_enddef(nc%id), "collision_io_netcdf_initialize_output nf90_enddef" ) + call netcdf_io_check( nf90_enddef(nc%id), "collision_io_netcdf_initialize_output nf90_enddef") ! Add in the space and stage dimension coordinates - call netcdf_io_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "collision_io_netcdf_initialize_output nf90_put_var space" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(1), start=[1,1], count=[len(nc%stage_coords(1)),1]), "collision_io_netcdf_initialize_output nf90_put_var stage 1" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(2), start=[1,2], count=[len(nc%stage_coords(2)),1]), "collision_io_netcdf_initialize_output nf90_put_var stage 2" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), & + "collision_io_netcdf_initialize_output nf90_put_var space" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(1), start=[1,1], & + count=[len(nc%stage_coords(1)),1]), "collision_io_netcdf_initialize_output nf90_put_var stage 1") + call netcdf_io_check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(2), start=[1,2], & + count=[len(nc%stage_coords(2)),1]), "collision_io_netcdf_initialize_output nf90_put_var stage 2" ) end associate end select @@ -306,44 +348,72 @@ module subroutine collision_io_netcdf_open(self, param, readonly) self%lfile_is_open = .true. ! Dimensions - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%collision_id_varname, nc%collision_id_dimid), "collision_io_netcdf_open nf90_inq_dimid collision_id_dimid" ) - call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%collision_id_dimid, nc%collision_id_varname, len=nc%max_idslot), "collision_io_find_idslot nf90_inquire_dimension max_idslot" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "collision_io_netcdf_open nf90_inq_dimid space_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "collision_io_netcdf_open nf90_inq_dimid name_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "collision_io_netcdf_open nf90_inq_dimid str_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%stage_dimname, nc%stage_dimid), "collision_io_netcdf_open nf90_inq_dimid stage_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%collision_id_varname, nc%collision_id_dimid), & + "collision_io_netcdf_open nf90_inq_dimid collision_id_dimid" ) + call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%collision_id_dimid, nc%collision_id_varname, len=nc%max_idslot),& + "collision_io_find_idslot nf90_inquire_dimension max_idslot" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), & + "collision_io_netcdf_open nf90_inq_dimid space_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), & + "collision_io_netcdf_open nf90_inq_dimid name_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), & + "collision_io_netcdf_open nf90_inq_dimid str_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%stage_dimname, nc%stage_dimid), & + "collision_io_netcdf_open nf90_inq_dimid stage_dimid" ) ! Dimension coordinates - call netcdf_io_check( nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid), "collision_io_netcdf_open nf90_inq_varid collision_id_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "collision_io_netcdf_open nf90_inq_varid space_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "collision_io_netcdf_open nf90_inq_varid name_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%stage_dimname, nc%stage_varid), "collision_io_netcdf_open nf90_inq_varid stage_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%collision_id_varname, nc%collision_id_varid), & + "collision_io_netcdf_open nf90_inq_varid collision_id_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), & + "collision_io_netcdf_open nf90_inq_varid space_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), & + "collision_io_netcdf_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%stage_dimname, nc%stage_varid), & + "collision_io_netcdf_open nf90_inq_varid stage_varid" ) ! Required Variables - call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "collision_io_netcdf_open nf90_inq_varid name_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "collision_io_netcdf_open nf90_inq_varid time_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%regime_varname, nc%regime_varid), "collision_io_netcdf_open nf90_inq_varid regime_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%Qloss_varname, nc%Qloss_varid), "collision_io_netcdf_open nf90_inq_varid Qloss_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid), "collision_io_netcdf_open nf90_inq_varid ptype_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "collision_io_netcdf_open nf90_inq_varid rh_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "collision_io_netcdf_open nf90_inq_varid vh_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), "collision_io_netcdf_open nf90_inq_varid Gmass_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "collision_io_netcdf_open nf90_inq_varid radius_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), & + "collision_io_netcdf_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), & + "collision_io_netcdf_open nf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%regime_varname, nc%regime_varid), & + "collision_io_netcdf_open nf90_inq_varid regime_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Qloss_varname, nc%Qloss_varid), & + "collision_io_netcdf_open nf90_inq_varid Qloss_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid), & + "collision_io_netcdf_open nf90_inq_varid ptype_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), & + "collision_io_netcdf_open nf90_inq_varid rh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), & + "collision_io_netcdf_open nf90_inq_varid vh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), & + "collision_io_netcdf_open nf90_inq_varid Gmass_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), & + "collision_io_netcdf_open nf90_inq_varid radius_varid" ) if (param%lrotation) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "collision_io_netcdf_open nf90_inq_varid Ip_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "collision_io_netcdf_open nf90_inq_varid rot_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), & + "collision_io_netcdf_open nf90_inq_varid Ip_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), & + "collision_io_netcdf_open nf90_inq_varid rot_varid" ) end if if (param%lenergy) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%ke_orb_varid), "collision_io_netcdf_open nf90_inq_varid ke_orb_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%pe_varname, nc%pe_varid), "collision_io_netcdf_open nf90_inq_varid pe_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%be_varname, nc%be_varid), "collision_io_netcdf_open nf90_inq_varid be_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%te_varname, nc%te_varid), "collision_io_netcdf_open nf90_inq_varid te_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid), "collision_io_netcdf_open nf90_inq_varid L_orbit_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%ke_orb_varid), & + "collision_io_netcdf_open nf90_inq_varid ke_orb_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%pe_varname, nc%pe_varid), & + "collision_io_netcdf_open nf90_inq_varid pe_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%be_varname, nc%be_varid), & + "collision_io_netcdf_open nf90_inq_varid be_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%te_varname, nc%te_varid), & + "collision_io_netcdf_open nf90_inq_varid te_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid), & + "collision_io_netcdf_open nf90_inq_varid L_orbit_varid" ) if (param%lrotation) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%ke_spin_varid), "collision_io_netcdf_open nf90_inq_varid ke_spin_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid), "collision_io_netcdf_open nf90_inq_varid L_spin_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%ke_spin_varid), & + "collision_io_netcdf_open nf90_inq_varid ke_spin_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid), & + "collision_io_netcdf_open nf90_inq_varid L_spin_varid" ) end if end if @@ -372,15 +442,21 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) select type(nc => history%nc) class is (collision_netcdf_parameters) - associate(collider => self%collider, impactors => self%collider%impactors, fragments => self%collider%fragments, eslot => self%collider%collision_id) - call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "collision_io_netcdf_write_frame_snapshot nf90_set_fill" ) + associate(collider => self%collider, impactors => self%collider%impactors, fragments => self%collider%fragments, & + eslot => self%collider%collision_id) + call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), & + "collision_io_netcdf_write_frame_snapshot nf90_set_fill" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, eslot, start=[eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var collision_id_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var time_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, eslot, start=[eslot]), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var collision_id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var time_varid" ) charstring = trim(adjustl(REGIME_NAMES(impactors%regime))) - call netcdf_io_check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var regime_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Qloss_varid, impactors%Qloss, start=[eslot] ), "collision_io_netcdf_write_frame_snapshot nf90_put_var Qloss_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], & + count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var regime_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Qloss_varid, impactors%Qloss, start=[eslot] ), & + "collision_io_netcdf_write_frame_snapshot nf90_put_var Qloss_varid" ) select type(before =>self%collider%before) class is (swiftest_nbody_system) @@ -401,48 +477,71 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) 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 + ! 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" ) + 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" ) + 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" ) - 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" ) - 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" ) - 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" ) - 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" ) + 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" ) - 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" ) + 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" ) + 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" ) + 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" ) - 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" ) - 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" ) + 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 do end select end select if (param%lenergy) then - call netcdf_io_check( nf90_put_var(nc%id, nc%ke_orb_varid, collider%ke_orbit(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_orb_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%ke_spin_varid, collider%ke_spin(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_spin_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%pe_varid, collider%pe(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%be_varid, collider%be(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%te_varid, collider%te(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid tefore" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_orbit_varid, collider%L_orbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_orbit_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_spin_varid, collider%L_spin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_spin_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ke_orb_varid, collider%ke_orbit(:), start=[ 1, eslot], & + count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_orb_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%ke_spin_varid, collider%ke_spin(:), start=[ 1, eslot], & + count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_spin_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%pe_varid, collider%pe(:), start=[ 1, eslot], & + count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%be_varid, collider%be(:), start=[ 1, eslot], & + count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%te_varid, collider%te(:), start=[ 1, eslot], & + count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid tefore" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_orbit_varid, collider%L_orbit(:,:), start=[1, 1, eslot], & + count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_orbit_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_spin_varid, collider%L_spin(:,:), start=[1, 1, eslot], & + count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_spin_varid before" ) end if call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) ) From 7d5a7df547d7bccc0cfaca4503d977b51778bc31 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 22 Feb 2024 14:43:45 -0500 Subject: [PATCH 4/6] Refactored to enforce line limit --- src/collision/collision_check.f90 | 56 ++++++++++++++++++------------- 1 file changed, 32 insertions(+), 24 deletions(-) diff --git a/src/collision/collision_check.f90 b/src/collision/collision_check.f90 index c95547a6e..02173f6a0 100644 --- a/src/collision/collision_check.f90 +++ b/src/collision/collision_check.f90 @@ -30,7 +30,8 @@ pure elemental subroutine collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot, real(DP), intent(in) :: dt !! Step size logical, intent(in) :: lvdotr !! Logical flag indicating that these two bodies are approaching in the current substep logical, intent(out) :: lcollision !! Logical flag indicating whether these two bodies will collide or not - logical, intent(out) :: lclosest !! Logical flag indicating that, while not a collision, this is the closest approach for this pair of bodies + logical, intent(out) :: lclosest !! Logical flag indicating that, while not a collision, this is the closest approach for + !! this pair of bodies ! Internals real(DP) :: r2, rlim2, a, e, q, vdotr, tcr2, dt2 @@ -39,7 +40,8 @@ pure elemental subroutine collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot, lclosest = .false. if (r2 <= rlim2) then ! checks if bodies are actively colliding in this time step lcollision = .true. - else ! if they are not actively colliding in this time step, checks if they are going to collide next time step based on velocities and q + else ! if they are not actively colliding in this time step, checks if they are going to collide next time step based on + ! velocities and q lcollision = .false. vdotr = xr * vxr + yr * vyr + zr * vzr if (lvdotr .and. (vdotr > 0.0_DP)) then @@ -67,13 +69,13 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l !! Adapted from Hal Levison's Swift routine symba5_merge.f implicit none ! Arguments - class(collision_list_plpl), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(base_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - logical, intent(out) :: lany_collision !! Returns true if cany pair of encounters resulted in a collision + class(collision_list_plpl), intent(inout) :: self !! SyMBA pl-tp encounter list object + class(base_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object + class(base_parameters),intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical, intent(out) :: lany_collision !! Returns true if cany pair of encounters resulted in a collision ! Internals logical, dimension(:), allocatable :: lcollision, lmask real(DP), dimension(NDIM) :: xr, vr @@ -113,14 +115,16 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l vr(:) = pl%vb(:, i) - pl%vb(:, j) rlim = pl%radius(i) + pl%radius(j) Gmtot = pl%Gmass(i) + pl%Gmass(j) - call collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k), lcollision(k), self%lclosest(k)) + call collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k), lcollision(k), & + self%lclosest(k)) end do lany_collision = any(lcollision(:)) lany_closest = (param%lenc_save_closest .and. any(self%lclosest(:))) if (lany_collision .or. lany_closest) then - call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary + call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into + ! barycentric coordinates do k = 1_I8B, nenc if (.not.lcollision(k) .and. .not. self%lclosest(k)) cycle i = self%index1(k) @@ -134,10 +138,12 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l self%r2(:,k) = pl%rh(:,j) + nbody_system%cb%rb(:) self%v2(:,k) = pl%vb(:,j) if (lcollision(k)) then - ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collider pair + ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a + ! collider pair if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_impactors([i,j]) - ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in + ! the step pl%lcollision([i, j]) = .true. pl%status([i, j]) = COLLIDED call pl%info(i)%set_value(status="COLLIDED") @@ -169,13 +175,13 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l !! Adapted from Hal Levison's Swift routine symba5_merge.f implicit none ! Arguments - class(collision_list_pltp), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(base_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - logical, intent(out) :: lany_collision !! Returns true if cany pair of encounters resulted in a collision + class(collision_list_pltp),intent(inout) :: self !! SyMBA pl-tp encounter list object + class(base_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical, intent(out) :: lany_collision !! Returns true if cany pair of encounters resulted in a collision ! Internals logical, dimension(:), allocatable :: lcollision, lmask real(DP), dimension(NDIM) :: xr, vr @@ -217,7 +223,8 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l j = self%index2(k) xr(:) = pl%rh(:, i) - tp%rh(:, j) vr(:) = pl%vb(:, i) - tp%vb(:, j) - call collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k), lcollision(k), self%lclosest(k)) + call collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k), & + lcollision(k), self%lclosest(k)) end do lany_collision = any(lcollision(:)) @@ -225,7 +232,8 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l if (lany_collision .or. lany_closest) then - call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary + call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into + ! barycentric coordiantes do k = 1, nenc if (.not.lcollision(k) .and. .not. self%lclosest(k)) cycle i = self%index1(k) @@ -247,8 +255,8 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l write(timestr, *) t call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_rh=tp%rh(:,j), discard_vh=tp%vh(:,j)) write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & - // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " at t = " // trim(adjustl(timestr)) + // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) & + // ")" // " at t = " // trim(adjustl(timestr)) call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) end if end do From 9eee64132f10511b3ca911ffbf039366b72cf42a Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 22 Feb 2024 16:00:59 -0500 Subject: [PATCH 5/6] 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) From 2736a1acc1c997dace5fd2559028e3dc31d3c29c Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 22 Feb 2024 16:36:16 -0500 Subject: [PATCH 6/6] Added stage and collision_id dimensions to the collision.nc output for name and id variables --- src/collision/collision_io.f90 | 44 ++++++++++++++++++++-------------- 1 file changed, 26 insertions(+), 18 deletions(-) diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 78d4b8d0f..c426310df 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -185,24 +185,32 @@ module subroutine collision_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, & [nc%collision_id_dimid], nc%collision_id_varid), & "collision_io_netcdf_initialize_output nf90_def_var collision_id_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), & + + call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, & + nc%space_dimid, nc%space_varid), & "collision_io_netcdf_initialize_output nf90_def_var space_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, & - [nc%str_dimid, nc%name_dimid], nc%name_varid), & - "collision_io_netcdf_initialize_output nf90_def_var name_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, & [nc%str_dimid, nc%stage_dimid], nc%stage_varid), & "collision_io_netcdf_initialize_output nf90_def_var stage_varid") ! Variables - call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), & - "collision_io_netcdf_initialize_output nf90_def_var id_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, & + [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%name_varid), & + "collision_io_netcdf_initialize_output nf90_def_var name_varid") + + call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, & + [nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%id_varid), & + "collision_io_netcdf_initialize_output nf90_def_var id_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & nc%collision_id_dimid, nc%time_varid), & - "collision_io_netcdf_initialize_output nf90_def_var time_varid" ) + "collision_io_netcdf_initialize_output nf90_def_var time_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & [nc%str_dimid, nc%collision_id_dimid], nc%regime_varid), & "collision_io_netcdf_initialize_output nf90_def_var regime_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & [nc%collision_id_dimid], nc%Qloss_varid), & "collision_io_netcdf_initialize_output nf90_def_var Qloss_varid") @@ -211,7 +219,6 @@ module subroutine collision_io_netcdf_initialize_output(self, param) [nc%str_dimid, nc%name_dimid,nc%stage_dimid, nc%collision_id_dimid], nc%ptype_varid), & "collision_io_netcdf_initialize_output nf90_def_var ptype_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, & [nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rh_varid), & "collision_io_netcdf_initialize_output nf90_def_var rh_varid") @@ -224,7 +231,6 @@ module subroutine collision_io_netcdf_initialize_output(self, param) [nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Gmass_varid), & "collision_io_netcdf_initialize_output nf90_def_var Gmass_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, & [nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%radius_varid), & "collision_io_netcdf_initialize_output nf90_def_var radius_varid") @@ -449,12 +455,14 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, eslot, start=[eslot]), & "collision_io_netcdf_write_frame_snapshot nf90_put_var collision_id_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), & "collision_io_netcdf_write_frame_snapshot nf90_put_var time_varid" ) charstring = trim(adjustl(REGIME_NAMES(impactors%regime))) call netcdf_io_check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], & count=[NAMELEN, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var regime_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%Qloss_varid, impactors%Qloss, start=[eslot] ), & "collision_io_netcdf_write_frame_snapshot nf90_put_var Qloss_varid" ) @@ -483,14 +491,14 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) 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 ]), & + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot, stage, eslot]), & "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") + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot, stage, eslot], & + count=[NAMELEN,1,1,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") + count=[NAMELEN,1,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], & @@ -500,7 +508,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) 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], & + 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]), & @@ -516,14 +524,14 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) 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 ]), & + call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[ idslot, stage, eslot]), & "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" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot, stage, eslot], & + count=[NAMELEN,1,1,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" ) + count=[NAMELEN,1,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], &