From 2ffdbe4b0f56fb4704f2a116444c848cd7af64af Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 7 Jan 2023 17:42:04 -0500 Subject: [PATCH] Made adjustments to initial position generator in Fraggle. Also moved NetCDF file I/O operations out of the nbody system constructor, as that is now used internally to make temporary systems that shouldn't be opening files. --- src/fraggle/fraggle_generate.f90 | 4 +-- src/swiftest/swiftest_util.f90 | 47 +++----------------------------- src/symba/symba_module.f90 | 1 - src/symba/symba_util.f90 | 41 ++++++++++++++++++++++++++++ 4 files changed, 47 insertions(+), 46 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index df6fb6d5f..28359ea6e 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -252,7 +252,7 @@ module subroutine fraggle_generate_pos_vec(collider) ! We will treat the first two fragments of the list as special cases. ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap loverlap(:) = .true. - rdistance(:) = rdistance_scale_factor * .mag.impactors%vc(:,:) + rdistance(:) = rdistance_scale_factor * impactors%radius(:) ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be. call random_number(mass_rscale) @@ -332,7 +332,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) ! Internals real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir real(DP) :: v_init, v_final, mass_init, mass_final, rotmag - real(DP), parameter :: random_scale_factor = 0.1_DP !! The relative scale factor to apply to the random component of the rotation vector + real(DP), parameter :: random_scale_factor = 0.01_DP !! The relative scale factor to apply to the random component of the rotation vector integer(I4B) :: i associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 30e3080ad..093756af3 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -2541,13 +2541,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) ! Arguments class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - type(encounter_storage) :: encounter_history - type(collision_storage) :: collision_history - allocate(swiftest_storage(param%dump_cadence) :: param%system_history) - allocate(swiftest_netcdf_parameters :: param%system_history%nc) - call param%system_history%reset() select case(param%integrator) case (INT_BS) @@ -2602,32 +2596,6 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) allocate(symba_list_plpl :: nbody_system%plpl_encounter) allocate(collision_list_plpl :: nbody_system%plpl_collision) - if (param%lenc_save_trajectory .or. param%lenc_save_closest) then - allocate(encounter_netcdf_parameters :: encounter_history%nc) - call encounter_history%reset() - select type(nc => encounter_history%nc) - class is (encounter_netcdf_parameters) - nc%file_name = ENCOUNTER_OUTFILE - if (.not.param%lrestart) then - call nc%initialize(param) - call nc%close() - end if - end select - allocate(nbody_system%encounter_history, source=encounter_history) - end if - - allocate(collision_netcdf_parameters :: collision_history%nc) - call collision_history%reset() - select type(nc => collision_history%nc) - class is (collision_netcdf_parameters) - nc%file_name = COLLISION_OUTFILE - if (.not.param%lrestart) then - call nc%initialize(param) - call nc%close() - end if - end select - allocate(nbody_system%collision_history, source=collision_history) - end select case (INT_RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' @@ -2638,15 +2606,6 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) allocate(swiftest_particle_info :: nbody_system%cb%info) - select case(param%collision_model) - case("MERGE") - allocate(collision_basic :: nbody_system%collider) - case("BOUNCE") - allocate(collision_bounce :: nbody_system%collider) - case("FRAGGLE") - allocate(collision_fraggle :: nbody_system%collider) - end select - call nbody_system%collider%setup(nbody_system) nbody_system%t = param%tstart @@ -2698,8 +2657,11 @@ module subroutine swiftest_util_setup_initialize_system(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - associate(nbody_system => self, cb => self%cb, pl => self%pl, tp => self%tp, nc => param%system_history%nc) + allocate(swiftest_storage(param%dump_cadence) :: param%system_history) + allocate(swiftest_netcdf_parameters :: param%system_history%nc) + call param%system_history%reset() + associate(nbody_system => self, cb => self%cb, pl => self%pl, tp => self%tp, nc => param%system_history%nc) call nbody_system%read_in(param) call nbody_system%validate_ids(param) call nbody_system%set_msys() @@ -2722,7 +2684,6 @@ module subroutine swiftest_util_setup_initialize_system(self, param) nc%file_name = param%outfile call nbody_system%write_frame(param) call nc%close() - end associate return diff --git a/src/symba/symba_module.f90 b/src/symba/symba_module.f90 index bf3ab8df7..bdc7011e3 100644 --- a/src/symba/symba_module.f90 +++ b/src/symba/symba_module.f90 @@ -22,7 +22,6 @@ module symba real(DP), private, parameter :: RHSCALE = 6.5_DP real(DP), private, parameter :: RSHELL = 0.48075_DP - !> SyMBA central body particle class type, extends(helio_cb) :: symba_cb end type symba_cb diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 0a78693b2..c3dc1be18 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -9,6 +9,7 @@ submodule(symba) s_symba_util use swiftest + use fraggle contains module subroutine symba_util_append_pl(self, source, lsource_mask) @@ -267,6 +268,9 @@ module subroutine symba_util_setup_initialize_system(self, param) ! Arguments class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + type(encounter_storage) :: encounter_history + type(collision_storage) :: collision_history ! Call parent method associate(nbody_system => self) @@ -274,6 +278,43 @@ module subroutine symba_util_setup_initialize_system(self, param) call nbody_system%pltp_encounter%setup(0_I8B) call nbody_system%plpl_encounter%setup(0_I8B) call nbody_system%plpl_collision%setup(0_I8B) + + if (param%lenc_save_trajectory .or. param%lenc_save_closest) then + allocate(encounter_netcdf_parameters :: encounter_history%nc) + call encounter_history%reset() + select type(nc => encounter_history%nc) + class is (encounter_netcdf_parameters) + nc%file_name = ENCOUNTER_OUTFILE + if (.not.param%lrestart) then + call nc%initialize(param) + call nc%close() + end if + end select + allocate(nbody_system%encounter_history, source=encounter_history) + end if + + allocate(collision_netcdf_parameters :: collision_history%nc) + call collision_history%reset() + select type(nc => collision_history%nc) + class is (collision_netcdf_parameters) + nc%file_name = COLLISION_OUTFILE + if (.not.param%lrestart) then + call nc%initialize(param) + call nc%close() + end if + end select + allocate(nbody_system%collision_history, source=collision_history) + + select case(param%collision_model) + case("MERGE") + allocate(collision_basic :: nbody_system%collider) + case("BOUNCE") + allocate(collision_bounce :: nbody_system%collider) + case("FRAGGLE") + allocate(collision_fraggle :: nbody_system%collider) + end select + call nbody_system%collider%setup(nbody_system) + end associate return