From 8bd3cdfb888f11321975f1d90b7775bddfe41d8c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 17:45:27 -0500 Subject: [PATCH 1/8] Started the process of writing the write method for encounter storage objects --- src/encounter/encounter_io.f90 | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 8ed849ba2..ac187d54f 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -105,6 +105,7 @@ module subroutine encounter_io_initialize_output(self, param) call check( nf90_def_var(nciu%id, nciu%Ip_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%Ip_varid), "encounter_io_initialize_output nf90_def_var Ip_varid" ) call check( nf90_def_var(nciu%id, nciu%rot_varname, nciu%out_type, [nciu%space_dimid, nciu%id_dimid, nciu%time_dimid], nciu%rot_varid), "encounter_io_initialize_output nf90_def_var rot_varid" ) end if + call check( nf90_inquire(nciu%id, nVariables=nvar), "encounter_io_initialize_output nf90_inquire nVariables" ) do varid = 1, nvar call check( nf90_inquire_variable(nciu%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize_output nf90_inquire_variable" ) @@ -139,16 +140,31 @@ module subroutine encounter_io_write_frame(self, iu, param) implicit none ! Arguments class(encounter_list), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: iu !! Parameters used to identify a particular encounter io NetCDF dataset + class(encounter_io_parameters), intent(inout) :: nciu !! Parameters used to identify a particular encounter io NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i,old_mode, n + integer(I4B) :: tslot,i,old_mode, n + character(len=NAMELEN) :: charstring + + tslot = nciu%ienc_frame + call check( nf90_set_fill(nciu%id, nf90_nofill, old_mode), "encounter_io_write_frame_base nf90_set_fill" ) + call check( nf90_put_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) + + ! charstring = trim(adjustl(self%info(j)%name)) + ! call check( nf90_put_var(nciu%id, nciu%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_info_base nf90_put_var name_varid" ) + + ! charstring = trim(adjustl(self%info(j)%particle_type)) + ! call check( nf90_put_var(nciu%id, nciu%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_write_info_base nf90_put_var particle_type_varid" ) - i = iu%ienc_frame - n = int(self%nenc, kind=I4B) - call check( nf90_set_fill(iu%id, nf90_nofill, old_mode), "encounter_io_write_frame_base nf90_set_fill" ) - call check( nf90_put_var(iu%id, iu%time_varid, self%t, start=[i]), "encounter_io_write_frame nf90_put_var time_varid" ) + ! call check( nf90_put_var(nciu%id, nciu%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var rh_varid" ) + ! call check( nf90_put_var(nciu%id, nciu%vh_varid, self%vh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var vh_varid" ) + ! call check( nf90_put_var(nciu%id, nciu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "encounter_io_write_frame_base nf90_put_var body Gmass_varid" ) + ! if (param%lclose) call check( nf90_put_var(nciu%id, nciu%radius_varid, self%radius(j), start=[idslot, tslot]), "encounter_io_write_frame_base nf90_put_var body radius_varid" ) + ! if (param%lrotation) then + ! call check( nf90_put_var(nciu%id, nciu%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var body Ip_varid" ) + ! call check( nf90_put_var(nciu%id, nciu%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "encounter_io_write_frame_base nf90_put_var body rotx_varid" ) + ! end if return end subroutine encounter_io_write_frame From c094acc315f854c820cd55a42fd48bb01709fc4f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 17:49:33 -0500 Subject: [PATCH 2/8] fixed argument list in encounter write interface and removed unnecessary encounter_list variables --- src/encounter/encounter_io.f90 | 2 +- src/encounter/encounter_setup.f90 | 20 -------------------- src/encounter/encounter_util.f90 | 24 ------------------------ src/modules/encounter_classes.f90 | 10 ++-------- src/symba/symba_util.f90 | 18 ------------------ 5 files changed, 3 insertions(+), 71 deletions(-) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index ac187d54f..62920db7e 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -133,7 +133,7 @@ module subroutine encounter_io_initialize_output(self, param) end subroutine encounter_io_initialize_output - module subroutine encounter_io_write_frame(self, iu, param) + module subroutine encounter_io_write_frame(self, nciu, param) !! author: David A. Minton !! !! Write a frame of output of an encounter list structure. diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 index c741cf0a7..22253b4e4 100644 --- a/src/encounter/encounter_setup.f90 +++ b/src/encounter/encounter_setup.f90 @@ -77,12 +77,6 @@ module subroutine encounter_setup_list(self, n) if (allocated(self%x2)) deallocate(self%x2) if (allocated(self%v1)) deallocate(self%v1) if (allocated(self%v2)) deallocate(self%v2) - if (allocated(self%Gmass1)) deallocate(self%Gmass1) - if (allocated(self%Gmass2)) deallocate(self%Gmass2) - if (allocated(self%radius1)) deallocate(self%radius1) - if (allocated(self%radius2)) deallocate(self%radius2) - if (allocated(self%name1)) deallocate(self%name1) - if (allocated(self%name2)) deallocate(self%name2) self%nenc = n if (n == 0_I8B) return @@ -98,12 +92,6 @@ module subroutine encounter_setup_list(self, n) allocate(self%x2(NDIM,n)) allocate(self%v1(NDIM,n)) allocate(self%v2(NDIM,n)) - allocate(self%Gmass1(n)) - allocate(self%Gmass2(n)) - allocate(self%radius1(n)) - allocate(self%radius2(n)) - allocate(self%name1(n)) - allocate(self%name2(n)) self%lvdotr(:) = .false. self%status(:) = INACTIVE @@ -115,14 +103,6 @@ module subroutine encounter_setup_list(self, n) self%x2(:,:) = 0.0_DP self%v1(:,:) = 0.0_DP self%v2(:,:) = 0.0_DP - self%Gmass1(:) = 0.0_DP - self%Gmass2(:) = 0.0_DP - self%radius1(:) = 0.0_DP - self%radius2(:) = 0.0_DP - do i = 1_I8B, n - self%name1(i) = "UNNAMED" - self%name2(i) = "UNNAMED" - end do return end subroutine encounter_setup_list diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index cf7dc3a71..d0c801ab4 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -36,12 +36,6 @@ module subroutine encounter_util_append_list(self, source, lsource_mask) call util_append(self%x2, source%x2, nold, nsrc, lsource_mask) call util_append(self%v1, source%v1, nold, nsrc, lsource_mask) call util_append(self%v2, source%v2, nold, nsrc, lsource_mask) - call util_append(self%Gmass1, source%Gmass1, nold, nsrc, lsource_mask) - call util_append(self%Gmass2, source%Gmass2, nold, nsrc, lsource_mask) - call util_append(self%radius1, source%radius1, nold, nsrc, lsource_mask) - call util_append(self%radius2, source%radius2, nold, nsrc, lsource_mask) - call util_append(self%name1, source%name1, nold, nsrc, lsource_mask) - call util_append(self%name2, source%name2, nold, nsrc, lsource_mask) self%nenc = nold + count(lsource_mask(1:nsrc)) return @@ -70,12 +64,6 @@ module subroutine encounter_util_copy_list(self, source) self%x2(:,1:n) = source%x2(:,1:n) self%v1(:,1:n) = source%v1(:,1:n) self%v2(:,1:n) = source%v2(:,1:n) - self%Gmass1(1:n) = source%Gmass1(1:n) - self%Gmass2(1:n) = source%Gmass2(1:n) - self%radius1(1:n) = source%radius1(1:n) - self%radius2(1:n) = source%radius2(1:n) - self%name1(1:n) = source%name1(1:n) - self%name2(1:n) = source%name2(1:n) end associate return @@ -115,12 +103,6 @@ module subroutine encounter_util_dealloc_list(self) if (allocated(self%x2)) deallocate(self%x2) if (allocated(self%v1)) deallocate(self%v1) if (allocated(self%v2)) deallocate(self%v2) - if (allocated(self%Gmass1)) deallocate(self%Gmass1) - if (allocated(self%Gmass2)) deallocate(self%Gmass2) - if (allocated(self%radius1)) deallocate(self%radius1) - if (allocated(self%radius2)) deallocate(self%radius2) - if (allocated(self%name1)) deallocate(self%name1) - if (allocated(self%name2)) deallocate(self%name2) return end subroutine encounter_util_dealloc_list @@ -216,12 +198,6 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive) call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive) call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive) - call util_spill(keeps%Gmass1, discards%Gmass1, lspill_list, ldestructive) - call util_spill(keeps%Gmass2, discards%Gmass2, lspill_list, ldestructive) - call util_spill(keeps%radius1, discards%radius1, lspill_list, ldestructive) - call util_spill(keeps%radius2, discards%radius2, lspill_list, ldestructive) - call util_spill(keeps%name1, discards%name1, lspill_list, ldestructive) - call util_spill(keeps%name2, discards%name2, lspill_list, ldestructive) nenc_old = keeps%nenc diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index ab8dd0f11..3b775370e 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -31,12 +31,6 @@ module encounter_classes real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter - real(DP), dimension(:), allocatable :: Gmass1 !! G*mass of body 1 in the encounter - real(DP), dimension(:), allocatable :: Gmass2 !! G*mass of body 2 in the encounter - real(DP), dimension(:), allocatable :: radius1 !! radius of body 1 in the encounter - real(DP), dimension(:), allocatable :: radius2 !! radius of body 2 in the encounter - character(NAMELEN), dimension(:), allocatable :: name1 !! name body 1 in the encounter - character(NAMELEN), dimension(:), allocatable :: name2 !! name of body 2 in the encounter contains procedure :: setup => encounter_setup_list !! A constructor that sets the number of encounters and allocates and initializes all arrays procedure :: append => encounter_util_append_list !! Appends elements from one structure to another @@ -217,10 +211,10 @@ module subroutine encounter_io_initialize_output(self, param) class(swiftest_parameters), intent(in) :: param end subroutine encounter_io_initialize_output - module subroutine encounter_io_write_frame(self, iu, param) + module subroutine encounter_io_write_frame(self, nciu, param) implicit none class(encounter_list), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: iu !! Parameters used to identify a particular encounter io NetCDF dataset + class(encounter_io_parameters), intent(inout) :: nciu !! Parameters used to identify a particular encounter io NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine encounter_io_write_frame diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 621287433..e01392599 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -726,12 +726,6 @@ module subroutine symba_util_rearray_pl(self, system, param) ! This is an encounter we already know about, so save the old information system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) system%plplenc_list%status(k) = plplenc_old%status(k) - system%plplenc_list%Gmass1(k) = plplenc_old%Gmass1(k) - system%plplenc_list%Gmass2(k) = plplenc_old%Gmass2(k) - system%plplenc_list%radius1(k) = plplenc_old%radius1(k) - system%plplenc_list%radius2(k) = plplenc_old%radius2(k) - system%plplenc_list%name1(k) = plplenc_old%name1(k) - system%plplenc_list%name2(k) = plplenc_old%name2(k) system%plplenc_list%x1(:,k) = plplenc_old%x1(:,k) system%plplenc_list%x2(:,k) = plplenc_old%x2(:,k) system%plplenc_list%v1(:,k) = plplenc_old%v1(:,k) @@ -742,12 +736,6 @@ module subroutine symba_util_rearray_pl(self, system, param) ! This is an encounter we already know about, but with the order reversed, so save the old information system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) system%plplenc_list%status(k) = plplenc_old%status(k) - system%plplenc_list%Gmass1(k) = plplenc_old%Gmass2(k) - system%plplenc_list%Gmass2(k) = plplenc_old%Gmass1(k) - system%plplenc_list%radius1(k) = plplenc_old%radius2(k) - system%plplenc_list%radius2(k) = plplenc_old%radius1(k) - system%plplenc_list%name1(k) = plplenc_old%name2(k) - system%plplenc_list%name2(k) = plplenc_old%name1(k) system%plplenc_list%x1(:,k) = plplenc_old%x2(:,k) system%plplenc_list%x2(:,k) = plplenc_old%x1(:,k) system%plplenc_list%v1(:,k) = plplenc_old%v2(:,k) @@ -775,12 +763,6 @@ module subroutine symba_util_rearray_pl(self, system, param) system%plplenc_list%id2(1:nencmin) = pack(system%plplenc_list%id2(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%lvdotr(1:nencmin) = pack(system%plplenc_list%lvdotr(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%status(1:nencmin) = pack(system%plplenc_list%status(1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%Gmass1(1:nencmin) = pack(system%plplenc_list%Gmass1(1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%Gmass2(1:nencmin) = pack(system%plplenc_list%Gmass2(1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%radius1(1:nencmin) = pack(system%plplenc_list%radius1(1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%radius2(1:nencmin) = pack(system%plplenc_list%radius2(1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%name1(1:nencmin) = pack(system%plplenc_list%name1(1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%name2(1:nencmin) = pack(system%plplenc_list%name2(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%tcollision(1:nencmin) = pack(system%plplenc_list%tcollision(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%level(1:nencmin) = pack(system%plplenc_list%level(1:nenc_old), lmask(1:nenc_old)) do i = 1, NDIM From c659ca1e7b4f17c1742a79abef5e827ea0090ad9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 18:12:03 -0500 Subject: [PATCH 3/8] Added new symba_system_snapshot object to track the system through encounters --- src/encounter/encounter_io.f90 | 2 +- src/modules/symba_classes.f90 | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 62920db7e..f41c25f85 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -140,7 +140,7 @@ module subroutine encounter_io_write_frame(self, nciu, param) implicit none ! Arguments class(encounter_list), intent(in) :: self !! Swiftest encounter structure - class(encounter_io_parameters), intent(inout) :: nciu !! Parameters used to identify a particular encounter io NetCDF dataset + class(encounter_io_parameters), intent(inout) :: nciu !! Parameters used to identify a particular encounter io NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: tslot,i,old_mode, n diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index cd97b74bd..d72097834 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -201,7 +201,14 @@ module symba_classes final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables end type symba_nbody_system + type, extends(helio_nbody_system) :: symba_system_snapshot + contains + procedure :: snapshot => symba_util_take_system_snapshot + final :: symba_util_final_snapshot + end type + interface + module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) use swiftest_classes, only : swiftest_parameters implicit none @@ -374,6 +381,15 @@ module subroutine symba_util_set_renc(self, scale) integer(I4B), intent(in) :: scale !! Current recursion depth end subroutine symba_util_set_renc + module subroutine symba_util_take_system_snapshot(self, system, param, t) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system snapshot object + class(symba_nbody_system), intent(in) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time + end subroutine symba_util_take_system_snapshot + module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none class(symba_parameters), intent(inout) :: self !! Current run configuration parameters with SyMBA additionss @@ -648,6 +664,11 @@ module subroutine symba_util_final_pl(self) type(symba_pl), intent(inout) :: self !! SyMBA massive body object end subroutine symba_util_final_pl + module subroutine symba_util_final_snapshot(self) + implicit none + type(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system object + end subroutine symba_util_final_snapshot + module subroutine symba_util_final_system(self) implicit none type(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object From f1502c059f188c0f88ad1c47c715b0a7d3ab3308 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 18:15:50 -0500 Subject: [PATCH 4/8] Added the symba_util_take_system_snapshot implementation interface --- src/symba/symba_util.f90 | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index e01392599..3fd00efc8 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -499,6 +499,20 @@ module subroutine symba_util_final_system(self) end subroutine symba_util_final_system + module subroutine symba_util_final_snapshot(self) + !! author: David A. Minton + !! + !! Finalize the SyMBA nbody system object - deallocates all allocatables + implicit none + ! Argument + type(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system object + + call self%dealloc() + + return + end subroutine symba_util_final_snapshot + + module subroutine symba_util_final_tp(self) !! author: David A. Minton !! @@ -1279,4 +1293,19 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) return end subroutine symba_util_spill_tp + + module subroutine symba_util_take_system_snapshot(self, system, param, t) + !! author: David A. Minton + !! + !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories + !! Can be played back through the encounter + implicit none + class(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system snapshot object + class(symba_nbody_system), intent(in) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time + + return + end subroutine symba_util_take_system_snapshot + end submodule s_symba_util From 71b7f9940fbf604d4adae0ddf6611f19f22cba7e Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 21:02:24 -0500 Subject: [PATCH 5/8] Lingering issues with THE SPACE DIMENSION --- python/swiftest/swiftest/init_cond.py | 5 +- python/swiftest/swiftest/io.py | 84 +++++++------------- python/swiftest/swiftest/simulation_class.py | 1 + 3 files changed, 33 insertions(+), 57 deletions(-) diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index b24eff0d8..8dcf08cf2 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -284,10 +284,9 @@ def vec2xr(param: Dict, **kwargs: Any): kwargs = {k:kwargs[k] for k,v in kwargs.items() if v is not None} if param['ROTATION']: if "rot" not in kwargs and "Gmass" in kwargs: - warnings.warn("Rotation vectors must be given when rotation is enabled for massive bodies",stacklevel=2) - return + kwargs['rot'] = np.zeros((len(kwargs['Gmass']),3)) if "Ip" not in kwargs and "rot" in kwargs: - kwargs['Ip'] = np.full_like(rot, 0.4) + kwargs['Ip'] = np.full_like(kwargs['rot'], 0.4) if "time" not in kwargs: kwargs["time"] = np.array([0.0]) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index dbcb7430d..40ce7ae33 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -529,12 +529,8 @@ def swifter_stream(f, param): tlab = [] if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL': - tlab.append('rhx') - tlab.append('rhy') - tlab.append('rhz') - tlab.append('vhx') - tlab.append('vhy') - tlab.append('vhz') + tlab.append('rh') + tlab.append('vh') if param['OUT_FORM'] == 'EL' or param['OUT_FORM'] == 'XVEL': tlab.append('a') tlab.append('e') @@ -577,12 +573,8 @@ def make_swiftest_labels(param): """ tlab = [] if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL': - tlab.append('rhx') - tlab.append('rhy') - tlab.append('rhz') - tlab.append('vhx') - tlab.append('vhy') - tlab.append('vhz') + tlab.append('rh') + tlab.append('vh') if param['OUT_FORM'] == 'EL' or param['OUT_FORM'] == 'XVEL': tlab.append('a') @@ -599,18 +591,10 @@ def make_swiftest_labels(param): plab.append('rhill') clab = ['Gmass', 'radius', 'j2rp2', 'j4rp4'] if param['ROTATION']: - clab.append('Ip1') - clab.append('Ip2') - clab.append('Ip3') - clab.append('rotx') - clab.append('roty') - clab.append('rotz') - plab.append('Ip1') - plab.append('Ip2') - plab.append('Ip3') - plab.append('rotx') - plab.append('roty') - plab.append('rotz') + clab.append('Ip') + clab.append('rot') + plab.append('Ip') + plab.append('rot') if param['TIDES']: clab.append('k2') clab.append('Q') @@ -619,19 +603,11 @@ def make_swiftest_labels(param): infolab_float = [ "origin_time", - "origin_rhx", - "origin_rhy", - "origin_rhz", - "origin_vhx", - "origin_vhy", - "origin_vhz", + "origin_rh", + "origin_vh", "discard_time", - "discard_rhx", - "discard_rhy", - "discard_rhz", - "discard_vhx", - "discard_vhy", - "discard_vhz", + "discard_rh", + "discard_vh", ] infolab_int = [ "collision_id", @@ -1042,7 +1018,7 @@ def swiftest_particle_2xr(param): ------- infoxr : xarray dataset """ - veclab = ['time_origin', 'rhx_origin', 'py_origin', 'pz_origin', 'vhx_origin', 'vhy_origin', 'vhz_origin'] + veclab = ['time_origin', 'rh_origin', 'vh_origin'] id_list = [] origin_type_list = [] origin_vec_list = [] @@ -1099,7 +1075,7 @@ def select_active_from_frame(ds, param, framenum=-1): # Select only the active particles at this time step # Remove the inactive particles if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL': - iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['rhx'])), drop=True)[count_dim] + iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['rh'].isel(space=0))), drop=True)[count_dim] else: iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['a'])), drop = True)[count_dim] if count_dim == "id": @@ -1164,12 +1140,12 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram J4 = np.double(cb['j4rp4']) cbname = cb['name'].values[0] if param['ROTATION']: - Ip1cb = np.double(cb['Ip1']) - Ip2cb = np.double(cb['Ip2']) - Ip3cb = np.double(cb['Ip3']) - rotxcb = np.double(cb['rotx']) - rotycb = np.double(cb['roty']) - rotzcb = np.double(cb['rotz']) + Ip1cb = np.double(cb['Ip'].values[0]) + Ip2cb = np.double(cb['Ip'].values[1]) + Ip3cb = np.double(cb['Ip'].values[2]) + rotxcb = np.double(cb['rot'].values[0]) + rotycb = np.double(cb['rot'].values[1]) + rotzcb = np.double(cb['rot'].values[2]) cbid = int(0) if in_type == 'ASCII': @@ -1196,16 +1172,16 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram if param['CHK_CLOSE']: print(pli['radius'].values[0], file=plfile) if param['IN_FORM'] == 'XV': - print(pli['rhx'].values[0], pli['rhy'].values[0], pli['rhz'].values[0], file=plfile) - print(pli['vhx'].values[0], pli['vhy'].values[0], pli['vhz'].values[0], file=plfile) + print(pli['rh'].values[0,0], pli['rh'].values[0,1], pli['rh'].values[0,2], file=plfile) + print(pli['vh'].values[0,0], pli['vh'].values[0,1], pli['vh'].values[0,2], file=plfile) elif param['IN_FORM'] == 'EL': print(pli['a'].values[0], pli['e'].values[0], pli['inc'].values[0], file=plfile) print(pli['capom'].values[0], pli['omega'].values[0], pli['capm'].values[0], file=plfile) else: print(f"{param['IN_FORM']} is not a valid input format type.") if param['ROTATION']: - print(pli['Ip1'].values[0], pli['Ip2'].values[0], pli['Ip3'].values[0], file=plfile) - print(pli['rotx'].values[0], pli['roty'].values[0], pli['rotz'].values[0], file=plfile) + print(pli['Ip'].values[0,0], pli['Ip'].values[0,1], pli['Ip'].values[0,2], file=plfile) + print(pli['rot'].values[0,0], pli['rot'].values[0,1], pli['rot'].values[0,2], file=plfile) plfile.close() # TP file @@ -1215,8 +1191,8 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram tpi = tp.sel(id=i) print(tpi['name'].values[0], file=tpfile) if param['IN_FORM'] == 'XV': - print(tpi['rhx'].values[0], tpi['rhy'].values[0], tpi['rhz'].values[0], file=tpfile) - print(tpi['vhx'].values[0], tpi['vhy'].values[0], tpi['vhz'].values[0], file=tpfile) + print(tpi['rh'].values[0,0], tpi['rh'].values[0,1], tpi['rh'].values[0,2], file=tpfile) + print(tpi['vh'].values[0,0], tpi['vh'].values[0,1], tpi['vh'].values[0,2], file=tpfile) elif param['IN_FORM'] == 'EL': print(tpi['a'].values[0], tpi['e'].values[0], tpi['inc'].values[0], file=tpfile) print(tpi['capom'].values[0], tpi['omega'].values[0], tpi['capm'].values[0], file=tpfile) @@ -1279,8 +1255,8 @@ def swifter_xr2infile(ds, param, framenum=-1): print(i.values, pli['Gmass'].values, file=plfile) if param['CHK_CLOSE']: print(pli['radius'].values, file=plfile) - print(pli['rhx'].values, pli['rhy'].values, pli['rhz'].values, file=plfile) - print(pli['vhx'].values, pli['vhy'].values, pli['vhz'].values, file=plfile) + print(pli['rh'].values[0,0], pli['ry'].values[0,1], pli['rh'].values[0,2], file=plfile) + print(pli['vh'].values[0,0], pli['vh'].values[0,1], pli['vh'].values[0,2], file=plfile) plfile.close() # TP file @@ -1289,8 +1265,8 @@ def swifter_xr2infile(ds, param, framenum=-1): for i in tp.id: tpi = tp.sel(id=i) print(i.values, file=tpfile) - print(tpi['rhx'].values, tpi['rhy'].values, tpi['rhz'].values, file=tpfile) - print(tpi['vhx'].values, tpi['vhy'].values, tpi['vhz'].values, file=tpfile) + print(tpi['rh'].values[0,0], tpi['ry'].values[0,1], tpi['rh'].values[0,2], file=tpfile) + print(tpi['vh'].values[0,0], tpi['vh'].values[0,1], tpi['vh'].values[0,2], file=tpfile) tpfile.close() else: # Now make Swiftest files diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 22f2001ac..6cd75a521 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -2811,6 +2811,7 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil A dataset containing the extracted initial condition data. """ + frame = None if codename != "Swiftest": self.save(new_param_file, framenum, codename) return From 086872160622bbddb94422dc1e425e5022684488 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 21:03:04 -0500 Subject: [PATCH 6/8] More updates for preparing to output encounter data --- examples/Fragmentation/Fragmentation_Movie.py | 13 ++++++------- src/modules/symba_classes.f90 | 12 +++++++----- src/setup/setup.f90 | 2 ++ src/symba/symba_step.f90 | 2 ++ src/symba/symba_util.f90 | 5 +++++ 5 files changed, 22 insertions(+), 12 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 24fc40e30..fc308687a 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -36,7 +36,7 @@ # Define the names and initial conditions of the various fragmentation simulation types # ---------------------------------------------------------------------------------------------------------------------- available_movie_styles = ["disruption_headon", "supercatastrophic_off_axis", "hitandrun"] -movie_title_list = ["Head-on Disrutption", "Off-axis Supercatastrophic", "Hit and Run"] +movie_title_list = ["Head-on Disruption", "Off-axis Supercatastrophic", "Hit and Run"] movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error @@ -135,8 +135,8 @@ def center(Gmass, x, y): y_com = np.sum(Gmass * y) / np.sum(Gmass) return x_com, y_com - Gmass, x, y, point_rad = next(self.data_stream(frame)) - x_com, y_com = center(Gmass, x, y) + Gmass, rh, point_rad = next(self.data_stream(frame)) + x_com, y_com = center(Gmass, rh[:,1], rh[:,2]) self.scatter_artist.set_offsets(np.c_[x - x_com, y - y_com]) self.scatter_artist.set_sizes(point_rad) return self.scatter_artist, @@ -147,10 +147,9 @@ def data_stream(self, frame=0): ds = ds.where(ds['name'] != "Sun", drop=True) radius = ds['radius'].values Gmass = ds['Gmass'].values - x = ds['xhx'].values - y = ds['xhy'].values + rh = ds['rh'].values point_rad = 2 * radius * self.ax_pt_size - yield Gmass, x, y, point_rad + yield Gmass, rh, point_rad if __name__ == "__main__": @@ -188,7 +187,7 @@ def data_stream(self, frame=0): if run_new: sim = swiftest.Simulation(param_file=param_file, rotation=True, init_cond_format = "XV", compute_conservation_values=True) sim.add_solar_system_body("Sun") - sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style], rot=rot_vectors[style]) + sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style]) #, rot=rot_vectors[style]) # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index d72097834..432126e4e 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -177,6 +177,12 @@ module symba_classes procedure :: resolve_collision => symba_collision_resolve_plplenc !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the c end type symba_plplenc + type, extends(helio_nbody_system) :: symba_system_snapshot + contains + procedure :: snapshot => symba_util_take_system_snapshot + final :: symba_util_final_snapshot + end type + !******************************************************************************************************************************** ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** @@ -188,6 +194,7 @@ module symba_classes integer(I4B) :: irec !! System recursion level type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file integer(I4B) :: ienc_frame = 0 !! Encounter history frame number + type(symba_system_snapshot) :: snapshot contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps @@ -201,11 +208,6 @@ module symba_classes final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables end type symba_nbody_system - type, extends(helio_nbody_system) :: symba_system_snapshot - contains - procedure :: snapshot => symba_util_take_system_snapshot - final :: symba_util_final_snapshot - end type interface diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 26aed237c..4a3d98aed 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -68,6 +68,8 @@ module subroutine setup_construct_system(system, param) allocate(symba_pltpenc :: system%pltpenc_list) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_plplenc :: system%plplcollision_list) + allocate(symba_pl :: system%snapshot%pl) + allocate(symba_tp :: system%snapshot%tp) end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index ecf1874f6..9e9d2de41 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -218,6 +218,8 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) call pl%drift(system, param, dtl) call tp%drift(system, param, dtl) + !call system% + if (lencounter) call system%recursive_step(param, t+dth,irecp) system%irec = ireci diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 3fd00efc8..63ba3224b 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -1300,10 +1300,15 @@ module subroutine symba_util_take_system_snapshot(self, system, param, t) !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories !! Can be played back through the encounter implicit none + ! Internals class(symba_system_snapshot), intent(inout) :: self !! SyMBA nbody system snapshot object class(symba_nbody_system), intent(in) :: system !! SyMBA nbody system object class(symba_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! current time + ! Arguments + logical, dimension(:), allocatable :: lmask + + !if (system%pl) return end subroutine symba_util_take_system_snapshot From f4fba0e972ebdc9e358f362a66e7e5f006d0a15d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 21:07:15 -0500 Subject: [PATCH 7/8] Fixed input file for fragmentation movie --- examples/Fragmentation/Fragmentation_Movie.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index fc308687a..5d660bcec 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -167,7 +167,6 @@ def data_stream(self, frame=0): movie_styles = available_movie_styles.copy() for style in movie_styles: - param_file = Path(style) / "param.in" bin_file = Path(style) / "bin.nc" if bin_file.exists(): user_selection = input(f"An older simulation of {movie_titles[style]} exists. Do you want to re-run it? [y/N] ") @@ -185,7 +184,7 @@ def data_stream(self, frame=0): # Pull in the Swiftest output data from the parameter file and store it as a Xarray dataset. if run_new: - sim = swiftest.Simulation(param_file=param_file, rotation=True, init_cond_format = "XV", compute_conservation_values=True) + sim = swiftest.Simulation(simdir=style, rotation=True, init_cond_format = "XV", compute_conservation_values=True) sim.add_solar_system_body("Sun") sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style]) #, rot=rot_vectors[style]) From e8f01d59101710e67cd61db1c3c84594c545eec0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 4 Dec 2022 22:46:00 -0500 Subject: [PATCH 8/8] Fixed some bad loop indexing --- src/io/io.f90 | 4 ++-- src/main/swiftest_driver.f90 | 2 +- src/modules/swiftest_classes.f90 | 2 +- src/netcdf/netcdf.f90 | 8 ++++---- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 81aca06d1..69c32f5dc 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -256,7 +256,7 @@ module subroutine io_dump_system(self, param) dump_param%out_form = "XV" dump_param%outfile = trim(adjustl(DUMP_NC_FILE(idx))) - dump_param%ioutput = 0 + dump_param%ioutput = 1 call dump_param%nciu%initialize(dump_param) call self%write_frame(dump_param%nciu, dump_param) call dump_param%nciu%close() @@ -285,7 +285,7 @@ module subroutine io_dump_storage(self, param) integer(I4B) :: i integer(I8B) :: iloop_start - iloop_start = param%iloop - int(param%istep_out * param%dump_cadence + 1, kind=I8B) + iloop_start = max(param%iloop - int(param%istep_out * param%dump_cadence, kind=I8B),1) do i = 1, param%dump_cadence param%ioutput = int(iloop_start / param%istep_out, kind=I4B) + i if (allocated(self%frame(i)%item)) then diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 846915444..528e0c73d 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -80,7 +80,7 @@ program swiftest_driver ! Set up loop and output cadence variables t = tstart nloops = ceiling((tstop - t0) / dt, kind=I8B) - istart = ceiling((tstart - t0) / dt + 1, kind=I8B) + istart = ceiling((tstart - t0) / dt + 1.0_DP, kind=I8B) ioutput = int(istart / istep_out, kind=I4B) ! Set up system storage for intermittent file dumps diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 0eb932079..4f87b7707 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -159,7 +159,7 @@ module swiftest_classes real(DP) :: tstop = -1.0_DP !! Integration stop time real(DP) :: dt = -1.0_DP !! Time step integer(I8B) :: iloop = 0_I8B !! Main loop counter - integer(I4B) :: ioutput = 0 !! Output counter + integer(I4B) :: ioutput = 1 !! Output counter character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index b790c7f75..261fa6355 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -456,7 +456,7 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr) call pl%setup(npl, param) call tp%setup(ntp, param) - tslot = param%ioutput + 1 + tslot = param%ioutput call check( nf90_inquire_dimension(nciu%id, nciu%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid" ) allocate(rtemp(idmax)) @@ -699,7 +699,7 @@ module subroutine netcdf_read_hdr_system(self, nciu, param) logical, dimension(:), allocatable :: plmask, tpmask, plmmask - tslot = param%ioutput + 1 + tslot = param%ioutput call check( nf90_inquire_dimension(nciu%id, nciu%id_dimid, len=idmax), "netcdf_read_hdr_system nf90_inquire_dimension id_dimid" ) call check( nf90_get_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" ) @@ -1031,7 +1031,7 @@ module subroutine netcdf_write_frame_base(self, nciu, param) call self%write_info(nciu, param) - tslot = param%ioutput + 1 + tslot = param%ioutput call check( nf90_set_fill(nciu%id, nf90_nofill, old_mode), "netcdf_write_frame_base nf90_set_fill" ) select type(self) @@ -1242,7 +1242,7 @@ module subroutine netcdf_write_hdr_system(self, nciu, param) ! Internals integer(I4B) :: tslot - tslot = param%ioutput + 1 + tslot = param%ioutput call check( nf90_put_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid" ) call check( nf90_put_var(nciu%id, nciu%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var npl_varid" )