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

Commit

Permalink
Fixed issues getting particle types when they aren't in the input file
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 4, 2022
1 parent a527bc3 commit 63877ff
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 27 deletions.
2 changes: 1 addition & 1 deletion examples/Basic_Simulation/initial_conditions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,4 +63,4 @@
sim.get_parameter()

# Run the simulation
#sim.run(tstart=0.0, tstop=1.0e2, dt=0.01, tstep_out=1.0e0, dump_cadence=0, )
sim.run(tstart=0.0, tstop=1.0e2, dt=0.01, tstep_out=1.0e0, dump_cadence=0, )
78 changes: 52 additions & 26 deletions src/netcdf/netcdf.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,7 +336,6 @@ module subroutine netcdf_open(self, param, readonly)

! Required Variables
call check( nf90_inq_varid(nciu%id, nciu%name_varname, nciu%name_varid), "netcdf_open nf90_inq_varid name_varid" )
call check( nf90_inq_varid(nciu%id, nciu%ptype_varname, nciu%ptype_varid), "netcdf_open nf90_inq_varid ptype_varid" )
call check( nf90_inq_varid(nciu%id, nciu%gmass_varname, nciu%Gmass_varid), "netcdf_open nf90_inq_varid Gmass_varid" )

if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then
Expand Down Expand Up @@ -389,6 +388,7 @@ module subroutine netcdf_open(self, param, readonly)
status = nf90_inq_varid(nciu%id, nciu%status_varname, nciu%status_varid)
status = nf90_inq_varid(nciu%id, nciu%j2rp2_varname, nciu%j2rp2_varid)
status = nf90_inq_varid(nciu%id, nciu%j4rp4_varname, nciu%j4rp4_varid)
status = nf90_inq_varid(nciu%id, nciu%ptype_varname, nciu%ptype_varid)

if (param%integrator == SYMBA) then
status = nf90_inq_varid(nciu%id, nciu%nplm_varname, nciu%nplm_varid)
Expand Down Expand Up @@ -472,8 +472,8 @@ module function netcdf_read_frame_system(self, nciu, param) result(ierr)
end if

! Next, filter only bodies that don't have mass (test particles)
call check( nf90_get_var(nciu%id, nciu%Gmass_varid, vectemp(:,:), start=[1, 1, tslot]), "netcdf_read_frame_system nf90_getvar Gmass_varid" )
plmask(:) = vectemp(1,:) == vectemp(1,:) .and. validmask(:)
call check( nf90_get_var(nciu%id, nciu%Gmass_varid, rtemp(:), start=[1, tslot]), "netcdf_read_frame_system nf90_getvar tp finder Gmass_varid" )
plmask(:) = rtemp(:) == rtemp(:) .and. validmask(:)
tpmask(:) = .not. plmask(:) .and. validmask(:)
plmask(1) = .false. ! This is the central body

Expand Down Expand Up @@ -694,15 +694,15 @@ module subroutine netcdf_read_hdr_system(self, nciu, param)


tslot = param%ioutput + 1
call check( nf90_inquire_dimension(nciu%id, nciu%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid" )
call check( nf90_inquire_dimension(nciu%id, nciu%id_dimid, len=idmax), "netcdf_read_hdr_system nf90_inquire_dimension id_dimid" )
call check( nf90_get_var(nciu%id, nciu%time_varid, self%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" )

allocate(gmtemp(idmax))
allocate(tpmask(idmax))
allocate(plmask(idmax))
allocate(plmmask(idmax))

call check( nf90_get_var(nciu%id, nciu%Gmass_varid, gmtemp, start=[1,1], count=[idmax,1]), "netcdf_read_frame_system nf90_getvar Gmass_varid" )
call check( nf90_get_var(nciu%id, nciu%Gmass_varid, gmtemp, start=[1,1], count=[idmax,1]), "netcdf_read_hdr_system nf90_getvar Gmass_varid" )

plmask(:) = gmtemp(:) == gmtemp(:)
tpmask(:) = .not. plmask(:)
Expand Down Expand Up @@ -827,14 +827,40 @@ module subroutine netcdf_read_particle_info_system(self, nciu, param, plmask, tp
call tp%info(i)%set_value(name=ctemp(tpind(i)))
end do

call check( nf90_get_var(nciu%id, nciu%ptype_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_read_particle_info_system nf90_getvar ptype_varid" )
call cb%info%set_value(particle_type=ctemp(1))
do i = 1, npl
call pl%info(i)%set_value(particle_type=ctemp(plind(i)))
end do
do i = 1, ntp
call tp%info(i)%set_value(particle_type=ctemp(tpind(i)))
end do
status = nf90_get_var(nciu%id, nciu%ptype_varid, ctemp, count=[NAMELEN, idmax])
if (status /= nf90_noerr) then ! Set default particle types
call cb%info%set_value(particle_type=CB_TYPE_NAME)

! Handle semi-interacting bodies in SyMBA
select type(pl)
class is (symba_pl)
select type (param)
class is (symba_parameters)
do i = 1, npl
if (pl%Gmass(i) < param%GMTINY) then
call pl%info(i)%set_value(particle_type=PL_TINY_TYPE_NAME)
else
call pl%info(i)%set_value(particle_type=PL_TYPE_NAME)
end if
end do
end select
class default ! Non-SyMBA massive bodies
do i = 1, npl
call pl%info(i)%set_value(particle_type=PL_TYPE_NAME)
end do
end select
do i = 1, ntp
call tp%info(i)%set_value(particle_type=TP_TYPE_NAME)
end do
else ! Use particle types defined in input file
call cb%info%set_value(particle_type=ctemp(1))
do i = 1, npl
call pl%info(i)%set_value(particle_type=ctemp(plind(i)))
end do
do i = 1, ntp
call tp%info(i)%set_value(particle_type=ctemp(tpind(i)))
end do
end if

status = nf90_inq_varid(nciu%id, nciu%status_varname, nciu%status_varid)
if (status == nf90_noerr) then
Expand Down Expand Up @@ -1049,28 +1075,28 @@ module subroutine netcdf_write_frame_base(self, nciu, param)
self%vh(1,j), self%vh(2,j), self%vh(3,j), &
a, e, inc, capom, omega, capm)
end if
call check( nf90_put_var(nciu%id, nciu%a_varid, a, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var a_varid" )
call check( nf90_put_var(nciu%id, nciu%e_varid, e, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var e_varid" )
call check( nf90_put_var(nciu%id, nciu%inc_varid, inc * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var inc_varid" )
call check( nf90_put_var(nciu%id, nciu%capom_varid, capom * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var capom_varid" )
call check( nf90_put_var(nciu%id, nciu%omega_varid, omega * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var omega_varid" )
call check( nf90_put_var(nciu%id, nciu%capm_varid, capm * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var capm_varid" )
call check( nf90_put_var(nciu%id, nciu%a_varid, a, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body a_varid" )
call check( nf90_put_var(nciu%id, nciu%e_varid, e, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body e_varid" )
call check( nf90_put_var(nciu%id, nciu%inc_varid, inc * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body inc_varid" )
call check( nf90_put_var(nciu%id, nciu%capom_varid, capom * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capom_varid" )
call check( nf90_put_var(nciu%id, nciu%omega_varid, omega * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body omega_varid" )
call check( nf90_put_var(nciu%id, nciu%capm_varid, capm * RAD2DEG, start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body capm_varid" )
end if

select type(self)
class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body
call check( nf90_put_var(nciu%id, nciu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var Gmass_varid" )
call check( nf90_put_var(nciu%id, nciu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body Gmass_varid" )
if (param%lrhill_present) then
call check( nf90_put_var(nciu%id, nciu%rhill_varid, self%rhill(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var rhill_varid" )
call check( nf90_put_var(nciu%id, nciu%rhill_varid, self%rhill(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body rhill_varid" )
end if
if (param%lclose) call check( nf90_put_var(nciu%id, nciu%radius_varid, self%radius(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var radius_varid" )
if (param%lclose) call check( nf90_put_var(nciu%id, nciu%radius_varid, self%radius(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body radius_varid" )
if (param%lrotation) then
call check( nf90_put_var(nciu%id, nciu%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var Ip_varid" )
call check( nf90_put_var(nciu%id, nciu%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var rotx_varid" )
call check( nf90_put_var(nciu%id, nciu%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var body Ip_varid" )
call check( nf90_put_var(nciu%id, nciu%rot_varid, self%rot(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]), "netcdf_write_frame_base nf90_put_var body rotx_varid" )
end if
! if (param%ltides) then
! call check( nf90_put_var(nciu%id, nciu%k2_varid, self%k2(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var k2_varid" )
! call check( nf90_put_var(nciu%id, nciu%Q_varid, self%Q(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var Q_varid" )
! call check( nf90_put_var(nciu%id, nciu%k2_varid, self%k2(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body k2_varid" )
! call check( nf90_put_var(nciu%id, nciu%Q_varid, self%Q(j), start=[idslot, tslot]), "netcdf_write_frame_base nf90_put_var body Q_varid" )
! end if

end select
Expand Down

0 comments on commit 63877ff

Please sign in to comment.