From e031a2e00a4bca04025d5c60d2bf90be5250a281 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 16 Dec 2022 15:07:02 -0500 Subject: [PATCH] Cleaning up after the restructuring --- src/CMakeLists.txt | 1 - src/collision/collision_io.f90 | 147 ++++++++++------- src/collision/collision_placeholder.f90 | 53 ------ src/collision/collision_setup.f90 | 35 +++- src/collision/collision_util.f90 | 207 ++++++++++++++---------- src/encounter/encounter_io.f90 | 98 ++++------- src/encounter/encounter_placeholder.f90 | 53 ------ src/fraggle/fraggle_generate.f90 | 25 ++- src/fraggle/fraggle_io.f90 | 4 +- src/fraggle/fraggle_regime.f90 | 4 +- src/fraggle/fraggle_set.f90 | 34 ++-- src/fraggle/fraggle_setup.f90 | 9 +- src/fraggle/fraggle_util.f90 | 54 +++++-- src/modules/collision_classes.f90 | 76 ++++----- src/modules/encounter_classes.f90 | 6 +- src/modules/fraggle_classes.f90 | 41 +++-- src/modules/swiftest_classes.f90 | 4 +- src/modules/symba_classes.f90 | 3 +- src/netcdf/netcdf.f90 | 12 +- src/symba/symba_collision.f90 | 30 ++-- 20 files changed, 436 insertions(+), 460 deletions(-) delete mode 100644 src/collision/collision_placeholder.f90 delete mode 100644 src/encounter/encounter_placeholder.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5ada15b9b..90ace7141 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,7 +27,6 @@ SET(FAST_MATH_FILES ${SRC}/modules/whm_classes.f90 ${SRC}/modules/swiftest.f90 ${SRC}/collision/collision_io.f90 - ${SRC}/collision/collision_placeholder.f90 ${SRC}/collision/collision_setup.f90 ${SRC}/collision/collision_util.f90 ${SRC}/discard/discard.f90 diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 88c37cbbc..e22937577 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -12,7 +12,51 @@ contains - module subroutine collision_io_initialize(self, param) + + +module subroutine collision_io_dump(self, param) + !! author: David A. Minton + !! + !! Dumps the time history of an encounter to file. + implicit none + ! Arguments + class(collision_storage(*)), intent(inout) :: self !! Encounter storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i + + select type(nc => self%nc) + class is (collision_io_parameters) + if (self%iframe > 0) then + nc%file_number = nc%file_number + 1 + call self%make_index_map() + nc%event_dimsize = self%nt + nc%name_dimsize = self%nid + + write(nc%file_name, '("collision_",I0.6,".nc")') nc%file_number + call nc%initialize(param) + + do i = 1, self%nframes + if (allocated(self%frame(i)%item)) then + select type(snapshot => self%frame(i)%item) + class is (collision_snapshot) + param%ioutput = i + call snapshot%write_frame(nc,param) + end select + else + exit + end if + end do + + call nc%close() + call self%reset() + end if + end select + + return +end subroutine collision_io_dump + + module subroutine collision_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. @@ -51,93 +95,93 @@ module subroutine collision_io_initialize(self, param) close(unit=LUN, status="delete") end if - call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "collision_io_initialize nf90_create" ) + call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "collision_io_initialize_output 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" + call check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "collision_io_initialize_output 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_output 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_output 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_output 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_output 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" ) + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_initialize_output 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_output 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_output 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%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "collision_io_initialize_output 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" ) + nc%event_dimid, nc%time_varid), "collision_io_initialize_output 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") + [nc%str_dimid, nc%event_dimid], nc%regime_varid), "collision_io_initialize_output 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") + [ nc%event_dimid], nc%Qloss_varid), "collision_io_initialize_output 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") + [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "collision_io_initialize_output 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") + [ nc%event_dimid], nc%loop_varid), "collision_io_initialize_output 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") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "collision_io_initialize_output 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") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "collision_io_initialize_output 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") + [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "collision_io_initialize_output 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") + [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "collision_io_initialize_output 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") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "collision_io_initialize_output 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") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "collision_io_initialize_output 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") + [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "collision_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), "collision_io_initialize nf90_def_var KE_spin_varid" ) + [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "collision_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), "collision_io_initialize nf90_def_var PE_varid" ) + [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "collision_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), "collision_io_initialize nf90_def_var L_orb_varid" ) + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "collision_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), "collision_io_initialize nf90_def_var L_spin_varid" ) + call check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type,& + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%Lspin_varid), "collision_io_initialize_output nf90_def_var Lspin_varid" ) end if - call check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_initialize nf90_inquire nVariables" ) + call check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_initialize_output 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" ) + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "collision_io_initialize_output 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" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "collision_io_initialize_output 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" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "collision_io_initialize_output 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" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "collision_io_initialize_output 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" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "collision_io_initialize_output 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" ) + call check( nf90_enddef(nc%id), "collision_io_initialize_output 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" ) + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "collision_io_initialize_output 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_output 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_output nf90_put_var stage 2" ) end associate end select @@ -147,7 +191,7 @@ module subroutine collision_io_initialize(self, param) 667 continue write(*,*) "Error creating fragmentation output file. " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine collision_io_initialize + end subroutine collision_io_initialize_output module subroutine collision_io_write_frame_snapshot(self, nc, param) @@ -169,23 +213,23 @@ module subroutine collision_io_write_frame_snapshot(self, nc, param) 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) + associate(system => self%collision_system, impactors => self%collision_system%impactors, fragments => self%collision_system%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))) + charstring = trim(adjustl(REGIME_NAMES(impactors%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" ) + call check( nf90_put_var(nc%id, nc%Qloss_varid, impactors%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) + allocate(pl, source=system%before%pl) case(2) - allocate(pl, source=fragments%pl) + allocate(pl, source=system%after%pl) end select npl = pl%nbody do i = 1, npl @@ -204,16 +248,11 @@ module subroutine collision_io_write_frame_snapshot(self, nc, param) 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" ) + call check( nf90_put_var(nc%id, nc%ke_orb_varid, system%ke_orbit(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var ke_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%ke_spin_varid, system%ke_spin(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var ke_spin_varid before" ) + call check( nf90_put_var(nc%id, nc%pe_varid, system%pe(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_write_frame_snapshot nf90_put_var pe_varid before" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, system%Lorbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_write_frame_snapshot nf90_put_var L_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%Lspin_varid, system%Lspin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_write_frame_snapshot nf90_put_var Lspin_varid before" ) end if call check( nf90_set_fill(nc%id, old_mode, old_mode) ) diff --git a/src/collision/collision_placeholder.f90 b/src/collision/collision_placeholder.f90 deleted file mode 100644 index 61ae6eebd..000000000 --- a/src/collision/collision_placeholder.f90 +++ /dev/null @@ -1,53 +0,0 @@ -!! 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/collision/collision_setup.f90 b/src/collision/collision_setup.f90 index 24d565472..287d74109 100644 --- a/src/collision/collision_setup.f90 +++ b/src/collision/collision_setup.f90 @@ -34,22 +34,41 @@ module subroutine collision_setup_fragments(self, n, param) ! Arguments class(collision_fragments), intent(inout) :: self integer(I4B), intent(in) :: n - class(swiftest_parameters), intent(in) :: param + 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() + self%mtot = 0.0_DP + allocate(self%status(n)) + allocate(self%rb(NDIM,n)) + allocate(self%vb(NDIM,n)) + allocate(self%mass(n)) + allocate(self%rot(NDIM,n)) + allocate(self%Ip(NDIM,n)) + + allocate(self%rc(NDIM,n)) + allocate(self%vc(NDIM,n)) + allocate(self%vmag(n)) + allocate(self%rmag(n)) + allocate(self%rotmag(n)) + allocate(self%radius(n)) + allocate(self%density(n)) + + self%status(:) = INACTIVE + self%rb(:,:) = 0.0_DP + self%vb(:,:) = 0.0_DP + self%rc(:,:) = 0.0_DP + self%vc(:,:) = 0.0_DP + self%vmag(:) = 0.0_DP + self%rmag(:) = 0.0_DP + self%rotmag(:) = 0.0_DP + self%radius(:) = 0.0_DP + self%density(:) = 0.0_DP return end subroutine collision_setup_fragments diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index fc0926c91..55a8b23a1 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -7,36 +7,43 @@ !! 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 +submodule (collision_classes) s_collision_util use swiftest contains - module subroutine collision_util_final_impactors(self) + module subroutine collision_util_dealloc_fragments(self) !! author: David A. Minton !! - !! Finalizer will deallocate all allocatables + !! Deallocates all allocatables implicit none ! Arguments - type(collision_impactors), intent(inout) :: self !! Collision impactors storage object + class(collision_fragments), intent(inout) :: self - call self%reset() + call util_dealloc_pl(self) - return - end subroutine collision_util_final_impactors + 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_dealloc_fragments - module subroutine collision_util_final_system(self) + module subroutine collision_util_final_impactors(self) !! author: David A. Minton !! !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(collision_system), intent(inout) :: self !! Collision impactors storage object + type(collision_impactors), intent(inout) :: self !! Collision impactors storage object call self%reset() return - end subroutine collision_util_final_system + end subroutine collision_util_final_impactors module subroutine collision_util_final_snapshot(self) @@ -78,7 +85,7 @@ module subroutine collision_util_final_system(self) call self%reset() return - end subroutine collision_util_final_storage + end subroutine collision_util_final_system module subroutine collision_util_get_idvalues_snapshot(self, idvals) @@ -87,28 +94,41 @@ module subroutine collision_util_get_idvalues_snapshot(self, idvals) !! Returns an array of all id values saved in this snapshot implicit none ! Arguments - class(collision_snapshot), intent(in) :: self !! Fraggle snapshot object + 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 + integer(I4B) :: npl_before, ntp_before, npl_after, ntp_after, ntot, nlo, nhi - if (allocated(self%impactors)) then - ncoll = self%impactors%pl%nbody - else - ncoll = 0 + npl_before = 0; ntp_before = 0; npl_after = 0; ntp_after = 0 + if (allocated(self%collision_system%before%pl)) then + npl_before = self%collision_system%before%pl%nbody + endif + + if (allocated(self%collision_system%before%tp)) then + ntp_before = self%collision_system%before%tp%nbody end if - if (allocated(self%fragments)) then - nfrag = self%fragments%pl%nbody - else - nfrag = 0 + if (allocated(self%collision_system%after%pl)) then + npl_after = self%collision_system%after%pl%nbody end if - if (ncoll + nfrag == 0) return - allocate(idvals(ncoll+nfrag)) + if (allocated(self%collision_system%after%tp)) then + ntp_after = self%collision_system%after%tp%nbody + end if + + ntot = npl_before + ntp_before + npl_after + ntp_after + if (ntot == 0) return + allocate(idvals(ntot)) - if (ncoll > 0) idvals(1:ncoll) = self%impactors%pl%id(:) - if (nfrag > 0) idvals(ncoll+1:ncoll+nfrag) = self%fragments%pl%id(:) + nlo = 1; nhi = npl_before + if (npl_before > 0) idvals(nlo:nhi) = self%collision_system%before%pl%id(1:npl_before) + nlo = nhi + 1; nhi = nhi + ntp_before + if (ntp_before > 0) idvals(nlo:nhi) = self%collision_system%before%tp%id(1:ntp_before) + + nlo = nhi + 1; nhi = nhi + npl_after + if (npl_after > 0) idvals(nlo:nhi) = self%collision_system%after%pl%id(1:npl_after) + nlo = nhi + 1; nhi = nhi + ntp_after + if (ntp_after > 0) idvals(nlo:nhi) = self%collision_system%after%tp%id(1:ntp_after) return @@ -124,14 +144,14 @@ module subroutine collision_util_get_energy_momentum(self, system, param, lbefo !! 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 + 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 + integer(I4B) :: npl_before, npl_after, stage associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody, pl => system%pl, cb => system%cb) @@ -166,22 +186,18 @@ module subroutine collision_util_get_energy_momentum(self, system, param, lbefo ! 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 + stage = 1 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. + stage = 2 end if + self%Lorbit(:,stage) = tmpsys%Lorbit(:) + self%Lspin(:,stage) = tmpsys%Lspin(:) + self%Ltot(:,stage) = tmpsys%Ltot(:) + self%ke_orbit(stage) = tmpsys%ke_orbit + self%ke_spin(stage) = tmpsys%ke_spin + self%pe(stage) = tmpsys%pe + self%Etot(stage) = tmpsys%te + if (stage == 2) self%Etot(stage) = self%Etot(stage) - (self%pe(2) - self%pe(1)) ! Gotta be careful with PE when number of bodies changes. end associate return @@ -211,27 +227,42 @@ module subroutine collision_util_index_map(self) return end subroutine collision_util_index_map - - module subroutine collision_util_reset_fragments(self) - !! author: David A. Minton - !! - !! Deallocates all allocatables + !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class + module subroutine collision_util_placeholder_accel(self, system, param, t, lbeg) implicit none - ! Arguments - class(collision_fragments), intent(inout) :: self + 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_util_placeholder_accel - call self%swiftest_pl%dealloc() + module subroutine collision_util_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_util_placeholder_kick - 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) + module subroutine collision_util_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_util_reset_fragments + end subroutine collision_util_placeholder_step module subroutine collision_util_reset_impactors(self) @@ -244,25 +275,25 @@ module subroutine collision_util_reset_impactors(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 + self%ncoll = 0 + self%rb(:,:) = 0.0_DP + self%vb(:,:) = 0.0_DP + self%rot(:,:) = 0.0_DP + self%Lspin(:,:) = 0.0_DP + self%Lorbit(:,:) = 0.0_DP + self%Ip(:,:) = 0.0_DP + self%mass(:) = 0.0_DP + self%radius(:) = 0.0_DP + self%Qloss = 0.0_DP + self%regime = 0 + + self%x_unit(:) = 0.0_DP + self%y_unit(:) = 0.0_DP + self%z_unit(:) = 0.0_DP + self%v_unit(:) = 0.0_DP + self%rbcom(:) = 0.0_DP + self%vbcom(:) = 0.0_DP + self%rbimp(:) = 0.0_DP return end subroutine collision_util_reset_impactors @@ -281,16 +312,16 @@ module subroutine collision_util_reset_system(self) 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 + self%Lorbit(:,:) = 0.0_DP + self%Lspin(:,:) = 0.0_DP + self%Ltot(:,:) = 0.0_DP + self%ke_orbit(:) = 0.0_DP + self%ke_spin(:) = 0.0_DP + self%pe(:) = 0.0_DP + self%Etot(:) = 0.0_DP return - end subroutine collision_util_reset_impactors + end subroutine collision_util_reset_system subroutine collision_util_save_snapshot(collision_history, snapshot) @@ -338,4 +369,4 @@ subroutine collision_util_save_snapshot(collision_history, snapshot) end subroutine collision_util_save_snapshot -end submodule s_encounter_util \ No newline at end of file +end submodule s_collision_util \ No newline at end of file diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index fb84067ca..70b782948 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -11,50 +11,6 @@ use swiftest contains - - module subroutine collision_io_dump(self, param) - !! author: David A. Minton - !! - !! Dumps the time history of an encounter to file. - implicit none - ! Arguments - class(collision_storage(*)), intent(inout) :: self !! Encounter storage object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i - - select type(nc => self%nc) - class is (collision_io_parameters) - if (self%iframe > 0) then - nc%file_number = nc%file_number + 1 - call self%make_index_map() - nc%event_dimsize = self%nt - nc%name_dimsize = self%nid - - write(nc%file_name, '("collision_",I0.6,".nc")') nc%file_number - call nc%initialize(param) - - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) then - select type(snapshot => self%frame(i)%item) - class is (collision_snapshot) - param%ioutput = i - call snapshot%write_frame(nc,param) - end select - else - exit - end if - end do - - call nc%close() - call self%reset() - end if - end select - - return - end subroutine collision_io_dump - - module subroutine encounter_io_dump(self, param) ! author: David A. Minton !! @@ -98,7 +54,7 @@ module subroutine encounter_io_dump(self, param) end subroutine encounter_io_dump - module subroutine encounter_io_initialize(self, param) + module subroutine encounter_io_initialize_output(self, param) !! author: David A. Minton !! !! Initialize a NetCDF encounter file system. This is a simplified version of the main simulation output NetCDF file, but with fewer variables. @@ -135,54 +91,54 @@ module subroutine encounter_io_initialize(self, param) close(unit=LUN, status="delete") end if - call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "encounter_io_initialize nf90_create" ) + call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "encounter_io_initialize_output nf90_create" ) ! Dimensions - call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize nf90_def_dim time_dimid" ) ! Simulation time dimension - call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "encounter_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), "encounter_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), "encounter_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "encounter_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "encounter_io_initialize_output 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), "encounter_io_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) ! Dimension coordinates - call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize nf90_def_var time_varid" ) - call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_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), "encounter_io_initialize nf90_def_var id_varid" ) + call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize_output nf90_def_var time_varid" ) + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize_output 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), "encounter_io_initialize_output nf90_def_var id_varid" ) ! Variables - call check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "encounter_io_initialize nf90_def_var id_varid" ) - call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%ptype_varid), "encounter_io_initialize nf90_def_var ptype_varid" ) - call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rh_varid), "encounter_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%time_dimid], nc%vh_varid), "encounter_io_initialize nf90_def_var vh_varid" ) - call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize nf90_def_var Gmass_varid" ) - call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_initialize nf90_def_var loop_varid" ) + call check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "encounter_io_initialize_output nf90_def_var id_varid" ) + call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%ptype_varid), "encounter_io_initialize_output nf90_def_var ptype_varid" ) + call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize_output 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%time_dimid], nc%vh_varid), "encounter_io_initialize_output nf90_def_var vh_varid" ) + call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize_output nf90_def_var Gmass_varid" ) + call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_initialize_output nf90_def_var loop_varid" ) if (param%lclose) then - call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize nf90_def_var radius_varid" ) + call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize_output nf90_def_var radius_varid" ) end if if (param%lrotation) then - call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%Ip_varid), "encounter_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%time_dimid], nc%rot_varid), "encounter_io_initialize nf90_def_var rot_varid" ) + call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize_output 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%time_dimid], nc%rot_varid), "encounter_io_initialize_output nf90_def_var rot_varid" ) end if - call check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize nf90_inquire nVariables" ) + call check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize_output nf90_inquire nVariables" ) do varid = 1, nvar - call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable" ) + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize_output nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "encounter_io_initialize_output nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "encounter_io_initialize_output nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "encounter_io_initialize_output nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "encounter_io_initialize_output nf90_def_var_fill NF90_CHAR" ) end select end do ! Take the file out of define mode - call check( nf90_enddef(nc%id), "encounter_io_initialize nf90_enddef" ) + call check( nf90_enddef(nc%id), "encounter_io_initialize_output nf90_enddef" ) ! Add in the space dimension coordinates - call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize nf90_put_var space" ) + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize_output nf90_put_var space" ) end associate @@ -191,7 +147,7 @@ module subroutine encounter_io_initialize(self, param) 667 continue write(*,*) "Error creating encounter output file. " // trim(adjustl(errmsg)) call util_exit(FAILURE) - end subroutine encounter_io_initialize + end subroutine encounter_io_initialize_output module subroutine encounter_io_write_frame_snapshot(self, nc, param) diff --git a/src/encounter/encounter_placeholder.f90 b/src/encounter/encounter_placeholder.f90 deleted file mode 100644 index 61ae6eebd..000000000 --- a/src/encounter/encounter_placeholder.f90 +++ /dev/null @@ -1,53 +0,0 @@ -!! 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/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 4ad5c75d1..7c9cfa358 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -17,15 +17,14 @@ contains - module subroutine collision_generate_fragments(self, impactors, system, param, lfailure) + module subroutine fraggle_generate_fragments(self, 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. use, intrinsic :: ieee_exceptions implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation - class(collision_impactors), intent(inout) :: impactors !! Fraggle impactors object containing the two-body equivalent values of the colliding bodies + class(fraggle_system), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation 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? @@ -125,9 +124,9 @@ module subroutine collision_generate_fragments(self, impactors, system, param, l dLmag = .mag. (fragments%Ltot_after(:) - fragments%Ltot_before(:)) exit - lfailure = ((abs(dEtot + fragments%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) + lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) if (lfailure) then - write(message, *) dEtot, abs(dEtot + fragments%Qloss) / FRAGGLE_ETOL + write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high energy error: " // & trim(adjustl(message))) cycle @@ -166,7 +165,7 @@ module subroutine collision_generate_fragments(self, impactors, system, param, l 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 collision_generate_fragments + end subroutine fraggle_generate_fragments subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) @@ -190,8 +189,8 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) associate(nfrag => fragments%nbody) allocate(loverlap(nfrag)) - lnoncat = (fragments%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point - lhitandrun = (fragments%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! Disruptive hit and runs have their own fragment distribution + lnoncat = (impactors%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point + lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! Disruptive hit and runs have their own fragment distribution ! 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) @@ -279,12 +278,12 @@ subroutine fraggle_generate_spins(fragments, impactors, 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(:) = impactors%L_spin(:,1) + f_spin * (impactors%L_orbit(:,2) + impactors%L_spin(:,2)) + L(:) = impactors%Lspin(:,1) + f_spin * (impactors%Lorbit(:,2) + impactors%Lspin(:,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(:) = impactors%L_spin(:,2) / (nfrag - 1) + L(:) = impactors%Lspin(:,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 @@ -425,8 +424,8 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, " ") call io_log_one_message(FRAGGLE_LOG_OUT, "Tangential velocity failure diagnostics") - call fragments%get_ang_mtm() - L_frag_tot = fragments%L_spin(:) + fragments%L_orbit(:) + call fragments%get_angular_momentum() + L_frag_tot = fragments%Lspin(:) + fragments%Lorbit(:) write(message, *) .mag.(fragments%L_budget(:) - L_frag_tot(:)) / (.mag.fragments%Ltot_before(:)) call io_log_one_message(FRAGGLE_LOG_OUT, "|L_remainder| : " // trim(adjustl(message))) write(message, *) fragments%ke_budget @@ -479,7 +478,7 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) end if end do b(1:3) = -L_lin_others(:) - b(4:6) = fragments%L_budget(:) - fragments%L_spin(:) - L_orb_others(:) + b(4:6) = fragments%L_budget(:) - fragments%Lspin(:) - L_orb_others(:) allocate(v_t_mag_output(nfrag)) v_t_mag_output(1:6) = util_solve_linear_system(A, b, 6, lfailure) if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 1ddc4b4e9..c55db4173 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -30,7 +30,7 @@ module subroutine fraggle_io_log_regime(impactors, fragments) write(LUN, *) "--------------------------------------------------------------------" write(LUN, *) "True number of impactors : ",impactors%ncoll write(LUN, *) "Index list of true impactors : ",impactors%idx(1:impactors%ncoll) - select case(fragments%regime) + select case(impactors%regime) case(COLLRESOLVE_REGIME_MERGE) write(LUN, *) "Regime: Merge" case(COLLRESOLVE_REGIME_DISRUPTION) @@ -42,7 +42,7 @@ module subroutine fraggle_io_log_regime(impactors, fragments) case(COLLRESOLVE_REGIME_HIT_AND_RUN) write(LUN, *) "Regime: Hit and run" end select - write(LUN, *) "Energy loss : ", fragments%Qloss + write(LUN, *) "Energy loss : ", impactors%Qloss write(LUN, *) "--------------------------------------------------------------------" close(LUN) diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index 1563a4ccf..2873d236d 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -60,7 +60,7 @@ module subroutine encounter_regime_impactors(self, fragments, system, param) !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime 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) + min_mfrag_si, impactors%regime, mlr, mslr, impactors%Qloss) fragments%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) fragments%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) @@ -78,7 +78,7 @@ module subroutine encounter_regime_impactors(self, fragments, system, param) ! 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 + impactors%Qloss = impactors%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG call fraggle_io_log_regime(impactors, fragments) end associate diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 122f593ab..fa0e28aeb 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -17,18 +17,18 @@ module subroutine fraggle_set_budgets_fragments(self) !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum implicit none ! Arguments - class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object + class(collision_system), intent(inout) :: self !! Fraggle fragment system object ! Internals real(DP) :: dEtot real(DP), dimension(NDIM) :: dL - associate(fragments => self) + associate(impactors => self%impactors, fragments => self%fragments) - dEtot = fragments%Etot_after - fragments%Etot_before - dL(:) = fragments%Ltot_after(:) - fragments%Ltot_before(:) + dEtot = self%Etot(2) - self%Etot(1) + dL(:) = self%Ltot(:,2) - self%Ltot(:,1) fragments%L_budget(:) = -dL(:) - fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(fragments%vbcom(:), fragments%vbcom(:))) - fragments%Qloss + fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(fragments%vbcom(:), fragments%vbcom(:))) - impactors%Qloss end associate return @@ -71,7 +71,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) jproj = 1 end if - select case(fragments%regime) + select case(impactors%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 encounter_regime. ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on @@ -110,7 +110,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) return case default - write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",fragments%regime + write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime end select ! Make the first two bins the same as the Mlr and Mslr values that came from encounter_regime @@ -139,7 +139,7 @@ module subroutine fraggle_set_mass_dist_fragments(self, impactors, param) fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining ! Compute physical properties of the new fragments - select case(fragments%regime) + select case(impactors%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) = impactors%radius(jtarg) fragments%density(1) = fragments%mass_dist(iMlr) / volume(jtarg) @@ -183,7 +183,7 @@ module subroutine encounter_set_coordinate_system(self, impactors) ! 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) + Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2) L_mag = .mag.Ltot(:) if (L_mag > sqrt(tiny(L_mag))) then @@ -249,24 +249,24 @@ module subroutine fraggle_set_natural_scale_factors(self, impactors) 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 + impactors%Lspin(:,:) = impactors%Lspin(:,:) / fragments%Lscale + impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / fragments%Lscale do i = 1, 2 - impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) + impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do fragments%mtot = fragments%mtot / fragments%mscale fragments%mass = fragments%mass / fragments%mscale fragments%radius = fragments%radius / fragments%dscale - fragments%Qloss = fragments%Qloss / fragments%Escale + impactors%Qloss = impactors%Qloss / fragments%Escale end associate return end subroutine fraggle_set_natural_scale_factors - module subroutine fraggle_set_original_scale_factors(self, impactors) + module subroutine fraggle_set_original_scale_factors(self) !! author: David A. Minton !! !! Restores dimenional quantities back to the system units @@ -293,9 +293,9 @@ module subroutine fraggle_set_original_scale_factors(self, impactors) 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 + impactors%Lspin = impactors%Lspin * fragments%Lscale do i = 1, 2 - impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) + impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do fragments%mtot = fragments%mtot * fragments%mscale @@ -310,7 +310,7 @@ module subroutine fraggle_set_original_scale_factors(self, impactors) fragments%vb(:, i) = fragments%v_coll(:, i) + fragments%vbcom(:) end do - fragments%Qloss = fragments%Qloss * fragments%Escale + impactors%Qloss = impactors%Qloss * fragments%Escale fragments%Lorbit_before(:) = fragments%Lorbit_before * fragments%Lscale fragments%Lspin_before(:) = fragments%Lspin_before * fragments%Lscale diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 index 3d8916967..7b6551d88 100644 --- a/src/fraggle/fraggle_setup.f90 +++ b/src/fraggle/fraggle_setup.f90 @@ -22,8 +22,6 @@ module subroutine fraggle_setup_reset_fragments(self) self%rb(:,:) = 0.0_DP self%vb(:,:) = 0.0_DP self%rot(:,:) = 0.0_DP - self%r_coll(:,:) = 0.0_DP - self%v_coll(:,:) = 0.0_DP self%v_r_unit(:,:) = 0.0_DP self%v_t_unit(:,:) = 0.0_DP self%v_n_unit(:,:) = 0.0_DP @@ -33,11 +31,6 @@ module subroutine fraggle_setup_reset_fragments(self) self%v_r_mag(:) = 0.0_DP self%v_t_mag(:) = 0.0_DP - self%ke_orbit = 0.0_DP - self%ke_spin = 0.0_DP - self%L_orbit(:) = 0.0_DP - self%L_spin(:) = 0.0_DP - return end subroutine fraggle_setup_reset_fragments @@ -52,7 +45,7 @@ module subroutine fraggle_setup_fragments(self, n, param) integer(I4B), intent(in) :: n class(swiftest_parameters), intent(in) :: param - call self%collision_fragments%setup(n, param) + call collision_util_setup_fragments(n, param) if (n < 0) return if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index bb32167b1..4d1a203fb 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -59,7 +59,7 @@ module subroutine fraggle_util_add_fragments_to_system(fragments, impactors, sys end subroutine fraggle_util_add_fragments_to_system - module subroutine fraggle_util_ang_mtm(self) + module subroutine fraggle_util_get_angular_momentum(self) !! Author: David A. Minton !! !! Calcualtes the current angular momentum of the fragments @@ -70,17 +70,17 @@ module subroutine fraggle_util_ang_mtm(self) integer(I4B) :: i associate(fragments => self, nfrag => self%nbody) - fragments%L_orbit(:) = 0.0_DP - fragments%L_spin(:) = 0.0_DP + fragments%Lorbit(:) = 0.0_DP + fragments%Lspin(:) = 0.0_DP do i = 1, nfrag - fragments%L_orbit(:) = fragments%L_orbit(:) + fragments%mass(i) * (fragments%r_coll(:, i) .cross. fragments%v_coll(:, i)) - fragments%L_spin(:) = fragments%L_spin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) + fragments%Lorbit(:) = fragments%Lorbit(:) + fragments%mass(i) * (fragments%rc(:, i) .cross. fragments%vc(:, i)) + fragments%Lspin(:) = fragments%Lspin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) end do end associate return - end subroutine fraggle_util_ang_mtm + end subroutine fraggle_util_get_angular_momentum module subroutine fraggle_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) @@ -131,6 +131,24 @@ module subroutine fraggle_util_construct_temporary_system(fragments, system, par end subroutine fraggle_util_construct_temporary_system + module subroutine fraggle_util_dealloc_fragments(self) + !! author: David A. Minton + !! + !! Deallocates all allocatables + implicit none + ! Arguments + class(fraggle_fragments), intent(inout) :: self + + call collision_util_deallocate_fragments(self) + + 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_n_mag) + + return + end subroutine fraggle_util_dealloc_fragments + + module subroutine fraggle_util_final_impactors(self) !! author: David A. Minton @@ -154,20 +172,26 @@ module subroutine fraggle_util_final_fragments(self) ! Arguments type(fraggle_fragments), intent(inout) :: self !! Fraggle encountar storage object - 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) + call self%dealloc() return end subroutine fraggle_util_final_fragments + module subroutine fraggle_util_final_system(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(fraggle_system), intent(inout) :: self !! Collision impactors storage object + + call self%reset() + + return + end subroutine fraggle_util_final_system + + module subroutine fraggle_util_restructure(self, impactors, try, f_spin, r_max_start) !! Author: David A. Minton !! diff --git a/src/modules/collision_classes.f90 b/src/modules/collision_classes.f90 index d3ede5fcb..d2d2bd9b1 100644 --- a/src/modules/collision_classes.f90 +++ b/src/modules/collision_classes.f90 @@ -36,8 +36,8 @@ module collision_classes 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) :: Lspin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: Lorbit !! 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 @@ -77,12 +77,12 @@ module collision_classes 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 + procedure :: accel => collision_util_placeholder_accel !! Placeholder subroutine to fulfill requirement for an accel method + procedure :: kick => collision_util_placeholder_kick !! Placeholder subroutine to fulfill requirement for a kick method + procedure :: step => collision_util_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 :: dealloc => collision_util_dealloc_fragments !! Deallocates all allocatable arrays end type collision_fragments type :: collision_system @@ -97,12 +97,11 @@ module collision_classes 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 + 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.) @@ -110,11 +109,25 @@ module collision_classes final :: collision_util_final_system !! Finalizer will deallocate all allocatables end type collision_system - - !> NetCDF dimension and variable names for the enounter save object + !! NetCDF dimension and variable names for the enounter save object type, extends(encounter_io_parameters) :: collision_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 => collision_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object + procedure :: initialize => collision_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object end type collision_io_parameters type, extends(encounter_snapshot) :: collision_snapshot @@ -136,26 +149,17 @@ module collision_classes 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) + module subroutine collision_io_initialize_output(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 + end subroutine collision_io_initialize_output module subroutine collision_io_write_frame_snapshot(self, nc, param) implicit none @@ -165,7 +169,7 @@ module subroutine collision_io_write_frame_snapshot(self, nc, param) 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) + module subroutine collision_util_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 @@ -173,9 +177,9 @@ module subroutine collision_placeholder_accel(self, system, param, t, lbeg) 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 + end subroutine collision_util_placeholder_accel - module subroutine collision_placeholder_kick(self, system, param, t, dt, lbeg) + module subroutine collision_util_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 @@ -184,9 +188,9 @@ module subroutine collision_placeholder_kick(self, system, param, t, dt, lbeg) 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 + end subroutine collision_util_placeholder_kick - module subroutine collision_placeholder_step(self, system, param, t, dt) + module subroutine collision_util_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 @@ -194,7 +198,7 @@ module subroutine collision_placeholder_step(self, system, param, t, dt) 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 + end subroutine collision_util_placeholder_step module subroutine collision_regime_system(self, system, param) implicit none @@ -235,10 +239,10 @@ module subroutine collision_setup_system(self, system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine collision_setup_system - module subroutine collision_util_reset_fragments(self) + module subroutine collision_util_dealloc_fragments(self) implicit none class(collision_fragments), intent(inout) :: self - end subroutine collision_util_reset_fragments + end subroutine collision_util_dealloc_fragments module subroutine collision_util_final_impactors(self) implicit none @@ -257,7 +261,7 @@ end subroutine collision_util_final_snapshot module subroutine collision_util_final_system(self) implicit none - type(collision_system), intent(inout) :: self !! Collision system object + type(collision_system), intent(inout) :: self !! Collision system object end subroutine collision_util_final_system module subroutine collision_util_get_idvalues_snapshot(self, idvals) diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 84a7fdf4c..e3aa2f9fe 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -72,7 +72,7 @@ module encounter_classes 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 + procedure :: initialize => encounter_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object end type encounter_io_parameters type encounter_bounding_box_1D @@ -213,11 +213,11 @@ module subroutine encounter_io_dump(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine encounter_io_dump - module subroutine encounter_io_initialize(self, param) + module subroutine encounter_io_initialize_output(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 + end subroutine encounter_io_initialize_output module subroutine encounter_io_write_frame_snapshot(self, nc, param) implicit none diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 6ef73ffb1..b8a558af8 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -27,10 +27,14 @@ module fraggle_classes !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates type, extends(collision_fragments) :: fraggle_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 - real(DP), dimension(:), allocatable :: v_n_mag !! Array of normal direction velocity 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 + real(DP), dimension(:), allocatable :: v_n_mag !! Array of normal direction velocity magnitudes of individual fragments + real(DP), dimension(NDIM) :: Lorbit !! Orbital angular momentum vector of all fragments + real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments + real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments + real(DP) :: ke_spin !! Spin kinetic energy of all fragments + real(DP) :: ke_budget !! Kinetic energy budget for computing quanities ! 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 @@ -40,32 +44,36 @@ module fraggle_classes real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) 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 :: set_budgets => fraggle_set_budgets_fragments !! Sets the energy and momentum budgets of the fragments based on the collider value 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_angular_momentum => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments + procedure :: dealloc => fraggle_util_dealloc_fragments !! Deallocates all allocatables 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 + type, extends(collision_system) :: fraggle_system + contains + procedure :: generate_fragments => fraggle_generate_fragments !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + final :: fraggle_util_final_system !! Finalizer will deallocate all allocatables + end type fraggle_system + interface - module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure) + module subroutine fraggle_generate_fragments(self, 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(collision_impactors), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies + class(fraggle_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 fraggle_generate_fragments - module subroutine fraggle_io_log_regime(impactors, fragments) implicit none class(collision_impactors), intent(in) :: impactors @@ -117,10 +125,10 @@ module subroutine fraggle_util_add_fragments_to_system(fragments, impactors, sys class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine fraggle_util_add_fragments_to_system - module subroutine fraggle_util_ang_mtm(self) + module subroutine fraggle_util_get_angular_momentum(self) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object - end subroutine fraggle_util_ang_mtm + end subroutine fraggle_util_get_angular_momentum module subroutine fraggle_util_construct_temporary_system(fragments, system, param, tmpsys, tmpparam) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters @@ -132,6 +140,11 @@ 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_dealloc_fragments(self) + implicit none + class(fraggle_fragments), intent(inout) :: self + end subroutine fraggle_util_dealloc_fragments + module subroutine fraggle_util_final_impactors(self) implicit none type(collision_impactors), intent(inout) :: self !! Fraggle impactors object @@ -142,6 +155,10 @@ module subroutine fraggle_util_final_fragments(self) type(fraggle_fragments), intent(inout) :: self !! Fraggle frgments object end subroutine fraggle_util_final_fragments + module subroutine fraggle_util_final_system(self) + implicit none + type(fraggle_system), intent(inout) :: self !! Collision system object + end subroutine fraggle_util_final_system module subroutine fraggle_util_restructure(self, impactors, try, f_spin, r_max_start) implicit none diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 090446da8..cf2ae97f4 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -105,8 +105,8 @@ module swiftest_classes integer(I4B) :: PE_varid !! ID for the system potential energy variable character(NAMELEN) :: L_orb_varname = "L_orb" !! name of the orbital angular momentum vector variable integer(I4B) :: L_orb_varid !! ID for the system orbital angular momentum vector variable - character(NAMELEN) :: L_spin_varname = "L_spin" !! name of the spin angular momentum vector variable - integer(I4B) :: L_spin_varid !! ID for the system spin angular momentum vector variable + character(NAMELEN) :: Lspin_varname = "Lspin" !! name of the spin angular momentum vector variable + integer(I4B) :: Lspin_varid !! ID for the system spin angular momentum vector variable character(NAMELEN) :: L_escape_varname = "L_escape" !! name of the escaped angular momentum vector variable integer(I4B) :: L_escape_varid !! ID for the escaped angular momentum vector variable character(NAMELEN) :: Ecollisions_varname = "Ecollisions" !! name of the escaped angular momentum y variable diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 9e710c8fe..98143ed08 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -16,7 +16,8 @@ module symba_classes 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 : collision_impactors, fraggle_fragments - use encounter_classes, only : encounter_list, encounter_storage, collision_storage + use encounter_classes, only : encounter_list, encounter_storage + use collision_classes, only : collision_storage implicit none public diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 4ce217124..eaf11ddd4 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -106,7 +106,7 @@ module function netcdf_get_old_t_final_system(self, param) result(old_t_final) self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + self%Ecollisions + self%Euntracked call check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_orb_varid" ) - call check( nf90_get_var(nc%id, nc%L_spin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_spin_varid" ) + call check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system Lspin_varid" ) call check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_get_old_t_final_system L_escape_varid" ) self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) @@ -265,7 +265,7 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, nc%time_dimid, nc%KE_spin_varid), "netcdf_initialize_output nf90_def_var KE_spin_varid" ) call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), "netcdf_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%time_dimid], nc%L_orb_varid), "netcdf_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%time_dimid], nc%L_spin_varid), "netcdf_initialize_output nf90_def_var L_spin_varid" ) + call check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%Lspin_varid), "netcdf_initialize_output nf90_def_var Lspin_varid" ) call check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_escape_varid), "netcdf_initialize_output nf90_def_var L_escape_varid" ) call check( nf90_def_var(nc%id, nc%Ecollisions_varname, nc%out_type, nc%time_dimid, nc%Ecollisions_varid), "netcdf_initialize_output nf90_def_var Ecollisions_varid" ) call check( nf90_def_var(nc%id, nc%Euntracked_varname, nc%out_type, nc%time_dimid, nc%Euntracked_varid), "netcdf_initialize_output nf90_def_var Euntracked_varid" ) @@ -429,7 +429,7 @@ module subroutine netcdf_open(self, param, readonly) status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) - status = nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid) + status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) @@ -768,8 +768,8 @@ module subroutine netcdf_read_hdr_system(self, nc, param) if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_read_hdr_system nf90_getvar PE_varid" ) status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar L_orb_varid" ) - status = nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid) - if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%L_spin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar L_spin_varid" ) + status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) + if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar Lspin_varid" ) status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) if (status == nf90_noerr) call check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_read_hdr_system nf90_getvar L_escape_varid" ) status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) @@ -1277,7 +1277,7 @@ module subroutine netcdf_write_hdr_system(self, nc, param) call check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_write_hdr_system nf90_put_var KE_spin_varid" ) call check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_write_hdr_system nf90_put_var PE_varid" ) call check( nf90_put_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_orb_varid" ) - call check( nf90_put_var(nc%id, nc%L_spin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_spin_varid" ) + call check( nf90_put_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var Lspin_varid" ) call check( nf90_put_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_write_hdr_system nf90_put_var L_escape_varid" ) call check( nf90_put_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_write_hdr_system nf90_put_var Ecollisions_varid" ) call check( nf90_put_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_write_hdr_system nf90_put_var Euntracked_varid" ) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 90d3847c4..16fd0d7bd 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -32,7 +32,7 @@ module function symba_collision_casedisruption(system, param, t) result(status) associate(impactors => system%impactors, fragments => system%fragments) - select case(fragments%regime) + select case(impactors%regime) case(COLLRESOLVE_REGIME_DISRUPTION) message = "Disruption between" case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) @@ -66,7 +66,7 @@ module function symba_collision_casedisruption(system, param, t) result(status) nfrag = fragments%nbody write(message, *) nfrag call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - select case(fragments%regime) + select case(impactors%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTION ibiggest = impactors%idx(maxloc(system%pl%Gmass(impactors%idx(:)), dim=1)) @@ -183,7 +183,7 @@ module function symba_collision_casemerge(system, param, t) result(status) integer(I4B) :: status !! Status flag assigned to this outcome ! Internals integer(I4B) :: i, j, k, ibiggest - real(DP), dimension(NDIM) :: L_spin_new + real(DP), dimension(NDIM) :: Lspin_new real(DP) :: dpe character(len=STRMAX) :: message @@ -207,12 +207,12 @@ module function symba_collision_casemerge(system, param, t) result(status) if (param%lrotation) then ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = impactors%L_orbit(:,1) + impactors%L_orbit(:,2) + impactors%L_spin(:,1) + impactors%L_spin(:,2) + Lspin_new(:) = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,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) + fragments%rot(:,1) = Lspin_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(:) + impactors%L_orbit(:,1) + impactors%L_orbit(:,2) + system%Lescape(:) = system%Lescape(:) + impactors%Lorbit(:,1) + impactors%Lorbit(:,2) end if ! Keep track of the component of potential energy due to the pre-impact impactors%idx for book-keeping @@ -471,7 +471,7 @@ function symba_collision_consolidate_impactors(pl, cb, param, idx_parent, impact !! author: David A. Minton !! !! 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. + !! and pairs of quantities (x and v vectors, mass, radius, Lspin, and Ip) that can be used to resolve the collisional outcome. implicit none ! Arguments class(symba_pl), intent(inout) :: pl !! SyMBA massive body object @@ -535,7 +535,7 @@ function symba_collision_consolidate_impactors(pl, cb, param, idx_parent, impact impactors%ncoll = count(pl%lcollision(impactors%idx(:))) impactors%idx = pack(impactors%idx(:), pl%lcollision(impactors%idx(:))) - impactors%L_spin(:,:) = 0.0_DP + impactors%Lspin(:,:) = 0.0_DP impactors%Ip(:,:) = 0.0_DP ! Find the barycenter of each body along with its children, if it has any @@ -545,7 +545,7 @@ function symba_collision_consolidate_impactors(pl, cb, param, idx_parent, impact ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then 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)) + impactors%Lspin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) end if if (nchild(j) > 0) then @@ -565,14 +565,14 @@ function symba_collision_consolidate_impactors(pl, cb, param, idx_parent, impact xc(:) = impactors%rb(:, j) - xcom(:) vc(:) = impactors%vb(:, j) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - impactors%L_spin(:, j) = impactors%L_spin(:, j) + impactors%mass(j) * xcrossv(:) + impactors%Lspin(:, j) = impactors%Lspin(:, j) + impactors%mass(j) * xcrossv(:) xc(:) = xchild(:) - xcom(:) vc(:) = vchild(:) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * xcrossv(:) + impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * xcrossv(:) - impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & + impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * pl%Ip(3, idx_child) & * pl%radius(idx_child)**2 & * pl%rot(:, idx_child) impactors%Ip(:, j) = impactors%Ip(:, j) + mchild * pl%Ip(:, idx_child) @@ -596,7 +596,7 @@ function symba_collision_consolidate_impactors(pl, cb, param, idx_parent, impact 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(:,:) + impactors%Lorbit(:,:) = mxc(:,:) .cross. vcc(:,:) ! Destroy the kinship relationships for all members of this impactors%idx call pl%reset_kinship(impactors%idx(:)) @@ -937,7 +937,7 @@ subroutine symba_resolve_collision(plplcollision_list , system, param, t) call system%impactors%regime(system%fragments, system, param) else associate(fragments => system%fragments, impactors => system%impactors) - fragments%regime = COLLRESOLVE_REGIME_MERGE + impactors%regime = COLLRESOLVE_REGIME_MERGE fragments%mtot = sum(impactors%mass(:)) fragments%mass_dist(1) = fragments%mtot fragments%mass_dist(2) = 0.0_DP @@ -948,7 +948,7 @@ subroutine symba_resolve_collision(plplcollision_list , system, param, t) end if call collision_history%take_snapshot(param,system, t, "before") - select case (system%fragments%regime) + select case (system%impactors%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) plplcollision_list%status(i) = symba_collision_casedisruption(system, param, t) case (COLLRESOLVE_REGIME_HIT_AND_RUN)