Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Moved rotation and tide variables to base swiftest_pl and swiftest_cb…
Browse files Browse the repository at this point in the history 
… classes
  • Loading branch information
@daminton
daminton committed Jul 13, 2021
1 parent 211fe0d commit 1c8d750
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 349 deletions.
97 changes: 83 additions & 14 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,13 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
case ("GR")
call io_toupper(param_value)
if (param_value == "YES" .or. param_value == 'T') self%lgr = .true.
case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "ROTATION", "TIDES", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
case ("ROTATION")
call io_toupper(param_value)
if (param_value == "YES" .or. param_value == 'T') self%lrotation = .true.
case ("TIDES")
call io_toupper(param_value)
if (param_value == "YES" .or. param_value == 'T') self%ltides = .true.
case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
case default
write(iomsg,*) "Unknown parameter -> ",param_name
iostat = -1
Expand Down Expand Up @@ -192,6 +198,11 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
return
end if
end if
if (self%ltides .and. .not. self%lrotation) then
write(iomsg,*) 'Tides require rotation to be turned on'
iostat = -1
return
end if

write(*,*) "T0 = ",self%t0
write(*,*) "TSTOP = ",self%tstop
Expand All @@ -217,6 +228,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
write(*,*) "EXTRA_FORCE = ",self%lextra_force
write(*,*) "BIG_DISCARD = ",self%lbig_discard
write(*,*) "RHILL_PRESENT = ",self%lrhill_present
write(*,*) "ROTATION = ", self%lrotation
write(*,*) "TIDES = ", self%ltides
write(*,*) "ENERGY = ",self%lenergy
write(*,*) "MU2KG = ",self%MU2KG
write(*,*) "TU2S = ",self%TU2S
Expand Down Expand Up @@ -325,6 +338,8 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
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)
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)
iostat = 0
iomsg = "UDIO not implemented"
end associate
Expand Down Expand Up @@ -610,6 +625,14 @@ 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%id(i)
end select
Expand Down Expand Up @@ -811,11 +834,23 @@ module subroutine io_read_frame_body(self, iu, param, form, ierr)
read(iu, iostat = ierr) self%vh(2, 1:n)
read(iu, iostat = ierr) self%vh(3, 1:n)
end select
select type(self)
select type(pl => self)
class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body
read(iu, iostat = ierr) self%Gmass(1:n)
self%mass(1:n) = self%Gmass / param%GU
read(iu, iostat = ierr) self%radius(1:n)
read(iu, iostat = ierr) pl%Gmass(1:n)
pl%mass(1:n) = pl%Gmass / param%GU
read(iu, iostat = ierr) pl%radius(1:n)
if (param%lrotation) then
read(iu, iostat = ierr) pl%rot(1, 1:n)
read(iu, iostat = ierr) pl%rot(2, 1:n)
read(iu, iostat = ierr) pl%rot(3, 1:n)
read(iu, iostat = ierr) pl%Ip(1, 1:n)
read(iu, iostat = ierr) pl%Ip(2, 1:n)
read(iu, iostat = ierr) pl%Ip(3, 1:n)
end if
if (param%ltides) then
read(iu, iostat = ierr) pl%k2(1:n)
read(iu, iostat = ierr) pl%Q(1:n)
end if
end select
end associate

Expand Down Expand Up @@ -849,6 +884,14 @@ module subroutine io_read_frame_cb(self, iu, param, form, ierr)
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)
Expand Down Expand Up @@ -1133,10 +1176,22 @@ module subroutine io_write_frame_body(self, iu, param)
write(iu) self%vh(2, 1:n)
write(iu) self%vh(3, 1:n)
end select
select type(self)
select type(pl => self)
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)
write(iu) pl%Gmass(1:n)
write(iu) pl%radius(1:n)
if (param%lrotation) then
write(iu) pl%rot(1, 1:n)
write(iu) pl%rot(2, 1:n)
write(iu) pl%rot(3, 1:n)
write(iu) pl%Ip(1, 1:n)
write(iu) pl%Ip(2, 1:n)
write(iu) pl%Ip(3, 1:n)
end if
if (param%ltides) then
write(iu) pl%k2(1:n)
write(iu) pl%Q(1:n)
end if
end select
end associate

Expand All @@ -1156,12 +1211,26 @@ 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
associate(cb => self)
!write(iu) cb%name
write(iu) cb%id
write(iu) cb%Gmass
write(iu) cb%radius
write(iu) cb%j2rp2
write(iu) cb%j4rp4
if (param%lrotation) then
write(iu) cb%rot(1)
write(iu) cb%rot(2)
write(iu) cb%rot(3)
write(iu) cb%Ip(1)
write(iu) cb%Ip(2)
write(iu) cb%Ip(3)
end if
if (param%ltides) then
write(iu) cb%k2
write(iu) cb%Q
end if
end associate
return
end subroutine io_write_frame_cb

Expand Down
12 changes: 12 additions & 0 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,8 @@ module swiftest_classes
logical :: lclose = .false. !! Turn on close encounters
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)
logical :: lrotation = .false. !! Include rotation states of big bodies
logical :: ltides = .false. !! Include tidal dissipation

! Future features not implemented or in development
logical :: lgr = .false. !! Turn on GR
Expand Down Expand Up @@ -114,6 +116,12 @@ module swiftest_classes
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) :: agr = 0.0_DP !! Acceleration due to post-Newtonian correction
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
contains
private
procedure, public :: initialize => io_read_cb_in !! I/O routine for reading in central body data
Expand Down Expand Up @@ -192,6 +200,10 @@ module swiftest_classes
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)
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
!! 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
Expand Down
69 changes: 0 additions & 69 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,6 @@ module symba_classes
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 :: reader => symba_io_param_reader
Expand All @@ -33,21 +31,12 @@ module symba_classes
!*******************************************************************************************************************************
!> 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

!********************************************************************************************************************************
Expand Down Expand Up @@ -83,11 +72,6 @@ module symba_classes
!*******************************************************************************************************************************
!> 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
Expand All @@ -103,9 +87,6 @@ module symba_classes
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
procedure, public :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle
end type symba_pl

Expand Down Expand Up @@ -234,23 +215,6 @@ module subroutine symba_io_initialize_particle_info(self, param)
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
Expand All @@ -261,31 +225,6 @@ module subroutine symba_io_read_frame_info(self, iu, param, form, ierr)
integer(I4B), intent(out) :: ierr !! Error code
end subroutine symba_io_read_frame_info

module subroutine symba_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
Expand All @@ -294,14 +233,6 @@ module subroutine symba_io_write_frame_info(self, iu, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine symba_io_write_frame_info

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_setup_pl(self,n)
implicit none
class(symba_pl), intent(inout) :: self !! SyMBA test particle object
Expand Down
Loading

0 comments on commit 1c8d750

Please sign in to comment.