From 19be7d8a10ea54e2e684fcddb91d59b73687f175 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 2 Aug 2021 16:43:52 -0400 Subject: [PATCH] Cleaned up discard io routines and started adding SyMBA-specific discard writer --- .../1pl_1pl_encounter/cb.swiftest.in | Bin 53 -> 80 bytes .../1pl_1pl_encounter/init_cond.py | 1 + .../1pl_1pl_encounter/param.swiftest.in | 3 +- .../1pl_1pl_encounter/pl.swifter.in | 2 +- .../1pl_1pl_encounter/pl.swiftest.in | Bin 228 -> 256 bytes .../1pl_1pl_encounter/tp.swiftest.in | Bin 2 -> 16 bytes python/swiftest/swiftest/io.py | 12 +++- src/io/io.f90 | 61 +++++++++------- src/modules/swiftest_classes.f90 | 11 +-- src/modules/symba_classes.f90 | 63 ++++++++++------- src/rmvs/rmvs_step.f90 | 4 +- src/symba/symba_io.f90 | 66 +++++++++++++++++- 12 files changed, 156 insertions(+), 67 deletions(-) diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/cb.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/cb.swiftest.in index 4c5d870405c9daad1d597cb1d3f2bd78a1b2227e..d0ae0ed15fe3ea8dd15557055a926fce3c60b59c 100644 GIT binary patch literal 80 ncmd;JKmZOP6NHU2HoW29>+AsI-}OJ>6US3*597mhVB-S-U7iOf literal 53 wcmW;Axe)*$3rbo3*74~+Eq5}gct*a?m%+)WyI1?iYw*UYD diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py index 7600320c2..20be5a433 100755 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py @@ -173,6 +173,7 @@ print(f'ENC_OUT {swiftest_enc}') print(f'EXTRA_FORCE no') print(f'BIG_DISCARD no') +print(f'DISCARD_OUT discard.swiftest.out') print(f'ROTATION no') print(f'GR no') print(f'MU2KG {MU2KG}') diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in index 1866557b2..d44f4df0e 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in @@ -5,7 +5,7 @@ DT 0.0006844626967830253 CB_IN cb.swiftest.in PL_IN pl.swiftest.in TP_IN tp.swiftest.in -IN_TYPE ASCII +IN_TYPE REAL8 ISTEP_OUT 1 ISTEP_DUMP 1 BIN_OUT bin.swiftest.dat @@ -22,6 +22,7 @@ CHK_QMIN_RANGE 0.004650467260962157 1000.0 ENC_OUT enc.swiftest.dat EXTRA_FORCE no BIG_DISCARD no +DISCARD_OUT discard.swiftest.out ROTATION no GR no MU2KG 1.988409870698051e+30 diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/pl.swifter.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/pl.swifter.in index 0eb21018b..9f0548fc1 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/pl.swifter.in +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/pl.swifter.in @@ -1,5 +1,5 @@ 3 ! Planet input file generated using init_cond.py -1 39.47692640889762629 +1 39.476926408897625196 0.0 0.0 0.0 0.0 0.0 0.0 2 0.00012002693582795244940133 0.010044724833237892 diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/pl.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/pl.swiftest.in index 19c6d6e3a2436162bb5984e0ecb1cb9352cb223f..d8da7a92a44b1e9caa3907ead959cdec31e066cc 100644 GIT binary patch literal 256 zcmd;JU|?VZVi4c}VgVqA@l!y8KmZa0VF>tOuNl*S=&QyDdsK0lJi2<~#U*rILVhbs zIyBlY7vp7;bRcBDlutC@{W5v`KZ2`e<&?MB!PKu_Vy_x9sl|Sa{`cIZU5Rja3YkwR SWDH@mKPFMv^xXC_SUmurogOa$ literal 228 zcmZ9Gxe)^~30`5m% zX0lVQ37Yye{^|yd{X~mu+c^*|K`` Ygvg?I8#T=l^4}PjYzRwv=DU9V0CM9iPyhe` diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/tp.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/tp.swiftest.in index 573541ac9702dd3969c9bc859d2b91ec1f7e6e56..64bf92f74a457d2f4bc42798493db15cc3ab1008 100644 GIT binary patch literal 16 Ncmd;JKmZOP6953P01*HH literal 2 JcmXru0ssJP06PEx diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 2dd4ef7b3..5782c6444 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -5,7 +5,7 @@ import sys import tempfile -newfeaturelist = ("FRAGMENTATION", "ROTATION", "TIDES", "ENERGY", "GR", "YARKOVSKY", "YORP" ) +newfeaturelist = ("FRAGMENTATION", "ROTATION", "TIDES", "ENERGY", "GR", "YARKOVSKY", "YORP") def real2float(realstr): """ @@ -279,6 +279,7 @@ def write_labeled_param(param, param_file_name): 'CB_IN', 'BIN_OUT', 'ENC_OUT', + 'DISCARD_OUT', 'CHK_QMIN', 'CHK_RMIN', 'CHK_RMAX', @@ -889,7 +890,7 @@ def swift2swifter(swift_param, plname="", tpname="", conversion_questions={}): swifter_param['ENC_OUT'] = input("ENC_OUT: Encounter file name: [enc.dat]> ") if swifter_param['ENC_OUT'] == '': swifter_param['ENC_OUT'] = "enc.dat" - + intxt = conversion_questions.get('EXTRA_FORCE', None) if not intxt: intxt = input("EXTRA_FORCE: Use additional user-specified force routines? (y/N)> ") @@ -1228,6 +1229,13 @@ def swifter2swiftest(swifter_param, plname="", tpname="", cbname="", conversion_ swiftest_param.pop('J2', None) swiftest_param.pop('J4', None) swiftest_param.pop('RHILL_PRESENT', None) + + swiftest_param['DISCARD_OUT'] = conversion_questions.get('DISCARD_OUT', '') + if not swiftest_param['DISCARD_OUT']: + swiftest_param['DISCARD_OUT'] = input("DISCARD_OUT: Discard file name: [discard.out]> ") + if swiftest_param['DISCARD_OUT'] == '': + swiftest_param['DISCARD_OUT'] = "discard.out" + swiftest_param['! VERSION'] = "Swiftest parameter file converted from Swifter" return swiftest_param diff --git a/src/io/io.f90 b/src/io/io.f90 index b424094eb..e3556b11c 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -99,7 +99,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param_value = io_get_token(line, ifirst, ilast, iostat) read(param_value, *) self%qmin_ahi case ("ENC_OUT") - self%encounter_file = param_value + self%enc_out = param_value + case ("DISCARD_OUT") + self%discard_out = param_value case ("EXTRA_FORCE") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') self%lextra_force = .true. @@ -225,9 +227,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) write(*,*) "CHK_QMIN = ",self%qmin write(*,*) "CHK_QMIN_COORD = ",trim(adjustl(self%qmin_coord)) write(*,*) "CHK_QMIN_RANGE = ",self%qmin_alo, self%qmin_ahi - write(*,*) "ENC_OUT = ",trim(adjustl(self%encounter_file)) write(*,*) "EXTRA_FORCE = ",self%lextra_force - write(*,*) "BIG_DISCARD = ",self%lbig_discard write(*,*) "RHILL_PRESENT = ",self%lrhill_present write(*,*) "ROTATION = ", self%lrotation write(*,*) "TIDES = ", self%ltides @@ -235,6 +235,18 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) write(*,*) "MU2KG = ",self%MU2KG write(*,*) "TU2S = ",self%TU2S write(*,*) "DU2M = ",self%DU2M + if (trim(adjustl(self%enc_out)) /= "") then + write(*,*) "ENC_OUT = ",trim(adjustl(self%enc_out)) + else + write(*,*) "! ENC_OUT not set: Encounters will not be recorded to file" + end if + if (trim(adjustl(self%discard_out)) /= "") then + write(*,*) "DISCARD_OUT = ",trim(adjustl(self%discard_out)) + write(*,*) "BIG_DISCARD = ",self%lbig_discard + else + write(*,*) "! DISCARD_OUT not set: Discards will not be recorded to file" + write(*,*) "! BIG_DISCARD = ",self%lbig_discard + end if if ((self%MU2KG < 0.0_DP) .or. (self%TU2S < 0.0_DP) .or. (self%DU2M < 0.0_DP)) then write(iomsg,*) 'Invalid unit conversion factor' @@ -317,7 +329,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "OUT_FORM"; write(param_value, Afmt) trim(adjustl(param%out_form)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "OUT_STAT"; write(param_value, Afmt) "APPEND"; write(unit, Afmt) adjustl(param_name), adjustl(param_value) end if - write(param_name, Afmt) "ENC_OUT"; write(param_value, Afmt) trim(adjustl(param%encounter_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "ENC_OUT"; write(param_value, Afmt) trim(adjustl(param%enc_out)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) if (param%istep_dump > 0) then write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt) adjustl(param_name), adjustl(param_value) end if @@ -761,7 +773,7 @@ end subroutine io_read_param_in function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, encounter_file, out_type) result(ierr) + xh1, xh2, vh1, vh2, enc_out, out_type) result(ierr) !! author: David A. Minton !! !! Read close encounter data from input binary files @@ -773,7 +785,7 @@ function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & integer(I4B), intent(out) :: name1, name2 real(DP), intent(out) :: t, mass1, mass2, radius1, radius2 real(DP), dimension(:), intent(out) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file, out_type + character(*), intent(in) :: enc_out, out_type ! Result integer(I4B) :: ierr ! Internals @@ -782,7 +794,7 @@ function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & integer(I4B), save :: iu = lun if (lfirst) then - open(unit = iu, file = encounter_file, status = 'OLD', form = 'UNFORMATTED', iostat = ierr) + open(unit = iu, file = enc_out, status = 'OLD', form = 'UNFORMATTED', iostat = ierr) if (ierr /= 0) then write(*, *) "Swiftest Error:" write(*, *) " unable to open binary encounter file" @@ -1046,31 +1058,27 @@ module subroutine io_write_discard(self, param) character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp - associate(t => param%t, discards => self%tp_discards, nsp => self%tp_discards%nbody, dxh => self%tp_discards%xh, dvh => self%tp_discards%vh, & - dname => self%tp_discards%id, dstatus => self%tp_discards%status) - + associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody) + if (nsp == 0) return select case(param%out_stat) case('APPEND') - open(unit = LUN, file = param%outfile, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = LUN, file = param%outfile, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr) + open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr) case default write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) call util_exit(FAILURE) end select lfirst = .false. - if (param%lgr) call discards%pv2v(param) + if (param%lgr) call tp_discards%pv2v(param) - write(LUN, HDRFMT) t, nsp, param%lbig_discard + write(LUN, HDRFMT) param%t, nsp, param%lbig_discard do i = 1, nsp - write(LUN, NAMEFMT) sub, dname(i), dstatus(i) - write(LUN, VECFMT) dxh(1, i), dxh(2, i), dxh(3, i) - write(LUN, VECFMT) dvh(1, i), dvh(2, i), dvh(3, i) + write(LUN, NAMEFMT) sub, tp_discards%id(i), tp_discards%status(i) + write(LUN, VECFMT) tp_discards%xh(1, i), tp_discards%xh(2, i), tp_discards%xh(3, i) + write(LUN, VECFMT) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i) end do if (param%lbig_discard) then - associate(npl => self%pl%nbody, pl => self%pl, GMpl => self%pl%Gmass, & - Rpl => self%pl%radius, name => self%pl%id, xh => self%pl%xh) - if (param%lgr) then allocate(pltemp, source = pl) call pltemp%pv2v(param) @@ -1082,12 +1090,11 @@ module subroutine io_write_discard(self, param) write(LUN, NPLFMT) npl do i = 1, npl - write(LUN, PLNAMEFMT) name(i), GMpl(i), Rpl(i) - write(LUN, VECFMT) xh(1, i), xh(2, i), xh(3, i) + write(LUN, PLNAMEFMT) pl%id(i), pl%Gmass(i), pl%radius(i) + write(LUN, VECFMT) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i) write(LUN, VECFMT) vh(1, i), vh(2, i), vh(3, i) end do deallocate(vh) - end associate end if close(LUN) end associate @@ -1097,7 +1104,7 @@ end subroutine io_write_discard module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, encounter_file, out_type) + xh1, xh2, vh1, vh2, enc_out, out_type) !! author: David A. Minton !! !! Write close encounter data to output binary files @@ -1110,16 +1117,16 @@ module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, rad integer(I4B), intent(in) :: name1, name2 real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file, out_type + character(*), intent(in) :: enc_out, out_type ! Internals logical , save :: lfirst = .true. integer(I4B), parameter :: lun = 30 integer(I4B) :: ierr integer(I4B), save :: iu = lun - open(unit = iu, file = encounter_file, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + open(unit = iu, file = enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) if ((ierr /= 0) .and. lfirst) then - open(unit = iu, file = encounter_file, status = 'NEW', form = 'UNFORMATTED', iostat = ierr) + open(unit = iu, file = enc_out, status = 'NEW', form = 'UNFORMATTED', iostat = ierr) end if if (ierr /= 0) then write(*, *) "Swiftest Error:" diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index c7a7939a1..83514a2ab 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -38,7 +38,8 @@ module swiftest_classes character(STRMAX) :: qmin_coord = 'HELIO' !! Coordinate frame to use for qmin real(DP) :: qmin_alo = -1.0_DP !! Minimum semimajor axis for qmin real(DP) :: qmin_ahi = -1.0_DP !! Maximum semimajor axis for qmin - character(STRMAX) :: encounter_file = ENC_OUTFILE !! Name of output file for encounters + character(STRMAX) :: enc_out = "" !! Name of output file for encounters + character(STRMAX) :: discard_out = "" !! Name of output file for discards real(QP) :: MU2KG = -1.0_QP !! Converts mass units to grams real(QP) :: TU2S = -1.0_QP !! Converts time units to seconds real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters @@ -277,8 +278,8 @@ module swiftest_classes ! Concrete classes that are common to the basic integrator (only test particles considered for discard) procedure :: discard => discard_system !! Perform a discard step on the system procedure :: dump => io_dump_system !! Dump the state of the system to a file - procedure :: read_frame => io_read_frame_system !! Append a frame of output data to file - procedure :: write_discard => io_write_discard !! Append a frame of output data to file + procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file + procedure :: write_discard => io_write_discard !! Write out information about discarded test particles procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file procedure :: initialize => setup_initialize_system !! Initialize the system from input files procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. @@ -584,12 +585,12 @@ module subroutine io_toupper(string) end subroutine io_toupper module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, encounter_file, out_type) + xh1, xh2, vh1, vh2, enc_out, out_type) implicit none integer(I4B), intent(in) :: name1, name2 real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file, out_type + character(*), intent(in) :: enc_out, out_type end subroutine io_write_encounter module subroutine io_write_frame_body(self, iu, param) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 01af9a48f..2094b533a 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -167,6 +167,7 @@ module symba_classes class(symba_pl), allocatable :: pl_discards !! Discarded test particle data structure integer(I4B) :: irec !! System recursion level 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 procedure :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step procedure :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system @@ -265,34 +266,12 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp - module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters - implicit none - class(symba_pl), intent(inout) :: self !! SyMBA massive body particle data structure - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step - end subroutine symba_kick_getacch_pl - - module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters - implicit none - class(symba_tp), intent(inout) :: self !! SyMBA test particle data structure - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current time - logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step - end subroutine symba_kick_getacch_tp - - module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) + module subroutine symba_io_write_discard(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration - end subroutine symba_kick_pltpenc + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine symba_io_write_discard module subroutine symba_io_dump_particle_info(self, param, msg) use swiftest_classes, only : swiftest_parameters @@ -341,6 +320,36 @@ module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) integer(I4B), intent(out) :: ierr !! Error code end subroutine symba_io_read_frame_info + module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(symba_pl), intent(inout) :: self !! SyMBA massive body particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current simulation time + logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step + end subroutine symba_kick_getacch_pl + + module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(symba_tp), intent(inout) :: self !! SyMBA test particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step + end subroutine symba_kick_getacch_tp + + module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) + implicit none + class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration + end subroutine symba_kick_pltpenc + + module subroutine symba_io_write_frame_info(self, iu, param) use swiftest_classes, only : swiftest_parameters implicit none diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 74e6958c9..0385aeecc 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -539,7 +539,7 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) call orbel_xv2aqt(mu, xpc(:, i), vpc(:, i), a, peri, capm, tperi) r2 = dot_product(xpc(:, i), xpc(:, i)) if ((abs(tperi) > FACQDT * dt) .or. (r2 > rhill2)) peri = sqrt(r2) - if (param%encounter_file /= "") then + if (param%enc_out /= "") then id1 = pl%id(ipleP) rpl = pl%radius(ipleP) xh1(:) = pl%inner(inner_index)%x(:, ipleP) @@ -548,7 +548,7 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) xh2(:) = xpc(:, i) + xh1(:) vh2(:) = xpc(:, i) + vh1(:) call io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, xh1(:), xh2(:), vh1(:), vh2(:), & - param%encounter_file, param%out_type) + param%enc_out, param%out_type) end if if (tp%lperi(i)) then if (peri < tp%peri(i)) then diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 403204017..97e4ab3f6 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -207,7 +207,70 @@ module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) ierr = 0 end subroutine symba_io_read_frame_info - + + + module subroutine symba_io_write_discard(self, param) + implicit none + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I4B), parameter :: LUN = 40 + integer(I4B) :: i, ierr + logical, save :: lfirst = .true. + real(DP), dimension(:,:), allocatable :: vh + character(*), parameter :: HDRFMT = '(E23.16, 1X, I8, 1X, L1)' + character(*), parameter :: NAMEFMT = '(A, 2(1X, I8))' + character(*), parameter :: VECFMT = '(3(E23.16, 1X))' + character(*), parameter :: NPLFMT = '(I8)' + character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' + class(swiftest_body), allocatable :: pltemp + + associate(pl_discards => self%pl_discards, nsppl => self%pl_discards%nbody, pl => self%pl, npl => self%pl%nbody) + if (self%tp_discards%nbody > 0) call io_write_discard(self, param) + + if (nsppl == 0) return + select case(param%out_stat) + case('APPEND') + open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) + case('NEW', 'REPLACE', 'UNKNOWN') + open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr) + case default + write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) + call util_exit(FAILURE) + end select + lfirst = .false. + if (param%lgr) call pl_discards%pv2v(param) + + ! write(LUN, HDRFMT) param%t, nsp, param%lbig_discard + ! do i = 1, nsp + ! write(LUN, NAMEFMT) sub, tp_discards%id(i), tp_discards%status(i) + ! write(LUN, VECFMT) tp_discards%xh(1, i), tp_discards%xh(2, i), tp_discards%xh(3, i) + ! write(LUN, VECFMT) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i) + ! end do + ! if (param%lbig_discard) then + ! if (param%lgr) then + ! allocate(pltemp, source = pl) + ! call pltemp%pv2v(param) + ! allocate(vh, source = pltemp%vh) + ! deallocate(pltemp) + ! else + ! allocate(vh, source = pl%vh) + ! end if + + ! write(LUN, NPLFMT) npl + ! do i = 1, npl + ! write(LUN, PLNAMEFMT) pl%id(i), pl%Gmass(i), pl%radius(i) + ! write(LUN, VECFMT) pl%xh(1, i), pl%xh(2, i), pl%xh(3, i) + ! write(LUN, VECFMT) vh(1, i), vh(2, i), vh(3, i) + ! end do + ! deallocate(vh) + ! end if + ! close(LUN) + end associate + + return + end subroutine symba_io_write_discard + module subroutine symba_io_write_frame_info(self, iu, param) implicit none @@ -216,6 +279,5 @@ module subroutine symba_io_write_frame_info(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_io_write_frame_info - end submodule s_symba_io