From 195fd1ac282dc7d2fbd6c83eec3b059435efa2cd Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 24 Sep 2021 17:17:54 -0400 Subject: [PATCH] Added methods and structures for reading in NetCDF files. It hasn't been tested, but it compiles --- src/io/io.f90 | 68 ++++++++++++++-- src/modules/swiftest_classes.f90 | 23 +++--- src/netcdf/netcdf.f90 | 131 +++++++++++++++++-------------- src/setup/setup.f90 | 11 +-- 4 files changed, 147 insertions(+), 86 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 89b0a5988..95b5abb7a 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -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 @@ -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') @@ -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 @@ -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 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 2dbcdd38a..5fe807071 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index d97f9169e..3ad813659 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -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 !! @@ -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]) ) @@ -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 !! @@ -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 diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 2b02fce28..208183e55 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -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)