diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c74ea07a0..5aa3c4f8f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,18 +13,23 @@ # Add the source files SET(FAST_MATH_FILES + ${SRC}/modules/swiftest_globals.f90 + ${SRC}/modules/lambda_function.f90 + ${SRC}/modules/swiftest_operators.f90 + ${SRC}/modules/walltime_classes.f90 + ${SRC}/modules/swiftest_classes.f90 ${SRC}/modules/encounter_classes.f90 + ${SRC}/modules/collision_classes.f90 ${SRC}/modules/fraggle_classes.f90 ${SRC}/modules/helio_classes.f90 - ${SRC}/modules/lambda_function.f90 ${SRC}/modules/rmvs_classes.f90 - ${SRC}/modules/swiftest_classes.f90 - ${SRC}/modules/swiftest_globals.f90 - ${SRC}/modules/swiftest_operators.f90 - ${SRC}/modules/swiftest.f90 ${SRC}/modules/symba_classes.f90 - ${SRC}/modules/walltime_classes.f90 ${SRC}/modules/whm_classes.f90 + ${SRC}/modules/swiftest.f90 + ${SRC}/fraggle/collision_io.f90 + ${SRC}/fraggle/collision_placeholder.f90 + ${SRC}/fraggle/collision_setup.f90 + ${SRC}/fraggle/collision_util.f90 ${SRC}/discard/discard.f90 ${SRC}/drift/drift.f90 ${SRC}/encounter/encounter_check.f90 @@ -33,7 +38,6 @@ SET(FAST_MATH_FILES ${SRC}/encounter/encounter_io.f90 ${SRC}/fraggle/fraggle_generate.f90 ${SRC}/fraggle/fraggle_io.f90 - ${SRC}/fraggle/fraggle_placeholder.f90 ${SRC}/fraggle/fraggle_regime.f90 ${SRC}/fraggle/fraggle_set.f90 ${SRC}/fraggle/fraggle_setup.f90 diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 new file mode 100644 index 000000000..88c37cbbc --- /dev/null +++ b/src/collision/collision_io.f90 @@ -0,0 +1,227 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule(collision_classes) s_collision_io + use swiftest + +contains + + module subroutine collision_io_initialize(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. + use, intrinsic :: ieee_arithmetic + use netcdf + implicit none + ! Arguments + class(collision_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param + ! Internals + integer(I4B) :: nvar, varid, vartype + real(DP) :: dfill + real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 + logical :: fileExists + character(len=STRMAX) :: errmsg + integer(I4B) :: ndims + + select type(param) + class is (symba_parameters) + associate(nc => self, collision_history => param%collision_history) + dfill = ieee_value(dfill, IEEE_QUIET_NAN) + sfill = ieee_value(sfill, IEEE_QUIET_NAN) + + select case (param%out_type) + case("NETCDF_FLOAT") + self%out_type = NF90_FLOAT + case("NETCDF_DOUBLE") + self%out_type = NF90_DOUBLE + end select + + ! Check if the file exists, and if it does, delete it + inquire(file=nc%file_name, exist=fileExists) + if (fileExists) then + open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) + close(unit=LUN, status="delete") + end if + + call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "collision_io_initialize nf90_create" ) + + ! Dimensions + call check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "collision_io_initialize nf90_def_dim event_dimid" ) ! Dimension to store individual collision events + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "collision_io_initialize nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers + call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_initialize nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" + + ! Dimension coordinates + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_initialize nf90_def_var space_varid" ) + call check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_initialize nf90_def_var name_varid") + call check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_initialize nf90_def_var stage_varid" ) + + ! Variables + call check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "collision_io_initialize nf90_def_var id_varid" ) + call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & + nc%event_dimid, nc%time_varid), "collision_io_initialize nf90_def_var time_varid" ) + call check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & + [nc%str_dimid, nc%event_dimid], nc%regime_varid), "collision_io_initialize nf90_def_var regime_varid") + call check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & + [ nc%event_dimid], nc%Qloss_varid), "collision_io_initialize nf90_def_var Qloss_varid") + + call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & + [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "collision_io_initialize nf90_def_var ptype_varid") + + call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & + [ nc%event_dimid], nc%loop_varid), "collision_io_initialize nf90_def_var loop_varid") + + call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "collision_io_initialize nf90_def_var rh_varid") + + call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "collision_io_initialize nf90_def_var vh_varid") + + call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& + [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "collision_io_initialize nf90_def_var Gmass_varid") + + + call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& + [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "collision_io_initialize nf90_def_var radius_varid") + + call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "collision_io_initialize nf90_def_var Ip_varid") + + 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), "collision_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), "collision_io_initialize 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), "collision_io_initialize 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), "collision_io_initialize 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), "collision_io_initialize 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), "collision_io_initialize nf90_def_var L_spin_varid" ) + end if + + call check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_initialize nf90_inquire nVariables" ) + do varid = 1, nvar + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "collision_io_initialize nf90_inquire_variable" ) + select case(vartype) + case(NF90_INT) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "collision_io_initialize nf90_def_var_fill NF90_INT" ) + case(NF90_FLOAT) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "collision_io_initialize nf90_def_var_fill NF90_FLOAT" ) + case(NF90_DOUBLE) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "collision_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + case(NF90_CHAR) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "collision_io_initialize nf90_def_var_fill NF90_CHAR" ) + end select + end do + ! Take the file out of define mode + call check( nf90_enddef(nc%id), "collision_io_initialize nf90_enddef" ) + + ! Add in the space and stage dimension coordinates + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "collision_io_initialize nf90_put_var space" ) + call 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_initialize nf90_put_var stage 1" ) + call 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_initialize nf90_put_var stage 2" ) + + end associate + end select + + return + + 667 continue + write(*,*) "Error creating fragmentation output file. " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine collision_io_initialize + + + module subroutine collision_io_write_frame_snapshot(self, nc, param) + !! author: David A. Minton + !! + !! Write a frame of output of a collision result + use netcdf + implicit none + ! Arguments + class(collision_snapshot), intent(in) :: self !! Swiftest encounter structure + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, idslot, old_mode, npl, stage + character(len=:), allocatable :: charstring + class(swiftest_pl), allocatable :: pl + + select type(nc) + class is (collision_io_parameters) + select type (param) + class is (symba_parameters) + associate(impactors => self%impactors, fragments => self%fragments, collision_history => param%collision_history, eslot => param%ioutput) + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "collision_io_write_frame_snapshot nf90_set_fill" ) + + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "collision_io_write_frame_snapshot nf90_put_var time_varid" ) + call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[eslot]), "collision_io_write_frame_snapshot nf90_put_varloop_varid" ) + + charstring = trim(adjustl(REGIME_NAMES(fragments%regime))) + call check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[len(charstring), 1]), "collision_io_write_frame_snapshot nf90_put_var regime_varid" ) + call check( nf90_put_var(nc%id, nc%Qloss_varid, fragments%Qloss, start=[eslot] ), "collision_io_write_frame_snapshot nf90_put_var Qloss_varid" ) + + do stage = 1,2 + if (allocated(pl)) deallocate(pl) + select case(stage) + case(1) + allocate(pl, source=impactors%pl) + case(2) + allocate(pl, source=fragments%pl) + end select + npl = pl%nbody + do i = 1, npl + idslot = findloc(collision_history%idvals,pl%id(i),dim=1) + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), "collision_io_write_frame_snapshot nf90_put_var id_varid" ) + charstring = trim(adjustl(pl%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], count=[len(charstring), 1]), "collision_io_write_frame_snapshot nf90_put_var name_varid" ) + charstring = trim(adjustl(pl%info(i)%particle_type)) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], count=[len(charstring), 1, 1]), "collision_io_write_frame_snapshot nf90_put_var particle_type_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var vh_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "collision_io_write_frame_snapshot nf90_put_var Gmass_varid" ) + call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "collision_io_write_frame_snapshot nf90_put_var radius_varid" ) + call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_write_frame_snapshot nf90_put_var rotx_varid" ) + end do + end do + if (param%lenergy) then + call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_before, start=[ 1, eslot]), "collision_io_write_frame_snapshot 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]), "collision_io_write_frame_snapshot 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]), "collision_io_write_frame_snapshot 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]), "collision_io_write_frame_snapshot nf90_put_var ke_spin_varid after" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_before, start=[ 1, eslot]), "collision_io_write_frame_snapshot nf90_put_var pe_varid before" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_after, start=[ 2, eslot]), "collision_io_write_frame_snapshot 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]), "collision_io_write_frame_snapshot 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]), "collision_io_write_frame_snapshot 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]), "collision_io_write_frame_snapshot 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]), "collision_io_write_frame_snapshot nf90_put_var L_spin_varid after" ) + end if + + call check( nf90_set_fill(nc%id, old_mode, old_mode) ) + end associate + end select + end select + return + end subroutine collision_io_write_frame_snapshot + + +end submodule s_collision_io \ No newline at end of file diff --git a/src/fraggle/fraggle_placeholder.f90 b/src/collision/collision_placeholder.f90 similarity index 74% rename from src/fraggle/fraggle_placeholder.f90 rename to src/collision/collision_placeholder.f90 index 35d5ea960..61ae6eebd 100644 --- a/src/fraggle/fraggle_placeholder.f90 +++ b/src/collision/collision_placeholder.f90 @@ -7,47 +7,47 @@ !! You should have received a copy of the GNU General Public License along with Swiftest. !! If not, see: https://www.gnu.org/licenses. -submodule(fraggle_classes) s_fraggle_placeholder +submodule(collision_classes) s_collision_placeholder use swiftest contains !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class - module subroutine fraggle_placeholder_accel(self, system, param, t, lbeg) + module subroutine collision_placeholder_accel(self, system, param, t, lbeg) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(collision_fragments), intent(inout) :: self !! Fraggle fragment system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step - write(*,*) "The type-bound procedure 'accel' is not defined for type fraggle_fragments" + write(*,*) "The type-bound procedure 'accel' is not defined for the collision_fragments class" return - end subroutine fraggle_placeholder_accel + end subroutine collision_placeholder_accel - module subroutine fraggle_placeholder_kick(self, system, param, t, dt, lbeg) + module subroutine collision_placeholder_kick(self, system, param, t, dt, lbeg) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(collision_fragments), intent(inout) :: self !! Fraggle fragment system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. - write(*,*) "The type-bound procedure 'kick' is not defined for type fraggle_fragments" + write(*,*) "The type-bound procedure 'kick' is not defined for the collision_fragments class" return - end subroutine fraggle_placeholder_kick + end subroutine collision_placeholder_kick - module subroutine fraggle_placeholder_step(self, system, param, t, dt) + module subroutine collision_placeholder_step(self, system, param, t, dt) implicit none - class(fraggle_fragments), intent(inout) :: self !! Swiftest body object + class(collision_fragments), intent(inout) :: self !! Swiftest body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Simulation time real(DP), intent(in) :: dt !! Current stepsize - write(*,*) "The type-bound procedure 'step' is not defined for type fraggle_fragments" + write(*,*) "The type-bound procedure 'step' is not defined for the collision_fragments class" return - end subroutine fraggle_placeholder_step + end subroutine collision_placeholder_step -end submodule s_fraggle_placeholder \ No newline at end of file +end submodule s_collision_placeholder \ No newline at end of file diff --git a/src/collision/collision_setup.f90 b/src/collision/collision_setup.f90 new file mode 100644 index 000000000..24d565472 --- /dev/null +++ b/src/collision/collision_setup.f90 @@ -0,0 +1,73 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule (collision_classes) s_collision_setup + use swiftest +contains + + module subroutine collision_setup_system(self, system, param) + !! author: David A. Minton + !! + !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots + implicit none + ! Arguments + class(collision_system), intent(inout) :: self !! Encounter collision system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + + + return + end subroutine collision_setup_system + + module subroutine collision_setup_fragments(self, n, param) + !! author: David A. Minton + !! + !! Allocates arrays for n fragments in a collision system. Passing n = 0 deallocates all arrays. + implicit none + ! Arguments + class(collision_fragments), intent(inout) :: self + integer(I4B), intent(in) :: n + class(swiftest_parameters), intent(in) :: param + + + call self%swiftest_pl%setup(n, param) + if (n < 0) return + + call self%dealloc() + + if (n == 0) return + + ! allocate(self%rotmag(n)) + ! allocate(self%v_r_mag(n)) + ! allocate(self%v_t_mag(n)) + ! allocate(self%v_n_mag(n)) + + call self%reset() + + return + end subroutine collision_setup_fragments + + + module subroutine collision_setup_impactors(self, system, param) + !! author: David A. Minton + !! + !! Initializes a collider object + implicit none + ! Arguments + class(collision_impactors), intent(inout) :: self !! Fragment system object + class(swiftest_nbody_system), intent(in) :: system + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + + return + end subroutine collision_setup_impactors + +end submodule s_collision_setup + + \ No newline at end of file diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 new file mode 100644 index 000000000..fc0926c91 --- /dev/null +++ b/src/collision/collision_util.f90 @@ -0,0 +1,341 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule (encounter_classes) s_encounter_util + use swiftest +contains + + module subroutine collision_util_final_impactors(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_impactors), intent(inout) :: self !! Collision impactors storage object + + call self%reset() + + return + end subroutine collision_util_final_impactors + + + module subroutine collision_util_final_system(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_system), intent(inout) :: self !! Collision impactors storage object + + call self%reset() + + return + end subroutine collision_util_final_system + + + module subroutine collision_util_final_snapshot(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_snapshot), intent(inout) :: self !! Fraggle encountar storage object + + call encounter_util_final_snapshot(self%encounter_snapshot) + + return + end subroutine collision_util_final_snapshot + + + module subroutine collision_util_final_storage(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_storage(*)), intent(inout) :: self !! Collision storage object + + call util_final_storage(self%swiftest_storage) + + return + end subroutine collision_util_final_storage + + + module subroutine collision_util_final_system(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_system), intent(inout) :: self !! Collision system object + + call self%reset() + + return + end subroutine collision_util_final_storage + + + module subroutine collision_util_get_idvalues_snapshot(self, idvals) + !! author: David A. Minton + !! + !! Returns an array of all id values saved in this snapshot + implicit none + ! Arguments + class(collision_snapshot), intent(in) :: self !! Fraggle snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + ! Internals + integer(I4B) :: ncoll, nfrag + + if (allocated(self%impactors)) then + ncoll = self%impactors%pl%nbody + else + ncoll = 0 + end if + + if (allocated(self%fragments)) then + nfrag = self%fragments%pl%nbody + else + nfrag = 0 + end if + + if (ncoll + nfrag == 0) return + allocate(idvals(ncoll+nfrag)) + + if (ncoll > 0) idvals(1:ncoll) = self%impactors%pl%id(:) + if (nfrag > 0) idvals(ncoll+1:ncoll+nfrag) = self%fragments%pl%id(:) + + return + + end subroutine collision_util_get_idvalues_snapshot + + + module subroutine collision_util_get_energy_momentum(self, system, param, lbefore) + !! Author: David A. Minton + !! + !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) + !! This subrourtine works by building a temporary internal massive body object out of the non-excluded bodies and optionally with fragments appended. + !! This will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it is in the main program. + !! This will temporarily expand the massive body object in a temporary system object called tmpsys to feed it into symba_energy + implicit none + ! Arguments + class(collision_system), intent(inout) :: self !! Encounter collision system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with impactors included and fragments excluded or vice versa + ! Internals + class(swiftest_nbody_system), allocatable, save :: tmpsys + class(swiftest_parameters), allocatable, save :: tmpparam + integer(I4B) :: npl_before, npl_after + + associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody, pl => system%pl, cb => system%cb) + + ! Because we're making a copy of the massive body object with the excludes/fragments appended, we need to deallocate the + ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by + ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this + ! subroutine should be called relatively infrequently, it shouldn't matter too much. + + npl_before = pl%nbody + npl_after = npl_before + nfrag + + if (lbefore) then + call encounter_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) + ! Build the exluded body logical mask for the *before* case: Only the original bodies are used to compute energy and momentum + tmpsys%pl%status(impactors%idx(1:impactors%ncoll)) = ACTIVE + tmpsys%pl%status(npl_before+1:npl_after) = INACTIVE + else + if (.not.allocated(tmpsys)) then + write(*,*) "Error in collision_util_get_energy_momentum. " // & + " This must be called with lbefore=.true. at least once before calling it with lbefore=.false." + call util_exit(FAILURE) + end if + ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum + call encounter_util_add_fragments_to_system(fragments, impactors, tmpsys, tmpparam) + tmpsys%pl%status(impactors%idx(1:impactors%ncoll)) = INACTIVE + tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE + end if + + if (param%lflatten_interactions) call tmpsys%pl%flatten(param) + + call tmpsys%get_energy_and_momentum(param) + + ! Calculate the current fragment energy and momentum balances + if (lbefore) then + fragments%Lorbit_before(:) = tmpsys%Lorbit(:) + fragments%Lspin_before(:) = tmpsys%Lspin(:) + fragments%Ltot_before(:) = tmpsys%Ltot(:) + fragments%ke_orbit_before = tmpsys%ke_orbit + fragments%ke_spin_before = tmpsys%ke_spin + fragments%pe_before = tmpsys%pe + fragments%Etot_before = tmpsys%te + else + fragments%Lorbit_after(:) = tmpsys%Lorbit(:) + fragments%Lspin_after(:) = tmpsys%Lspin(:) + fragments%Ltot_after(:) = tmpsys%Ltot(:) + fragments%ke_orbit_after = tmpsys%ke_orbit + fragments%ke_spin_after = tmpsys%ke_spin + fragments%pe_after = tmpsys%pe + fragments%Etot_after = tmpsys%te - (fragments%pe_after - fragments%pe_before) ! Gotta be careful with PE when number of bodies changes. + end if + end associate + + return + end subroutine collision_util_get_energy_momentum + + + module subroutine collision_util_index_map(self) + !! author: David A. Minton + !! + !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id + implicit none + ! Arguments + class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + ! Internals + integer(I4B), dimension(:), allocatable :: idvals + real(DP), dimension(:), allocatable :: tvals + + call encounter_util_get_vals_storage(self, idvals, tvals) + + ! Consolidate ids to only unique values + call util_unique(idvals,self%idvals,self%idmap) + self%nid = size(self%idvals) + + ! Don't consolidate time values (multiple collisions can happen in a single time step) + self%nt = size(self%tvals) + + return + end subroutine collision_util_index_map + + + module subroutine collision_util_reset_fragments(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables + implicit none + ! Arguments + class(collision_fragments), intent(inout) :: self + + call self%swiftest_pl%dealloc() + + if (allocated(self%rc)) deallocate(self%rc) + if (allocated(self%vc)) deallocate(self%vc) + if (allocated(self%rmag)) deallocate(self%rmag) + if (allocated(self%rotmag)) deallocate(self%rotmag) + if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) + if (allocated(self%v_t_unit)) deallocate(self%v_t_unit) + if (allocated(self%v_n_unit)) deallocate(self%v_n_unit) + + return + end subroutine collision_util_reset_fragments + + + module subroutine collision_util_reset_impactors(self) + !! author: David A. Minton + !! + !! Resets the collider object variables to 0 and deallocates the index and mass distributions + implicit none + ! Arguments + class(collision_impactors), intent(inout) :: self + + if (allocated(self%idx)) deallocate(self%idx) + if (allocated(self%mass_dist)) deallocate(self%mass_dist) + ncoll = 0 + rb(:,:) = 0.0_DP + vb(:,:) = 0.0_DP + rot(:,:) = 0.0_DP + L_spin(:,:) = 0.0_DP + L_orbit(:,:) = 0.0_DP + Ip(:,:) = 0.0_DP + mass(:) = 0.0_DP + radius(:) = 0.0_DP + Qloss = 0.0_DP + regime = 0 + + x_unit(:) = 0.0_DP + y_unit(:) = 0.0_DP + z_unit(:) = 0.0_DP + v_unit(:) = 0.0_DP + rbcom(:) = 0.0_DP + vbcom(:) = 0.0_DP + rbimp(:) = 0.0_DP + + return + end subroutine collision_util_reset_impactors + + + module subroutine collision_util_reset_system(self) + !! author: David A. Minton + !! + !! Resets the collider system and deallocates all allocatables + implicit none + ! Arguments + class(collision_system), intent(inout) :: self + + if (allocated(self%impactors)) deallocate(self%impactors) + if (allocated(self%fragments)) deallocate(self%fragments) + if (allocated(self%before)) deallocate(self%before) + if (allocated(self%after)) deallocate(self%after) + + Lorbit(:,:) = 0.0_DP + Lspin(:,:) = 0.0_DP + Ltot(:,:) = 0.0_DP + ke_orbit(:) = 0.0_DP + ke_spin(:) = 0.0_DP + pe(:) = 0.0_DP + Etot(:) = 0.0_DP + + return + end subroutine collision_util_reset_impactors + + + subroutine collision_util_save_snapshot(collision_history, snapshot) + !! author: David A. Minton + !! + !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. + !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an + !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment + !! Memory usage grows by a factor of 2 each time it fills up, but no more. + implicit none + ! Arguments + type(collision_storage(*)), allocatable, intent(inout) :: collision_history !! Collision history object + class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object + ! Internals + type(collision_storage(nframes=:)), allocatable :: tmp + integer(I4B) :: i, nnew, nold, nbig + + ! Advance the snapshot frame counter + collision_history%iframe = collision_history%iframe + 1 + + ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 + nnew = collision_history%iframe + nold = collision_history%nframes + + if (nnew > nold) then + nbig = nold + do while (nbig < nnew) + nbig = nbig * 2 + end do + allocate(collision_storage(nbig) :: tmp) + tmp%iframe = collision_history%iframe + call move_alloc(collision_history%nc, tmp%nc) + + do i = 1, nold + if (allocated(collision_history%frame(i)%item)) call move_alloc(collision_history%frame(i)%item, tmp%frame(i)%item) + end do + deallocate(collision_history) + call move_alloc(tmp,collision_history) + nnew = nbig + end if + + collision_history%frame(nnew) = snapshot + + return + end subroutine collision_util_save_snapshot + + +end submodule s_encounter_util \ No newline at end of file diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 60748c5b0..fb84067ca 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -12,7 +12,7 @@ contains - module subroutine encounter_io_dump_collision(self, param) + module subroutine collision_io_dump(self, param) !! author: David A. Minton !! !! Dumps the time history of an encounter to file. @@ -24,7 +24,7 @@ module subroutine encounter_io_dump_collision(self, param) integer(I4B) :: i select type(nc => self%nc) - class is (fraggle_io_parameters) + class is (collision_io_parameters) if (self%iframe > 0) then nc%file_number = nc%file_number + 1 call self%make_index_map() @@ -37,7 +37,7 @@ module subroutine encounter_io_dump_collision(self, param) do i = 1, self%nframes if (allocated(self%frame(i)%item)) then select type(snapshot => self%frame(i)%item) - class is (fraggle_snapshot) + class is (collision_snapshot) param%ioutput = i call snapshot%write_frame(nc,param) end select @@ -52,10 +52,10 @@ module subroutine encounter_io_dump_collision(self, param) end select return - end subroutine encounter_io_dump_collision + end subroutine collision_io_dump - module subroutine encounter_io_dump_encounter(self, param) + module subroutine encounter_io_dump(self, param) ! author: David A. Minton !! !! Dumps the time history of an encounter to file. @@ -95,7 +95,7 @@ module subroutine encounter_io_dump_encounter(self, param) end select return - end subroutine encounter_io_dump_encounter + end subroutine encounter_io_dump module subroutine encounter_io_initialize(self, param) @@ -194,7 +194,7 @@ module subroutine encounter_io_initialize(self, param) end subroutine encounter_io_initialize - module subroutine encounter_io_write_frame(self, nc, param) + module subroutine encounter_io_write_frame_snapshot(self, nc, param) !! author: David A. Minton !! !! Write a frame of output of an encounter trajectory. @@ -213,43 +213,43 @@ module subroutine encounter_io_write_frame(self, nc, param) select type (param) class is (symba_parameters) associate(pl => self%pl, tp => self%tp, encounter_history => param%encounter_history, tslot => param%ioutput) - call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill" ) + call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame_snapshot nf90_set_fill" ) - call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) - call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame nf90_put_var pl loop_varid" ) + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame_snapshot nf90_put_var time_varid" ) + call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl loop_varid" ) npl = pl%nbody do i = 1, npl idslot = findloc(encounter_history%idvals,pl%id(i),dim=1) - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl vh_varid" ) - call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame_snapshot nf90_put_var pl id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl vh_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl Gmass_varid" ) - if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid" ) + if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame_snapshot nf90_put_var pl radius_varid" ) if (param%lrotation) then - call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl Ip_varid" ) - call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rotx_varid" ) + call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var pl rotx_varid" ) end if charstring = trim(adjustl(pl%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var pl name_varid" ) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var pl name_varid" ) charstring = trim(adjustl(pl%info(i)%particle_type)) - call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var pl particle_type_varid" ) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var pl particle_type_varid" ) end do ntp = tp%nbody do i = 1, ntp idslot = findloc(param%encounter_history%idvals,tp%id(i),dim=1) - call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var tp id_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var tp vh_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_write_frame_snapshot nf90_put_var tp id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var tp rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame_snapshot nf90_put_var tp vh_varid" ) charstring = trim(adjustl(tp%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var tp name_varid" ) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var tp name_varid" ) charstring = trim(adjustl(tp%info(i)%particle_type)) - call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame nf90_put_var tp particle_type_varid" ) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[len(charstring), 1]), "encounter_io_write_frame_snapshot nf90_put_var tp particle_type_varid" ) end do call check( nf90_set_fill(nc%id, old_mode, old_mode) ) @@ -258,7 +258,7 @@ module subroutine encounter_io_write_frame(self, nc, param) end select return - end subroutine encounter_io_write_frame + end subroutine encounter_io_write_frame_snapshot diff --git a/src/encounter/encounter_placeholder.f90 b/src/encounter/encounter_placeholder.f90 new file mode 100644 index 000000000..61ae6eebd --- /dev/null +++ b/src/encounter/encounter_placeholder.f90 @@ -0,0 +1,53 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule(collision_classes) s_collision_placeholder + use swiftest + +contains + + !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class + module subroutine collision_placeholder_accel(self, system, param, t, lbeg) + implicit none + class(collision_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + write(*,*) "The type-bound procedure 'accel' is not defined for the collision_fragments class" + return + end subroutine collision_placeholder_accel + + module subroutine collision_placeholder_kick(self, system, param, t, dt, lbeg) + implicit none + class(collision_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + + write(*,*) "The type-bound procedure 'kick' is not defined for the collision_fragments class" + return + end subroutine collision_placeholder_kick + + module subroutine collision_placeholder_step(self, system, param, t, dt) + implicit none + class(collision_fragments), intent(inout) :: self !! Swiftest body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize + + write(*,*) "The type-bound procedure 'step' is not defined for the collision_fragments class" + return + end subroutine collision_placeholder_step + + +end submodule s_collision_placeholder \ No newline at end of file diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 9d9c97f11..c4acc1df9 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -156,20 +156,6 @@ module subroutine encounter_util_final_snapshot(self) end subroutine encounter_util_final_snapshot - module subroutine encounter_util_final_collision_storage(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_storage(*)), intent(inout) :: self !! SyMBA nbody system object - - call util_final_storage(self%swiftest_storage) - - return - end subroutine encounter_util_final_collision_storage - - module subroutine encounter_util_final_storage(self) !! author: David A. Minton !! @@ -275,7 +261,7 @@ subroutine encounter_util_get_vals_storage(storage, idvals, tvals) end subroutine encounter_util_get_vals_storage - module subroutine encounter_util_index_map_encounter(self) + module subroutine encounter_util_index_map(self) !! author: David A. Minton !! !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id. @@ -298,32 +284,7 @@ module subroutine encounter_util_index_map_encounter(self) self%nt = size(self%tvals) return - end subroutine encounter_util_index_map_encounter - - - - module subroutine encounter_util_index_map_collision(self) - !! author: David A. Minton - !! - !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - implicit none - ! Arguments - class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object - ! Internals - integer(I4B), dimension(:), allocatable :: idvals - real(DP), dimension(:), allocatable :: tvals - - call encounter_util_get_vals_storage(self, idvals, tvals) - - ! Consolidate ids to only unique values - call util_unique(idvals,self%idvals,self%idmap) - self%nid = size(self%idvals) - - ! Don't consolidate time values (multiple collisions can happen in a single time step) - self%nt = size(self%tvals) - - return - end subroutine encounter_util_index_map_collision + end subroutine encounter_util_index_map module subroutine encounter_util_resize_list(self, nnew) @@ -402,53 +363,7 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru end subroutine encounter_util_spill_list - - subroutine encounter_util_save_collision(collision_history, snapshot) - !! author: David A. Minton - !! - !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. - !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an - !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment - !! Memory usage grows by a factor of 2 each time it fills up, but no more. - implicit none - ! Arguments - type(collision_storage(*)), allocatable, intent(inout) :: collision_history !! Collision history object - class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object - ! Internals - type(collision_storage(nframes=:)), allocatable :: tmp - integer(I4B) :: i, nnew, nold, nbig - - ! Advance the snapshot frame counter - collision_history%iframe = collision_history%iframe + 1 - - ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 - nnew = collision_history%iframe - nold = collision_history%nframes - - if (nnew > nold) then - nbig = nold - do while (nbig < nnew) - nbig = nbig * 2 - end do - allocate(collision_storage(nbig) :: tmp) - tmp%iframe = collision_history%iframe - call move_alloc(collision_history%nc, tmp%nc) - - do i = 1, nold - if (allocated(collision_history%frame(i)%item)) call move_alloc(collision_history%frame(i)%item, tmp%frame(i)%item) - end do - deallocate(collision_history) - call move_alloc(tmp,collision_history) - nnew = nbig - end if - - collision_history%frame(nnew) = snapshot - - return - end subroutine encounter_util_save_collision - - - subroutine encounter_util_save_encounter(encounter_history, snapshot) + subroutine encounter_util_save_snapshot(encounter_history, snapshot) !! author: David A. Minton !! !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. @@ -492,261 +407,7 @@ subroutine encounter_util_save_encounter(encounter_history, snapshot) encounter_history%frame(nnew) = snapshot return - end subroutine encounter_util_save_encounter - - - module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) - !! author: David A. Minton - !! - !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories - !! can be played back through the encounter - implicit none - ! Internals - class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store - real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. - ! Arguments - class(fraggle_snapshot), allocatable :: snapshot - type(symba_pl) :: pl - character(len=:), allocatable :: stage - - if (present(arg)) then - stage = arg - else - stage = "" - end if - - select type (system) - class is (symba_nbody_system) - - select case(stage) - case("before") - ! Saves the states of the bodies involved in the collision before the collision is resolved - associate (idx => system%colliders%idx, ncoll => system%colliders%ncoll) - call pl%setup(ncoll, param) - pl%id(:) = system%pl%id(idx(:)) - pl%Gmass(:) = system%pl%Gmass(idx(:)) - pl%radius(:) = system%pl%radius(idx(:)) - pl%rot(:,:) = system%pl%rot(:,idx(:)) - pl%Ip(:,:) = system%pl%Ip(:,idx(:)) - pl%rh(:,:) = system%pl%rh(:,idx(:)) - pl%vh(:,:) = system%pl%vh(:,idx(:)) - pl%info(:) = system%pl%info(idx(:)) - !end select - allocate(system%colliders%pl, source=pl) - end associate - case("after") - allocate(fraggle_snapshot :: snapshot) - allocate(snapshot%colliders, source=system%colliders) - allocate(snapshot%fragments, source=system%fragments) - snapshot%t = t - select type(param) - class is (symba_parameters) - call encounter_util_save_collision(param%collision_history,snapshot) - end select - case default - write(*,*) "encounter_util_snapshot_collision requies either 'before' or 'after' passed to 'arg'" - end select - - end select - - return - end subroutine encounter_util_snapshot_collision - - - module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) - !! author: David A. Minton - !! - !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories - !! can be played back through the encounter - implicit none - ! Internals - class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store - real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) - ! Arguments - class(encounter_snapshot), allocatable :: snapshot - integer(I4B) :: i, pi, pj, k, npl_snap, ntp_snap, iflag - real(DP), dimension(NDIM) :: rrel, vrel, rcom, vcom - real(DP) :: Gmtot, a, q, capm, tperi - real(DP), dimension(NDIM,2) :: rb,vb - - if (.not.present(t)) then - write(*,*) "encounter_util_snapshot_encounter requires `t` to be passed" - return - end if + end subroutine encounter_util_save_snapshot - if (.not.present(arg)) then - write(*,*) "encounter_util_snapshot_encounter requires `arg` to be passed" - return - end if - - select type(param) - class is (symba_parameters) - select type (system) - class is (symba_nbody_system) - select type(pl => system%pl) - class is (symba_pl) - select type (tp => system%tp) - class is (symba_tp) - associate(npl => pl%nbody, ntp => tp%nbody) - if (npl + ntp == 0) return - allocate(encounter_snapshot :: snapshot) - allocate(snapshot%pl, mold=pl) - allocate(snapshot%tp, mold=tp) - snapshot%iloop = param%iloop - - select type(pl_snap => snapshot%pl) - class is (symba_pl) - select type(tp_snap => snapshot%tp) - class is (symba_tp) - - select case(arg) - case("trajectory") - snapshot%t = t - - npl_snap = npl - ntp_snap = ntp - - if (npl > 0) then - pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == system%irec - npl_snap = count(pl%lmask(1:npl)) - end if - if (ntp > 0) then - tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == system%irec - ntp_snap = count(tp%lmask(1:ntp)) - end if - - if (npl_snap + ntp_snap == 0) return ! Nothing to snapshot - - pl_snap%nbody = npl_snap - - ! Take snapshot of the currently encountering massive bodies - if (npl_snap > 0) then - call pl_snap%setup(npl_snap, param) - pl_snap%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) - pl_snap%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) - pl_snap%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) - pl_snap%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) - do i = 1, NDIM - pl_snap%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) - pl_snap%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) - end do - if (param%lclose) then - pl_snap%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) - end if - - if (param%lrotation) then - do i = 1, NDIM - pl_snap%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) - pl_snap%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) - end do - end if - call pl_snap%sort("id", ascending=.true.) - end if - - ! Take snapshot of the currently encountering test particles - tp_snap%nbody = ntp_snap - if (ntp_snap > 0) then - call tp_snap%setup(ntp_snap, param) - tp_snap%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) - tp_snap%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) - do i = 1, NDIM - tp_snap%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) - tp_snap%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) - end do - end if - - ! Save the snapshot - param%encounter_history%nid = param%encounter_history%nid + ntp_snap + npl_snap - call encounter_util_save_encounter(param%encounter_history,snapshot) - case("closest") - associate(plplenc_list => system%plplenc_list, pltpenc_list => system%pltpenc_list) - if (any(plplenc_list%lclosest(:))) then - call pl_snap%setup(2, param) - do k = 1, plplenc_list%nenc - if (plplenc_list%lclosest(k)) then - pi = plplenc_list%index1(k) - pj = plplenc_list%index2(k) - pl_snap%levelg(:) = pl%levelg([pi,pj]) - pl_snap%id(:) = pl%id([pi,pj]) - pl_snap%info(:) = pl%info([pi,pj]) - pl_snap%Gmass(:) = pl%Gmass([pi,pj]) - Gmtot = sum(pl_snap%Gmass(:)) - if (param%lclose) pl_snap%radius(:) = pl%radius([pi,pj]) - if (param%lrotation) then - do i = 1, NDIM - pl_snap%Ip(i,:) = pl%Ip(i,[pi,pj]) - pl_snap%rot(i,:) = pl%rot(i,[pi,pj]) - end do - end if - - ! Compute pericenter passage time to get the closest approach parameters - rrel(:) = plplenc_list%r2(:,k) - plplenc_list%r1(:,k) - vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k) - call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) - snapshot%t = t + tperi - if ((snapshot%t < maxval(pl_snap%info(:)%origin_time)) .or. & - (snapshot%t > minval(pl_snap%info(:)%discard_time))) cycle - - ! Computer the center mass of the pair - rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot - vcom(:) = (plplenc_list%v1(:,k) * pl_snap%Gmass(1) + plplenc_list%v2(:,k) * pl_snap%Gmass(2)) / Gmtot - rb(:,1) = plplenc_list%r1(:,k) - rcom(:) - rb(:,2) = plplenc_list%r2(:,k) - rcom(:) - vb(:,1) = plplenc_list%v1(:,k) - vcom(:) - vb(:,2) = plplenc_list%v2(:,k) - vcom(:) - - ! Drift the relative orbit to get the new relative position and velocity - call drift_one(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), tperi, iflag) - if (iflag /= 0) write(*,*) "Danby error in encounter_util_snapshot_encounter. Closest approach positions and vectors may not be accurate." - - ! Get the new position and velocity vectors - rb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * rrel(:) - rb(:,2) = (pl_snap%Gmass(1)) / Gmtot * rrel(:) - - vb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * vrel(:) - vb(:,2) = (pl_snap%Gmass(1)) / Gmtot * vrel(:) - - ! Move the CoM assuming constant velocity over the time it takes to reach periapsis - rcom(:) = rcom(:) + vcom(:) * tperi - - ! Compute the heliocentric position and velocity vector at periapsis - pl_snap%rh(:,1) = rb(:,1) + rcom(:) - pl_snap%rh(:,2) = rb(:,2) + rcom(:) - pl_snap%vh(:,1) = vb(:,1) + vcom(:) - pl_snap%vh(:,2) = vb(:,2) + vcom(:) - - call pl_snap%sort("id", ascending=.true.) - call encounter_util_save_encounter(param%encounter_history,snapshot) - end if - end do - - plplenc_list%lclosest(:) = .false. - end if - - if (any(pltpenc_list%lclosest(:))) then - do k = 1, pltpenc_list%nenc - end do - pltpenc_list%lclosest(:) = .false. - end if - end associate - case default - write(*,*) "encounter_util_snapshot_encounter requires `arg` to be either `trajectory` or `closest`" - end select - end select - end select - end associate - end select - end select - end select - end select - - return - end subroutine encounter_util_snapshot_encounter end submodule s_encounter_util \ No newline at end of file diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index ac5f166ba..4ad5c75d1 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -17,7 +17,7 @@ contains - module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure) + module subroutine collision_generate_fragments(self, impactors, system, param, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. @@ -25,7 +25,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa implicit none ! Arguments class(fraggle_fragments), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies + class(collision_impactors), intent(inout) :: impactors !! Fraggle impactors object containing the two-body equivalent values of the colliding bodies class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? @@ -68,22 +68,22 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lk_plpl = .false. end if - call fragments%set_natural_scale(colliders) + call fragments%set_natural_scale(impactors) call fragments%reset() ! Calculate the initial energy of the system without the collisional family - call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.true.) + call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.true.) ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution - r_max_start = 1.2_DP * (norm2(colliders%rb(:,2) - colliders%rb(:,1))) + r_max_start = 1.2_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) lfailure = .false. try = 1 do while (try < MAXTRY) write(message,*) try call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle try " // trim(adjustl(message))) if (lfailure) then - call fragments%restructure(colliders, try, f_spin, r_max_start) + call fragments%restructure(impactors, try, f_spin, r_max_start) call fragments%reset() try = try + 1 end if @@ -91,36 +91,36 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = .false. call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - call fraggle_generate_pos_vec(fragments, colliders, r_max_start) - call fragments%set_coordinate_system(colliders) + call fraggle_generate_pos_vec(fragments, impactors, r_max_start) + call fragments%set_coordinate_system(impactors) ! Initial velocity guess will be the barycentric velocity of the colliding system so that the budgets are based on the much smaller collisional-frame velocities do concurrent (i = 1:nfrag) fragments%vb(:, i) = fragments%vbcom(:) end do - call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.false.) call fragments%set_budgets() - call fraggle_generate_spins(fragments, colliders, f_spin, lfailure) + call fraggle_generate_spins(fragments, impactors, f_spin, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find spins") cycle end if - call fraggle_generate_tan_vel(fragments, colliders, lfailure) + call fraggle_generate_tan_vel(fragments, impactors, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find tangential velocities") cycle end if - call fraggle_generate_rad_vel(fragments, colliders, lfailure) + call fraggle_generate_rad_vel(fragments, impactors, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find radial velocities") cycle end if - call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.false.) dEtot = fragments%Etot_after - fragments%Etot_before dLmag = .mag. (fragments%Ltot_after(:) - fragments%Ltot_before(:)) exit @@ -158,7 +158,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa trim(adjustl(message)) // " tries") end if - call fragments%set_original_scale(colliders) + call fragments%set_original_scale(impactors) ! Restore the big array if (lk_plpl) call pl%flatten(param) @@ -166,10 +166,10 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily return - end subroutine fraggle_generate_fragments + end subroutine collision_generate_fragments - subroutine fraggle_generate_pos_vec(fragments, colliders, r_max_start) + subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the @@ -178,7 +178,7 @@ subroutine fraggle_generate_pos_vec(fragments, colliders, r_max_start) implicit none ! Arguments class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object real(DP), intent(in) :: r_max_start !! Initial guess for the starting maximum radial distance of fragments ! Internals real(DP) :: dis, rad, r_max, fdistort @@ -196,14 +196,14 @@ subroutine fraggle_generate_pos_vec(fragments, colliders, r_max_start) ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) r_max = r_max_start - rad = sum(colliders%radius(:)) + rad = sum(impactors%radius(:)) ! Get the unit vectors for the relative position and velocity vectors. These are used to shift the fragment cloud depending on the - runit(:) = colliders%rb(:,2) - colliders%rb(:,1) + runit(:) = impactors%rb(:,2) - impactors%rb(:,1) runit(:) = runit(:) / (.mag. runit(:)) - vunit(:) = colliders%vb(:,2) - colliders%vb(:,1) + vunit(:) = impactors%vb(:,2) - impactors%vb(:,1) vunit(:) = vunit(:) / (.mag. vunit(:)) ! This is a factor that will "distort" the shape of the frgment cloud in the direction of the impact velocity @@ -214,8 +214,8 @@ subroutine fraggle_generate_pos_vec(fragments, colliders, r_max_start) call random_number(fragments%r_coll(:,3:nfrag)) loverlap(:) = .true. do while (any(loverlap(3:nfrag))) - fragments%r_coll(:, 1) = colliders%rb(:, 1) - fragments%rbcom(:) - fragments%r_coll(:, 2) = colliders%rb(:, 2) - fragments%rbcom(:) + fragments%r_coll(:, 1) = impactors%rb(:, 1) - fragments%rbcom(:) + fragments%r_coll(:, 2) = impactors%rb(:, 2) - fragments%rbcom(:) r_max = r_max + 0.1_DP * rad do i = 3, nfrag if (loverlap(i)) then @@ -230,13 +230,13 @@ subroutine fraggle_generate_pos_vec(fragments, colliders, r_max_start) loverlap(:) = .false. do j = 1, nfrag do i = j + 1, nfrag - dis = norm2(fragments%r_coll(:,j) - fragments%r_coll(:,i)) + dis = .mag.(fragments%r_coll(:,j) - fragments%r_coll(:,i)) loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) end do end do end do call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%r_coll) - call fragments%set_coordinate_system(colliders) + call fragments%set_coordinate_system(impactors) do concurrent(i = 1:nfrag) fragments%rb(:,i) = fragments%r_coll(:,i) + fragments%rbcom(:) @@ -253,7 +253,7 @@ subroutine fraggle_generate_pos_vec(fragments, colliders, r_max_start) end subroutine fraggle_generate_pos_vec - subroutine fraggle_generate_spins(fragments, colliders, f_spin, lfailure) + subroutine fraggle_generate_spins(fragments, impactors, f_spin, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Calculates the spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. @@ -262,7 +262,7 @@ subroutine fraggle_generate_spins(fragments, colliders, f_spin, lfailure) implicit none ! Arguments class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object real(DP), intent(in) :: f_spin !! Fraction of energy or momentum that goes into spin (whichever gives the lowest kinetic energy) logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! ! Internals @@ -279,12 +279,12 @@ subroutine fraggle_generate_spins(fragments, colliders, f_spin, lfailure) ke_remainder = fragments%ke_budget ! Add a fraction of the orbit angular momentum of the second body to the spin of the first body to kick things off - L(:) = colliders%L_spin(:,1) + f_spin * (colliders%L_orbit(:,2) + colliders%L_spin(:,2)) + L(:) = impactors%L_spin(:,1) + f_spin * (impactors%L_orbit(:,2) + impactors%L_spin(:,2)) fragments%rot(:,1) = L(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) L_remainder(:) = L_remainder(:) - L(:) ! Partition the spin momentum of the second body into the spin of the fragments, with some random variation - L(:) = colliders%L_spin(:,2) / (nfrag - 1) + L(:) = impactors%L_spin(:,2) / (nfrag - 1) call random_number(frot_rand(:,2:nfrag)) frot_rand(:,2:nfrag) = 2 * (frot_rand(:,2:nfrag) - 0.5_DP) * frot_rand_mag @@ -297,17 +297,17 @@ subroutine fraggle_generate_spins(fragments, colliders, f_spin, lfailure) ! Make sure we didn't blow our kinetic energy budget do i = 1, nfrag - ke_remainder = ke_remainder - 0.5_DP * fragments%mass(i) * fragments%Ip(3, i) * fragments%radius(i)**2 * norm2(fragments%rot(:, i)) + ke_remainder = ke_remainder - 0.5_DP * fragments%mass(i) * fragments%Ip(3, i) * fragments%radius(i)**2 * .mag.(fragments%rot(:, i)) end do ! Distributed most of the remaining angular momentum amongst all the particles fragments%ke_spin = 0.0_DP - if (norm2(L_remainder(:)) > FRAGGLE_LTOL) then + if (.mag.(L_remainder(:)) > FRAGGLE_LTOL) then do i = nfrag, 1, -1 ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets - rot_ke(:) = sqrt(2 * f_spin * ke_remainder / (i * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i))) * L_remainder(:) / norm2(L_remainder(:)) + rot_ke(:) = sqrt(2 * f_spin * ke_remainder / (i * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i))) * L_remainder(:) / .mag.(L_remainder(:)) rot_L(:) = f_spin * L_remainder(:) / (i * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3, i)) - if (norm2(rot_ke) < norm2(rot_L)) then + if (.mag.(rot_ke) < .mag.(rot_L)) then fragments%rot(:,i) = fragments%rot(:,i) + rot_ke(:) else fragments%rot(:, i) = fragments%rot(:,i) + rot_L(:) @@ -341,7 +341,7 @@ subroutine fraggle_generate_spins(fragments, colliders, f_spin, lfailure) end subroutine fraggle_generate_spins - subroutine fraggle_generate_tan_vel(fragments, colliders, lfailure) + subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Adjusts the tangential velocities and spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. @@ -357,7 +357,7 @@ subroutine fraggle_generate_tan_vel(fragments, colliders, lfailure) implicit none ! Arguments class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds ! Internals integer(I4B) :: i, try @@ -379,9 +379,9 @@ subroutine fraggle_generate_tan_vel(fragments, colliders, lfailure) allocate(v_t_initial, mold=fragments%v_t_mag) do try = 1, MAXTRY - v_t_initial(1) = dot_product(colliders%vb(:,1),fragments%v_t_unit(:,1)) + v_t_initial(1) = dot_product(impactors%vb(:,1),fragments%v_t_unit(:,1)) do i = 2, nfrag - v_t_initial(i) = dot_product(colliders%vb(:,2), fragments%v_t_unit(:,i)) + v_t_initial(i) = dot_product(impactors%vb(:,2), fragments%v_t_unit(:,i)) end do fragments%v_t_mag(:) = v_t_initial @@ -525,7 +525,7 @@ end function tangential_objective_function end subroutine fraggle_generate_tan_vel - subroutine fraggle_generate_rad_vel(fragments, colliders, lfailure) + subroutine fraggle_generate_rad_vel(fragments, impactors, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! @@ -533,7 +533,7 @@ subroutine fraggle_generate_rad_vel(fragments, colliders, lfailure) implicit none ! Arguments class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! ! Internals real(DP), parameter :: TOL_MIN = FRAGGLE_ETOL ! This needs to be more accurate than the tangential step, as we are trying to minimize the total residual energy @@ -552,7 +552,7 @@ subroutine fraggle_generate_rad_vel(fragments, colliders, lfailure) allocate(v_r_initial, source=fragments%v_r_mag) ! Initialize radial velocity magnitudes with a random value that related to equipartition of kinetic energy with some noise and scaled with respect to the initial distance - v_r_initial(1) = dot_product(colliders%vb(:,1),fragments%v_r_unit(:,1)) + v_r_initial(1) = dot_product(impactors%vb(:,1),fragments%v_r_unit(:,1)) fragments%ke_orbit = 0.5_DP * fragments%mass(1) * v_r_initial(1)**2 ke_radial = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 3253bcdb2..1ddc4b4e9 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -12,225 +12,13 @@ contains - module subroutine fraggle_io_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. - use, intrinsic :: ieee_arithmetic - use netcdf - implicit none - ! Arguments - class(fraggle_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param - ! Internals - integer(I4B) :: nvar, varid, vartype - real(DP) :: dfill - real(SP) :: sfill - integer(I4B), parameter :: NO_FILL = 0 - logical :: fileExists - character(len=STRMAX) :: errmsg - integer(I4B) :: ndims - - select type(param) - class is (symba_parameters) - associate(nc => self, collision_history => param%collision_history) - dfill = ieee_value(dfill, IEEE_QUIET_NAN) - sfill = ieee_value(sfill, IEEE_QUIET_NAN) - - select case (param%out_type) - case("NETCDF_FLOAT") - self%out_type = NF90_FLOAT - case("NETCDF_DOUBLE") - self%out_type = NF90_DOUBLE - end select - - ! Check if the file exists, and if it does, delete it - inquire(file=nc%file_name, exist=fileExists) - if (fileExists) then - open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) - close(unit=LUN, status="delete") - end if - - call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "fraggle_io_initialize nf90_create" ) - - ! Dimensions - call check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "fraggle_io_initialize nf90_def_dim event_dimid" ) ! Dimension to store individual collision events - call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "fraggle_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension - call check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "fraggle_io_initialize nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers - call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "fraggle_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - call check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "fraggle_io_initialize nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" - - ! Dimension coordinates - call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "fraggle_io_initialize nf90_def_var space_varid" ) - call check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "fraggle_io_initialize nf90_def_var name_varid") - call check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "fraggle_io_initialize nf90_def_var stage_varid" ) - - ! Variables - call check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "fraggle_io_initialize nf90_def_var id_varid" ) - call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & - nc%event_dimid, nc%time_varid), "fraggle_io_initialize nf90_def_var time_varid" ) - call check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & - [nc%str_dimid, nc%event_dimid], nc%regime_varid), "fraggle_io_initialize nf90_def_var regime_varid") - call check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & - [ nc%event_dimid], nc%Qloss_varid), "fraggle_io_initialize nf90_def_var Qloss_varid") - - call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & - [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "fraggle_io_initialize nf90_def_var ptype_varid") - - call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & - [ nc%event_dimid], nc%loop_varid), "fraggle_io_initialize nf90_def_var loop_varid") - - call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "fraggle_io_initialize nf90_def_var rh_varid") - - call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "fraggle_io_initialize nf90_def_var vh_varid") - - call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "fraggle_io_initialize nf90_def_var Gmass_varid") - - - call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "fraggle_io_initialize nf90_def_var radius_varid") - - call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "fraggle_io_initialize nf90_def_var Ip_varid") - - 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%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 - call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "fraggle_io_initialize nf90_inquire_variable" ) - select case(vartype) - case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "fraggle_io_initialize nf90_def_var_fill NF90_INT" ) - case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "fraggle_io_initialize nf90_def_var_fill NF90_FLOAT" ) - case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "fraggle_io_initialize nf90_def_var_fill NF90_DOUBLE" ) - case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "fraggle_io_initialize nf90_def_var_fill NF90_CHAR" ) - end select - end do - ! Take the file out of define mode - call check( nf90_enddef(nc%id), "fraggle_io_initialize nf90_enddef" ) - - ! Add in the space and stage dimension coordinates - call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "fraggle_io_initialize nf90_put_var space" ) - call check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(1), start=[1,1], count=[len(nc%stage_coords(1)),1]), "fraggle_io_initialize nf90_put_var stage 1" ) - call check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(2), start=[1,2], count=[len(nc%stage_coords(2)),1]), "fraggle_io_initialize nf90_put_var stage 2" ) - - end associate - end select - - return - - 667 continue - write(*,*) "Error creating fragmentation output file. " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine fraggle_io_initialize_output - - - module subroutine fraggle_io_write_frame(self, nc, param) - !! author: David A. Minton - !! - !! Write a frame of output of a collision result - use netcdf - implicit none - ! Arguments - class(fraggle_snapshot), intent(in) :: self !! Swiftest encounter structure - class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i, idslot, old_mode, npl, stage - character(len=:), allocatable :: charstring - class(swiftest_pl), allocatable :: pl - - select type(nc) - class is (fraggle_io_parameters) - select type (param) - class is (symba_parameters) - associate(colliders => self%colliders, fragments => self%fragments, collision_history => param%collision_history, eslot => param%ioutput) - call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "fraggle_io_write_frame nf90_set_fill" ) - - call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "fraggle_io_write_frame nf90_put_var time_varid" ) - call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[eslot]), "fraggle_io_write_frame nf90_put_varloop_varid" ) - - charstring = trim(adjustl(REGIME_NAMES(fragments%regime))) - call check( nf90_put_var(nc%id, nc%regime_varid, charstring, start=[1, eslot], count=[len(charstring), 1]), "fraggle_io_write_frame nf90_put_var regime_varid" ) - call check( nf90_put_var(nc%id, nc%Qloss_varid, fragments%Qloss, start=[eslot] ), "fraggle_io_write_frame nf90_put_var Qloss_varid" ) - - do stage = 1,2 - if (allocated(pl)) deallocate(pl) - select case(stage) - case(1) - allocate(pl, source=colliders%pl) - case(2) - allocate(pl, source=fragments%pl) - end select - npl = pl%nbody - do i = 1, npl - idslot = findloc(collision_history%idvals,pl%id(i),dim=1) - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), "fraggle_io_write_frame nf90_put_var id_varid" ) - charstring = trim(adjustl(pl%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], count=[len(charstring), 1]), "fraggle_io_write_frame nf90_put_var name_varid" ) - charstring = trim(adjustl(pl%info(i)%particle_type)) - call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot, stage, eslot], count=[len(charstring), 1, 1]), "fraggle_io_write_frame nf90_put_var particle_type_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var vh_varid" ) - call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "fraggle_io_write_frame nf90_put_var Gmass_varid" ) - call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "fraggle_io_write_frame nf90_put_var radius_varid" ) - call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var Ip_varid" ) - 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 - 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 - end select - end select - return - end subroutine fraggle_io_write_frame - - - module subroutine fraggle_io_log_regime(colliders, fragments) + module subroutine fraggle_io_log_regime(impactors, fragments) !! author: David A. Minton !! !! Writes a log of the results of the collisional regime determination implicit none ! Arguments - class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment object ! Internals character(STRMAX) :: errmsg @@ -240,8 +28,8 @@ module subroutine fraggle_io_log_regime(colliders, fragments) write(LUN, *) "--------------------------------------------------------------------" write(LUN, *) " Fraggle collisional regime determination results" write(LUN, *) "--------------------------------------------------------------------" - write(LUN, *) "True number of colliders : ",colliders%ncoll - write(LUN, *) "Index list of true colliders : ",colliders%idx(1:colliders%ncoll) + write(LUN, *) "True number of impactors : ",impactors%ncoll + write(LUN, *) "Index list of true impactors : ",impactors%idx(1:impactors%ncoll) select case(fragments%regime) case(COLLRESOLVE_REGIME_MERGE) write(LUN, *) "Regime: Merge" diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index dc5b07f58..1563a4ccf 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -7,19 +7,19 @@ !! You should have received a copy of the GNU General Public License along with Swiftest. !! If not, see: https://www.gnu.org/licenses. -submodule(fraggle_classes) s_fraggle_regime +submodule(fraggle_classes) s_encounter_regime use swiftest contains - module subroutine fraggle_regime_colliders(self, fragments, system, param) + module subroutine encounter_regime_impactors(self, fragments, system, param) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! - !! Determine which fragmentation regime the set of colliders will be. This subroutine is a wrapper for the non-polymorphic raggle_regime_collresolve subroutine. + !! Determine which fragmentation regime the set of impactors will be. This subroutine is a wrapper for the non-polymorphic raggle_regime_collresolve subroutine. !! It converts to SI units prior to calling implicit none ! Arguments - class(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object + class(collision_impactors), intent(inout) :: self !! Fraggle impactors object class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters @@ -30,22 +30,22 @@ module subroutine fraggle_regime_colliders(self, fragments, system, param) real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si, runit real(DP) :: mlr, mslr, mtot, dentot - associate(colliders => self) + associate(impactors => self) ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them to the regime determination subroutine - if (colliders%mass(1) > colliders%mass(2)) then + if (impactors%mass(1) > impactors%mass(2)) then jtarg = 1 jproj = 2 else jtarg = 2 jproj = 1 end if - mass_si(:) = colliders%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system - radius_si(:) = colliders%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system + mass_si(:) = impactors%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system + radius_si(:) = impactors%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The two-body equivalent density of the collider system - x1_si(:) = colliders%rb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system - v1_si(:) = colliders%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system - x2_si(:) = colliders%rb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system - v2_si(:) = colliders%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system + x1_si(:) = impactors%rb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system + v1_si(:) = impactors%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system + x2_si(:) = impactors%rb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system + v2_si(:) = impactors%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system Mcb_si = system%cb%mass * param%MU2KG !! The central body mass of the system select type(param) class is (symba_parameters) @@ -58,7 +58,7 @@ module subroutine fraggle_regime_colliders(self, fragments, system, param) dentot = sum(mass_si(:) * density_si(:)) / mtot !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime - call fraggle_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & + call encounter_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & x1_si(:), x2_si(:), v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), & min_mfrag_si, fragments%regime, mlr, mslr, fragments%Qloss) @@ -67,27 +67,27 @@ module subroutine fraggle_regime_colliders(self, fragments, system, param) fragments%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) ! Find the center of mass of the collisional system - fragments%mtot = sum(colliders%mass(:)) - 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 + fragments%mtot = sum(impactors%mass(:)) + fragments%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot + fragments%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot ! Find the point of impact between the two bodies - runit(:) = colliders%rb(:,2) - colliders%rb(:,1) + runit(:) = impactors%rb(:,2) - impactors%rb(:,1) runit(:) = runit(:) / (.mag. runit(:)) - fragments%rbimp(:) = colliders%rb(:,1) + colliders%radius(1) * runit(:) + fragments%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * runit(:) ! Convert quantities back to the system units and save them into the fragment system fragments%mass_dist(:) = (fragments%mass_dist(:) / param%MU2KG) fragments%Qloss = fragments%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG - call fraggle_io_log_regime(colliders, fragments) + call fraggle_io_log_regime(impactors, fragments) end associate return - end subroutine fraggle_regime_colliders + end subroutine encounter_regime_impactors - subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & + subroutine encounter_regime_collresolve(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & regime, Mlr, Mslr, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -379,6 +379,6 @@ function calc_c_star(Rc1) result(c_star) return end function calc_c_star - end subroutine fraggle_regime_collresolve + end subroutine encounter_regime_collresolve -end submodule s_fraggle_regime \ No newline at end of file +end submodule s_encounter_regime \ No newline at end of file diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 3fbe50d2f..122f593ab 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -35,7 +35,7 @@ module subroutine fraggle_set_budgets_fragments(self) end subroutine fraggle_set_budgets_fragments - module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) + module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) !! author: David A. Minton !! !! Sets the mass of fragments based on the mass distribution returned by the regime calculation. @@ -44,7 +44,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) implicit none ! Arguments class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals integer(I4B) :: i, jproj, jtarg, nfrag, istart @@ -61,9 +61,9 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) associate(fragments => self) ! Get mass weighted mean of Ip and density - volume(1:2) = 4._DP / 3._DP * PI * colliders%radius(1:2)**3 - Ip_avg(:) = (colliders%mass(1) * colliders%Ip(:,1) + colliders%mass(2) * colliders%Ip(:,2)) / fragments%mtot - if (colliders%mass(1) > colliders%mass(2)) then + volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 + Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / fragments%mtot + if (impactors%mass(1) > impactors%mass(2)) then jtarg = 1 jproj = 2 else @@ -73,7 +73,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) select case(fragments%regime) case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN) - ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of fraggle_regime. + ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of encounter_regime. ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on ! the limits bracketed above and the model size distribution of fragments. ! Check to see if our size distribution would give us a smaller number of fragments than the maximum number @@ -105,15 +105,15 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) call fragments%setup(1, param) fragments%mass(1) = fragments%mass_dist(1) - fragments%radius(1) = colliders%radius(jtarg) + fragments%radius(1) = impactors%radius(jtarg) fragments%density(1) = fragments%mass_dist(1) / volume(jtarg) - if (param%lrotation) fragments%Ip(:, 1) = colliders%Ip(:,1) + if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) return case default write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",fragments%regime end select - ! Make the first two bins the same as the Mlr and Mslr values that came from fraggle_regime + ! Make the first two bins the same as the Mlr and Mslr values that came from encounter_regime fragments%mass(1) = fragments%mass_dist(iMlr) fragments%mass(2) = fragments%mass_dist(iMslr) @@ -141,9 +141,9 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) ! Compute physical properties of the new fragments select case(fragments%regime) case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment - fragments%radius(1) = colliders%radius(jtarg) + fragments%radius(1) = impactors%radius(jtarg) fragments%density(1) = fragments%mass_dist(iMlr) / volume(jtarg) - fragments%Ip(:, 1) = colliders%Ip(:,1) + fragments%Ip(:, 1) = impactors%Ip(:,1) istart = 2 case default istart = 1 @@ -160,61 +160,66 @@ module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) end subroutine fraggle_set_mass_dist_fragments - module subroutine fraggle_set_coordinate_system(self, colliders) + module subroutine encounter_set_coordinate_system(self, impactors) !! author: David A. Minton !! !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. implicit none ! Arguments class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object ! Internals integer(I4B) :: i real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot - real(DP) :: r_col_norm, v_col_norm, L_mag + real(DP) :: L_mag real(DP), dimension(NDIM, self%nbody) :: L_sigma associate(fragments => self, nfrag => self%nbody) - delta_v(:) = colliders%vb(:, 2) - colliders%vb(:, 1) - v_col_norm = .mag. delta_v(:) - delta_r(:) = colliders%rb(:, 2) - colliders%rb(:, 1) - r_col_norm = .mag. delta_r(:) + delta_v(:) = impactors%vb(:, 2) - impactors%vb(:, 1) + delta_r(:) = impactors%rb(:, 2) - impactors%rb(:, 1) ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector ! and the y-axis aligned with the pre-impact distance vector. - Ltot = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) - fragments%y_coll_unit(:) = delta_r(:) / r_col_norm + + ! y-axis is the separation distance + fragments%y_coll_unit(:) = .unit.delta_r(:) + Ltot = impactors%L_orbit(:,1) + impactors%L_orbit(:,2) + impactors%L_spin(:,1) + impactors%L_spin(:,2) + L_mag = .mag.Ltot(:) if (L_mag > sqrt(tiny(L_mag))) then - fragments%z_coll_unit(:) = Ltot(:) / L_mag - else + fragments%z_coll_unit(:) = .unit.Ltot(:) + else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction call random_number(fragments%z_coll_unit(:)) - fragments%z_coll_unit(:) = fragments%z_coll_unit(:) / (.mag.fragments%z_coll_unit(:)) + fragments%z_coll_unit(:) = .unit.fragments%z_coll_unit(:) end if + ! The cross product of the y- by z-axis will give us the x-axis fragments%x_coll_unit(:) = fragments%y_coll_unit(:) .cross. fragments%z_coll_unit(:) - + + fragments%v_coll_unit(:) = .unit.delta_v(:) + if (.not.any(fragments%r_coll(:,:) > 0.0_DP)) return fragments%rmag(:) = .mag. fragments%r_coll(:,:) - - call random_number(L_sigma(:,:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, - ! otherwise we can get an ill-conditioned system - + + ! Randomize the tangential velocity direction. + ! This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, otherwise we can get an ill-conditioned system + call random_number(L_sigma(:,:)) do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP) - fragments%v_r_unit(:, i) = fragments%r_coll(:, i) / fragments%rmag(i) fragments%v_n_unit(:, i) = fragments%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) - fragments%v_n_unit(:, i) = fragments%v_n_unit(:, i) / (.mag. fragments%v_n_unit(:, i)) - fragments%v_t_unit(:, i) = fragments%v_n_unit(:, i) .cross. fragments%v_r_unit(:, i) - fragments%v_t_unit(:, i) = fragments%v_t_unit(:, i) / (.mag. fragments%v_t_unit(:, i)) end do + ! Define the radial, normal, and tangential unit vectors for each individual fragment + fragments%v_r_unit(:,:) = .unit. fragments%r_coll(:,:) + fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:) + fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:)) + end associate return - end subroutine fraggle_set_coordinate_system + end subroutine encounter_set_coordinate_system - module subroutine fraggle_set_natural_scale_factors(self, colliders) + module subroutine fraggle_set_natural_scale_factors(self, impactors) !! author: David A. Minton !! !! Scales dimenional quantities to ~O(1) with respect to the collisional system. @@ -222,33 +227,33 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) implicit none ! Arguments class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object ! Internals integer(I4B) :: i associate(fragments => self) ! Set scale factors - fragments%Escale = 0.5_DP * (colliders%mass(1) * dot_product(colliders%vb(:,1), colliders%vb(:,1)) & - + colliders%mass(2) * dot_product(colliders%vb(:,2), colliders%vb(:,2))) - fragments%dscale = sum(colliders%radius(:)) + fragments%Escale = 0.5_DP * (impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & + + impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2))) + fragments%dscale = sum(impactors%radius(:)) fragments%mscale = fragments%mtot fragments%vscale = sqrt(fragments%Escale / fragments%mscale) fragments%tscale = fragments%dscale / fragments%vscale fragments%Lscale = fragments%mscale * fragments%dscale * fragments%vscale - ! Scale all dimensioned quantities of colliders and fragments + ! Scale all dimensioned quantities of impactors and fragments fragments%rbcom(:) = fragments%rbcom(:) / fragments%dscale fragments%vbcom(:) = fragments%vbcom(:) / fragments%vscale fragments%rbimp(:) = fragments%rbimp(:) / fragments%dscale - colliders%rb(:,:) = colliders%rb(:,:) / fragments%dscale - colliders%vb(:,:) = colliders%vb(:,:) / fragments%vscale - colliders%mass(:) = colliders%mass(:) / fragments%mscale - colliders%radius(:) = colliders%radius(:) / fragments%dscale - colliders%L_spin(:,:) = colliders%L_spin(:,:) / fragments%Lscale - colliders%L_orbit(:,:) = colliders%L_orbit(:,:) / fragments%Lscale + impactors%rb(:,:) = impactors%rb(:,:) / fragments%dscale + impactors%vb(:,:) = impactors%vb(:,:) / fragments%vscale + impactors%mass(:) = impactors%mass(:) / fragments%mscale + impactors%radius(:) = impactors%radius(:) / fragments%dscale + impactors%L_spin(:,:) = impactors%L_spin(:,:) / fragments%Lscale + impactors%L_orbit(:,:) = impactors%L_orbit(:,:) / fragments%Lscale do i = 1, 2 - colliders%rot(:,i) = colliders%L_spin(:,i) / (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i)) + impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do fragments%mtot = fragments%mtot / fragments%mscale @@ -261,7 +266,7 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) end subroutine fraggle_set_natural_scale_factors - module subroutine fraggle_set_original_scale_factors(self, colliders) + module subroutine fraggle_set_original_scale_factors(self, impactors) !! author: David A. Minton !! !! Restores dimenional quantities back to the system units @@ -269,7 +274,7 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) implicit none ! Arguments class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object ! Internals integer(I4B) :: i logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes @@ -284,13 +289,13 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) fragments%vbcom(:) = fragments%vbcom(:) * fragments%vscale fragments%rbimp(:) = fragments%rbimp(:) * fragments%dscale - colliders%mass = colliders%mass * fragments%mscale - colliders%radius = colliders%radius * fragments%dscale - colliders%rb = colliders%rb * fragments%dscale - colliders%vb = colliders%vb * fragments%vscale - colliders%L_spin = colliders%L_spin * fragments%Lscale + impactors%mass = impactors%mass * fragments%mscale + impactors%radius = impactors%radius * fragments%dscale + impactors%rb = impactors%rb * fragments%dscale + impactors%vb = impactors%vb * fragments%vscale + impactors%L_spin = impactors%L_spin * fragments%Lscale do i = 1, 2 - colliders%rot(:,i) = colliders%L_spin(:,i) * (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i)) + impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do fragments%mtot = fragments%mtot * fragments%mscale diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 index bd660dc98..3d8916967 100644 --- a/src/fraggle/fraggle_setup.f90 +++ b/src/fraggle/fraggle_setup.f90 @@ -52,30 +52,19 @@ module subroutine fraggle_setup_fragments(self, n, param) integer(I4B), intent(in) :: n class(swiftest_parameters), intent(in) :: param - call setup_pl(self, n, param) + call self%collision_fragments%setup(n, param) if (n < 0) return - if (allocated(self%r_coll)) deallocate(self%r_coll) - if (allocated(self%v_coll)) deallocate(self%v_coll) - if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) - if (allocated(self%v_t_unit)) deallocate(self%v_t_unit) - if (allocated(self%v_n_unit)) deallocate(self%v_n_unit) - if (allocated(self%rmag)) deallocate(self%rmag) - if (allocated(self%rotmag)) deallocate(self%rotmag) if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) if (allocated(self%v_t_mag)) deallocate(self%v_t_mag) + if (allocated(self%v_n_mag)) deallocate(self%v_t_mag) if (n == 0) return - allocate(self%r_coll(NDIM,n)) - allocate(self%v_coll(NDIM,n)) - allocate(self%v_r_unit(NDIM,n)) - allocate(self%v_t_unit(NDIM,n)) - allocate(self%v_n_unit(NDIM,n)) - allocate(self%rmag(n)) allocate(self%rotmag(n)) allocate(self%v_r_mag(n)) allocate(self%v_t_mag(n)) + allocate(self%v_n_mag(n)) call self%reset() diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index c4fffe4a6..bb32167b1 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -11,14 +11,14 @@ use swiftest contains - module subroutine fraggle_util_add_fragments_to_system(fragments, colliders, system, param) + module subroutine fraggle_util_add_fragments_to_system(fragments, impactors, system, param) !! Author: David A. Minton !! !! Adds fragments to the temporary system pl object implicit none ! Arguments class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment system object - class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters ! Internals @@ -44,9 +44,9 @@ module subroutine fraggle_util_add_fragments_to_system(fragments, colliders, sys pl%Ip(:,npl_before+1:npl_after) = fragments%Ip(:,1:nfrag) pl%rot(:,npl_before+1:npl_after) = fragments%rot(:,1:nfrag) end if - ! This will remove the colliders from the system since we've replaced them with fragments + ! This will remove the impactors from the system since we've replaced them with fragments lexclude(1:npl_after) = .false. - lexclude(colliders%idx(1:colliders%ncoll)) = .true. + lexclude(impactors%idx(1:impactors%ncoll)) = .true. where(lexclude(1:npl_after)) pl%status(1:npl_after) = INACTIVE elsewhere @@ -132,18 +132,18 @@ end subroutine fraggle_util_construct_temporary_system - module subroutine fraggle_util_final_colliders(self) + module subroutine fraggle_util_final_impactors(self) !! author: David A. Minton !! !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(fraggle_colliders), intent(inout) :: self !! Fraggle encountar storage object + type(collision_impactors), intent(inout) :: self !! Fraggle encountar storage object if (allocated(self%idx)) deallocate(self%idx) return - end subroutine fraggle_util_final_colliders + end subroutine fraggle_util_final_impactors module subroutine fraggle_util_final_fragments(self) @@ -168,137 +168,14 @@ module subroutine fraggle_util_final_fragments(self) end subroutine fraggle_util_final_fragments - module subroutine fraggle_util_final_snapshot(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(fraggle_snapshot), intent(inout) :: self !! Fraggle encountar storage object - - call encounter_util_final_snapshot(self%encounter_snapshot) - - return - end subroutine fraggle_util_final_snapshot - - - module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore) - !! Author: David A. Minton - !! - !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) - !! This subrourtine works by building a temporary internal massive body object out of the non-excluded bodies and optionally with fragments appended. - !! This will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it is in the main program. - !! This will temporarily expand the massive body object in a temporary system object called tmpsys to feed it into symba_energy - implicit none - ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters - logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa - ! Internals - class(swiftest_nbody_system), allocatable, save :: tmpsys - class(swiftest_parameters), allocatable, save :: tmpparam - integer(I4B) :: npl_before, npl_after - - associate(fragments => self, nfrag => self%nbody, pl => system%pl, cb => system%cb) - - ! Because we're making a copy of the massive body object with the excludes/fragments appended, we need to deallocate the - ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by - ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this - ! subroutine should be called relatively infrequently, it shouldn't matter too much. - - npl_before = pl%nbody - npl_after = npl_before + nfrag - - if (lbefore) then - call fraggle_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) - ! Build the exluded body logical mask for the *before* case: Only the original bodies are used to compute energy and momentum - tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = ACTIVE - tmpsys%pl%status(npl_before+1:npl_after) = INACTIVE - else - if (.not.allocated(tmpsys)) then - write(*,*) "Error in fraggle_util_get_energy_momentum. " // & - " This must be called with lbefore=.true. at least once before calling it with lbefore=.false." - call util_exit(FAILURE) - end if - ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum - call fraggle_util_add_fragments_to_system(fragments, colliders, tmpsys, tmpparam) - tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = INACTIVE - tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE - end if - - if (param%lflatten_interactions) call tmpsys%pl%flatten(param) - - call tmpsys%get_energy_and_momentum(param) - - - ! Calculate the current fragment energy and momentum balances - if (lbefore) then - fragments%Lorbit_before(:) = tmpsys%Lorbit(:) - fragments%Lspin_before(:) = tmpsys%Lspin(:) - fragments%Ltot_before(:) = tmpsys%Ltot(:) - fragments%ke_orbit_before = tmpsys%ke_orbit - fragments%ke_spin_before = tmpsys%ke_spin - fragments%pe_before = tmpsys%pe - fragments%Etot_before = tmpsys%te - else - fragments%Lorbit_after(:) = tmpsys%Lorbit(:) - fragments%Lspin_after(:) = tmpsys%Lspin(:) - fragments%Ltot_after(:) = tmpsys%Ltot(:) - fragments%ke_orbit_after = tmpsys%ke_orbit - fragments%ke_spin_after = tmpsys%ke_spin - fragments%pe_after = tmpsys%pe - fragments%Etot_after = tmpsys%te - (fragments%pe_after - fragments%pe_before) ! Gotta be careful with PE when number of bodies changes. - end if - end associate - - return - end subroutine fraggle_util_get_energy_momentum - - - module subroutine fraggle_util_get_idvalues_snapshot(self, idvals) - !! author: David A. Minton - !! - !! Returns an array of all id values saved in this snapshot - implicit none - ! Arguments - class(fraggle_snapshot), intent(in) :: self !! Fraggle snapshot object - integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot - ! Internals - integer(I4B) :: ncoll, nfrag - - if (allocated(self%colliders)) then - ncoll = self%colliders%pl%nbody - else - ncoll = 0 - end if - - if (allocated(self%fragments)) then - nfrag = self%fragments%pl%nbody - else - nfrag = 0 - end if - - if (ncoll + nfrag == 0) return - allocate(idvals(ncoll+nfrag)) - - if (ncoll > 0) idvals(1:ncoll) = self%colliders%pl%id(:) - if (nfrag > 0) idvals(ncoll+1:ncoll+nfrag) = self%fragments%pl%id(:) - - return - - end subroutine fraggle_util_get_idvalues_snapshot - - - module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_start) + module subroutine fraggle_util_restructure(self, impactors, try, f_spin, r_max_start) !! Author: David A. Minton !! !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints implicit none ! Arguments class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object integer(I4B), intent(in) :: try !! The current number of times Fraggle has tried to find a solution real(DP), intent(inout) :: f_spin !! Fraction of energy/momentum that goes into spin. This decreases ater a failed attempt real(DP), intent(inout) :: r_max_start !! The maximum radial distance that the position calculation starts with. This increases after a failed attempt @@ -310,7 +187,7 @@ module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_s ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions associate(fragments => self) call random_number(delta_r_max) - delta_r_max = sum(colliders%radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) + delta_r_max = sum(impactors%radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) if (try == 1) then ke_tot_deficit = - (fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin) ke_avg_deficit = ke_tot_deficit diff --git a/src/modules/collision_classes.f90 b/src/modules/collision_classes.f90 new file mode 100644 index 000000000..d3ede5fcb --- /dev/null +++ b/src/modules/collision_classes.f90 @@ -0,0 +1,305 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +module collision_classes + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Definition of classes and methods used to determine close encounters + use swiftest_globals + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system, swiftest_pl, swiftest_storage, netcdf_parameters + use encounter_classes, only : encounter_snapshot, encounter_io_parameters, encounter_storage, encounter_io_parameters + implicit none + public + + + !>Symbolic names for collisional outcomes from collresolve_resolve: + integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE = 1 + integer(I4B), parameter :: COLLRESOLVE_REGIME_DISRUPTION = 2 + integer(I4B), parameter :: COLLRESOLVE_REGIME_SUPERCATASTROPHIC = 3 + integer(I4B), parameter :: COLLRESOLVE_REGIME_GRAZE_AND_MERGE = 4 + integer(I4B), parameter :: COLLRESOLVE_REGIME_HIT_AND_RUN = 5 + character(len=*),dimension(5), parameter :: REGIME_NAMES = ["Merge", "Disruption", "Supercatastrophic", "Graze and Merge", "Hit and Run"] + + !******************************************************************************************************************************** + ! collision_impactors class definitions and method interfaces + !******************************************************************************************************************************* + !> Class definition for the variables that describe the bodies involved in the collision + type :: collision_impactors + integer(I4B) :: ncoll !! Number of bodies involved in the collision + integer(I4B), dimension(:), allocatable :: idx !! Index of bodies involved in the collision + real(DP), dimension(NDIM,2) :: rb !! Two-body equivalent position vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision + real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision + real(DP), dimension(2) :: radius !! Two-body equivalent radii of the collider bodies prior to the collision + real(DP) :: Qloss !! Energy lost during the collision + integer(I4B) :: regime !! Collresolve regime code for this collision + real(DP), dimension(:), allocatable :: mass_dist !! Distribution of fragment mass determined by the regime calculation (largest fragment, second largest, and remainder) + + ! Values in a coordinate frame centered on the collider barycenter and collisional system unit vectors + real(DP), dimension(NDIM) :: x_unit !! x-direction unit vector of collisional system + real(DP), dimension(NDIM) :: y_unit !! y-direction unit vector of collisional system + real(DP), dimension(NDIM) :: z_unit !! z-direction unit vector of collisional system + real(DP), dimension(NDIM) :: v_unit !! z-direction unit vector of collisional system + real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider system in system barycentric coordinates + real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider system in system barycentric coordinates + real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider system in system barycentric coordinates + + contains + procedure :: set_coordinate_system => collision_set_coordinate_impactors !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + procedure :: setup => collision_setup_impactors !! Allocates arrays for n fragments in a fragment system. Passing n = 0 deallocates all arrays. + procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions + final :: collision_util_final_impactors !! Finalizer will deallocate all allocatables + end type collision_impactors + + !******************************************************************************************************************************** + ! collision_fragments class definitions and method interfaces + !******************************************************************************************************************************* + !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates + type, abstract, extends(swiftest_pl) :: collision_fragments + real(DP) :: mtot !! Total mass of fragments + real(DP), dimension(:,:), allocatable :: rc !! Position vectors in the collision coordinate frame + real(DP), dimension(:,:), allocatable :: vc !! Velocity vectors in the collision coordinate frame + real(DP), dimension(:), allocatable :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame + real(DP), dimension(:), allocatable :: vmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame + real(DP), dimension(:), allocatable :: rotmag !! Array of rotation magnitudes of individual fragments + real(DP), dimension(:,:), allocatable :: v_r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(:,:), allocatable :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(:,:), allocatable :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame + + contains + procedure :: accel => collision_placeholder_accel !! Placeholder subroutine to fulfill requirement for an accel method + procedure :: kick => collision_placeholder_kick !! Placeholder subroutine to fulfill requirement for a kick method + procedure :: step => collision_placeholder_step !! Placeholder subroutine to fulfill requirement for a step method + procedure :: set_coordinate_system => collision_set_coordinate_fragments !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + procedure :: setup => collision_setup_fragments !! Allocates arrays for n fragments in a Fraggle system. Passing n = 0 deallocates all arrays. + procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays + end type collision_fragments + + type :: collision_system + !! This class defines a collisional system that stores impactors and fragments. This is written so that various collision models (i.e. Fraggle) could potentially be used + !! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments + class(collision_impactors), allocatable :: impactors !! Object containing information on the pre-collision system + class(collision_fragments), allocatable :: fragments !! Object containing information on the post-collision system + class(swiftest_nbody_system), allocatable :: before !! A snapshot of the subset of the system involved in the collision + class(swiftest_nbody_system), allocatable :: after !! A snapshot of the subset of the system containing products of the collision + + ! For the following variables, index 1 refers to the *entire* n-body system in its pre-collisional state and index 2 refers to the system in its post-collisional state + real(DP), dimension(NDIM,2) :: Lorbit !! Before/after orbital angular momentum + real(DP), dimension(NDIM,2) :: Lspin !! Before/after spin angular momentum + real(DP), dimension(NDIM,2) :: Ltot !! Before/after total system angular momentum + real(DP), dimension(2) :: ke_orbit !! Before/after orbital kinetic energy + real(DP), dimension(2) :: ke_spin !! Before/after spin kinetic energy + real(DP), dimension(2) :: pe !! Before/after potential energy + real(DP), dimension(2) :: Etot !! Before/after total system energy + contains + procedure :: generate_fragments => collision_generate_fragment_system !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + procedure :: regime => collision_regime_system !! Determine which fragmentation regime the set of impactors will be + procedure :: setup => collision_setup_system !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots + procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) + procedure :: reset => collision_util_reset_system !! Deallocates all allocatables + final :: collision_util_final_system !! Finalizer will deallocate all allocatables + end type collision_system + + + !> NetCDF dimension and variable names for the enounter save object + type, extends(encounter_io_parameters) :: collision_io_parameters + contains + procedure :: initialize => collision_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object + end type collision_io_parameters + + type, extends(encounter_snapshot) :: collision_snapshot + logical :: lcollision !! Indicates that this snapshot contains at least one collision + class(collision_system), allocatable :: collision_system !! impactors object at this snapshot + contains + procedure :: write_frame => collision_io_write_frame_snapshot !! Writes a frame of encounter data to file + procedure :: get_idvals => collision_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot + final :: collision_util_final_snapshot !! Finalizer deallocates all allocatables + end type collision_snapshot + + !> A class that that is used to store simulation history data between file output + type, extends(swiftest_storage) :: collision_storage + contains + procedure :: dump => collision_io_dump !! Dumps contents of encounter history to file + procedure :: take_snapshot => collision_util_snapshot !! Take a minimal snapshot of the system through an encounter + procedure :: make_index_map => collision_util_index_map !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id + final :: collision_util_final_storage !! Finalizer deallocates all allocatables + end type collision_storage + + interface + module subroutine collision_generate_fragment_system(self, system, param, lfailure) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(collision_system), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? + end subroutine collision_generate_fragment_system + + module subroutine collision_io_dump(self, param) + implicit none + class(collision_storage(*)), intent(inout) :: self !! Collision storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine collision_io_dump + + module subroutine collision_io_initialize(self, param) + implicit none + class(collision_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine collision_io_initialize + + module subroutine collision_io_write_frame_snapshot(self, nc, param) + implicit none + class(collision_snapshot), intent(in) :: self !! Swiftest encounter structure + class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine collision_io_write_frame_snapshot + + !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class + module subroutine collision_placeholder_accel(self, system, param, t, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(collision_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + end subroutine collision_placeholder_accel + + module subroutine collision_placeholder_kick(self, system, param, t, dt, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(collision_fragments), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + end subroutine collision_placeholder_kick + + module subroutine collision_placeholder_step(self, system, param, t, dt) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(collision_fragments), intent(inout) :: self !! Helio massive body particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Stepsiz + end subroutine collision_placeholder_step + + module subroutine collision_regime_system(self, system, param) + implicit none + class(collision_system), intent(inout) :: self !! Collision system object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + end subroutine collision_regime_system + + module subroutine collision_set_coordinate_impactors(self) + implicit none + class(collision_impactors), intent(inout) :: self !! Collider system object + end subroutine collision_set_coordinate_impactors + + module subroutine collision_set_coordinate_fragments(self) + implicit none + class(collision_fragments), intent(inout) :: self !! Fragment system object + end subroutine collision_set_coordinate_fragments + + module subroutine collision_setup_fragments(self, n, param) + implicit none + class(collision_fragments), intent(inout) :: self !! Fragment system object + integer(I4B), intent(in) :: n !! Number of fragments + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + end subroutine collision_setup_fragments + + module subroutine collision_setup_impactors(self, system, param) + implicit none + class(collision_impactors), intent(inout) :: self !! Fragment system object + class(swiftest_nbody_system), intent(in) :: system + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + end subroutine collision_setup_impactors + + module subroutine collision_setup_system(self, system, param) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(collision_system), intent(inout) :: self !! Encounter collision system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine collision_setup_system + + module subroutine collision_util_reset_fragments(self) + implicit none + class(collision_fragments), intent(inout) :: self + end subroutine collision_util_reset_fragments + + module subroutine collision_util_final_impactors(self) + implicit none + type(collision_impactors), intent(inout) :: self !! Collision impactors storage object + end subroutine collision_util_final_impactors + + module subroutine collision_util_final_storage(self) + implicit none + type(collision_storage(*)), intent(inout) :: self !! SyMBA nbody system object + end subroutine collision_util_final_storage + + module subroutine collision_util_final_snapshot(self) + implicit none + type(collision_snapshot), intent(inout) :: self !! Fraggle storage snapshot object + end subroutine collision_util_final_snapshot + + module subroutine collision_util_final_system(self) + implicit none + type(collision_system), intent(inout) :: self !! Collision system object + end subroutine collision_util_final_system + + module subroutine collision_util_get_idvalues_snapshot(self, idvals) + implicit none + class(collision_snapshot), intent(in) :: self !! Fraggle snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + end subroutine collision_util_get_idvalues_snapshot + + module subroutine collision_util_get_energy_momentum(self, system, param, lbefore) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(collision_system), intent(inout) :: self !! Encounter collision system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with impactors included and fragments excluded or vice versa + end subroutine collision_util_get_energy_momentum + + module subroutine collision_util_index_map(self) + implicit none + class(collision_storage(*)), intent(inout) :: self !! Collision storage object + end subroutine collision_util_index_map + + module subroutine collision_util_reset_impactors(self) + implicit none + class(collision_impactors), intent(inout) :: self + end subroutine collision_util_reset_impactors + + module subroutine collision_util_reset_system(self) + implicit none + class(collision_system), intent(inout) :: self + end subroutine collision_util_reset_system + + module subroutine collision_util_snapshot(self, param, system, t, arg) + implicit none + class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store + real(DP), intent(in), optional :: t !! Time of snapshot if different from system time + character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. + end subroutine collision_util_snapshot + + end interface + +end module collision_classes + diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 2b313eeb5..0eb0dc3f7 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -12,7 +12,7 @@ module encounter_classes !! !! Definition of classes and methods used to determine close encounters use swiftest_globals - use swiftest_classes + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system, swiftest_cb, swiftest_tp, swiftest_pl, swiftest_storage, netcdf_parameters implicit none public @@ -50,37 +50,17 @@ module encounter_classes real(DP) :: t !! Simulation time when snapshot was taken integer(I8B) :: iloop !! Loop number at time of snapshot contains - procedure :: write_frame => encounter_io_write_frame !! Writes a frame of encounter data to file - procedure :: get_idvals => encounter_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot + procedure :: write_frame => encounter_io_write_frame_snapshot !! Writes a frame of encounter data to file + procedure :: get_idvals => encounter_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot final :: encounter_util_final_snapshot end type encounter_snapshot - !> NetCDF dimension and variable names for the enounter save object - type, extends(netcdf_parameters) :: encounter_io_parameters - character(NAMELEN) :: loop_varname = "loopnum" !! Loop number for encounter - integer(I4B) :: loop_varid !! ID for the recursion level variable - integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot - integer(I4B) :: name_dimsize = 0 !! Number of potential id values in snapshot - integer(I4B) :: file_number = 1 !! The number to append on the output file - contains - procedure :: initialize => encounter_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object - end type encounter_io_parameters - - !> A class that that is used to store simulation history data between file output - type, extends(swiftest_storage) :: collision_storage - contains - procedure :: dump => encounter_io_dump_collision !! Dumps contents of encounter history to file - procedure :: take_snapshot => encounter_util_snapshot_collision !! Take a minimal snapshot of the system through an encounter - procedure :: make_index_map => encounter_util_index_map_collision !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - final :: encounter_util_final_collision_storage - end type collision_storage - !> A class that that is used to store simulation history data between file output type, extends(swiftest_storage) :: encounter_storage contains - procedure :: dump => encounter_io_dump_encounter !! Dumps contents of encounter history to file - procedure :: make_index_map => encounter_util_index_map_encounter !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - procedure :: take_snapshot => encounter_util_snapshot_encounter !! Take a minimal snapshot of the system through an encounter + procedure :: dump => encounter_io_dump !! Dumps contents of encounter history to file + procedure :: make_index_map => encounter_util_index_map !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id + procedure :: take_snapshot => encounter_util_snapshot !! Take a minimal snapshot of the system through an encounter final :: encounter_util_final_storage end type encounter_storage @@ -216,30 +196,18 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_single_list - module subroutine encounter_io_dump_collision(self, param) - implicit none - class(collision_storage(*)), intent(inout) :: self !! Collision storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine encounter_io_dump_collision - - module subroutine encounter_io_dump_encounter(self, param) - implicit none - class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine encounter_io_dump_encounter - - module subroutine encounter_io_initialize(self, param) + module subroutine encounter_io_dump(self, param) implicit none - class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param - end subroutine encounter_io_initialize + class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine encounter_io_dump - module subroutine encounter_io_write_frame(self, nc, param) + module subroutine encounter_io_write_frame_snapshot(self, nc, param) implicit none class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine encounter_io_write_frame + end subroutine encounter_io_write_frame_snapshot module subroutine encounter_setup_aabb(self, n, n_last) implicit none @@ -282,11 +250,6 @@ module subroutine encounter_util_final_aabb(self) type(encounter_bounding_box_1D), intent(inout) :: self !!Bounding box structure along a single dimension end subroutine encounter_util_final_aabb - module subroutine encounter_util_final_collision_storage(self) - implicit none - type(collision_storage(*)), intent(inout) :: self !! SyMBA nbody system object - end subroutine encounter_util_final_collision_storage - module subroutine encounter_util_final_list(self) implicit none type(encounter_list), intent(inout) :: self !! Swiftest encounter list object @@ -308,15 +271,10 @@ module subroutine encounter_util_get_idvalues_snapshot(self, idvals) integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot end subroutine encounter_util_get_idvalues_snapshot - module subroutine encounter_util_index_map_collision(self) - implicit none - class(collision_storage(*)), intent(inout) :: self !! Collision storage object - end subroutine encounter_util_index_map_collision - - module subroutine encounter_util_index_map_encounter(self) + module subroutine encounter_util_index_map(self) implicit none class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object - end subroutine encounter_util_index_map_encounter + end subroutine encounter_util_index_map module subroutine encounter_util_resize_list(self, nnew) implicit none @@ -324,23 +282,14 @@ module subroutine encounter_util_resize_list(self, nnew) integer(I8B), intent(in) :: nnew !! New size of list needed end subroutine encounter_util_resize_list - module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) - implicit none - class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store - real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. - end subroutine encounter_util_snapshot_collision - - module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) + module subroutine encounter_util_snapshot(self, param, system, t, arg) implicit none class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store real(DP), intent(in), optional :: t !! Time of snapshot if different from system time character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) - end subroutine encounter_util_snapshot_encounter + end subroutine encounter_util_snapshot module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 235451379..6ef73ffb1 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -14,77 +14,23 @@ module fraggle_classes use swiftest_globals use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system, swiftest_cb, swiftest_pl, swiftest_storage, netcdf_parameters use encounter_classes, only : encounter_snapshot, encounter_io_parameters, encounter_storage + use collision_classes, only : collision_impactors, collision_fragments, collision_system implicit none public integer(I4B), parameter :: FRAGGLE_NMASS_DIST = 3 !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder) character(len=*), parameter :: FRAGGLE_LOG_OUT = "fraggle.log" !! Name of log file for Fraggle diagnostic information - !******************************************************************************************************************************** - ! fraggle_colliders class definitions and method interfaces - !******************************************************************************************************************************* - !> Class definition for the variables that describe the bodies involved in the collision - type :: fraggle_colliders - integer(I4B) :: ncoll !! Number of bodies involved in the collision - integer(I4B), dimension(:), allocatable :: idx !! Index of bodies involved in the collision - real(DP), dimension(NDIM,2) :: rb !! Two-body equivalent position vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision - real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision - real(DP), dimension(2) :: radius !! Two-body equivalent radii of the collider bodies prior to the collision - class(swiftest_pl), allocatable :: pl !! A snapshot of the planets involved in the collision - contains - procedure :: regime => fraggle_regime_colliders !! Determine which fragmentation regime the set of colliders will be - final :: fraggle_util_final_colliders !! Finalizer will deallocate all allocatables - end type fraggle_colliders - !******************************************************************************************************************************** ! fraggle_fragments class definitions and method interfaces !******************************************************************************************************************************* !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates - type, extends(swiftest_pl) :: fraggle_fragments - real(DP) :: mtot !! Total mass of fragments - real(DP) :: Qloss !! Energy lost during the collision - real(DP), dimension(FRAGGLE_NMASS_DIST) :: mass_dist !! Distribution of fragment mass determined by the regime calculation (largest fragment, second largest, and remainder) - integer(I4B) :: regime !! Collresolve regime code for this collision - - ! Values in a coordinate frame centered on the collider barycenter and collisional system unit vectors (these are used internally by the fragment generation subroutine) - real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider system in system barycentric coordinates - real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider system in system barycentric coordinates - real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider system in system barycentric coordinates - real(DP), dimension(NDIM) :: x_coll_unit !! x-direction unit vector of collisional system - real(DP), dimension(NDIM) :: y_coll_unit !! y-direction unit vector of collisional system - real(DP), dimension(NDIM) :: z_coll_unit !! z-direction unit vector of collisional system - real(DP), dimension(:,:), allocatable :: r_coll !! Array of fragment position vectors in the collisional coordinate frame - real(DP), dimension(:,:), allocatable :: v_coll !! Array of fragment velocity vectors in the collisional coordinate frame - real(DP), dimension(:,:), allocatable :: v_r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(:,:), allocatable :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(:,:), allocatable :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(:), allocatable :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame - real(DP), dimension(:), allocatable :: rotmag !! Array of rotation magnitudes of individual fragments - real(DP), dimension(:), allocatable :: v_r_mag !! Array of radial direction velocity magnitudes of individual fragments - real(DP), dimension(:), allocatable :: v_t_mag !! Array of tangential direction velocity magnitudes of individual fragments - class(swiftest_pl), allocatable :: pl !! A snapshot of the fragments created in the collision + type, extends(collision_fragments) :: fraggle_fragments - ! Energy and momentum book-keeping variables that characterize the whole system of fragments - real(DP) :: ke_orbit !! Current orbital kinetic energy of the system of fragments in the collisional frame - real(DP) :: ke_spin !! Current spin kinetic energy of the system of fragments in the collisional frame - real(DP), dimension(NDIM) :: L_orbit !! Current orbital angular momentum of the system of fragments in the collisional frame - real(DP), dimension(NDIM) :: L_spin !! Current spin angular momentum of the system of fragments in the collisional frame - real(DP) :: ke_budget !! Total kinetic energy budget for the system of fragmens in the collisional frame - real(DP), dimension(NDIM) :: L_budget !! Total angular momentum budget for the system of fragmens in the collisional frame + real(DP), dimension(:), allocatable :: v_r_mag !! Array of radial direction velocity magnitudes of individual fragments + real(DP), dimension(:), allocatable :: v_t_mag !! Array of tangential direction velocity magnitudes of individual fragments + real(DP), dimension(:), allocatable :: v_n_mag !! Array of normal direction velocity magnitudes of individual fragments - ! For the following variables, "before" refers to the *entire* n-body system in its pre-collisional state and "after" refers to the system in its post-collisional state - real(DP), dimension(NDIM) :: Lorbit_before, Lorbit_after !! Before/after orbital angular momentum - real(DP), dimension(NDIM) :: Lspin_before, Lspin_after !! Before/after spin angular momentum - real(DP), dimension(NDIM) :: Ltot_before, Ltot_after !! Before/after total system angular momentum - real(DP) :: ke_orbit_before, ke_orbit_after !! Before/after orbital kinetic energy - real(DP) :: ke_spin_before, ke_spin_after !! Before/after spin kinetic energy - real(DP) :: pe_before, pe_after !! Before/after potential energy - real(DP) :: Etot_before, Etot_after !! Before/after total system energy ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor @@ -95,159 +41,66 @@ module fraggle_classes real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains procedure :: generate_fragments => fraggle_generate_fragments !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. - procedure :: accel => fraggle_placeholder_accel !! Placeholder subroutine to fulfill requirement for an accel method - procedure :: kick => fraggle_placeholder_kick !! Placeholder subroutine to fulfill requirement for a kick method - procedure :: step => fraggle_placeholder_step !! Placeholder subroutine to fulfill requirement for a step method procedure :: set_budgets => fraggle_set_budgets_fragments !! Sets the energy and momentum budgets of the fragments based on the collider value - procedure :: set_coordinate_system => fraggle_set_coordinate_system !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. procedure :: set_mass_dist => fraggle_set_mass_dist_fragments !! Sets the distribution of mass among the fragments depending on the regime type procedure :: set_natural_scale => fraggle_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. procedure :: set_original_scale => fraggle_set_original_scale_factors !! Restores dimenional quantities back to the original system units procedure :: setup => fraggle_setup_fragments !! Allocates arrays for n fragments in a Fraggle system. Passing n = 0 deallocates all arrays. procedure :: reset => fraggle_setup_reset_fragments !! 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) procedure :: get_ang_mtm => fraggle_util_ang_mtm !! Calcualtes the current angular momentum of the fragments - procedure :: get_energy_and_momentum => fraggle_util_get_energy_momentum !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: restructure => fraggle_util_restructure !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints final :: fraggle_util_final_fragments !! Finalizer will deallocate all allocatables end type fraggle_fragments - !! NetCDF dimension and variable names for the enounter save object - type, extends(encounter_io_parameters) :: fraggle_io_parameters - integer(I4B) :: stage_dimid !! ID for the stage dimension - integer(I4B) :: stage_varid !! ID for the stage variable - character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after) - character(len=6), dimension(2) :: stage_coords = ["before", "after"] !! The stage coordinate labels - - character(NAMELEN) :: event_dimname = "collision" !! Name of collision event dimension - integer(I4B) :: event_dimid !! ID for the collision event dimension - integer(I4B) :: event_varid !! ID for the collision event variable - integer(I4B) :: event_dimsize = 0 !! Number of events - - character(NAMELEN) :: Qloss_varname = "Qloss" !! name of the energy loss variable - integer(I4B) :: Qloss_varid !! ID for the energy loss variable - character(NAMELEN) :: regime_varname = "regime" !! name of the collision regime variable - integer(I4B) :: regime_varid !! ID for the collision regime variable - - contains - procedure :: initialize => fraggle_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object - end type fraggle_io_parameters - - type, extends(encounter_snapshot) :: fraggle_snapshot - logical :: lcollision !! Indicates that this snapshot contains at least one collision - class(fraggle_colliders), allocatable :: colliders !! Colliders object at this snapshot - class(fraggle_fragments), allocatable :: fragments !! Fragments object at this snapshot - contains - procedure :: write_frame => fraggle_io_write_frame !! Writes a frame of encounter data to file - procedure :: get_idvals => fraggle_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot - final :: fraggle_util_final_snapshot - end type fraggle_snapshot interface module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies + class(collision_impactors), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fraggle_generate_fragments - module subroutine fraggle_io_initialize_output(self, param) - implicit none - class(fraggle_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine fraggle_io_initialize_output - - module subroutine fraggle_io_write_frame(self, nc, param) - implicit none - class(fraggle_snapshot), intent(in) :: self !! Swiftest encounter structure - class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine fraggle_io_write_frame - module subroutine fraggle_io_log_regime(colliders, fragments) + module subroutine fraggle_io_log_regime(impactors, fragments) implicit none - class(fraggle_colliders), intent(in) :: colliders + class(collision_impactors), intent(in) :: impactors class(fraggle_fragments), intent(in) :: fragments end subroutine fraggle_io_log_regime - !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class - module subroutine fraggle_placeholder_accel(self, system, param, t, lbeg) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters - implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step - end subroutine fraggle_placeholder_accel - - module subroutine fraggle_placeholder_kick(self, system, param, t, dt, lbeg) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters - implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Stepsize - logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. - end subroutine fraggle_placeholder_kick - - module subroutine fraggle_placeholder_step(self, system, param, t, dt) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters - implicit none - class(fraggle_fragments), intent(inout) :: self !! Helio massive body particle object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Stepsiz - end subroutine fraggle_placeholder_step - - module subroutine fraggle_regime_colliders(self, fragments, system, param) - implicit none - class(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object - class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object - class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - end subroutine fraggle_regime_colliders - module subroutine fraggle_set_budgets_fragments(self) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object end subroutine fraggle_set_budgets_fragments - module subroutine fraggle_set_coordinate_system(self, colliders) - implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object - end subroutine fraggle_set_coordinate_system - - module subroutine fraggle_set_mass_dist_fragments(self, colliders, param) + module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine fraggle_set_mass_dist_fragments - module subroutine fraggle_set_natural_scale_factors(self, colliders) + module subroutine fraggle_set_natural_scale_factors(self, impactors) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object end subroutine fraggle_set_natural_scale_factors - module subroutine fraggle_set_original_scale_factors(self, colliders) + module subroutine fraggle_set_original_scale_factors(self, impactors) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(inout) :: impactors !! Fraggle collider system object end subroutine fraggle_set_original_scale_factors module subroutine fraggle_setup_fragments(self, n, param) implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - integer(I4B), intent(in) :: n !! Number of fragments - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(fraggle_fragments), intent(inout) :: self + integer(I4B), intent(in) :: n + class(swiftest_parameters), intent(in) :: param end subroutine fraggle_setup_fragments module subroutine fraggle_setup_reset_fragments(self) @@ -255,11 +108,11 @@ module subroutine fraggle_setup_reset_fragments(self) class(fraggle_fragments), intent(inout) :: self end subroutine fraggle_setup_reset_fragments - module subroutine fraggle_util_add_fragments_to_system(fragments, colliders, system, param) + module subroutine fraggle_util_add_fragments_to_system(fragments, impactors, system, param) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none class(fraggle_fragments), intent(in) :: fragments !! Fraggle fragment system object - class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine fraggle_util_add_fragments_to_system @@ -279,41 +132,21 @@ module subroutine fraggle_util_construct_temporary_system(fragments, system, par class(swiftest_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters end subroutine fraggle_util_construct_temporary_system - module subroutine fraggle_util_final_colliders(self) + module subroutine fraggle_util_final_impactors(self) implicit none - type(fraggle_colliders), intent(inout) :: self !! Fraggle colliders object - end subroutine fraggle_util_final_colliders + type(collision_impactors), intent(inout) :: self !! Fraggle impactors object + end subroutine fraggle_util_final_impactors module subroutine fraggle_util_final_fragments(self) implicit none type(fraggle_fragments), intent(inout) :: self !! Fraggle frgments object end subroutine fraggle_util_final_fragments - module subroutine fraggle_util_final_snapshot(self) - implicit none - type(fraggle_snapshot), intent(inout) :: self !! Fraggle storage snapshot object - end subroutine fraggle_util_final_snapshot - - module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters - implicit none - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters - logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa - end subroutine fraggle_util_get_energy_momentum - - module subroutine fraggle_util_get_idvalues_snapshot(self, idvals) - implicit none - class(fraggle_snapshot), intent(in) :: self !! Fraggle snapshot object - integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot - end subroutine fraggle_util_get_idvalues_snapshot - module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_start) + module subroutine fraggle_util_restructure(self, impactors, try, f_spin, r_max_start) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - class(fraggle_colliders), intent(in) :: colliders !! Fraggle collider system object + class(collision_impactors), intent(in) :: impactors !! Fraggle collider system object integer(I4B), intent(in) :: try !! The current number of times Fraggle has tried to find a solution real(DP), intent(inout) :: f_spin !! Fraction of energy/momentum that goes into spin. This decreases ater a failed attempt real(DP), intent(inout) :: r_max_start !! The maximum radial distance that the position calculation starts with. This increases after a failed attempt diff --git a/src/modules/swiftest.f90 b/src/modules/swiftest.f90 index edc41f134..8ca51ffdb 100644 --- a/src/modules/swiftest.f90 +++ b/src/modules/swiftest.f90 @@ -14,15 +14,16 @@ module swiftest !! 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. use swiftest_globals use swiftest_operators + use lambda_function use swiftest_classes use whm_classes use rmvs_classes use helio_classes use symba_classes + use encounter_classes + use collision_classes use fraggle_classes - use lambda_function use walltime_classes - use encounter_classes use io_progress_bar !use advisor_annotate !$ use omp_lib diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a3bf66ad7..090446da8 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -347,8 +347,8 @@ module swiftest_classes logical, dimension(:), allocatable :: ldiscard !! Body should be discarded logical, dimension(:), allocatable :: lmask !! Logical mask used to select a subset of bodies when performing certain operations (drift, kick, accel, etc.) real(DP), dimension(:), allocatable :: mu !! G * (Mcb + [m]) - real(DP), dimension(:,:), allocatable :: rh !! Swiftestcentric position - real(DP), dimension(:,:), allocatable :: vh !! Swiftestcentric velocity + real(DP), dimension(:,:), allocatable :: rh !! Heliocentric position + real(DP), dimension(:,:), allocatable :: vh !! Heliocentric velocity real(DP), dimension(:,:), allocatable :: rb !! Barycentric position real(DP), dimension(:,:), allocatable :: vb !! Barycentric velocity real(DP), dimension(:,:), allocatable :: ah !! Total heliocentric acceleration @@ -550,7 +550,6 @@ module swiftest_classes procedure :: read_in => io_read_in_system !! Reads the initial conditions for an nbody system procedure :: read_particle_info => netcdf_read_particle_info_system !! Read in particle metadata from file procedure :: obl_pot => obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body - procedure :: finalize => setup_finalize_system !! Runs any finalization subroutines when ending the simulation. procedure :: initialize => setup_initialize_system !! Initialize the system from input files procedure :: init_particle_info => setup_initialize_particle_info_system !! Initialize the system from input files ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. @@ -1173,12 +1172,6 @@ module subroutine setup_construct_system(system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine setup_construct_system - module subroutine setup_finalize_system(self, param) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine setup_finalize_system - module subroutine setup_initialize_particle_info_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index d9590b59e..fb669b559 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -95,14 +95,6 @@ module swiftest_globals integer(I4B), parameter :: NEW_PARTICLE = -15 integer(I4B), parameter :: OLD_PARTICLE = -16 - !>Symbolic names for collisional outcomes from collresolve_resolve: - integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE = 1 - integer(I4B), parameter :: COLLRESOLVE_REGIME_DISRUPTION = 2 - integer(I4B), parameter :: COLLRESOLVE_REGIME_SUPERCATASTROPHIC = 3 - integer(I4B), parameter :: COLLRESOLVE_REGIME_GRAZE_AND_MERGE = 4 - integer(I4B), parameter :: COLLRESOLVE_REGIME_HIT_AND_RUN = 5 - character(len=*),dimension(5), parameter :: REGIME_NAMES = ["Merge", "Disruption", "Supercatastrophic", "Graze and Merge", "Hit and Run"] - !> String labels for body/particle addition/subtraction in discard file character(*), parameter :: ADD = '+1' character(*), parameter :: SUB = '-1' diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 507b16745..9e710c8fe 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -15,7 +15,7 @@ module symba_classes use swiftest_globals use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_particle_info, swiftest_storage, netcdf_parameters use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system - use fraggle_classes, only : fraggle_colliders, fraggle_fragments + use fraggle_classes, only : collision_impactors, fraggle_fragments use encounter_classes, only : encounter_list, encounter_storage, collision_storage implicit none public @@ -84,7 +84,7 @@ module symba_classes real(DP), dimension(:), allocatable :: atp !! semimajor axis following perihelion passage type(symba_kinship), dimension(:), allocatable :: kin !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step contains - procedure :: make_colliders => symba_collision_make_colliders_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family + procedure :: make_impactors => symba_collision_make_impactors_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family procedure :: flatten => symba_util_flatten_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: discard => symba_discard_pl !! Process massive body discards procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level @@ -189,7 +189,7 @@ module symba_classes class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step integer(I4B) :: irec !! System recursion level - class(fraggle_colliders), allocatable :: colliders !! Fraggle colliders object + class(collision_impactors), allocatable :: impactors !! Fraggle impactors object class(fraggle_fragments), allocatable :: fragments !! Fraggle fragmentation system object contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA @@ -224,11 +224,11 @@ module subroutine symba_collision_extract_collisions_from_encounters(self, syste class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine - module subroutine symba_collision_make_colliders_pl(self,idx) + module subroutine symba_collision_make_impactors_pl(self,idx) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object 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 + end subroutine symba_collision_make_impactors_pl module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec) implicit none diff --git a/src/operators/operator_unit.f90 b/src/operators/operator_unit.f90 index e9c68f28f..8ba5d89e5 100644 --- a/src/operators/operator_unit.f90 +++ b/src/operators/operator_unit.f90 @@ -20,12 +20,16 @@ pure module function operator_unit_sp(A) result(B) implicit none ! Arguments real(SP), dimension(:), intent(in) :: A - real(SP), dimension(NDIM) :: B + real(SP), dimension(NDIM) :: B ! Internals real(SP) :: Amag Amag = norm2(A(:)) - B(:) = A(:) / Amag + if (Amag > tiny(1._SP)) then + B(:) = A(:) / Amag + else + B(:) = 0.0_SP + end if return end function operator_unit_sp @@ -40,7 +44,11 @@ pure module function operator_unit_dp(A) result(B) real(DP) :: Amag Amag = norm2(A(:)) - B(:) = A(:) / Amag + if (Amag > tiny(1._DP)) then + B(:) = A(:) / Amag + else + B(:) = 0.0_DP + end if return end function operator_unit_dp @@ -55,7 +63,11 @@ pure module function operator_unit_qp(A) result(B) real(QP) :: Amag Amag = norm2(A(:)) - B(:) = A(:) / Amag + if (Amag > tiny(1._QP)) then + B(:) = A(:) / Amag + else + B(:) = 0.0_QP + end if return end function operator_unit_qp @@ -67,7 +79,6 @@ pure module function operator_unit_el_sp(A) result(B) real(SP), dimension(:,:), intent(in) :: A real(SP), dimension(:,:), allocatable :: B ! Internals - real(SP) :: Amag integer(I4B) :: i,n n = size(A, 2) @@ -75,8 +86,7 @@ pure module function operator_unit_el_sp(A) result(B) allocate(B(NDIM,n)) do concurrent (i=1:n) - Amag = norm2(A(:, i)) - B(:,i) = A(:,i) / Amag + B(:,i) = operator_unit_sp(A(:,i)) end do return @@ -89,7 +99,6 @@ pure module function operator_unit_el_dp(A) result(B) real(DP), dimension(:,:), intent(in) :: A real(DP), dimension(:,:), allocatable :: B ! Internals - real(DP) :: Amag integer(I4B) :: i,n n = size(A, 2) @@ -97,8 +106,7 @@ pure module function operator_unit_el_dp(A) result(B) allocate(B(NDIM,n)) do concurrent (i=1:n) - Amag = norm2(A(:, i)) - B(:,i) = A(:,i) / Amag + B(:,i) = operator_unit_dp(A(:,i)) end do return @@ -110,7 +118,6 @@ pure module function operator_unit_el_qp(A) result(B) real(QP), dimension(:,:), intent(in) :: A real(QP), dimension(:,:), allocatable :: B ! Internals - real(QP) :: Amag integer(I4B) :: i,n n = size(A, 2) @@ -118,8 +125,7 @@ pure module function operator_unit_el_qp(A) result(B) allocate(B(NDIM,n)) do concurrent (i=1:n) - Amag = norm2(A(:, i)) - B(:,i) = A(:,i) / Amag + B(:,i) = operator_unit_qp(A(:,i)) end do return diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index c9ff0dc7d..353f60fb2 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -89,10 +89,10 @@ module subroutine setup_construct_system(system, param) allocate(collision_storage :: param%collision_history) associate (collision_history => param%collision_history) - allocate(fraggle_io_parameters :: collision_history%nc) + allocate(collision_io_parameters :: collision_history%nc) call collision_history%reset() select type(nc => collision_history%nc) - class is (fraggle_io_parameters) + class is (collision_io_parameters) nc%file_number = param%iloop / param%dump_cadence end select end associate @@ -106,32 +106,10 @@ module subroutine setup_construct_system(system, param) call util_exit(FAILURE) end select - - - - return end subroutine setup_construct_system - module subroutine setup_finalize_system(self, param) - !! author: David A. Minton - !! - !! Runs any finalization subroutines when ending the simulation. - !! - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - associate(system => self) - call param%system_history%nc%close() - end associate - - return - end subroutine setup_finalize_system - - module subroutine setup_initialize_particle_info_system(self, param) !! author: David A. Minton !! diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index e9c619d02..90d3847c4 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -30,7 +30,7 @@ module function symba_collision_casedisruption(system, param, t) result(status) character(len=STRMAX) :: message real(DP) :: dpe - associate(colliders => system%colliders, fragments => system%fragments) + associate(impactors => system%impactors, fragments => system%fragments) select case(fragments%regime) case(COLLRESOLVE_REGIME_DISRUPTION) @@ -38,14 +38,14 @@ module function symba_collision_casedisruption(system, param, t) result(status) case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) message = "Supercatastrophic disruption between" end select - call symba_collision_collider_message(system%pl, colliders%idx, message) + call symba_collision_collider_message(system%pl, impactors%idx, message) call io_log_one_message(FRAGGLE_LOG_OUT, message) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - call fragments%set_mass_dist(colliders, param) + call fragments%set_mass_dist(impactors, param) ! Generate the position and velocity distributions of the fragments - call fragments%generate_fragments(colliders, system, param, lfailure) + call fragments%generate_fragments(impactors, system, param, lfailure) dpe = fragments%pe_after - fragments%pe_before system%Ecollisions = system%Ecollisions - dpe @@ -57,9 +57,9 @@ module function symba_collision_casedisruption(system, param, t) result(status) nfrag = 0 select type(pl => system%pl) class is (symba_pl) - pl%status(colliders%idx(:)) = status - pl%ldiscard(colliders%idx(:)) = .false. - pl%lcollision(colliders%idx(:)) = .false. + pl%status(impactors%idx(:)) = status + pl%ldiscard(impactors%idx(:)) = .false. + pl%lcollision(impactors%idx(:)) = .false. end select else ! Populate the list of new bodies @@ -69,7 +69,7 @@ module function symba_collision_casedisruption(system, param, t) result(status) select case(fragments%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTION - ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) + ibiggest = impactors%idx(maxloc(system%pl%Gmass(impactors%idx(:)), dim=1)) fragments%id(1) = system%pl%id(ibiggest) fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] param%maxid = fragments%id(nfrag) @@ -105,12 +105,12 @@ module function symba_collision_casehitandrun(system, param, t) result(status) character(len=STRMAX) :: message real(DP) :: dpe - associate(colliders => system%colliders, fragments => system%fragments) + associate(impactors => system%impactors, fragments => system%fragments) message = "Hit and run between" - call symba_collision_collider_message(system%pl, colliders%idx, message) + call symba_collision_collider_message(system%pl, impactors%idx, message) call io_log_one_message(FRAGGLE_LOG_OUT, trim(adjustl(message))) - if (colliders%mass(1) > colliders%mass(2)) then + if (impactors%mass(1) > impactors%mass(2)) then jtarg = 1 jproj = 2 else @@ -118,16 +118,16 @@ module function symba_collision_casehitandrun(system, param, t) result(status) jproj = 1 end if - if (fragments%mass_dist(2) > 0.9_DP * colliders%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + if (fragments%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched call io_log_one_message(FRAGGLE_LOG_OUT, "Pure hit and run. No new fragments generated.") nfrag = 0 lpure = .true. else ! Imperfect hit and run, so we'll keep the largest body and destroy the other lpure = .false. - call fragments%set_mass_dist(colliders, param) + call fragments%set_mass_dist(impactors, param) ! Generate the position and velocity distributions of the fragments - call fragments%generate_fragments(colliders, system, param, lpure) + call fragments%generate_fragments(impactors, system, param, lpure) dpe = fragments%pe_after - fragments%pe_before system%Ecollisions = system%Ecollisions - dpe @@ -146,13 +146,13 @@ module function symba_collision_casehitandrun(system, param, t) result(status) status = HIT_AND_RUN_PURE select type(pl => system%pl) class is (symba_pl) - pl%status(colliders%idx(:)) = ACTIVE - pl%ldiscard(colliders%idx(:)) = .false. - pl%lcollision(colliders%idx(:)) = .false. + pl%status(impactors%idx(:)) = ACTIVE + pl%ldiscard(impactors%idx(:)) = .false. + pl%lcollision(impactors%idx(:)) = .false. end select - allocate(system%fragments%pl, source=system%colliders%pl) ! Be sure to save the pl so that snapshots still work + allocate(system%fragments%pl, source=system%impactors%pl) ! Be sure to save the pl so that snapshots still work else - ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) + ibiggest = impactors%idx(maxloc(system%pl%Gmass(impactors%idx(:)), dim=1)) fragments%id(1) = system%pl%id(ibiggest) fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] param%maxid = fragments%id(nfrag) @@ -187,37 +187,37 @@ module function symba_collision_casemerge(system, param, t) result(status) real(DP) :: dpe character(len=STRMAX) :: message - associate(colliders => system%colliders, fragments => system%fragments) + associate(impactors => system%impactors, fragments => system%fragments) message = "Merging" - call symba_collision_collider_message(system%pl, colliders%idx, message) + call symba_collision_collider_message(system%pl, impactors%idx, message) call io_log_one_message(FRAGGLE_LOG_OUT, message) select type(pl => system%pl) class is (symba_pl) - call fragments%set_mass_dist(colliders, param) + call fragments%set_mass_dist(impactors, param) ! Calculate the initial energy of the system without the collisional family - call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.true.) + call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.true.) - ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) + ibiggest = impactors%idx(maxloc(pl%Gmass(impactors%idx(:)), dim=1)) fragments%id(1) = pl%id(ibiggest) fragments%rb(:,1) = fragments%rbcom(:) fragments%vb(:,1) = fragments%vbcom(:) if (param%lrotation) then ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) + L_spin_new(:) = impactors%L_orbit(:,1) + impactors%L_orbit(:,2) + impactors%L_spin(:,1) + impactors%L_spin(:,2) ! Assume prinicpal axis rotation on 3rd Ip axis fragments%rot(:,1) = L_spin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2) else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - system%Lescape(:) = system%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + system%Lescape(:) = system%Lescape(:) + impactors%L_orbit(:,1) + impactors%L_orbit(:,2) end if - ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping + ! Keep track of the component of potential energy due to the pre-impact impactors%idx for book-keeping ! Get the energy of the system after the collision - call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.false.) dpe = fragments%pe_after - fragments%pe_before system%Ecollisions = system%Ecollisions - dpe system%Euntracked = system%Euntracked + dpe @@ -225,8 +225,8 @@ module function symba_collision_casemerge(system, param, t) result(status) ! Update any encounter lists that have the removed bodies in them so that they instead point to the new do k = 1, system%plplenc_list%nenc - do j = 1, colliders%ncoll - i = colliders%idx(j) + do j = 1, impactors%ncoll + i = impactors%idx(j) if (i == ibiggest) cycle if (system%plplenc_list%id1(k) == pl%id(i)) then system%plplenc_list%id1(k) = pl%id(ibiggest) @@ -258,7 +258,7 @@ subroutine symba_collision_collider_message(pl, collidx, collider_message) implicit none ! Arguments class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - integer(I4B), dimension(:), intent(in) :: collidx !! Index of collisional colliders%idx members + integer(I4B), dimension(:), intent(in) :: collidx !! Index of collisional impactors%idx members character(*), intent(inout) :: collider_message !! The message to print to the screen. ! Internals integer(I4B) :: i, n @@ -375,7 +375,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir 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 - if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) + 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 pl%lcollision([i, j]) = .true. @@ -467,10 +467,10 @@ pure elemental subroutine symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, G end subroutine symba_collision_check_one - function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) result(lflag) + function symba_collision_consolidate_impactors(pl, cb, param, idx_parent, impactors) result(lflag) !! author: David A. Minton !! - !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all colliders%idx members, + !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all impactors%idx 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 @@ -478,9 +478,9 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid class(symba_cb), intent(inout) :: cb !! SyMBA central body object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - class(fraggle_colliders), intent(out) :: colliders + class(collision_impactors), intent(out) :: impactors ! Result - logical :: lflag !! Logical flag indicating whether a colliders%idx was successfully created or not + logical :: lflag !! Logical flag indicating whether a impactors%idx was successfully created or not ! Internals type collidx_array integer(I4B), dimension(:), allocatable :: id @@ -488,7 +488,7 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid end type collidx_array type(collidx_array), dimension(2) :: parent_child_index_array integer(I4B), dimension(2) :: nchild - integer(I4B) :: i, j, ncolliders, idx_child + integer(I4B) :: i, j, nimpactors, idx_child real(DP), dimension(2) :: volume, density real(DP) :: mchild, volchild real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv @@ -510,9 +510,9 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid pl%kin(idx_parent(2))%parent = idx_parent(1) end if - colliders%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si - colliders%radius(:) = pl%radius(idx_parent(:)) - volume(:) = (4.0_DP / 3.0_DP) * PI * colliders%radius(:)**3 + impactors%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si + impactors%radius(:) = pl%radius(idx_parent(:)) + volume(:) = (4.0_DP / 3.0_DP) * PI * impactors%radius(:)**3 ! Group together the ids and indexes of each collisional parent and its children do j = 1, 2 @@ -528,24 +528,24 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid end associate end do - ! Consolidate the groups of collsional parents with any children they may have into a single "colliders%idx" index array - ncolliders = 2 + sum(nchild(:)) - allocate(colliders%idx(ncolliders)) - colliders%idx = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] + ! Consolidate the groups of collsional parents with any children they may have into a single "impactors%idx" index array + nimpactors = 2 + sum(nchild(:)) + allocate(impactors%idx(nimpactors)) + impactors%idx = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] - colliders%ncoll = count(pl%lcollision(colliders%idx(:))) - colliders%idx = pack(colliders%idx(:), pl%lcollision(colliders%idx(:))) - colliders%L_spin(:,:) = 0.0_DP - colliders%Ip(:,:) = 0.0_DP + impactors%ncoll = count(pl%lcollision(impactors%idx(:))) + impactors%idx = pack(impactors%idx(:), pl%lcollision(impactors%idx(:))) + impactors%L_spin(:,:) = 0.0_DP + impactors%Ip(:,:) = 0.0_DP ! Find the barycenter of each body along with its children, if it has any do j = 1, 2 - colliders%rb(:, j) = pl%rh(:, idx_parent(j)) + cb%rb(:) - colliders%vb(:, j) = pl%vb(:, idx_parent(j)) + impactors%rb(:, j) = pl%rh(:, idx_parent(j)) + cb%rb(:) + impactors%vb(:, j) = pl%vb(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then - colliders%Ip(:, j) = colliders%mass(j) * pl%Ip(:, idx_parent(j)) - colliders%L_spin(:, j) = colliders%Ip(3, j) * colliders%radius(j)**2 * pl%rot(:, idx_parent(j)) + impactors%Ip(:, j) = impactors%mass(j) * pl%Ip(:, idx_parent(j)) + impactors%L_spin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) end if if (nchild(j) > 0) then @@ -560,49 +560,49 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid ! Get angular momentum of the child-parent pair and add that to the spin ! Add the child's spin if (param%lrotation) then - xcom(:) = (colliders%mass(j) * colliders%rb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) - vcom(:) = (colliders%mass(j) * colliders%vb(:,j) + mchild * vchild(:)) / (colliders%mass(j) + mchild) - xc(:) = colliders%rb(:, j) - xcom(:) - vc(:) = colliders%vb(:, j) - vcom(:) + xcom(:) = (impactors%mass(j) * impactors%rb(:,j) + mchild * xchild(:)) / (impactors%mass(j) + mchild) + vcom(:) = (impactors%mass(j) * impactors%vb(:,j) + mchild * vchild(:)) / (impactors%mass(j) + mchild) + xc(:) = impactors%rb(:, j) - xcom(:) + vc(:) = impactors%vb(:, j) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + colliders%mass(j) * xcrossv(:) + impactors%L_spin(:, j) = impactors%L_spin(:, j) + impactors%mass(j) * xcrossv(:) xc(:) = xchild(:) - xcom(:) vc(:) = vchild(:) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) + impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * xcrossv(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & + impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & * pl%radius(idx_child)**2 & * pl%rot(:, idx_child) - colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) + impactors%Ip(:, j) = impactors%Ip(:, j) + mchild * pl%Ip(:, idx_child) end if ! Merge the child and parent - colliders%mass(j) = colliders%mass(j) + mchild - colliders%rb(:, j) = xcom(:) - colliders%vb(:, j) = vcom(:) + impactors%mass(j) = impactors%mass(j) + mchild + impactors%rb(:, j) = xcom(:) + impactors%vb(:, j) = vcom(:) end do end if - density(j) = colliders%mass(j) / volume(j) - colliders%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) - if (param%lrotation) colliders%Ip(:, j) = colliders%Ip(:, j) / colliders%mass(j) + density(j) = impactors%mass(j) / volume(j) + impactors%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) + if (param%lrotation) impactors%Ip(:, j) = impactors%Ip(:, j) / impactors%mass(j) end do lflag = .true. - xcom(:) = (colliders%mass(1) * colliders%rb(:, 1) + colliders%mass(2) * colliders%rb(:, 2)) / sum(colliders%mass(:)) - vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) - mxc(:, 1) = colliders%mass(1) * (colliders%rb(:, 1) - xcom(:)) - mxc(:, 2) = colliders%mass(2) * (colliders%rb(:, 2) - xcom(:)) - vcc(:, 1) = colliders%vb(:, 1) - vcom(:) - vcc(:, 2) = colliders%vb(:, 2) - vcom(:) - colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) + xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / sum(impactors%mass(:)) + vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / sum(impactors%mass(:)) + mxc(:, 1) = impactors%mass(1) * (impactors%rb(:, 1) - xcom(:)) + mxc(:, 2) = impactors%mass(2) * (impactors%rb(:, 2) - xcom(:)) + vcc(:, 1) = impactors%vb(:, 1) - vcom(:) + vcc(:, 2) = impactors%vb(:, 2) - vcom(:) + impactors%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) - ! Destroy the kinship relationships for all members of this colliders%idx - call pl%reset_kinship(colliders%idx(:)) + ! Destroy the kinship relationships for all members of this impactors%idx + call pl%reset_kinship(impactors%idx(:)) return - end function symba_collision_consolidate_colliders + end function symba_collision_consolidate_impactors module subroutine symba_collision_extract_collisions_from_encounters(self, system, param) @@ -674,10 +674,10 @@ module subroutine symba_collision_extract_collisions_from_encounters(self, syste end subroutine symba_collision_extract_collisions_from_encounters - module subroutine symba_collision_make_colliders_pl(self, idx) + module subroutine symba_collision_make_impactors_pl(self, idx) !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton !! - !! When a single body is involved in more than one collision in a single step, it becomes part of a colliders%idx. + !! When a single body is involved in more than one collision in a single step, it becomes part of a impactors%idx. !! The largest body involved in a multi-body collision is the "parent" and all bodies that collide with it are its "children," !! including those that collide with the children. !! @@ -732,7 +732,7 @@ module subroutine symba_collision_make_colliders_pl(self, idx) end associate return - end subroutine symba_collision_make_colliders_pl + end subroutine symba_collision_make_impactors_pl subroutine symba_collision_mergeaddsub(system, param, t, status) @@ -747,7 +747,7 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) real(DP), intent(in) :: t !! Time of collision integer(I4B), intent(in) :: status !! Status flag to assign to adds ! Internals - integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag + integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, nimpactors, nfrag logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew, plsub character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' @@ -757,9 +757,9 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) class is (symba_pl) select type(pl_discards => system%pl_discards) class is (symba_merger) - associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody, colliders => system%colliders, fragments => system%fragments) - ! Add the colliders%idx bodies to the subtraction list - ncolliders = colliders%ncoll + associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody, impactors => system%impactors, fragments => system%fragments) + ! Add the impactors%idx bodies to the subtraction list + nimpactors = impactors%ncoll nfrag = fragments%nbody param%maxid_collision = max(param%maxid_collision, maxval(system%pl%info(:)%collision_id)) @@ -768,8 +768,8 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) ! Setup new bodies allocate(plnew, mold=pl) call plnew%setup(nfrag, param) - ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) - ismallest = colliders%idx(minloc(pl%Gmass(colliders%idx(:)), dim=1)) + ibiggest = impactors%idx(maxloc(pl%Gmass(impactors%idx(:)), dim=1)) + ismallest = impactors%idx(minloc(pl%Gmass(impactors%idx(:)), dim=1)) ! Copy over identification, information, and physical properties of the new bodies from the fragment list plnew%id(1:nfrag) = fragments%id(1:nfrag) @@ -796,13 +796,13 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) then + do i = 1, nimpactors + if (impactors%idx(i) == ibiggest) then iother = ismallest else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=t, & + call pl%info(impactors%idx(i))%set_value(status="Supercatastrophic", discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -820,21 +820,21 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) cycle + do i = 1, nimpactors + if (impactors%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=t, & + call pl%info(impactors%idx(i))%set_value(status=origin_type, discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do case(MERGED) call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE - do i = 1, ncolliders - if (colliders%idx(i) == ibiggest) cycle + do i = 1, nimpactors + if (impactors%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=t, discard_rh=pl%rh(:,i), & + call pl%info(impactors%idx(i))%set_value(status="MERGED", discard_time=t, discard_rh=pl%rh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select @@ -874,11 +874,11 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) pl_adds%ncomp(nstart:nend) = plnew%nbody ! Add the discarded bodies to the discard list - pl%status(colliders%idx(:)) = MERGED - pl%ldiscard(colliders%idx(:)) = .true. - pl%lcollision(colliders%idx(:)) = .true. + pl%status(impactors%idx(:)) = MERGED + pl%ldiscard(impactors%idx(:)) = .true. + pl%lcollision(impactors%idx(:)) = .true. lmask(:) = .false. - lmask(colliders%idx(:)) = .true. + lmask(impactors%idx(:)) = .true. call plnew%setup(0, param) deallocate(plnew) @@ -887,11 +887,11 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) call pl%spill(plsub, lmask, ldestructive=.false.) nstart = pl_discards%nbody + 1 - nend = pl_discards%nbody + ncolliders - call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, ncolliders)]) + nend = pl_discards%nbody + nimpactors + call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)]) ! Record how many bodies were subtracted in this event - pl_discards%ncomp(nstart:nend) = ncolliders + pl_discards%ncomp(nstart:nend) = nimpactors call plsub%setup(0, param) deallocate(plsub) @@ -926,24 +926,24 @@ subroutine symba_resolve_collision(plplcollision_list , system, param, t) select type (cb => system%cb) class is (symba_cb) do i = 1, ncollisions - allocate(fraggle_colliders :: system%colliders) + allocate(collision_impactors :: system%impactors) allocate(fraggle_fragments :: system%fragments) 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, system%colliders) + lgoodcollision = symba_collision_consolidate_impactors(pl, cb, param, idx_parent, system%impactors) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle if (param%lfragmentation) then - call system%colliders%regime(system%fragments, system, param) + call system%impactors%regime(system%fragments, system, param) else - associate(fragments => system%fragments, colliders => system%colliders) + associate(fragments => system%fragments, impactors => system%impactors) fragments%regime = COLLRESOLVE_REGIME_MERGE - fragments%mtot = sum(colliders%mass(:)) + fragments%mtot = sum(impactors%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 + fragments%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot + fragments%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot end associate end if @@ -960,7 +960,7 @@ subroutine symba_resolve_collision(plplcollision_list , system, param, t) call util_exit(FAILURE) end select call collision_history%take_snapshot(param,system, t, "after") - deallocate(system%colliders,system%fragments) + deallocate(system%impactors,system%fragments) end do end select end select