From 632d6dc5a02eb8728b0fd9a1270f4d9ef23c464c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 14 Dec 2022 10:45:03 -0500 Subject: [PATCH] Fixed issues with time-keeping and indexing of storage frames --- src/encounter/encounter_io.f90 | 11 ++++++----- src/encounter/encounter_util.f90 | 3 ++- src/fraggle/fraggle_io.f90 | 11 ++++++----- src/io/io.f90 | 6 +++--- src/netcdf/netcdf.f90 | 26 ++++++++++++++++++++------ src/setup/setup.f90 | 2 +- src/symba/symba_collision.f90 | 4 ++-- src/symba/symba_util.f90 | 15 ++++++--------- 8 files changed, 46 insertions(+), 32 deletions(-) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 22827f611..60748c5b0 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -112,9 +112,10 @@ module subroutine encounter_io_initialize(self, param) integer(I4B) :: nvar, varid, vartype real(DP) :: dfill real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 logical :: fileExists character(len=STRMAX) :: errmsg - integer(I4B) :: ndims, i + integer(I4B) :: ndims associate(nc => self) dfill = ieee_value(dfill, IEEE_QUIET_NAN) @@ -167,13 +168,13 @@ module subroutine encounter_io_initialize(self, param) call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) end select end do diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 0da74cae7..9a1be2901 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -690,7 +690,8 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k) call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) snapshot%t = t + tperi - if (snapshot%t < maxval(pl_snap%info(:)%origin_time)) cycle + if ((snapshot%t < maxval(pl_snap%info(:)%origin_time)) .or. & + (snapshot%t > minval(pl_snap%info(:)%discard_time))) cycle ! Computer the center mass of the pair rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index f47a64047..6f75c7d01 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -26,9 +26,10 @@ module subroutine fraggle_io_initialize_output(self, param) integer(I4B) :: nvar, varid, vartype real(DP) :: dfill real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 logical :: fileExists character(len=STRMAX) :: errmsg - integer(I4B) :: i, ndims + integer(I4B) :: ndims select type(param) class is (symba_parameters) @@ -120,13 +121,13 @@ module subroutine fraggle_io_initialize_output(self, param) call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "fraggle_io_initialize nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "fraggle_io_initialize nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "fraggle_io_initialize nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "fraggle_io_initialize nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "fraggle_io_initialize nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "fraggle_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "fraggle_io_initialize nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "fraggle_io_initialize nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "fraggle_io_initialize nf90_def_var_fill NF90_CHAR" ) end select end do ! Take the file out of define mode diff --git a/src/io/io.f90 b/src/io/io.f90 index f159e6ac7..54f99e42b 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -294,11 +294,11 @@ module subroutine io_dump_storage(self, param) integer(I8B) :: iloop_start if (self%iframe == 0) return - iloop_start = param%iloop - int(param%istep_out * param%dump_cadence, kind=I8B) + 1 + iloop_start = max(param%iloop - int(param%istep_out * param%dump_cadence, kind=I8B) + 1_I8B,0_I8B) call self%make_index_map() - do i = 1, param%dump_cadence - param%ioutput = iloop_start + self%tmap(i) + do i = 1, self%iframe if (allocated(self%frame(i)%item)) then + param%ioutput = iloop_start + self%tmap(i) select type(system => self%frame(i)%item) class is (swiftest_nbody_system) call system%write_frame(param) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 588a138d6..4ce217124 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -157,6 +157,7 @@ module subroutine netcdf_initialize_output(self, param) integer(I4B) :: nvar, varid, vartype real(DP) :: dfill real(SP) :: sfill + integer(I4B), parameter :: NO_FILL = 0 logical :: fileExists character(len=STRMAX) :: errmsg integer(I4B) :: ndims @@ -281,16 +282,24 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "netcdf_initialize_output nf90_inquire_variable" ) select case(vartype) case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "netcdf_initialize_output nf90_def_var_fill NF90_INT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), "netcdf_initialize_output nf90_def_var_fill NF90_INT" ) case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), "netcdf_initialize_output nf90_def_var_fill NF90_FLOAT" ) case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), "netcdf_initialize_output nf90_def_var_fill NF90_DOUBLE" ) case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) + call check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), "netcdf_initialize_output nf90_def_var_fill NF90_CHAR" ) end select end do + ! Set special fill mode for discard time so that we can make use of it for non-discarded bodies. + select case (vartype) + case(NF90_FLOAT) + call check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_SP)), "netcdf_initialize_output nf90_def_var_fill discard_time NF90_FLOAT" ) + case(NF90_DOUBLE) + call check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_DP)), "netcdf_initialize_output nf90_def_var_fill discard_time NF90_DOUBLE" ) + end select + ! Take the file out of define mode call check( nf90_enddef(nc%id), "netcdf_initialize_output nf90_enddef" ) @@ -944,7 +953,7 @@ module subroutine netcdf_read_particle_info_system(self, nc, param, plmask, tpma if (status == nf90_noerr) then call check( nf90_get_var(nc%id, nc%collision_id_varid, itemp), "netcdf_read_particle_info_system nf90_getvar collision_id_varid" ) else - itemp = 0.0_DP + itemp = 0 end if do i = 1, npl @@ -958,7 +967,12 @@ module subroutine netcdf_read_particle_info_system(self, nc, param, plmask, tpma if (status == nf90_noerr) then call check( nf90_get_var(nc%id, nc%discard_time_varid, rtemp), "netcdf_read_particle_info_system nf90_getvar discard_time_varid" ) else - rtemp = 0.0_DP + select case (param%out_type) + case("NETCDF_FLOAT") + rtemp(:) = huge(0.0_SP) + case("NETCDF_DOUBLE") + rtemp(:) = huge(0.0_DP) + end select end if call cb%info%set_value(discard_time=rtemp(1)) diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index ef8558aef..57679f622 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -248,7 +248,7 @@ module subroutine setup_body(self, n, param) origin_time = -huge(1.0_DP), & origin_rh = [0.0_DP, 0.0_DP, 0.0_DP], & origin_vh = [0.0_DP, 0.0_DP, 0.0_DP], & - discard_time = -huge(1.0_DP), & + discard_time = huge(1.0_DP), & discard_rh = [0.0_DP, 0.0_DP, 0.0_DP], & discard_vh = [0.0_DP, 0.0_DP, 0.0_DP], & discard_body_id = -1 & diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 8d554f018..0353cc7d5 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -367,8 +367,8 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step pl%lcollision([i, j]) = .true. pl%status([i, j]) = COLLISION - call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i)) - call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j)) + call pl%info(i)%set_value(status="COLLISION") + call pl%info(j)%set_value(status="COLLISION") end if else self%r2(:,k) = tp%rh(:,j) + system%cb%rb(:) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 157c3f5af..cf817622e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -659,15 +659,12 @@ module subroutine symba_util_rearray_pl(self, system, param) pl%lcollision(1:npl) = .false. pl%lmask(1:npl) = .true. - select type(param) - class is (symba_parameters) - pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY - where(pl%lmtiny(1:npl)) - pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME - elsewhere - pl%info(1:npl)%particle_type = PL_TYPE_NAME - end where - end select + pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY + where(pl%lmtiny(1:npl)) + pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME + elsewhere + pl%info(1:npl)%particle_type = PL_TYPE_NAME + end where call pl%write_info(param%system_history%nc, param) deallocate(ldump_mask)