diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index e37adf13d..61650e700 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -98,23 +98,24 @@ module subroutine fraggle_io_initialize_output(self, param) call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "fraggle_io_initialize nf90_def_var rot_varid") + + if (param%lenergy) then - call check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "fraggle_io_initialize_output nf90_def_var KE_orb_varid") - - call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "fraggle_io_initialize_output nf90_def_var KE_spin_varid" ) - - call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "fraggle_io_initialize_output nf90_def_var PE_varid" ) + call check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "fraggle_io_initialize_output nf90_def_var KE_orb_varid") - call check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "fraggle_io_initialize_output nf90_def_var L_orb_varid" ) + call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "fraggle_io_initialize_output nf90_def_var KE_spin_varid" ) - call check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_spin_varid), "fraggle_io_initialize_output nf90_def_var L_spin_varid" ) + call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "fraggle_io_initialize_output nf90_def_var PE_varid" ) + call check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "fraggle_io_initialize_output nf90_def_var L_orb_varid" ) + call check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_spin_varid), "fraggle_io_initialize_output nf90_def_var L_spin_varid" ) + end if call check( nf90_inquire(nc%id, nVariables=nvar), "fraggle_io_initialize nf90_inquire nVariables" ) do varid = 1, nvar @@ -202,16 +203,18 @@ module subroutine fraggle_io_write_frame(self, nc, param) call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var rotx_varid" ) end do end do - call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid before" ) - call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid after" ) - call check( nf90_put_var(nc%id, nc%ke_spin_varid, fragments%ke_spin_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_spin_varid before" ) - call check( nf90_put_var(nc%id, nc%ke_spin_varid, fragments%ke_spin_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_spin_varid after" ) - call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid before" ) - call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid after" ) - call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid before" ) - call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid after" ) - call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid before" ) - call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid after" ) + if (param%lenergy) then + call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid after" ) + call check( nf90_put_var(nc%id, nc%ke_spin_varid, fragments%ke_spin_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_spin_varid before" ) + call check( nf90_put_var(nc%id, nc%ke_spin_varid, fragments%ke_spin_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_spin_varid after" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid before" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid after" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid after" ) + call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid before" ) + call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid after" ) + end if call check( nf90_set_fill(nc%id, old_mode, old_mode) ) end associate diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 1a2c2ef8e..507b16745 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -176,8 +176,6 @@ module symba_classes type, extends(symba_encounter) :: symba_plplenc contains procedure :: extract_collisions => symba_collision_extract_collisions_from_encounters !! Processes the pl-pl encounter list remove only those encounters that led to a collision - procedure :: resolve_fragmentations => symba_resolve_collision_fragmentations !! Process list of collisions, determine the collisional regime, and then create fragments - procedure :: resolve_mergers => symba_resolve_collision_mergers !! Process list of collisions and merge colliding bodies together procedure :: resolve_collision => symba_resolve_collision_plplenc !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the c end type symba_plplenc @@ -232,22 +230,6 @@ module subroutine symba_collision_make_colliders_pl(self,idx) integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision end subroutine symba_collision_make_colliders_pl - module subroutine symba_resolve_collision_fragmentations(self, system, param, t) - implicit none - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - end subroutine symba_resolve_collision_fragmentations - - module subroutine symba_resolve_collision_mergers(self, system, param, t) - implicit none - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - end subroutine symba_resolve_collision_mergers - module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 52324c5a3..d9808f5cd 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -890,24 +890,24 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) end subroutine symba_collision_mergeaddsub - module subroutine symba_resolve_collision_fragmentations(self, system, param, t) + subroutine symba_resolve_collision(plplcollision_list , system, param, t) !! author: David A. Minton !! !! Process list of collisions, determine the collisional regime, and then create fragments. !! implicit none ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision + class(symba_plplenc), intent(inout) :: plplcollision_list !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Internals ! Internals integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: i - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, collision_history => param%collision_history) + associate(ncollisions => plplcollision_list%nenc, idx1 => plplcollision_list%index1, idx2 => plplcollision_list%index2, collision_history => param%collision_history) select type(pl => system%pl) class is (symba_pl) select type (cb => system%cb) @@ -920,7 +920,19 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param, t) lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, system%colliders) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle - call system%colliders%regime(system%fragments, system, param) + if (param%lfragmentation) then + call system%colliders%regime(system%fragments, system, param) + else + associate(fragments => system%fragments, colliders => system%colliders) + fragments%regime = COLLRESOLVE_REGIME_MERGE + fragments%mtot = sum(colliders%mass(:)) + fragments%mass_dist(1) = fragments%mtot + fragments%mass_dist(2) = 0.0_DP + fragments%mass_dist(3) = 0.0_DP + fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot + fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot + end associate + end if if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "before") select case (system%fragments%regime) @@ -942,52 +954,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param, t) end associate return - end subroutine symba_resolve_collision_fragmentations - - - module subroutine symba_resolve_collision_mergers(self, system, param, t) - !! author: David A. Minton - !! - !! Process list of collisions and merge colliding bodies together. - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - ! Internals - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical :: lgoodcollision - integer(I4B) :: i - - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, fragments => system%fragments, colliders => system%colliders) - select type(pl => system%pl) - class is (symba_pl) - select type(cb => system%cb) - class is (symba_cb) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) - if (.not. lgoodcollision) cycle - if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - - fragments%regime = COLLRESOLVE_REGIME_MERGE - fragments%mtot = sum(colliders%mass(:)) - fragments%mass_dist(1) = fragments%mtot - fragments%mass_dist(2) = 0.0_DP - fragments%mass_dist(3) = 0.0_DP - fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot - fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot - plplcollision_list%status(i) = symba_collision_casemerge(system, param, t) - end do - end select - end select - end associate - - return - end subroutine symba_resolve_collision_mergers + end subroutine symba_resolve_collision module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec) @@ -1035,11 +1002,8 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & "***********************************************************") allocate(tmp_param, source=param) - if (param%lfragmentation) then - call plplcollision_list%resolve_fragmentations(system, param, t) - else - call plplcollision_list%resolve_mergers(system, param, t) - end if + + call symba_resolve_collision(plplcollision_list, system, param, t) ! Destroy the collision list now that the collisions are resolved call plplcollision_list%setup(0_I8B)