From e1861566148329798780d8640b584fc19de45487 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 9 Jul 2021 18:00:36 -0400 Subject: [PATCH] Put together basic SyMBA infrastructure --- Makefile.Defines | 8 +- src/discard/discard.f90 | 10 +- src/helio/helio_drift.f90 | 4 +- src/io/io.f90 | 287 ++++++------------- src/modules/helio_classes.f90 | 4 +- src/modules/rmvs_classes.f90 | 2 +- src/modules/swiftest_classes.f90 | 160 ++++------- src/modules/swiftest_globals.f90 | 6 +- src/modules/symba_classes.f90 | 282 +++++++++++++++++-- src/rmvs/rmvs_discard.f90 | 2 +- src/rmvs/rmvs_step.f90 | 14 +- src/setup/setup.f90 | 12 +- src/symba/symba_io.f90 | 470 +++++++++++++++++++++++++++++++ src/util/util_copy.f90 | 154 ---------- src/util/util_spill_and_fill.f90 | 38 +-- src/util/util_valid.f90 | 4 +- src/whm/whm_drift.f90 | 4 +- 17 files changed, 904 insertions(+), 557 deletions(-) create mode 100644 src/symba/symba_io.f90 delete mode 100644 src/util/util_copy.f90 diff --git a/Makefile.Defines b/Makefile.Defines index 07126f842..1346f4b09 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -65,13 +65,13 @@ GPAR = -fopenmp -ftree-parallelize-loops=4 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries -FFLAGS = $(IDEBUG) $(HEAPARR) +#FFLAGS = $(IDEBUG) $(HEAPARR) #FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) -FORTRAN = ifort +#FORTRAN = ifort #AR = xiar -#FORTRAN = gfortran -#FFLAGS = -ffree-line-length-none $(GDEBUG) #$(GMEM) +FORTRAN = gfortran +FFLAGS = -ffree-line-length-none $(GDEBUG) #$(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 6609ed802..1a1aecea6 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -88,11 +88,11 @@ subroutine discard_sun_tp(tp, system, param) rh2 = dot_product(tp%xh(:, i), tp%xh(:, i)) if ((param%rmax >= 0.0_DP) .and. (rh2 > rmax2)) then tp%status(i) = DISCARDED_RMAX - write(*, *) "Particle ", tp%name(i), " too far from sun at t = ", t + write(*, *) "Particle ", tp%id(i), " too far from sun at t = ", t tp%ldiscard(i) = .true. else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then tp%status(i) = DISCARDED_RMIN - write(*, *) "Particle ", tp%name(i), " too close to sun at t = ", t + write(*, *) "Particle ", tp%id(i), " too close to sun at t = ", t tp%ldiscard(i) = .true. else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) @@ -100,7 +100,7 @@ subroutine discard_sun_tp(tp, system, param) energy = 0.5_DP * vb2 - msys / sqrt(rb2) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then tp%status(i) = DISCARDED_RMAXU - write(*, *) "Particle ", tp%name(i), " is unbound and too far from barycenter at t = ", t + write(*, *) "Particle ", tp%id(i), " is unbound and too far from barycenter at t = ", t tp%ldiscard(i) = .true. end if end if @@ -150,7 +150,7 @@ subroutine discard_peri_tp(tp, system, param) (tp%atp(i) <= param%qmin_ahi) .and. & (tp%peri(i) <= param%qmin)) then tp%status(i) = DISCARDED_PERI - write(*, *) "Particle ", tp%name(i), " perihelion distance too small at t = ", t + write(*, *) "Particle ", tp%id(i), " perihelion distance too small at t = ", t tp%ldiscard(i) = .true. end if end if @@ -191,7 +191,7 @@ subroutine discard_pl_tp(tp, system, param) if (isp /= 0) then tp%status(i) = DISCARDED_PLR pl%ldiscard(j) = .true. - write(*, *) "Particle ", tp%name(i), " too close to massive body ", pl%name(j), " at t = ", t + write(*, *) "Particle ", tp%id(i), " too close to massive body ", pl%id(j), " at t = ", t tp%ldiscard(i) = .true. exit end if diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index ce55797bf..40da379ee 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -47,7 +47,7 @@ module subroutine helio_drift_pl(self, system, param, dt) dtp(1:npl), iflag(1:npl)) if (any(iflag(1:npl) /= 0)) then do i = 1, npl - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" write(*, *) pl%xh(:,i) write(*, *) pl%vb(:,i) write(*, *) " stopping " @@ -136,7 +136,7 @@ module subroutine helio_drift_tp(self, system, param, dt) if (any(iflag(1:ntp) /= 0)) then tp%status = DISCARDED_DRIFTERR do i = 1, ntp - if (iflag(i) /= 0) write(*, *) "Particle ", tp%name(i), " lost due to error in Danby drift" + if (iflag(i) /= 0) write(*, *) "Particle ", tp%id(i), " lost due to error in Danby drift" end do end if end associate diff --git a/src/io/io.f90 b/src/io/io.f90 index 55789417b..57ee2176a 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -14,12 +14,12 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none ! Arguments class(swiftest_parameters), intent(inout) :: self !! Collection of parameters - integer, intent(in) :: unit !! File unit number - character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. + integer, intent(in) :: unit !! File unit number + character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader - integer, intent(out) :: iostat !! IO status code - character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 + integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader + integer, intent(out) :: iostat !! IO status code + character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals logical :: t0_set = .false. !! Is the initial time set in the input file? logical :: tstop_set = .false. !! Is the final time set in the input file? @@ -46,10 +46,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ifirst = ilast + 1 param_value = io_get_token(line_trim, ifirst, ilast, iostat) select case (param_name) - case ("NPLMAX") - read(param_value, *) self%nplmax - case ("NTPMAX") - read(param_value, *) self%ntpmax case ("T0") read(param_value, *) self%t0 t0_set = .true. @@ -110,36 +106,19 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("BIG_DISCARD") call util_toupper(param_value) if (param_value == "YES" .or. param_value == 'T' ) self%lbig_discard = .true. - case ("FRAGMENTATION") - call util_toupper(param_value) - if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. case ("MU2KG") read(param_value, *) self%MU2KG case ("TU2S") read(param_value, *) self%TU2S case ("DU2M") read(param_value, *) self%DU2M - case ("MTINY") - read(param_value, *) self%mtiny - mtiny_set = .true. case ("ENERGY") call util_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') self%lenergy = .true. - case ("ROTATION") - call util_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lrotation = .true. - case ("TIDES") - call util_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%ltides = .true. case ("GR") call util_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') self%lgr = .true. - case ("YARKOVSKY") - call util_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lyarkovsky = .true. - case ("YORP") - call util_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lyorp = .true. + case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "ROTATION", "TIDES", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default write(iomsg,*) "Unknown parameter -> ",param_name iostat = -1 @@ -212,8 +191,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) end if end if - write(*,*) "NPLMAX = ",self%nplmax - write(*,*) "NTPMAX = ",self%ntpmax write(*,*) "T0 = ",self%t0 write(*,*) "TSTOP = ",self%tstop write(*,*) "DT = ",self%dt @@ -255,24 +232,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) self%inv_c2 = einstinC * self%TU2S / self%DU2M self%inv_c2 = (self%inv_c2)**(-2) - ! The fragmentation model requires the user to set the unit system explicitly. - if ((integrator == SYMBA) .or. (integrator == RINGMOONS)) then - write(*,*) "FRAGMENTATION = ",self%lfragmentation - if (.not.mtiny_set) then - write(iomsg,*) 'SyMBA requres an MTINY value' - iostat = -1 - end if - else - if (self%lfragmentation) then - write(iomsg,*) 'This integrator does not support fragmentation. This parameter will be ignored.' - end if - if (mtiny_set) then - write(iomsg,*) 'This integrator does not support MTINY. This parameter will be ignored.' - return - end if - end if - - if ((integrator == SYMBA) .or. (integrator == RINGMOONS) .or. (integrator == RMVS)) then + if (integrator == RMVS) then if (.not.self%lclose) then write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.' iostat = -1 @@ -280,22 +240,12 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) end if end if - if (mtiny_set) then - if (self%mtiny < 0.0_DP) then - write(iomsg,*) "Invalid MTINY: ", self%mtiny - iostat = -1 - return - else - write(*,*) "MTINY = ", self%mtiny - end if - end if - ! Determine if the GR flag is set correctly for this integrator select case(integrator) case(WHM, RMVS) write(*,*) "GR = ", self%lgr case default - write(iomsg, *) 'GR is implemented compatible with this integrator. This parameter will be ignored.' + write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.' end select iostat = 0 @@ -313,83 +263,65 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) implicit none ! Arguments class(swiftest_parameters),intent(in) :: self !! Collection of parameters - integer, intent(in) :: unit !! File unit number - character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. + integer, intent(in) :: unit !! File unit number + character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. !! If you do not include a char-literal-constant, the iotype argument contains only DT. - integer, intent(in) :: v_list(:) !! Not used in this procedure - integer, intent(out) :: iostat !! IO status code - character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 + integer, intent(in) :: v_list(:) !! Not used in this procedure + integer, intent(out) :: iostat !! IO status code + character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals - !! In user-defined derived-type output, we need newline characters at the end of each format statement - !character(*),parameter :: Ifmt = '(A20,1X,I0/)' !! Format label for integer values - !character(*),parameter :: Rfmt = '(A20,1X,ES25.17/)' !! Format label for real values - !character(*),parameter :: R2fmt = '(A20,2(1X,ES25.17)/)' !! Format label for 2x real values - !character(*),parameter :: Sfmt = '(A20,1X,A/)' !! Format label for string values - !character(*),parameter :: Lfmt = '(A20,1X,L1/)' !! Format label for logical values - !character(*),parameter :: Pfmt = '(A20/)' !! Format label for single parameter string - character(*),parameter :: Ifmt = '(A20,1X,I0)' !! Format label for integer values - character(*),parameter :: Rfmt = '(A20,1X,ES25.17)' !! Format label for real values - character(*),parameter :: R2fmt = '(A20,2(1X,ES25.17))' !! Format label for 2x real values - character(*),parameter :: Sfmt = '(A20,1X,A)' !! Format label for string values - character(*),parameter :: Lfmt = '(A20,1X,L1)' !! Format label for logical values - character(*),parameter :: Pfmt = '(A20)' !! Format label for single parameter string - - write(unit, Ifmt) "NPLMAX", self%nplmax - write(unit, Ifmt) "NTPMAX", self%ntpmax - write(unit, Rfmt) "T0", self%t0 - write(unit, Rfmt) "TSTOP", self%tstop - write(unit, Rfmt) "DT", self%dt - write(unit, Sfmt) "CB_IN", trim(adjustl(self%incbfile)) - write(unit, Sfmt) "PL_IN", trim(adjustl(self%inplfile)) - write(unit, Sfmt) "TP_IN", trim(adjustl(self%intpfile)) - write(unit, Sfmt) "IN_TYPE", trim(adjustl(self%out_type)) - if (self%istep_out > 0) then - write(unit, Ifmt) "ISTEP_OUT", self%istep_out - write(unit, Sfmt) "BIN_OUT", trim(adjustl(self%outfile)) - write(unit, Sfmt) "OUT_TYPE", trim(adjustl(self%out_type)) - write(unit, Sfmt) "OUT_FORM", trim(adjustl(self%out_form)) - write(unit, Sfmt) "OUT_STAT", "APPEND" - else - write(unit, Pfmt) "!ISTEP_OUT " - write(unit, Pfmt) "!BIN_OUT" - write(unit, Pfmt) "!OUT_TYPE" - write(unit, Pfmt) "!OUT_FORM" - write(unit, Pfmt) "!OUT_STAT" - end if - write(unit, Sfmt) "ENC_OUT", trim(adjustl(self%encounter_file)) - if (self%istep_dump > 0) then - write(unit, Ifmt) "ISTEP_DUMP", self%istep_dump - else - write(unit, Pfmt) "!ISTEP_DUMP" - end if - write(unit, Rfmt) "CHK_RMIN", self%rmin - write(unit, Rfmt) "CHK_RMAX", self%rmax - write(unit, Rfmt) "CHK_EJECT", self%rmaxu - write(unit, Rfmt) "CHK_QMIN", self%qmin - if (self%qmin >= 0.0_DP) then - write(unit, Sfmt) "CHK_QMIN_COORD", trim(adjustl(self%qmin_coord)) - write(unit, R2fmt) "CHK_QMIN_RANGE", self%qmin_alo, self%qmin_ahi - else - write(unit, Pfmt) "!CHK_QMIN_COORD" - write(unit, Pfmt) "!CHK_QMIN_RANGE" - end if - if (self%lmtiny) write(unit, Rfmt) "MTINY", self%mtiny - write(unit, Rfmt) "MU2KG", self%MU2KG - write(unit, Rfmt) "TU2S", self%TU2S - write(unit, Rfmt) "DU2M", self%DU2M - - write(unit, Lfmt) "EXTRA_FORCE", self%lextra_force - write(unit, Lfmt) "BIG_DISCARD", self%lbig_discard - write(unit, Lfmt) "CHK_CLOSE", self%lclose - write(unit, Lfmt) "FRAGMENTATION", self%lfragmentation - write(unit, Lfmt) "ROTATION", self%lrotation - write(unit, Lfmt) "TIDES", self%ltides - write(unit, Lfmt) "GR", self%lgr - write(unit, Lfmt) "ENERGY", self%lenergy - !write(unit, Lfmt) "YARKOVSKY", self%lyarkovsky - !write(unit, Lfmt) "YORP", self%lyorp - iostat = 0 - iomsg = "UDIO not implemented" + character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values + character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values + character(*),parameter :: Rarrfmt = '(3(ES25.17,1X))' !! Format label for real values + character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values + character(len=*), parameter :: Afmt = '(A25,1X,64(:,A25,1X))' + character(256) :: param_name, param_value + type character_array + character(25) :: value + end type character_array + type(character_array), dimension(:), allocatable :: param_array + integer(I4B) :: i + + associate(param => self) + write(param_name, Afmt) "T0"; write(param_value,Rfmt) param%t0; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "TSTOP"; write(param_value, Rfmt) param%tstop; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "DT"; write(param_value, Rfmt) param%dt; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "PL_IN"; write(param_value, Afmt) trim(adjustl(param%inplfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "TP_in"; write(param_value, Afmt) trim(adjustl(param%intpfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "IN_TYPE"; write(param_value, Afmt) trim(adjustl(param%in_type)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + if (param%istep_out > 0) then + write(param_name, Afmt) "ISTEP_OUT"; write(param_value, Ifmt) param%istep_out; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "BIN_OUT"; write(param_value, Afmt) trim(adjustl(param%outfile)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "OUT_TYPE"; write(param_value, Afmt) trim(adjustl(param%out_type)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + 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) + 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 + write(param_name, Afmt) "CHK_RMIN"; write(param_value, Rfmt) param%rmin; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "CHK_RMAX"; write(param_value, Rfmt) param%rmax; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "CHK_EJECT"; write(param_value, Rfmt) param%rmaxu; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "CHK_QMIN"; write(param_value, Rfmt) param%qmin; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + if (param%qmin >= 0.0_DP) then + write(param_name, Afmt) "CHK_QMIN_COORD"; write(param_value, Afmt) trim(adjustl(param%qmin_coord)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + allocate(param_array(2)) + write(param_array(1)%value, Rfmt) param%qmin_alo + write(param_array(2)%value, Rfmt) param%qmin_ahi + write(param_name, Afmt) "CHK_QMIN_RANGE"; write(unit, Afmt) adjustl(param_name), adjustl(param_array(1)%value), adjustl(param_array(2)%value) + end if + write(param_name, Afmt) "MU2KG"; write(param_value, Rfmt) param%MU2KG; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "TU2S"; write(param_value, Rfmt) param%TU2S ; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "DU2M"; write(param_value, Rfmt) param%DU2M; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "EXTRA_FORCE"; write(param_value, Lfmt) param%lextra_force; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "BIG_DISCARD"; write(param_value, Lfmt) param%lbig_discard; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "CHK_CLOSE"; write(param_value, Lfmt) param%lclose; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "ENERGY"; write(param_value, Lfmt) param%lenergy; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "GR"; write(param_value, Lfmt) param%lgr; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + iostat = 0 + iomsg = "UDIO not implemented" + end associate return end subroutine io_param_writer @@ -570,7 +502,7 @@ module function io_get_args(integrator, param_file_name) result(ierr) if (ierr /= 0) call util_exit(USAGE) end function io_get_args - function io_get_token(buffer, ifirst, ilast, ierr) result(token) + module function io_get_token(buffer, ifirst, ilast, ierr) result(token) !! author: David A. Minton !! !! Retrieves a character token from an input string. Here a token is defined as any set of contiguous non-blank characters not @@ -658,7 +590,7 @@ module subroutine io_read_body_in(self, param) do i = 1, nbody select type(self) class is (swiftest_pl) - read(iu, *, iostat = ierr) self%name(i), val + read(iu, *, iostat = ierr) self%id(i), val self%mass(i) = real(val / param%GU, kind=DP) self%Gmass(i) = real(val, kind=DP) if (param%lclose) then @@ -667,16 +599,8 @@ module subroutine io_read_body_in(self, param) else self%radius(i) = 0.0_DP end if - if (param%lrotation) then - read(iu, iostat = ierr) self%Ip(:, i) - read(iu, iostat = ierr) self%rot(:, i) - end if - if (param%ltides) then - read(iu, iostat = ierr) self%k2(i) - read(iu, iostat = ierr) self%Q(i) - end if class is (swiftest_tp) - read(iu, *, iostat = ierr) self%name(i) + read(iu, *, iostat = ierr) self%id(i) end select if (ierr /= 0 ) exit read(iu, *, iostat = ierr) self%xh(1, i), self%xh(2, i), self%xh(3, i) @@ -735,15 +659,6 @@ module subroutine io_read_cb_in(self, param) read(iu, *, iostat = ierr) self%radius read(iu, *, iostat = ierr) self%j2rp2 read(iu, *, iostat = ierr) self%j4rp4 - if (param%lrotation) then - read(iu, *, iostat = ierr) self%Ip(:) - read(iu, *, iostat = ierr) self%rot(:) - end if - if (param%ltides) then - read(iu, *, iostat = ierr) self%k2 - read(iu, *, iostat = ierr) self%Q - end if - else open(unit = iu, file = param%incbfile, status = 'old', form = 'UNFORMATTED', iostat = ierr) call self%read_frame(iu, param, XV, ierr) @@ -852,20 +767,20 @@ end function io_read_encounter module subroutine io_read_frame_body(self, iu, param, form, ierr) !! author: David A. Minton !! - !! Reads a frame of output of either test particle or massive body data to the binary output file - !! Note: If outputting to orbital elements, but sure that the conversion is done prior to calling this method + !! Reads a frame of output of either test particle or massive body data from a binary output file !! !! Adapted from David E. Kaufmann's Swifter routine io_read_frame.f90 !! Adapted from Hal Levison's Swift routine io_read_frame.F implicit none ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest particle object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_body), intent(inout) :: self !! Swiftest particle object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - integer(I4B), intent(out) :: ierr !! Error code + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code associate(n => self%nbody) + read(iu, iostat = ierr) self%id(1:n) read(iu, iostat = ierr) self%name(1:n) select case (form) case (EL) @@ -888,18 +803,6 @@ module subroutine io_read_frame_body(self, iu, param, form, ierr) read(iu, iostat = ierr) self%Gmass(1:n) self%mass(1:n) = self%Gmass / param%GU read(iu, iostat = ierr) self%radius(1:n) - if (param%lrotation) then - read(iu, iostat = ierr) self%Ip(1, 1:n) - read(iu, iostat = ierr) self%Ip(2, 1:n) - read(iu, iostat = ierr) self%Ip(3, 1:n) - read(iu, iostat = ierr) self%rot(1, 1:n) - read(iu, iostat = ierr) self%rot(2, 1:n) - read(iu, iostat = ierr) self%rot(3, 1:n) - end if - if (param%ltides) then - read(iu, iostat = ierr) self%k2(1:n) - read(iu, iostat = ierr) self%Q(1:n) - end if end select end associate @@ -926,19 +829,13 @@ module subroutine io_read_frame_cb(self, iu, param, form, ierr) character(*), intent(in) :: form !! Input format code ("XV" or "EL") integer(I4B), intent(out) :: ierr !! Error cod + read(iu, iostat = ierr) self%id + read(iu, iostat = ierr) self%name read(iu, iostat = ierr) self%Gmass self%mass = self%Gmass / param%GU read(iu, iostat = ierr) self%radius read(iu, iostat = ierr) self%j2rp2 read(iu, iostat = ierr) self%j4rp4 - if (param%lrotation) then - read(iu, iostat = ierr) self%Ip(:) - read(iu, iostat = ierr) self%rot(:) - end if - if (param%ltides) then - read(iu, iostat = ierr) self%k2 - read(iu, iostat = ierr) self%Q - end if if (ierr /=0) then write(*,*) 'Error reading central body data' call util_exit(FAILURE) @@ -1065,7 +962,7 @@ module subroutine io_write_discard(self, param) 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%name, dstatus => self%tp_discards%status) + dname => self%tp_discards%id, dstatus => self%tp_discards%status) select case(param%out_stat) case('APPEND') @@ -1091,7 +988,7 @@ module subroutine io_write_discard(self, param) 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%name, xh => self%pl%xh) + Rpl => self%pl%radius, name => self%pl%id, xh => self%pl%xh) if (param%lgr) then allocate(pltemp, source = pl) @@ -1187,6 +1084,7 @@ module subroutine io_write_frame_body(self, iu, param) associate(n => self%nbody) if (n == 0) return + write(iu) self%id(1:n) write(iu) self%name(1:n) select case (param%out_form) case (EL) @@ -1208,18 +1106,6 @@ module subroutine io_write_frame_body(self, iu, param) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body write(iu) self%Gmass(1:n) write(iu) self%radius(1:n) - if (param%lrotation) then - write(iu) self%Ip(1, 1:n) - write(iu) self%Ip(2, 1:n) - write(iu) self%Ip(3, 1:n) - write(iu) self%rot(1, 1:n) - write(iu) self%rot(2, 1:n) - write(iu) self%rot(3, 1:n) - end if - if (param%ltides) then - write(iu) self%k2(1:n) - write(iu) self%Q(1:n) - end if end select end associate @@ -1239,23 +1125,12 @@ module subroutine io_write_frame_cb(self, iu, param) integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + write(iu) self%id + write(iu) self%name write(iu) self%Gmass write(iu) self%radius write(iu) self%j2rp2 write(iu) self%j4rp4 - if (param%lrotation) then - write(iu) self%Ip(1) - write(iu) self%Ip(2) - write(iu) self%Ip(3) - write(iu) self%rot(1) - write(iu) self%rot(2) - write(iu) self%rot(3) - end if - if (param%ltides) then - write(iu) self%k2 - write(iu) self%Q - end if - return end subroutine io_write_frame_cb diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index b7671883a..b7bdf826c 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -4,7 +4,7 @@ module helio_classes !! Definition of classes and methods specific to the Democratic Heliocentric Method !! Adapted from David E. Kaufmann's Swifter routine: helio.f90 use swiftest_globals - use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp + use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp, swiftest_nbody_system use whm_classes, only : whm_nbody_system implicit none @@ -171,7 +171,7 @@ module subroutine helio_step_pl(self, system, param, t, dt) end subroutine helio_step_pl module subroutine helio_step_tp(self, system, param, t, dt) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + use swiftest_classes, only : swiftest_cb, swiftest_parameters, swiftest_nbody_system implicit none class(helio_tp), intent(inout) :: self !! Helio test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 5e1d69437..7b4ecc36d 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -81,7 +81,7 @@ module rmvs_classes !******************************************************************************************************************************* !> RMVS massive body particle class - type, private, extends(whm_pl) :: rmvs_pl + type, public, extends(whm_pl) :: rmvs_pl integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet integer(I4B), dimension(:), allocatable :: plind ! Connects the planetocentric indices back to the heliocentric planet list diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 3d1cba23d..d8a5e46d5 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -9,7 +9,7 @@ module swiftest_classes public :: discard_pl, discard_system, discard_tp public :: drift_one public :: eucl_dist_index_plpl, eucl_dist_index_pltp, eucl_irij3_plpl - public :: io_dump_param, io_dump_swiftest, io_dump_system, io_get_args, io_param_reader, io_param_writer, io_read_body_in, & + public :: io_dump_param, io_dump_swiftest, io_dump_system, io_get_args, io_get_token, io_param_reader, io_param_writer, io_read_body_in, & io_read_cb_in, io_read_param_in, io_read_frame_body, io_read_frame_cb, io_read_frame_system, io_read_initialize_system, & io_write_discard, io_write_encounter, io_write_frame_body, io_write_frame_cb, io_write_frame_system public :: kickvh_body @@ -17,10 +17,9 @@ module swiftest_classes public :: orbel_el2xv_vec, orbel_xv2el_vec, orbel_scget, orbel_xv2aeq, orbel_xv2aqt public :: setup_body, setup_construct_system, setup_pl, setup_tp public :: user_getacch_body - public :: util_coord_b2h_pl, util_coord_b2h_tp, util_coord_h2b_pl, util_coord_h2b_tp, util_copy_body, util_copy_cb, util_copy_pl, & - util_copy_tp, util_copy_system, util_fill_body, util_fill_pl, util_fill_tp, util_reverse_status, util_set_beg_end_cb, & - util_set_beg_end_pl, util_set_ir3h, util_set_msys, util_set_mu_pl, util_set_mu_tp, util_set_rhill, & - util_spill_body, util_spill_pl, util_spill_tp + public :: util_coord_b2h_pl, util_coord_b2h_tp, util_coord_h2b_pl, util_coord_h2b_tp, util_fill_body, util_fill_pl, util_fill_tp, & + util_reverse_status, util_set_beg_end_cb, util_set_beg_end_pl, util_set_ir3h, util_set_msys, util_set_mu_pl, & + util_set_mu_tp, util_set_rhill, util_spill_body, util_spill_pl, util_spill_tp !******************************************************************************************************************************** ! swiftest_parameters class definitions @@ -54,7 +53,6 @@ module swiftest_classes 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 - real(DP) :: MTINY = -1.0_DP !! Smallest mass that is fully gravitating 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 @@ -65,10 +63,6 @@ module swiftest_classes logical :: lextra_force = .false. !! User defined force function turned on logical :: lbig_discard = .false. !! Save big bodies on every discard logical :: lclose = .false. !! Turn on close encounters - logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. - logical :: lmtiny = .false. !! Use the MTINY variable (Automatically set if running SyMBA) - logical :: lrotation = .false. !! Include rotation states of big bodies - logical :: ltides = .false. !! Include tidal dissipation logical :: lenergy = .false. !! Track the total energy of the system logical :: loblatecb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input) @@ -96,11 +90,10 @@ module swiftest_classes contains !! The minimal methods that all systems must have private - procedure :: dump => io_dump_swiftest + procedure :: dump => io_dump_swiftest procedure(abstract_initialize), public, deferred :: initialize - procedure(abstract_write_frame), public, deferred :: write_frame procedure(abstract_read_frame), public, deferred :: read_frame - procedure(abstract_copy), public, deferred :: copy + procedure(abstract_write_frame), public, deferred :: write_frame end type swiftest_base !******************************************************************************************************************************** @@ -108,6 +101,8 @@ module swiftest_classes !******************************************************************************************************************************** !> A concrete lass for the central body in a Swiftest simulation type, abstract, public, extends(swiftest_base) :: swiftest_cb + character(len=STRMAX) :: name !! Non-unique name + integer(I4B) :: id !! External identifier (unique) real(DP) :: mass = 0.0_DP !! Central body mass (units MU) real(DP) :: Gmass = 0.0_DP !! Central mass gravitational term G * mass (units GU * MU) real(DP) :: radius = 0.0_DP !! Central body radius (units DU) @@ -119,16 +114,11 @@ module swiftest_classes real(DP), dimension(NDIM) :: aoblend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step real(DP), dimension(NDIM) :: xb = 0.0_DP !! Barycentric position (units DU) real(DP), dimension(NDIM) :: vb = 0.0_DP !! Barycentric velocity (units DU / TU) - real(DP), dimension(NDIM) :: Ip = 0.0_DP !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. - real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) - real(DP) :: k2 = 0.0_DP !! Tidal Love number - real(DP) :: Q = 0.0_DP !! Tidal quality factor contains private procedure, public :: initialize => io_read_cb_in !! I/O routine for reading in central body data procedure, public :: write_frame => io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body procedure, public :: read_frame => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body - procedure, public :: copy => util_copy_cb !! Copies elements of one object to another. procedure, public :: set_beg_end => util_set_beg_end_cb !! Sets the beginning and ending oblateness acceleration term end type swiftest_cb @@ -138,25 +128,28 @@ module swiftest_classes !> An abstract class for a generic collection of Swiftest bodies type, abstract, public, extends(swiftest_base) :: swiftest_body !! Superclass that defines the generic elements of a Swiftest particle - logical :: lfirst = .true. !! Run the current step as a first - integer(I4B) :: nbody = 0 !! Number of bodies - integer(I4B), dimension(:), allocatable :: name !! External identifier - integer(I4B), dimension(:), allocatable :: status !! An integrator-specific status indicator - logical, dimension(:), allocatable :: ldiscard !! Body should be discarded - real(DP), dimension(:,:), allocatable :: xh !! Heliocentric position - real(DP), dimension(:,:), allocatable :: vh !! Heliocentric velocity - real(DP), dimension(:,:), allocatable :: xb !! Barycentric position - real(DP), dimension(:,:), allocatable :: vb !! Barycentric velocity - real(DP), dimension(:,:), allocatable :: ah !! Total heliocentric acceleration - real(DP), dimension(:,:), allocatable :: aobl !! Barycentric accelerations of bodies due to central body oblatenes - real(DP), dimension(:), allocatable :: ir3h !! Inverse heliocentric radius term (1/rh**3) - real(DP), dimension(:), allocatable :: a !! Semimajor axis (pericentric distance for a parabolic orbit) - real(DP), dimension(:), allocatable :: e !! Eccentricity - real(DP), dimension(:), allocatable :: inc !! Inclination - real(DP), dimension(:), allocatable :: capom !! Longitude of ascending node - real(DP), dimension(:), allocatable :: omega !! Argument of pericenter - real(DP), dimension(:), allocatable :: capm !! Mean anomaly - real(DP), dimension(:), allocatable :: mu !! G * (Mcb + [m]) + logical :: lfirst = .true. !! Run the current step as a first + integer(I4B) :: nbody = 0 !! Number of bodies + character(len=STRMAX), dimension(:), allocatable :: name !! Non-unique name + integer(I4B), dimension(:), allocatable :: id !! External identifier (unique) + integer(I4B), dimension(:), allocatable :: status !! An integrator-specific status indicator + logical, dimension(:), allocatable :: ldiscard !! Body should be discarded + real(DP), dimension(:,:), allocatable :: xh !! Heliocentric position + real(DP), dimension(:,:), allocatable :: vh !! Heliocentric velocity + real(DP), dimension(:,:), allocatable :: xb !! Barycentric position + real(DP), dimension(:,:), allocatable :: vb !! Barycentric velocity + real(DP), dimension(:,:), allocatable :: ah !! Total heliocentric acceleration + real(DP), dimension(:,:), allocatable :: aobl !! Barycentric accelerations of bodies due to central body oblatenes + real(DP), dimension(:), allocatable :: ir3h !! Inverse heliocentric radius term (1/rh**3) + real(DP), dimension(:), allocatable :: a !! Semimajor axis (pericentric distance for a parabolic orbit) + real(DP), dimension(:), allocatable :: e !! Eccentricity + real(DP), dimension(:), allocatable :: inc !! Inclination + real(DP), dimension(:), allocatable :: capom !! Longitude of ascending node + real(DP), dimension(:), allocatable :: omega !! Argument of pericenter + real(DP), dimension(:), allocatable :: capm !! Mean anomaly + real(DP), dimension(:), allocatable :: mu !! G * (Mcb + [m]) + integer(I4B), dimension(:,:), allocatable :: k_eucl !! Index array used to convert flattened the body-body comparison upper triangular matrix + integer(I8B) :: num_comparisons !! Number of body-body comparisons in the flattened upper triangular matrix !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the !! component list, such as setup_body and util_spill contains @@ -176,7 +169,6 @@ module swiftest_classes procedure, public :: set_ir3 => util_set_ir3h !! Sets the inverse heliocentric radius term (1/rh**3) procedure, public :: setup => setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays procedure, public :: accel_user => user_getacch_body !! Add user-supplied heliocentric accelerations to planets - procedure, public :: copy => util_copy_body !! Copies elements of one object to another. procedure, public :: fill => util_fill_body !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) procedure, public :: spill => util_spill_body !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) procedure, public :: reverse_status => util_reverse_status !! Reverses the active/inactive status of all particles in a structure @@ -192,19 +184,11 @@ module swiftest_classes real(DP), dimension(:), allocatable :: Gmass !! Mass gravitational term G * mass (units GU * MU) real(DP), dimension(:), allocatable :: rhill !! Body mass (units MU) real(DP), dimension(:), allocatable :: radius !! Body radius (units DU) - real(DP), dimension(:), allocatable :: density !! Body mass density - calculated internally (units MU / DU**3) - real(DP), dimension(:,:), allocatable :: Ip !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). - !! Principal axis rotation assumed. - real(DP), dimension(:,:), allocatable :: rot !! Body rotation vector in inertial coordinate frame (units rad / TU) - real(DP), dimension(:), allocatable :: k2 !! Tidal Love number - real(DP), dimension(:), allocatable :: Q !! Tidal quality factor - integer(I4B) :: num_comparisons !! Number of pl-pl Euclidean distance comparisons - integer(I4B), dimension(:,:), allocatable :: k_eucl !! Index array that converts i, j array indices into k index for use in - !! the Euclidean distance matrix real(DP), dimension(:), allocatable :: irij3 !! 1.0_DP / (rji2 * sqrt(rji2)) where rji2 is the square of the Euclidean distance real(DP), dimension(:,:), allocatable :: xbeg !! Position at beginning of step real(DP), dimension(:,:), allocatable :: xend !! Position at end of step real(DP), dimension(:,:), allocatable :: vbeg !! Velocity at beginning of step + real(DP), dimension(:), allocatable :: density !! Body mass density - calculated internally (units MU / DU**3) !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the !! component list, such as setup_pl and util_spill_pl contains @@ -220,7 +204,6 @@ module swiftest_classes procedure, public :: set_rhill => util_set_rhill !! Calculates the Hill's radii for each body procedure, public :: h2b => util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) procedure, public :: b2h => util_coord_b2h_pl !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity) - procedure, public :: copy => util_copy_pl !! Copies elements of one object to another. procedure, public :: fill => util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) procedure, public :: set_beg_end => util_set_beg_end_pl !! Sets the beginning and ending positions and velocities of planets. procedure, public :: spill => util_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) @@ -249,7 +232,6 @@ module swiftest_classes procedure, public :: set_mu => util_set_mu_tp !! Method used to construct the vectorized form of the central body mass procedure, public :: h2b => util_coord_h2b_tp !! Convert test particles from heliocentric to barycentric coordinates (position and velocity) procedure, public :: b2h => util_coord_b2h_tp !! Convert test particles from barycentric to heliocentric coordinates (position and velocity) - procedure, public :: copy => util_copy_tp !! Copies elements of one object to another. procedure, public :: fill => util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) procedure, public :: spill => util_spill_tp !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type swiftest_tp @@ -260,18 +242,22 @@ module swiftest_classes !> An abstract class for a basic Swiftest nbody system type, abstract, public, extends(swiftest_base) :: swiftest_nbody_system !! This superclass contains a minimial system of a set of test particles (tp), massive bodies (pl), and a central body (cb) - class(swiftest_cb), allocatable :: cb !! Central body data structure - class(swiftest_pl), allocatable :: pl !! Massive body data structure - class(swiftest_tp), allocatable :: tp !! Test particle data structure - class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure - real(DP) :: msys = 0.0_DP !! Total system mass - used for barycentric coordinate conversion - real(DP) :: ke = 0.0_DP !! System kinetic energy - real(DP) :: pe = 0.0_DP !! System potential energy - real(DP) :: te = 0.0_DP !! System total energy - real(DP), dimension(NDIM) :: htot = 0.0_DP !! System angular momentum vector - logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated - !! separately from massive bodies. Massive body variables are saved at half steps, and passed to - !! the test particles + class(swiftest_cb), allocatable :: cb !! Central body data structure + class(swiftest_pl), allocatable :: pl !! Massive body data structure + class(swiftest_tp), allocatable :: tp !! Test particle data structure + class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure + real(DP) :: msys = 0.0_DP !! Total system mass - used for barycentric coordinate conversion + real(DP) :: ke = 0.0_DP !! System kinetic energy + real(DP) :: pe = 0.0_DP !! System potential energy + real(DP) :: te = 0.0_DP !! System total energy + real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector + real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) + real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) + real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions + real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies + logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated + !! separately from massive bodies. Massive body variables are saved at half steps, and passed to + !! the test particles contains private !> Each integrator will have its own version of the step @@ -285,17 +271,9 @@ module swiftest_classes procedure, public :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. procedure, public :: write_discard => io_write_discard !! Append a frame of output data to file procedure, public :: write_frame => io_write_frame_system !! Append a frame of output data to file - procedure, public :: copy => util_copy_system !! Copies elements of one object to another. end type swiftest_nbody_system abstract interface - subroutine abstract_copy(self, src, mask) - import swiftest_base - class(swiftest_base), intent(inout) :: self - class(swiftest_base), intent(in) :: src - logical, dimension(:), intent(in) :: mask - end subroutine abstract_copy - subroutine abstract_discard_body(self, system, param) import swiftest_body, swiftest_nbody_system, swiftest_parameters class(swiftest_body), intent(inout) :: self !! Swiftest body object @@ -432,6 +410,15 @@ module function io_get_args(integrator, param_file_name) result(ierr) integer(I4B) :: ierr !! I/O error code end function io_get_args + module function io_get_token(buffer, ifirst, ilast, ierr) result(token) + implicit none + character(len=*), intent(in) :: buffer !! Input string buffer + integer(I4B), intent(inout) :: ifirst !! Index of the buffer at which to start the search for a token + integer(I4B), intent(out) :: ilast !! Index of the buffer at the end of the returned token + integer(I4B), intent(out) :: ierr !! Error code + character(len=:), allocatable :: token !! Returned token string + end function io_get_token + module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none class(swiftest_parameters), intent(inout) :: self !! Collection of parameters @@ -682,41 +669,6 @@ module subroutine util_coord_h2b_tp(self, cb) class(swiftest_cb), intent(in) :: cb !! Swiftest central body object end subroutine util_coord_h2b_tp - module subroutine util_copy_body(self, src, mask) - implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest body object to copy into - class(swiftest_base), intent(in) :: src !! Swiftest base object to copy from - logical, dimension(:), intent(in) :: mask !! Mask of elements in src object to copy into self - end subroutine util_copy_body - - module subroutine util_copy_cb(self, src, mask) - implicit none - class(swiftest_cb), intent(inout) :: self !! Swiftest central body object to copy into - class(swiftest_base), intent(in) :: src !! Swiftest base object to copy from - logical, dimension(:), intent(in) :: mask !! Mask of elements in src object to copy into selfk - end subroutine util_copy_cb - - module subroutine util_copy_pl(self, src, mask) - implicit none - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object to copy into - class(swiftest_base), intent(in) :: src !! Swiftest base object to copy from - logical, dimension(:), intent(in) :: mask !! Mask of elements in src object to copy into self - end subroutine util_copy_pl - - module subroutine util_copy_tp(self, src, mask) - implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object to copy into - class(swiftest_base), intent(in) :: src !! Swiftest base object to copy from - logical, dimension(:), intent(in) :: mask !! Mask of elements in src object to copy into self - end subroutine util_copy_tp - - module subroutine util_copy_system(self, src, mask) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object to copy into - class(swiftest_base), intent(in) :: src !! Swiftest base object to copy from - logical, dimension(:), intent(in) :: mask !! Mask of elements in src object to copy into self - end subroutine util_copy_system - module subroutine util_fill_body(self, inserts, lfill_list) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index f738dbad0..91db0adf3 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -109,9 +109,9 @@ module swiftest_globals character(*), parameter :: ENC_OUTFILE = 'encounter.out' character(*), parameter :: DISCARD_FILE = 'discard.out' character(*), parameter :: ENERGY_FILE = 'energy.out' - character(*), parameter :: CB_INFILE = 'cb_out.dat' - character(*), parameter :: PL_INFILE = 'pl_out.dat' - character(*), parameter :: TP_INFILE = 'tp_out.dat' + character(*), parameter :: CB_INFILE = 'cb.in' + character(*), parameter :: PL_INFILE = 'pl.in' + character(*), parameter :: TP_INFILE = 'tp.in' character(*), parameter :: BIN_OUTFILE = 'bin.dat' integer(I4B), parameter :: BINUNIT = 20 !! File unit number for the binary output file diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 14cf7c7b2..05c2838b6 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -1,61 +1,193 @@ module symba_classes !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! - !! Definition of classes and methods specific to the Democratic Heliocentric Method + !! Definition of classes and methods specific to the Democratic SyMBAcentric Method !! Adapted from David E. Kaufmann's Swifter routine: helio.f90 use swiftest_globals - use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system + use swiftest_classes, only : swiftest_parameters, swiftest_base + use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system implicit none + integer(I4B), parameter :: NENMAX = 32767 + integer(I4B), parameter :: NTENC = 3 + real(DP), parameter :: RHSCALE = 6.5_DP + real(DP), parameter :: RSHELL = 0.48075_DP + character(*), parameter :: PARTICLE_OUTFILE = 'particle.dat' + integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file - !******************************************************************************************************************************** - ! symba_nbody_system class definitions and method interfaces - !******************************************************************************************************************************** - type, public, extends(helio_nbody_system) :: symba_nbody_system + type, public, extends(swiftest_parameters) :: symba_parameters + character(STRMAX) :: particle_file = PARTICLE_OUTFILE !! Name of output particle information file + real(DP) :: MTINY = -1.0_DP !! Smallest mass that is fully gravitating + integer(I4B), dimension(:), allocatable :: seed !! Random seeds + logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. + logical :: lrotation = .false. !! Include rotation states of big bodies + logical :: ltides = .false. !! Include tidal dissipation contains private - procedure, public :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step - procedure, public :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system - end type symba_nbody_system + procedure, public :: reader => symba_io_param_reader + procedure, public :: writer => symba_io_param_writer + end type symba_parameters !******************************************************************************************************************************** ! symba_cb class definitions and method interfaces !******************************************************************************************************************************* - !> Helio central body particle class + !> SyMBA central body particle class type, public, extends(helio_cb) :: symba_cb + real(DP), dimension(NDIM) :: Ip = 0.0_DP !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. + real(DP), dimension(NDIM) :: rot = 0.0_DP !! Body rotation vector in inertial coordinate frame (units rad / TU) + real(DP) :: k2 = 0.0_DP !! Tidal Love number + real(DP) :: Q = 0.0_DP !! Tidal quality factor + real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body + real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body + real(DP) :: M0 = 0.0_DP !! Initial mass of the central body + real(DP) :: dM = 0.0_DP !! Change in mass of the central body + real(DP) :: R0 = 0.0_DP !! Initial radius of the central body + real(DP) :: dR = 0.0_DP !! Change in the radius of the central body contains private + procedure, public :: initialize => symba_io_read_cb_in !! I/O routine for reading in particle info data + procedure, public :: read_frame => symba_io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central bod + procedure, public :: write_frame => symba_io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body end type symba_cb !******************************************************************************************************************************** - ! symba_pl class definitions and method interfaces + ! symba_particle_info class definitions and method interfaces !******************************************************************************************************************************* + !> Class definition for the particle origin information object. This object is used to track time, location, and collisional regime + !> of fragments produced in collisional events. + type, public, extends(swiftest_base) :: symba_particle_info + character(len=32) :: origin_type !! String containing a description of the origin of the particle (e.g. Initial Conditions, Supercatastrophic, Disruption, etc.) + real(DP) :: origin_time !! The time of the particle's formation + real(DP), dimension(NDIM) :: origin_xh !! The heliocentric distance vector at the time of the particle's formation + real(DP), dimension(NDIM) :: origin_vh !! The heliocentric velocity vector at the time of the particle's formation + contains + private + procedure, public :: dump => symba_io_dump_particle_info !! I/O routine for dumping particle info to file + procedure, public :: initialize => symba_io_initialize_particle_info !! I/O routine for reading in particle info data + procedure, public :: read_frame => symba_io_read_frame_info !! I/O routine for reading in a single frame of particle info + procedure, public :: write_frame => symba_io_write_frame_info !! I/O routine for writing out a single frame of particle info + end type symba_particle_info + + !******************************************************************************************************************************** + ! symba_kinship class definitions and method interfaces + !******************************************************************************************************************************* + !> Class definition for the kinship relationships used in bookkeeping multiple collisions bodies in a single time step. + type symba_kinship + integer(I4B) :: parent ! Index of parent particle + integer(I4B) :: nchild ! number of children in merger list + integer(I4B), dimension(:), allocatable :: child ! Index of children particles + end type symba_kinship - !! Helio massive body particle class + !******************************************************************************************************************************** + ! symba_pl class definitions and method interfaces + !******************************************************************************************************************************* + !> SyMBA massive body class type, public, extends(helio_pl) :: symba_pl + real(DP), dimension(:,:), allocatable :: Ip !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). + !! Principal axis rotation assumed. + real(DP), dimension(:,:), allocatable :: rot !! Body rotation vector in inertial coordinate frame (units rad / TU) + real(DP), dimension(:), allocatable :: k2 !! Tidal Love number + real(DP), dimension(:), allocatable :: Q !! Tidal quality factor + logical, dimension(:), allocatable :: lcollision !! flag indicating whether body has merged with another this time step + logical, dimension(:), allocatable :: lencounter !! flag indicating whether body is part of an encounter this time step + integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with other planets this time step + integer(I4B), dimension(:), allocatable :: ntpenc !! number of encounters with test particles this time step + integer(I4B), dimension(:), allocatable :: levelg !! level at which this body should be moved + integer(I4B), dimension(:), allocatable :: levelm !! deepest encounter level achieved this time step + integer(I4B), dimension(:), allocatable :: isperi !! perihelion passage flag + real(DP), dimension(:), allocatable :: peri !! perihelion distance + real(DP), dimension(:), allocatable :: atp !! semimajor axis following perihelion passage + type(symba_kinship), dimension(:), allocatable :: kin !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step + type(symba_particle_info), dimension(:), allocatable :: info contains private procedure, public :: discard => symba_discard_pl !! Process massive body discards procedure, public :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other + procedure, public :: read_frame => symba_io_read_frame_pl !! I/O routine for reading out a single frame of time-series data for a massive body + procedure, public :: initialize => symba_io_read_pl_in !! I/O routine for reading in a massive body structure from file with SyMBA-specific parameters + procedure, public :: write_frame => symba_io_write_frame_pl !! I/O routine for writing out a single frame of time-series data for a massive body end type symba_pl !******************************************************************************************************************************** ! symba_tp class definitions and method interfaces !******************************************************************************************************************************* - - !! Helio test particle class + !> SyMBA test particle class type, public, extends(helio_tp) :: symba_tp + integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with planets this time step + integer(I4B), dimension(:), allocatable :: levelg !! level at which this particle should be moved + integer(I4B), dimension(:), allocatable :: levelm !! deepest encounter level achieved this time step contains private procedure, public :: discard => symba_discard_tp !! process test particle discards procedure, public :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body end type symba_tp + !******************************************************************************************************************************** + ! symba_plplenc class definitions and method interfaces + !******************************************************************************************************************************* + !> SyMBA class for tracking pl-pl close encounters in a step + type symba_plplenc + integer(I4B) :: nplplenc !! Total number of pl-pl encounters + logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag + integer(I4B), dimension(:), allocatable :: status !! status of the interaction + integer(I4B), dimension(:), allocatable :: level !! encounter recursion level + integer(I4B), dimension(:), allocatable :: index1 !! position of the first planet in encounter + integer(I4B), dimension(:), allocatable :: index2 !! position of the second planet in encounter + integer(I4B), dimension(:), allocatable :: enc_child !! the child of the encounter + integer(I4B), dimension(:), allocatable :: enc_parent !! the child of the encounter + real(DP), dimension(:,:), allocatable :: xh1 !! the heliocentric position of parent 1 in encounter + real(DP), dimension(:,:), allocatable :: xh2 !! the heliocentric position of parent 2 in encounter + real(DP), dimension(:,:), allocatable :: vb1 !! the barycentric velocity of parent 1 in encounter + real(DP), dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter + end type symba_plplenc + + !******************************************************************************************************************************** + ! symba_pltpenc class definitions and method interfaces + !******************************************************************************************************************************* + !> SyMBA class for tracking pl-tp close encounters in a step + type symba_pltpenc + integer(I4B) :: npltpenc !! Total number of pl-tp encounters + logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag + integer(I4B), dimension(:), allocatable :: status !! status of the interaction + integer(I4B), dimension(:), allocatable :: level !! encounter recursion level + integer(I4B), dimension(:), allocatable :: indexpl !! position of the planet in encounter + integer(I4B), dimension(:), allocatable :: indextp !! position of the test particle in encounter + end type symba_pltpenc + + !******************************************************************************************************************************** + ! symba_merger class definitions and method interfaces + !******************************************************************************************************************************* + !> SyMBA class for tracking pl-pl mergers in a step + type symba_merger + integer(I4B), dimension(:), allocatable :: id !! identifier + integer(I4B), dimension(:), allocatable :: index_ps !! position of the particle + integer(I4B), dimension(:), allocatable :: status !! status + integer(I4B), dimension(:), allocatable :: nadded !! number of resultant bodies from this collisional event aka number of fragments + real(DP), dimension(:,:), allocatable :: xb !! barycentric position + real(DP), dimension(:,:), allocatable :: vb !! barycentric velocity + real(DP), dimension(:), allocatable :: mass !! mass + real(DP), dimension(:), allocatable :: radius ! radius + real(DP), dimension(:,:), allocatable :: IP ! moment of intertia + real(DP), dimension(:,:), allocatable :: rot ! rotation + type(symba_particle_info), dimension(:), allocatable :: info + end type symba_merger + + !******************************************************************************************************************************** + ! symba_nbody_system class definitions and method interfaces + !******************************************************************************************************************************** + type, public, extends(helio_nbody_system) :: symba_nbody_system + contains + private + procedure, public :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step + procedure, public :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system + end type symba_nbody_system + + interface module subroutine symba_discard_pl(self, system, param) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(symba_pl), intent(inout) :: self !! RMVS test particle object + class(symba_pl), intent(inout) :: self !! SyMBA test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_discard_pl @@ -63,31 +195,137 @@ end subroutine symba_discard_pl module subroutine symba_discard_tp(self, system, param) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(symba_tp), intent(inout) :: self !! RMVS test particle object + class(symba_tp), intent(inout) :: self !! SyMBA test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_discard_tp module function symba_encounter_check_pl(self, system, dt) result(lencounter) implicit none - class(symba_pl), intent(inout) :: self !! RMVS test particle object - class(symba_nbody_system), intent(inout) :: system !! RMVS nbody system object + class(symba_pl), intent(inout) :: self !! SyMBA test particle object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size logical :: lencounter !! Returns true if there is at least one close encounter end function symba_encounter_check_pl module function symba_encounter_check_tp(self, system, dt) result(lencounter) implicit none - class(symba_tp), intent(inout) :: self !! RMVS test particle object - class(symba_nbody_system), intent(inout) :: system !! RMVS nbody system object + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size logical :: lencounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp + module subroutine symba_io_dump_particle_info(self, param, msg) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_particle_info), intent(inout) :: self !! Swiftest base object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + character(*), optional, intent(in) :: msg !! Message to display with dump operation + end subroutine symba_io_dump_particle_info + + module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) + implicit none + class(symba_parameters), intent(inout) :: self !! Collection of parameters + integer, intent(in) :: unit !! File unit number + character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. + !! If you do not include a char-literal-constant, the iotype argument contains only DT. + integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader + integer, intent(out) :: iostat !! IO status code + character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 + end subroutine symba_io_param_reader + + module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, iomsg) + implicit none + class(symba_parameters),intent(in) :: self !! Collection of SyMBA parameters + integer, intent(in) :: unit !! File unit number + character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. + !! If you do not include a char-literal-constant, the iotype argument contains only DT. + integer, intent(in) :: v_list(:) !! Not used in this procedure + integer, intent(out) :: iostat !! IO status code + character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 + end subroutine symba_io_param_writer + + module subroutine symba_io_initialize_particle_info(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_io_initialize_particle_info + + module subroutine symba_io_read_cb_in(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_cb), intent(inout) :: self + class(swiftest_parameters), intent(inout) :: param + end subroutine symba_io_read_cb_in + + module subroutine symba_io_read_frame_cb(self, iu, param, form, ierr) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_cb), intent(inout) :: self !! Swiftest central body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code + end subroutine symba_io_read_frame_cb + + module subroutine symba_io_read_frame_info(self, iu, param, form, ierr) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code + end subroutine symba_io_read_frame_info + + module subroutine symba_io_read_frame_pl(self, iu, param, form, ierr) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_pl), intent(inout) :: self !! Swiftest particle object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code + end subroutine symba_io_read_frame_pl + + + module subroutine symba_io_read_pl_in(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_io_read_pl_in + + module subroutine symba_io_write_frame_cb(self, iu, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_cb), intent(in) :: self !! SyMBA massive body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine symba_io_write_frame_cb + + module subroutine symba_io_write_frame_info(self, iu, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_particle_info), intent(in) :: self !! SyMBA particle info object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine symba_io_write_frame_info + + module subroutine symba_io_write_frame_pl(self, iu, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_pl), intent(in) :: self !! SyMBA massive body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine symba_io_write_frame_pl + module subroutine symba_step_system(self, param, t, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(symba_nbody_system), intent(inout) :: self !! RMVS nbody system object + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody 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 @@ -96,7 +334,7 @@ end subroutine symba_step_system module subroutine symba_step_interp_system(self, param, t, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(symba_nbody_system), intent(inout) :: self !! RMVS nbody system object + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody 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 diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 6dacd7969..14613724e 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -22,7 +22,7 @@ module subroutine rmvs_discard_tp(self, system, param) if ((tp%status(i) == ACTIVE) .and. (tp%lperi(i))) then if ((tp%peri(i) < pl%radius(iplperP))) then tp%status(i) = DISCARDED_PLQ - write(*, *) "Particle ",tp%name(i)," q with respect to Planet ",pl%name(iplperP)," is too small at t = ",t + write(*, *) "Particle ",tp%id(i)," q with respect to Planet ",pl%id(iplperP)," is too small at t = ",t tp%ldiscard(i) = .true. end if end if diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 71091f6dd..035783901 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -100,7 +100,7 @@ subroutine rmvs_interp_out(cb, pl, dt) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" write(*, *) GMcb(i), dto(i) write(*, *) xtmp(:,i) write(*, *) vtmp(:,i) @@ -122,7 +122,7 @@ subroutine rmvs_interp_out(cb, pl, dt) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" write(*, *) GMcb(i), -dto(i) write(*, *) xtmp(:,i) write(*, *) vtmp(:,i) @@ -254,7 +254,7 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /=0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" write(*, *) GMcb(i), dti(i) write(*, *) xtmp(:,i) write(*, *) vtmp(:,i) @@ -278,7 +278,7 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /=0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) " Planet ", pl%id(i), " is lost!!!!!!!!!!" write(*, *) GMcb(i), -dti(i) write(*, *) xtmp(:,i) write(*, *) vtmp(:,i) @@ -413,7 +413,7 @@ subroutine rmvs_make_planetocentric(cb, pl, tp) tpenci%ipleP = i tpenci%status(:) = ACTIVE ! Grab all the encountering test particles and convert them to a planetocentric frame - tpenci%name(:) = pack(tp%name(:), encmask(:)) + tpenci%id(:) = pack(tp%id(:), encmask(:)) do j = 1, NDIM tpenci%xheliocentric(j, :) = pack(tp%xh(j,:), encmask(:)) tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i) @@ -498,11 +498,11 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) r2 = dot_product(xpc(:, i), xpc(:, i)) if ((abs(tperi) > FACQDT * dt) .or. (r2 > rhill2)) peri = sqrt(r2) if (param%encounter_file /= "") then - id1 = pl%name(ipleP) + id1 = pl%id(ipleP) rpl = pl%radius(ipleP) xh1(:) = pl%inner(inner_index)%x(:, ipleP) vh1(:) = pl%inner(inner_index)%v(:, ipleP) - id2 = tp%name(i) + id2 = tp%id(i) 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(:), & diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 0f16e7c5b..a64f64a15 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -72,7 +72,7 @@ module subroutine setup_body(self,n) self%lfirst = .true. !write(*,*) 'Allocating the basic Swiftest particle' - allocate(self%name(n)) + allocate(self%id(n)) allocate(self%status(n)) allocate(self%ldiscard(n)) allocate(self%xh(NDIM, n)) @@ -90,7 +90,7 @@ module subroutine setup_body(self,n) allocate(self%capm(n)) allocate(self%mu(n)) - self%name(:) = 0 + self%id(:) = 0 self%status(:) = INACTIVE self%ldiscard(:) = .false. self%xh(:,:) = 0.0_DP @@ -131,20 +131,12 @@ module subroutine setup_pl(self,n) allocate(self%rhill(n)) allocate(self%radius(n)) allocate(self%density(n)) - allocate(self%Ip(NDIM, n)) - allocate(self%rot(NDIM, n)) - allocate(self%k2(n)) - allocate(self%Q(n)) self%mass(:) = 0.0_DP self%Gmass(:) = 0.0_DP self%rhill(:) = 0.0_DP self%radius(:) = 0.0_DP self%density(:) = 0.0_DP - self%Ip(:,:) = 0.0_DP - self%rot(:,:) = 0.0_DP - self%k2(:) = 0.0_DP - self%Q(:) = 0.0_DP self%num_comparisons = 0 return end subroutine setup_pl diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 new file mode 100644 index 000000000..70123d5bd --- /dev/null +++ b/src/symba/symba_io.f90 @@ -0,0 +1,470 @@ +submodule (symba_classes) s_symba_io + use swiftest +contains + module subroutine symba_io_dump_particle_info(self, param, msg) + !! author: David A. Minton + !! + !! Dumps the particle information data to a file + implicit none + class(symba_particle_info), intent(inout) :: self !! Swiftest base object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + character(*), optional, intent(in) :: msg !! Message to display with dump operation + end subroutine symba_io_dump_particle_info + + module subroutine symba_io_initialize_particle_info(self, param) + !! author: David A. Minton + !! + !! Initializes a particle info data structure, either starting a new one or reading one in + !! from a file if it is a restarted run + implicit none + class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_io_initialize_particle_info + + module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Read in parameters specific to the SyMBA integrator, then calls the base io_param_reader. + !! + !! Adapted from David E. Kaufmann's Swifter routine io_init_param.f90 + !! Adapted from Martin Duncan's Swift routine io_init_param.f + implicit none + ! Arguments + class(symba_parameters), intent(inout) :: self !! Collection of parameters + integer, intent(in) :: unit !! File unit number + character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. + !! If you do not include a char-literal-constant, the iotype argument contains only DT. + integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader + integer, intent(out) :: iostat !! IO status code + character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 + ! internals + integer(I4B) :: ilength, ifirst, ilast !! Variables used to parse input file + character(STRMAX) :: line !! Line of the input file + character (len=:), allocatable :: line_trim,param_name, param_value !! Strings used to parse the param file + integer(I4B) :: nseeds, nseeds_from_file, i + logical :: seed_set = .false. !! Is the random seed set in the input file? + character(len=*),parameter :: linefmt = '(A)' + + associate(param => self) + call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) + + call random_seed(size = nseeds) + if (allocated(param%seed)) deallocate(param%seed) + allocate(param%seed(nseeds)) + do + read(unit = unit, fmt = linefmt, iostat = iostat, end = 1) line + line_trim = trim(adjustl(line)) + ilength = len(line_trim) + if ((ilength /= 0)) then + ifirst = 1 + ! Read the pair of tokens. The first one is the parameter name, the second is the value. + param_name = io_get_token(line_trim, ifirst, ilast, iostat) + if (param_name == '') cycle ! No parameter name (usually because this line is commented out) + call util_toupper(param_name) + ifirst = ilast + 1 + param_value = io_get_token(line_trim, ifirst, ilast, iostat) + select case (param_name) + case ("FRAGMENTATION") + call util_toupper(param_value) + if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. + case ("ROTATION") + call util_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') self%lrotation = .true. + case ("TIDES") + call util_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') self%ltides = .true. + case ("MTINY") + read(param_value, *) param%mtiny + case("SEED") + read(param_value, *) nseeds_from_file + ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different + ! number of seeds than the current system. If the number of seeds in the file is smaller than required, we will use them as a source to fill in the missing elements. + ! If the number of seeds in the file is larger than required, we will truncate the seed array. + if (nseeds_from_file > nseeds) then + nseeds = nseeds_from_file + deallocate(param%seed) + allocate(param%seed(nseeds)) + do i = 1, nseeds + ifirst = ilast + 1 + param_value = io_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%seed(i) + end do + else ! Seed array in file is too small + do i = 1, nseeds_from_file + ifirst = ilast + 1 + param_value = io_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%seed(i) + end do + param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, i=nseeds_from_file+1, nseeds)] + end if + seed_set = .true. + end select + end if + end do + 1 continue + + write(*,*) "FRAGMENTATION = ", param%lfragmentation + if (param%lfragmentation) then + if (seed_set) then + call random_seed(put = param%seed) + else + call random_seed(get = param%seed) + end if + write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) + end if + write(*,*) "ROTATION = ", param%lrotation + write(*,*) "TIDES = ", param%ltides + + if (self%mtiny < 0.0_DP) then + write(iomsg,*) "MTINY invalid or not set: ", self%mtiny + iostat = -1 + return + else + write(*,*) "MTINY = ", self%mtiny + end if + + if (.not.self%lclose) then + write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.' + iostat = -1 + return + end if + end associate + return + + end subroutine symba_io_param_reader + + module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, iomsg) + !! author: David A. Minton + !! + !! Dump integration parameters specific to SyMBA to file and then call the base io_param_writer method. + !! + !! Adapted from David E. Kaufmann's Swifter routine io_dump_param.f90 + !! Adapted from Martin Duncan's Swift routine io_dump_param.f + implicit none + ! Arguments + class(symba_parameters),intent(in) :: self !! Collection of SyMBA parameters + integer, intent(in) :: unit !! File unit number + character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. + !! If you do not include a char-literal-constant, the iotype argument contains only DT. + integer, intent(in) :: v_list(:) !! Not used in this procedure + integer, intent(out) :: iostat !! IO status code + character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 + ! Internals + character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values + character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values + character(*),parameter :: Rarrfmt = '(3(ES25.17,1X))' !! Format label for real values + character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values + character(len=*), parameter :: Afmt = '(A25,1X,64(:,A25,1X))' + character(256) :: param_name, param_value + type character_array + character(25) :: value + end type character_array + type(character_array), dimension(:), allocatable :: param_array + integer(I4B) :: i + + associate(param => self) + call io_param_writer(param, unit, iotype, v_list, iostat, iomsg) + + ! Special handling is required for writing the random number seed array as its size is not known until runtime + ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements + write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "MTINY"; write(param_value, Rfmt) param%mtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + if (param%lfragmentation) then + write(param_name, Afmt) "SEED" + if (allocated(param_array)) deallocate(param_array) + allocate(param_array(0:size(param%seed))) + write(param_array(0)%value, Ifmt) size(param%seed) + do i = 1, size(param%seed) + write(param_array(i)%value, Ifmt) param%seed(i) + end do + write(unit, Afmt, advance='no') adjustl(param_name), adjustl(param_array(0)%value) + do i = 1, size(param%seed) + if (i < size(param%seed)) then + write(unit, Afmt, advance='no') adjustl(param_array(i)%value) + else + write(unit, Afmt) adjustl(param_array(i)%value) + end if + end do + end if + + iostat = 0 + end associate + + return + + end subroutine symba_io_param_writer + + module subroutine symba_io_read_frame_cb(self, iu, param, form, ierr) + !! author: David A. Minton + !! + !! Reads a frame of output of central body data to the binary output file + !! + !! Adapted from David E. Kaufmann's Swifter routine io_read_frame.f90 + !! Adapted from Hal Levison's Swift routine io_read_frame.F + implicit none + ! Arguments + class(symba_cb), intent(inout) :: self !! Swiftest central body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code + + call io_read_frame_cb(self, iu, param, form, ierr) + select type(param) + class is (symba_parameters) + if (param%lrotation) then + read(iu, iostat = ierr) self%Ip(:) + read(iu, iostat = ierr) self%rot(:) + end if + if (param%ltides) then + read(iu, iostat = ierr) self%k2 + read(iu, iostat = ierr) self%Q + end if + end select + if (ierr /=0) then + write(*,*) 'Error reading SyMBA central body data' + call util_exit(FAILURE) + end if + return + end subroutine symba_io_read_frame_cb + + module subroutine symba_io_read_frame_pl(self, iu, param, form, ierr) + !! author: David A. Minton + !! + !! Reads a frame of output of a SyMBA massive body object + !! + !! Adapted from David E. Kaufmann's Swifter routine io_read_frame.f90 + !! Adapted from Hal Levison's Swift routine io_read_frame.F + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! Swiftest particle object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code + + call io_read_frame_body(self, iu, param, form, ierr) + select type(param) + class is (symba_parameters) + associate(pl => self, npl => self%nbody) + if (param%lrotation) then + read(iu, iostat = ierr) pl%rot(1, 1:npl) + read(iu, iostat = ierr) pl%rot(2, 1:npl) + read(iu, iostat = ierr) pl%rot(3, 1:npl) + read(iu, iostat = ierr) pl%Ip(1, 1:npl) + read(iu, iostat = ierr) pl%Ip(2, 1:npl) + read(iu, iostat = ierr) pl%Ip(3, 1:npl) + end if + if (param%ltides) then + read(iu, iostat = ierr) pl%k2(1:npl) + read(iu, iostat = ierr) pl%Q(1:npl) + end if + end associate + end select + + if (ierr /=0) then + write(*,*) 'Error reading SyMBA massive body body data' + call util_exit(FAILURE) + end if + return + end subroutine symba_io_read_frame_pl + + subroutine symba_io_read_frame_info(self, iu, param, form, ierr) + !! author: David A. Minton + !! + !! Reads a single frame of a particle info data from a file. + implicit none + class(symba_particle_info), intent(inout) :: self !! SyMBA particle info object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code + end subroutine symba_io_read_frame_info + + module subroutine symba_io_read_cb_in(self, param) + !! author: David A. Minton + !! + !! Reads in central body data + !! + !! Adapted from David E. Kaufmann's Swifter routine swiftest_init_pl.f90 + !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f + implicit none + ! Arguments + class(symba_cb), intent(inout) :: self + class(swiftest_parameters), intent(inout) :: param + ! Internals + integer(I4B), parameter :: LUN = 7 !! Unit number of input file + integer(I4B) :: iu = LUN + integer(I4B) :: ierr + logical :: is_ascii + real(DP) :: t + real(QP) :: val + + select type(param) + class is (symba_parameters) + ierr = 0 + is_ascii = (param%in_type == 'ASCII') + if (is_ascii) then + open(unit = iu, file = param%incbfile, status = 'old', form = 'FORMATTED', iostat = ierr) + read(iu, *, iostat = ierr) val + self%Gmass = real(val, kind=DP) + self%mass = real(val / param%GU, kind=DP) + read(iu, *, iostat = ierr) self%radius + read(iu, *, iostat = ierr) self%j2rp2 + read(iu, *, iostat = ierr) self%j4rp4 + if (param%lrotation) then + read(iu, *, iostat = ierr) self%Ip(:) + read(iu, *, iostat = ierr) self%rot(:) + end if + if (param%ltides) then + read(iu, *, iostat = ierr) self%k2 + read(iu, *, iostat = ierr) self%Q + end if + else + open(unit = iu, file = param%incbfile, status = 'old', form = 'UNFORMATTED', iostat = ierr) + call self%read_frame(iu, param, XV, ierr) + end if + close(iu) + if (ierr /= 0) then + write(*,*) 'Error opening massive body initial conditions file ',trim(adjustl(param%incbfile)) + call util_exit(FAILURE) + end if + if (self%j2rp2 /= 0.0_DP) param%loblatecb = .true. + end select + + return + end subroutine symba_io_read_cb_in + + module subroutine symba_io_read_pl_in(self, param) + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Read in either test particle or massive body data + !! + !! Adapted from David E. Kaufmann's Swifter routine swiftest_init_pl.f90 and swiftest_init_tp.f90 + !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f and swiftest_init_tp.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B), parameter :: LUN = 7 !! Unit number of input file + integer(I4B) :: iu = LUN + integer(I4B) :: i, ierr, nbody + logical :: is_ascii + character(len=:), allocatable :: infile + real(DP) :: t + real(QP) :: val + + select type(param) + class is (symba_parameters) + ierr = 0 + is_ascii = (param%in_type == 'ASCII') + select case(param%in_type) + case(ASCII_TYPE) + open(unit = iu, file = infile, status = 'old', form = 'FORMATTED', iostat = ierr) + read(iu, *, iostat = ierr) nbody + call self%setup(nbody) + if (nbody > 0) then + do i = 1, nbody + read(iu, *, iostat = ierr) self%id(i), val + self%mass(i) = real(val / param%GU, kind=DP) + self%Gmass(i) = real(val, kind=DP) + read(iu, *, iostat = ierr) self%radius(i) + if (param%lrotation) then + read(iu, iostat = ierr) self%Ip(:, i) + read(iu, iostat = ierr) self%rot(:, i) + end if + if (param%ltides) then + read(iu, iostat = ierr) self%k2(i) + read(iu, iostat = ierr) self%Q(i) + end if + if (ierr /= 0 ) exit + read(iu, *, iostat = ierr) self%xh(1, i), self%xh(2, i), self%xh(3, i) + read(iu, *, iostat = ierr) self%vh(1, i), self%vh(2, i), self%vh(3, i) + if (ierr /= 0 ) exit + self%status(i) = ACTIVE + end do + end if + case (REAL4_TYPE, REAL8_TYPE) + open(unit = iu, file = infile, status = 'old', form = 'UNFORMATTED', iostat = ierr) + read(iu, iostat = ierr) nbody + call self%setup(nbody) + if (nbody > 0) then + call self%read_frame(iu, param, XV, ierr) + self%status(:) = ACTIVE + end if + case default + write(*,*) trim(adjustl(param%in_type)) // ' is an unrecognized file type' + ierr = -1 + end select + close(iu) + if (ierr /= 0 ) then + write(*,*) 'Error reading in initial conditions from ',trim(adjustl(infile)) + call util_exit(FAILURE) + end if + end select + + return + end subroutine symba_io_read_pl_in + + subroutine symba_io_write_frame_cb(self, iu, param) + !! author: David A. Minton + !! + !! Writes a single frame of a SyMBA pl file + implicit none + class(symba_cb), intent(in) :: self !! SyMBA massive body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + + call io_write_frame_cb(self, iu, param) + select type(param) + class is (symba_parameters) + if (param%lrotation) then + write(iu) self%rot(1) + write(iu) self%rot(2) + write(iu) self%rot(3) + write(iu) self%Ip(1) + write(iu) self%Ip(2) + write(iu) self%Ip(3) + end if + if (param%ltides) then + write(iu) self%k2 + write(iu) self%Q + end if + end select + end subroutine symba_io_write_frame_cb + + subroutine symba_io_write_frame_pl(self, iu, param) + !! author: David A. Minton + !! + !! Writes a single frame of a SyMBA pl file + implicit none + class(symba_pl), intent(in) :: self !! SyMBA massive body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + + call io_write_frame_body(self, iu, param) + select type(param) + class is (symba_parameters) + associate(pl => self, npl => self%nbody) + if (param%lrotation) then + write(iu) pl%rot(1, 1:npl) + write(iu) pl%rot(2, 1:npl) + write(iu) pl%rot(3, 1:npl) + write(iu) pl%Ip(1, 1:npl) + write(iu) pl%Ip(2, 1:npl) + write(iu) pl%Ip(3, 1:npl) + end if + if (param%ltides) then + write(iu) pl%k2(1:npl) + write(iu) pl%Q(1:npl) + end if + end associate + end select + end subroutine symba_io_write_frame_pl + +end submodule s_symba_io + diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 deleted file mode 100644 index 9731e9eff..000000000 --- a/src/util/util_copy.f90 +++ /dev/null @@ -1,154 +0,0 @@ -submodule (swiftest_classes) s_util_copy - use swiftest -contains - module procedure util_copy_cb - !! author: David A. Minton - !! - !! Copies elements of one Swiftest central body object to another. The mask is ignored for this case - !! - implicit none - select type(src) - class is(swiftest_cb) - self%mass = src%mass - self%Gmass = src%Gmass - self%radius = src%radius - self%density = src%density - self%j2rp2 = src%j2rp2 - self%j4rp4 = src%j4rp4 - self%aobl = src%aobl - self%xb = src%xb - self%vb = src%vb - self%Ip = src%Ip - self%rot = src%rot - self%k2 = src%k2 - self%Q = src%Q - class default - write(*,*) 'Error copying swiftest_cb object. Incompatible type.' - call util_exit(FAILURE) - end select - return - end procedure util_copy_cb - - module procedure util_copy_body - !! author: David A. Minton - !! - !! Copies elements of one Swiftest generic body object to another. The elements are determined by a mask. - !! - implicit none - - select type(src) - class is(swiftest_body) - where (mask(:)) - self%name(:) = src%name(:) - self%status(:) = src%status(:) - self%ldiscard(:) = src%ldiscard(:) - self%xh(1,:) = src%xh(1,:) - self%xh(2,:) = src%xh(2,:) - self%xh(3,:) = src%xh(3,:) - self%vh(1,:) = src%vh(1,:) - self%vh(2,:) = src%vh(2,:) - self%vh(3,:) = src%vh(3,:) - self%xb(1,:) = src%xb(1,:) - self%xb(2,:) = src%xb(2,:) - self%xb(3,:) = src%xb(3,:) - self%vb(1,:) = src%vb(1,:) - self%vb(2,:) = src%vb(2,:) - self%vb(3,:) = src%vb(3,:) - self%ah(1,:) = src%ah(1,:) - self%ah(2,:) = src%ah(2,:) - self%ah(3,:) = src%ah(3,:) - self%aobl(1,:) = src%aobl(1,:) - self%aobl(3,:) = src%aobl(2,:) - self%aobl(2,:) = src%aobl(3,:) - self%ir3h(:) = src%ir3h(:) - self%a(:) = src%a(:) - self%e(:) = src%e(:) - self%inc(:) = src%inc(:) - self%capom(:) = src%capom(:) - self%omega(:) = src%omega(:) - self%capm(:) = src%capm(:) - self%mu(:) = src%mu(:) - end where - class default - write(*,*) 'Error copying swiftest_body object. Incompatible type.' - call util_exit(FAILURE) - end select - - end procedure util_copy_body - - module procedure util_copy_pl - !! author: David A. Minton - !! - !! Copies elements of one Swiftest massive body object to another. The elements are determined by a mask. - !! - implicit none - select type(src) - class is(swiftest_pl) - where (mask(:)) - self%mass(:) = src%mass(:) - self%Gmass(:) = src%Gmass(:) - self%rhill(:) = src%rhill(:) - self%radius(:) = src%radius(:) - self%density(:) = src%density(:) - self%Ip(1,:) = src%Ip(1,:) - self%Ip(2,:) = src%Ip(2,:) - self%Ip(3,:) = src%Ip(3,:) - self%rot(1,:) = src%rot(1,:) - self%rot(2,:) = src%rot(2,:) - self%rot(3,:) = src%rot(3,:) - self%k2(:) = src%k2(:) - self%Q(:) = src%Q(:) - end where - call util_copy_body(self, src, mask) - class default - write(*,*) 'Error copying swiftest_pl object. Incompatible type.' - call util_exit(FAILURE) - end select - - end procedure util_copy_pl - - module procedure util_copy_tp - !! author: David A. Minton - !! - !! Copies elements of one swiftest test particle object to another. The elements are determined by a mask. - implicit none - - select type(src) - class is(swiftest_tp) - where (mask(:)) - self%isperi(:) = src%isperi(:) - self%peri(:) = src%peri(:) - self%atp(:) = src%atp(:) - end where - call util_copy_body(self, src, mask) - class default - write(*,*) 'Error copying swiftest_tp object. Incompatible type.' - call util_exit(FAILURE) - end select - - return - end procedure util_copy_tp - - module procedure util_copy_system - !! author: David A. Minton - !! - !! Copies elements of one Swiftest nbody system object to another. The elements are determined by a mask. - !! In this case the mask contains the pl elements followed by the tp elements. - !! - implicit none - - associate(pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) - - select type(src) - class is(swiftest_nbody_system) - call pl%copy(src%pl, mask(1:npl)) - call tp%copy(src%tp, mask(npl+1:npl+ntp)) - class default - write(*,*) 'Error copying swiftest_nbody_system object. Incompatible type.' - call util_exit(FAILURE) - end select - - - end associate - end procedure util_copy_system -end submodule s_util_copy \ No newline at end of file diff --git a/src/util/util_spill_and_fill.f90 b/src/util/util_spill_and_fill.f90 index da00ae72a..4513e5b24 100644 --- a/src/util/util_spill_and_fill.f90 +++ b/src/util/util_spill_and_fill.f90 @@ -17,7 +17,7 @@ module subroutine util_spill_body(self, discards, lspill_list) ! For each component, pack the discarded bodies into the discard object and do the inverse with the keeps !> Spill all the common components associate(keeps => self) - discards%name(:) = pack(keeps%name(:), lspill_list(:)) + discards%id(:) = pack(keeps%id(:), lspill_list(:)) discards%status(:) = pack(keeps%status(:), lspill_list(:)) discards%a(:) = pack(keeps%a(:), lspill_list(:)) discards%e(:) = pack(keeps%e(:), lspill_list(:)) @@ -34,7 +34,7 @@ module subroutine util_spill_body(self, discards, lspill_list) discards%aobl(i, :) = pack(keeps%aobl(i, :), lspill_list(:)) end do if (count(.not.lspill_list(:)) > 0) then - keeps%name(:) = pack(keeps%name(:), .not. lspill_list(:)) + keeps%id(:) = pack(keeps%id(:), .not. lspill_list(:)) keeps%status(:) = pack(keeps%status(:), .not. lspill_list(:)) keeps%a(:) = pack(keeps%a(:), .not. lspill_list(:)) keeps%e(:) = pack(keeps%e(:), .not. lspill_list(:)) @@ -82,10 +82,10 @@ module subroutine util_fill_body(self, inserts, lfill_list) ! For each component, pack the discarded bodies into the discard object and do the inverse with the keeps !> Spill all the common components - associate(keeps => self, insname => inserts%name, keepname => self%name) + associate(keeps => self, insname => inserts%id, keepname => self%id) - keeps%name(:) = unpack(keeps%name(:), .not.lfill_list(:), keeps%name(:)) - keeps%name(:) = unpack(inserts%name(:), lfill_list(:), keeps%name(:)) + keeps%id(:) = unpack(keeps%id(:), .not.lfill_list(:), keeps%id(:)) + keeps%id(:) = unpack(inserts%id(:), lfill_list(:), keeps%id(:)) keeps%status(:) = unpack(keeps%status(:), .not.lfill_list(:), keeps%status(:)) keeps%status(:) = unpack(inserts%status(:), lfill_list(:), keeps%status(:)) @@ -137,7 +137,7 @@ module subroutine util_fill_body(self, inserts, lfill_list) ! This is the base class, so will be the last to be called in the cascade. - keeps%nbody = size(keeps%name(:)) + keeps%nbody = size(keeps%id(:)) end associate end subroutine util_fill_body @@ -161,24 +161,12 @@ end subroutine util_fill_body discards%rhill(:) = pack(keeps%rhill(:), lspill_list(:)) discards%radius(:) = pack(keeps%radius(:), lspill_list(:)) discards%density(:) = pack(keeps%density(:), lspill_list(:)) - discards%k2(:) = pack(keeps%k2(:), lspill_list(:)) - discards%Q(:) = pack(keeps%Q(:), lspill_list(:)) - do i = 1, NDIM - discards%Ip(i, :) = pack(keeps%Ip(i, :), lspill_list(:)) - discards%rot(i, :) = pack(keeps%rot(i, :), lspill_list(:)) - end do if (count(.not.lspill_list(:)) > 0) then keeps%mass(:) = pack(keeps%mass(:), .not. lspill_list(:)) keeps%Gmass(:) = pack(keeps%Gmass(:), .not. lspill_list(:)) keeps%rhill(:) = pack(keeps%rhill(:), .not. lspill_list(:)) keeps%radius(:) = pack(keeps%radius(:), .not. lspill_list(:)) keeps%density(:) = pack(keeps%density(:), .not. lspill_list(:)) - keeps%k2(:) = pack(keeps%k2(:), .not. lspill_list(:)) - keeps%Q(:) = pack(keeps%Q(:), .not. lspill_list(:)) - do i = 1, NDIM - keeps%Ip(i, :) = pack(keeps%Ip(i, :), .not. lspill_list(:)) - keeps%rot(i, :) = pack(keeps%rot(i, :), .not. lspill_list(:)) - end do end if call util_spill_body(keeps, discards, lspill_list) @@ -219,20 +207,6 @@ end subroutine util_fill_body keeps%density(:) = unpack(keeps%density(:),.not.lfill_list(:), keeps%density(:)) keeps%density(:) = unpack(inserts%density(:),lfill_list(:), keeps%density(:)) - do i = 1, NDIM - keeps%Ip(i, :) = unpack(keeps%Ip(i, :), .not.lfill_list(:), keeps%Ip(i, :)) - keeps%Ip(i, :) = unpack(inserts%Ip(i, :), lfill_list(:), keeps%Ip(i, :)) - - keeps%rot(i, :) = unpack(keeps%rot(i, :), .not.lfill_list(:), keeps%rot(i, :)) - keeps%rot(i, :) = unpack(inserts%rot(i, :), lfill_list(:), keeps%rot(i, :)) - - end do - keeps%k2(:) = unpack(keeps%k2(:), .not.lfill_list(:), keeps%k2(:)) - keeps%k2(:) = unpack(inserts%k2(:), lfill_list(:), keeps%k2(:)) - - keeps%Q(:) = unpack(keeps%Q(:), .not.lfill_list(:), keeps%Q(:)) - keeps%Q(:) = unpack(inserts%Q(:), lfill_list(:), keeps%Q(:)) - call util_fill_body(keeps, inserts, lfill_list) class default write(*,*) 'Error! fill method called for incompatible return type on swiftest_pl' diff --git a/src/util/util_valid.f90 b/src/util/util_valid.f90 index 7dd755646..ec1110025 100644 --- a/src/util/util_valid.f90 +++ b/src/util/util_valid.f90 @@ -15,10 +15,10 @@ associate(npl => pl%nbody, ntp => tp%nbody) allocate(idarr(npl+ntp)) do i = 1, npl - idarr(i) = pl%name(i) + idarr(i) = pl%id(i) end do do i = 1, ntp - idarr(npl+i) = tp%name(i) + idarr(npl+i) = tp%id(i) end do call util_sort(idarr) do i = 1, npl + ntp - 1 diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index 0ac170991..a27897cfa 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -45,7 +45,7 @@ module subroutine whm_drift_pl(self, system, param, dt) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then - write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" + write(*, *) " Planet ", self%id(i), " is lost!!!!!!!!!!" write(*, *) pl%xj(:,i) write(*, *) pl%vj(:,i) write(*, *) " stopping " @@ -102,7 +102,7 @@ module subroutine whm_drift_tp(self, system, param, dt) if (any(iflag(1:ntp) /= 0)) then where(iflag(1:ntp) /= 0) tp%status(1:ntp) = DISCARDED_DRIFTERR do i = 1, ntp - if (iflag(i) /= 0) write(*, *) "Particle ", self%name(i), " lost due to error in Danby drift" + if (iflag(i) /= 0) write(*, *) "Particle ", self%id(i), " lost due to error in Danby drift" end do end if end associate