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

Commit

Permalink
Added methods and structures for reading in NetCDF files. It hasn't b…
Browse files Browse the repository at this point in the history 
…een tested, but it compiles
  • Loading branch information
@daminton
daminton committed Sep 24, 2021
1 parent d497a30 commit 195fd1a
Show file tree
Hide file tree
Showing 4 changed files with 147 additions and 86 deletions.
68 changes: 63 additions & 5 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1154,7 +1154,42 @@ module subroutine io_param_writer_one_QP(param_name, param_value, unit)
end subroutine io_param_writer_one_QP


module subroutine io_read_in_body(self, param)
module subroutine io_read_in_base(self,param)
!! author: Carlisle A. Wishard and David A. Minton
!!
!! Reads in either a central body, test particle, or massive body object. For the swiftest_body types (non-central body), it allocates array space for them
implicit none
class(swiftest_base), intent(inout) :: self !! Swiftest base object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: ierr !! Error code: returns 0 if the read is successful

select case(param%in_type)
case(NETCDF_DOUBLE_TYPE, NETCDF_FLOAT_TYPE)
select type(self)
class is (swiftest_body)
if (self%nbody == 0) return
call self%setup(self%nbody, param)
end select

ierr = self%read_frame(param%nciu, param)
if (ierr == 0) return
case default
select type(self)
class is (swiftest_body)
call io_read_in_body(self, param)
class is (swiftest_cb)
call io_read_in_cb(self, param)
end select
return
end select

667 continue
write(*,*) "Error reading body in io_read_in_base"
end subroutine io_read_in_base


subroutine io_read_in_body(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
Expand All @@ -1168,20 +1203,19 @@ module subroutine io_read_in_body(self, param)
! Internals
integer(I4B) :: iu = LUN
integer(I4B) :: i, nbody
logical :: is_ascii, is_pl
logical :: is_ascii
character(len=:), allocatable :: infile
real(DP) :: t
character(STRMAX) :: errmsg
integer(I4B) :: ierr

! Select the appropriate polymorphic class (test particle or massive body)

select type(self)
class is (swiftest_pl)
infile = param%inplfile
is_pl = .true.
class is (swiftest_tp)
infile = param%intpfile
is_pl = .false.
end select

is_ascii = (param%in_type == 'ASCII')
Expand Down Expand Up @@ -1218,7 +1252,7 @@ module subroutine io_read_in_body(self, param)
end subroutine io_read_in_body


module subroutine io_read_in_cb(self, param)
subroutine io_read_in_cb(self, param)
!! author: David A. Minton
!!
!! Reads in central body data
Expand Down Expand Up @@ -1279,6 +1313,30 @@ module subroutine io_read_in_cb(self, param)
end subroutine io_read_in_cb


module subroutine io_read_in_system(self, param)
!! author: David A. Minton and Carlisle A. Wishard
!!
!! Reads in the system from input files
implicit none
! Arguments
class(swiftest_nbody_system), intent(inout) :: self
class(swiftest_parameters), intent(inout) :: param
! Internals
integer(I4B) :: ierr

if ((param%in_type == NETCDF_DOUBLE_TYPE) .or. (param%in_type == NETCDF_FLOAT_TYPE)) then
ierr = self%read_frame(param%nciu, param)
if (ierr /=0) call util_exit(FAILURE)
else
call self%cb%read_in(param)
call self%pl%read_in(param)
call self%tp%read_in(param)
end if

return
end subroutine io_read_in_system


function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, &
xh1, xh2, vh1, vh2, enc_out, out_type) result(ierr)
!! author: David A. Minton
Expand Down
23 changes: 12 additions & 11 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,6 +200,7 @@ module swiftest_classes
!! The minimal methods that all systems must have
procedure :: dump => io_dump_base !! Dump contents to file
procedure :: dump_particle_info => io_dump_particle_info_base !! Dump contents of particle information metadata to file
procedure :: read_in => io_read_in_base !! Read in body initial conditions from a file
procedure :: write_frame_netcdf => netcdf_write_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format
procedure :: read_frame_netcdf => netcdf_read_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format
procedure :: write_particle_info_netcdf => netcdf_write_particle_info_base !! Writes out the particle information metadata to NetCDF file
Expand Down Expand Up @@ -238,7 +239,6 @@ module swiftest_classes
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
procedure :: read_in => io_read_in_cb !! I/O routine for reading in central body data
procedure :: read_frame_bin => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body
procedure :: write_frame_bin => io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body
generic :: write_frame => write_frame_bin !! Write a frame (either binary or NetCDF, using generic procedures)
Expand Down Expand Up @@ -286,7 +286,6 @@ module swiftest_classes
procedure :: drift => drift_body !! Loop through bodies and call Danby drift routine on heliocentric variables
procedure :: v2pv => gr_vh2pv_body !! Converts from velocity to psudeovelocity for GR calculations using symplectic integrators
procedure :: pv2v => gr_pv2vh_body !! Converts from psudeovelocity to velocity for GR calculations using symplectic integrators
procedure :: read_in => io_read_in_body !! Read in body initial conditions from a file
procedure :: read_frame_bin => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body
procedure :: write_frame_bin => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body
procedure :: accel_obl => obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
Expand Down Expand Up @@ -429,6 +428,7 @@ module swiftest_classes
!procedure :: read_hdr_bin => io_read_hdr !! Read a header for an output frame in Fortran binary format
procedure :: read_hdr_netcdf => netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format
procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format
procedure :: read_in => io_read_in_system !! Reads the initial conditions for an nbody system
procedure :: read_particle_info => io_read_particle_info_system !! Read in particle metadata from file
procedure :: write_discard => io_write_discard !! Write out information about discarded test particles
procedure :: obl_pot => obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body
Expand Down Expand Up @@ -819,17 +819,12 @@ end subroutine io_param_writer_one_QP
end interface io_param_writer_one

interface
module subroutine io_read_in_body(self, param)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
end subroutine io_read_in_body

module subroutine io_read_in_cb(self, param)
module subroutine io_read_in_base(self,param)
implicit none
class(swiftest_cb), intent(inout) :: self !! Swiftest central body object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
end subroutine io_read_in_cb
class(swiftest_base), intent(inout) :: self !! Swiftest base object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
end subroutine io_read_in_base

module subroutine io_read_in_param(self, param_file_name)
implicit none
Expand All @@ -843,6 +838,12 @@ module subroutine io_read_in_particle_info(self, iu)
integer(I4B), intent(in) :: iu !! Open file unit number
end subroutine io_read_in_particle_info

module subroutine io_read_in_system(self, param)
implicit none
class(swiftest_nbody_system), intent(inout) :: self
class(swiftest_parameters), intent(inout) :: param
end subroutine io_read_in_system

module function io_read_frame_body(self, iu, param) result(ierr)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
Expand Down
131 changes: 71 additions & 60 deletions src/netcdf/netcdf.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,7 @@ module subroutine netcdf_open(self, param)
return
end subroutine netcdf_open


module function netcdf_read_frame_base(self, iu, param) result(ierr)
!! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
!!
Expand All @@ -311,51 +312,68 @@ module function netcdf_read_frame_base(self, iu, param) result(ierr)
class(netcdf_parameters), intent(inout) :: iu !! Parameters used to for writing a NetCDF dataset to file
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
! Internals
integer(I4B) :: i, j, tslot, strlen, idslot
integer(I4B) :: i, j, tslot, strlen, idslot, idmax, n
integer(I4B), dimension(:), allocatable :: ind
character(len=:), allocatable :: charstring
integer(I4B) :: ierr !! Error code: returns 0 if the read is successful

call self%write_particle_info(iu)
real(DP), dimension(:), allocatable :: real_temp
logical, dimension(:), allocatable :: validmask, tpmask

tslot = int(param%ioutput, kind=I4B) + 1

select type(self)
class is (swiftest_body)
associate(n => self%nbody)
if (n == 0) return

allocate(ind(n))
call util_sort(self%id(1:n), ind)

do i = 1, n
j = ind(i)
idslot = self%id(j) + 1

if ((param%out_form == XV) .or. (param%out_form == XVEL)) then
if (self%nbody == 0) return
call check( nf90_inquire_dimension(iu%ncid, iu%id_dimid, len=idmax) )
allocate(ind(idmax))
allocate(real_temp(idmax))
allocate(validmask(idmax))
allocate(tpmask(idmax))
ind(:) = [(i, i = 1, idmax)]

! First filter out only the id slots that contain valid bodies
if (param%in_form == XV) then
call check( nf90_get_var(iu%ncid, iu%xhx_varid, real_temp(:)) )
else
call check( nf90_get_var(iu%ncid, iu%a_varid, real_temp(:)) )
end if
validmask(:) = real_temp(:) == real_temp(:)

! Filter out the central body, which is always in id dimension array position 1
validmask(1) = .false.

! Next, filter only bodies that don't have mass (test particles)
call check( nf90_get_var(iu%ncid, iu%Gmass_varid, real_temp(:)) )
tpmask(:) = real_temp(:) /= real_temp(:)

select type(self)
class is (swiftest_pl)
ind(:) = pack(ind(:), (.not.tpmask(:) .and. validmask(:)))
n = count(.not.tpmask(:))
class is (swiftest_tp)
ind(:) = pack(ind(:), (tpmask(:) .and. validmask(:)))
n = count(tpmask(:))
end select

do i = j, n
idslot = ind(j)

if (param%in_form == XV) then
call check( nf90_get_var(iu%ncid, iu%xhx_varid, self%xh(1, j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%xhy_varid, self%xh(2, j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%xhz_varid, self%xh(3, j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%vhx_varid, self%vh(1, j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%vhy_varid, self%vh(2, j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%vhz_varid, self%vh(3, j), start=[idslot, tslot]) )
end if

if ((param%out_form == EL) .or. (param%out_form == XVEL)) then
else if (param%in_form == EL) then
call check( nf90_get_var(iu%ncid, iu%a_varid, self%a(j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%e_varid, self%e(j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%inc_varid, self%inc(j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%capom_varid, self%capom(j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%omega_varid, self%omega(j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%capm_varid, self%capm(j), start=[idslot, tslot]) )

self%inc(j) = self%inc(j) * RAD2DEG
self%capom(j) = self%capom(j) * RAD2DEG
self%omega(j) = self%omega(j) * RAD2DEG
self%capm(j) = self%capm(j) * RAD2DEG

end if

select type(self)
class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body
call check( nf90_get_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) )
Expand All @@ -377,35 +395,40 @@ module function netcdf_read_frame_base(self, iu, param) result(ierr)
call check( nf90_get_var(iu%ncid, iu%k2_varid, self%k2(j), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%Q_varid, self%Q(j), start=[idslot, tslot]) )
end if

end select
end do
end associate
class is (swiftest_cb)
idslot = self%id + 1
call check( nf90_get_var(iu%ncid, iu%id_varid, self%id, start=[idslot]) )

call check( nf90_get_var(iu%ncid, iu%Gmass_varid, self%Gmass, start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%radius_varid, self%radius, start=[idslot, tslot]) )
if (param%lrotation) then
call check( nf90_get_var(iu%ncid, iu%Ip1_varid, self%Ip(1), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%Ip2_varid, self%Ip(2), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%Ip3_varid, self%Ip(3), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%rotx_varid, self%rot(1), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%roty_varid, self%rot(2), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%rotz_varid, self%rot(3), start=[idslot, tslot]) )
end if
if (param%ltides) then
call check( nf90_get_var(iu%ncid, iu%k2_varid, self%k2, start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%Q_varid, self%Q, start=[idslot, tslot]) )
end if
class is (swiftest_cb)

end select
idslot = 1
call check( nf90_get_var(iu%ncid, iu%id_varid, self%id, start=[idslot]) )

call check( nf90_get_var(iu%ncid, iu%Gmass_varid, self%Gmass, start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%radius_varid, self%radius, start=[idslot, tslot]) )
if (param%lrotation) then
call check( nf90_get_var(iu%ncid, iu%Ip1_varid, self%Ip(1), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%Ip2_varid, self%Ip(2), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%Ip3_varid, self%Ip(3), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%rotx_varid, self%rot(1), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%roty_varid, self%rot(2), start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%rotz_varid, self%rot(3), start=[idslot, tslot]) )
end if
if (param%ltides) then
call check( nf90_get_var(iu%ncid, iu%k2_varid, self%k2, start=[idslot, tslot]) )
call check( nf90_get_var(iu%ncid, iu%Q_varid, self%Q, start=[idslot, tslot]) )
end if


end select

!call self%read_particle_info(iu) ! THIS NEEDS TO BE IMPLEMENTED

return

end function netcdf_read_frame_base


module function netcdf_read_frame_system(self, iu, param) result(ierr)
!! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
!!
Expand All @@ -420,24 +443,12 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr)
call iu%open(param)

call self%read_hdr(iu, param)
ierr = self%cb%read_frame(iu, param)
if (ierr /= 0) then
write(*, *) "Cannot read central body frame."
goto 667
end if
ierr = self%pl%read_frame(iu, param)
if (ierr /= 0) then
write(*, *) "Cannot read massive body frame."
goto 667
end if
ierr = self%tp%read_frame(iu, param)
if (ierr /= 0) then
write(*, *) "Cannot read test particle frame."
goto 667
end if

call self%cb%read_in(param)
call self%pl%read_in(param)
call self%tp%read_in(param)
call iu%close(param)

ierr = 0
return

667 continue
Expand Down
11 changes: 1 addition & 10 deletions src/setup/setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,16 +167,7 @@ module subroutine setup_initialize_system(self, param)

associate(system => self, cb => self%cb, pl => self%pl, tp => self%tp)

if ((param%in_type == REAL8_TYPE) .or. (param%in_type == REAL4_TYPE)) then
!
else if ((param%in_type == NETCDF_DOUBLE_TYPE) .or. (param%in_type == NETCDF_FLOAT_TYPE)) then
!
else
call cb%read_in(param)
call pl%read_in(param)
call tp%read_in(param)
end if

call system%read_in(param)
call system%validate_ids(param)
call system%set_msys()
call pl%set_mu(cb)
Expand Down

0 comments on commit 195fd1a

Please sign in to comment.