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

Commit

Permalink
Testing improved restart capability by reading in more state variable…
Browse files Browse the repository at this point in the history 
…s from the paramter input file
  • Loading branch information
@daminton
daminton committed Jun 19, 2021
1 parent 9a68e0b commit 6f8359b
Show file tree
Hide file tree
Showing 12 changed files with 466 additions and 61 deletions.
4 changes: 2 additions & 2 deletions Makefile.Defines
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vec -vecabi=cmdtarget -simd -shared-intel
VTUNE_FLAGS = -g -O2 -vec -simd -shared-intel -qopenmp -debug inline-debug-info -parallel-source-info=2 -parallel -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -qopenmp -fp-model no-except -mp1 -xhost -traceback
#Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune

IDEBUG = -O0 -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin
IDEBUG = -O0 -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -init=snan,arrays

STRICTREAL = -mp1 -fp-model strict -prec-div -prec-sqrt -assume protect-parens
SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -prec-div -prec-sqrt -assume protect-parens
Expand All @@ -70,7 +70,7 @@ GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporari

#FFLAGS = -init=snan,arrays -traceback -no-wrap-margin -O3 -g -CB -nogen-interfaces -no-pie -fp-speculation=safe $(SIMDVEC) $(PAR) #$(HEAPARR)
#FFLAGS = $(IDEBUG) $(HEAPARR)
FFLAGS = -init=snan,arrays -no-wrap-margin -fp-model strict -fp-model no-except -traceback -g $(OPTIMIZE) -O3 $(SIMDVEC) $(PAR) $(HEAPARR)
FFLAGS = -init=snan,arrays -no-wrap-margin -fp-model strict -fp-model no-except -traceback -g -O3 $(SIMDVEC) $(PAR) $(HEAPARR)
FORTRAN = ifort
#AR = xiar

Expand Down
249 changes: 249 additions & 0 deletions ...son/9pl_18tp_encounters/.ipynb_checkpoints/swiftest_rmvs_vs_swifter_rmvs-checkpoint.ipynb
View file

Large diffs are not rendered by default.

35 changes: 35 additions & 0 deletions examples/symba_mars_disk/param.full.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
!Parameter file for the SyMBA-RINGMOONS test
T0 0.0
TSTOP 2400.0
DT 600.0
PL_IN mars.in
TP_IN tp.in
IN_TYPE ASCII
ISTEP_OUT 1
ISTEP_DUMP 1
BIN_OUT bin.dat
PARTICLE_FILE particle.dat
OUT_TYPE REAL8
OUT_FORM EL
OUT_STAT REPLACE
J2 0.0
J4 0.0
CHK_CLOSE yes
CHK_RMIN 3389500.0
CHK_RMAX 3389500000.0
CHK_EJECT 3389500000.0
CHK_QMIN 3389500.0
CHK_QMIN_COORD HELIO
CHK_QMIN_RANGE 3389500.0 338950000000.0
ENC_OUT /dev/null
EXTRA_FORCE no
BIG_DISCARD no
RHILL_PRESENT yes
MTINY 1000.0
ENERGY yes
FRAGMENTATION yes
ROTATION yes
MU2KG 1.0
DU2M 1.0
TU2S 1.0
SEED 2 3080983 2220830
2 changes: 1 addition & 1 deletion examples/symba_mars_disk/param.in
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
!Parameter file for the SyMBA-RINGMOONS test
T0 0.0
TSTOP 6e8
TSTOP 6000.0
DT 600.0
PL_IN mars.in
TP_IN tp.in
Expand Down
35 changes: 35 additions & 0 deletions examples/symba_mars_disk/param.partial.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
!Parameter file for the SyMBA-RINGMOONS test
T0 0.0
TSTOP 1200.0
DT 600.0
PL_IN mars.in
TP_IN tp.in
IN_TYPE ASCII
ISTEP_OUT 1
ISTEP_DUMP 1
BIN_OUT bin.dat
PARTICLE_FILE particle.dat
OUT_TYPE REAL8
OUT_FORM EL
OUT_STAT REPLACE
J2 0.0
J4 0.0
CHK_CLOSE yes
CHK_RMIN 3389500.0
CHK_RMAX 3389500000.0
CHK_EJECT 3389500000.0
CHK_QMIN 3389500.0
CHK_QMIN_COORD HELIO
CHK_QMIN_RANGE 3389500.0 338950000000.0
ENC_OUT /dev/null
EXTRA_FORCE no
BIG_DISCARD no
RHILL_PRESENT yes
MTINY 1000.0
ENERGY yes
FRAGMENTATION yes
ROTATION yes
MU2KG 1.0
DU2M 1.0
TU2S 1.0
SEED 2 3080983 2220830
37 changes: 37 additions & 0 deletions examples/symba_mars_disk/param.restart.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
TSTOP 2.40000000000000000E+03
T0 1.20000000000000000E+03
DT 6.00000000000000000E+02
PL_IN dump_pl2.bin
TP_in dump_tp2.bin
IN_TYPE REAL8
ISTEP_OUT 1
BIN_OUT bin.dat
PARTICLE_FILE particle.dat
OUT_TYPE REAL8
OUT_FORM EL
OUT_STAT APPEND
ENC_OUT /dev/null
ISTEP_DUMP 1
CHK_RMIN 3.38950000000000000E+06
CHK_RMAX 3.38950000000000000E+09
CHK_EJECT 3.38950000000000000E+09
CHK_QMIN 3.38950000000000000E+06
CHK_QMIN_COORD HELIO
CHK_QMIN_RANGE 3.38950000000000000E+06 3.38950000000000000E+11
MTINY 1.00000000000000000E+03
MU2KG 1.00000000000000000E+00
TU2S 1.00000000000000000E+00
DU2M 1.00000000000000000E+00
SEED 2 3080983 2220830
EXTRA_FORCE F
BIG_DISCARD F
RHILL_PRESENT T
CHK_CLOSE T
FRAGMENTATION T
ENERGY T
ROTATION T
TIDES F
GR F
YARKOVSKY F
YORP F
RINGMOONS F
9 changes: 4 additions & 5 deletions src/io/io_conservation_report.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,19 +12,16 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
type(symba_pl), intent(inout) :: symba_plA !! Swiftest planet data structure
integer(I4B), intent(in) :: npl !! Number of massive bodies
real(DP), intent(in) :: j2rp2, j4rp4 !! Central body oblateness terms
type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters
type(user_input_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters
logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen
! Internals
real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now
real(DP), dimension(NDIM), save :: Ltot_orig, Lorbit_orig, Lspin_orig
real(DP), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last
real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig
real(DP), save :: ke_orbit_last, ke_spin_last, pe_last, Eorbit_last
real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now
real(DP) :: Eorbit_error, Etotal_error, Ecoll_error
real(DP) :: Mtot_now, Merror
real(DP) :: Lmag_now, Lerror
logical, save :: lfirst = .true.
character(len=*), parameter :: egyfmt = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output
character(len=*), parameter :: egyheader = '("t,Eorbit,Ecollisions,Lx,Ly,Lz,Mtot")'
integer(I4B), parameter :: egyiu = 72
Expand All @@ -35,7 +32,9 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,

associate(Ecollisions => symba_plA%helio%swiftest%Ecollisions, Lescape => symba_plA%helio%swiftest%Lescape, Mescape => symba_plA%helio%swiftest%Mescape, &
Euntracked => symba_plA%helio%swiftest%Euntracked, mass => symba_plA%helio%swiftest%mass, dMcb => symba_plA%helio%swiftest%dMcb, &
Mcb_initial => symba_plA%helio%swiftest%Mcb_initial)
Mcb_initial => symba_plA%helio%swiftest%Mcb_initial, Eorbit_orig => param%Eorbit_orig, Mtot_orig => param%Mtot_orig , &
Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), &
lfirst => param%lfirstenergy)
if (lfirst) then
if (param%out_stat == "OLD") then
open(unit = egyiu, file = energy_file, form = "formatted", status = "old", action = "write", position = "append")
Expand Down
2 changes: 1 addition & 1 deletion src/modules/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
type(symba_pl), intent(inout) :: symba_plA !! Swiftest planet data structure
integer(I4B), intent(in) :: npl !! Number of massive bodies
real(DP), intent(in) :: j2rp2, j4rp4 !! Central body oblateness terms
type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters
type(user_input_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters
logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen
end subroutine io_conservation_report

Expand Down
8 changes: 8 additions & 0 deletions src/modules/user.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,14 @@ module user
logical :: ltides = .false. !! Include tidal dissipation
logical :: lringmoons = .false. !! Turn on the ringmoons code
logical :: lenergy = .false. !! Track the total energy of the system
logical :: lfirstenergy = .true.
real(DP) :: Eorbit_orig = 0.0_DP
real(DP) :: Mtot_orig = 0.0_DP
real(DP) :: Lmag_orig = 0.0_DP
real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP
real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP
real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP
logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step

! Future features not implemented or in development
logical :: lgr = .false. !! Turn on GR
Expand Down
103 changes: 52 additions & 51 deletions src/symba/symba_step_eucl.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,73 +87,74 @@ SUBROUTINE symba_step_eucl(t,dt,param,npl, ntp,symba_plA, symba_tpA, &
INTEGER(I4B) :: irec, nplm, ipl, ipl1, ipl2
INTEGER(I8B) :: i, k, counter
INTEGER(I8B), DIMENSION(:), ALLOCATABLE :: pltp_encounters_indices
LOGICAL, SAVE :: lfirst = .true.

LOGICAL(LGT), ALLOCATABLE, DIMENSION(:) :: pltp_lencounters
LOGICAL(lgt), ALLOCATABLE, DIMENSION(:) :: plpl_lvdotr, pltp_lvdotr

! Executable code
call symba_step_reset(npl, symba_plA, symba_tpA, plplenc_list, pltpenc_list, mergeadd_list, mergesub_list)
nplplenc = 0
npltpenc = 0
irec = 0
nplm = count(symba_plA%helio%swiftest%mass(1:npl) > param%mtiny)
associate(lfirst => param%lfirstkick)
call symba_step_reset(npl, symba_plA, symba_tpA, plplenc_list, pltpenc_list, mergeadd_list, mergesub_list)
nplplenc = 0
npltpenc = 0
irec = 0
nplm = count(symba_plA%helio%swiftest%mass(1:npl) > param%mtiny)

CALL symba_chk_eucl(npl, irec, symba_plA, dt, plplenc_list, nplplenc)

if(ntp>0)then
allocate(pltp_lencounters(symba_tpA%helio%swiftest%num_pltp_comparisons))
allocate(pltp_lvdotr(symba_tpA%helio%swiftest%num_pltp_comparisons))
pltp_lencounters = .false.
pltp_lvdotr = .false.
CALL symba_chk_eucl(npl, irec, symba_plA, dt, plplenc_list, nplplenc)
if(ntp>0)then
allocate(pltp_lencounters(symba_tpA%helio%swiftest%num_pltp_comparisons))
allocate(pltp_lvdotr(symba_tpA%helio%swiftest%num_pltp_comparisons))
pltp_lencounters = .false.
pltp_lvdotr = .false.

CALL symba_chk_eucl_pltp(symba_plA, symba_tpA, dt, pltp_lencounters, pltp_lvdotr, npltpenc)

! npltpenc = count(pltp_encounters > 0)
! print *,'step npltpenc: ',npltpenc
if(npltpenc>0)then
CALL symba_chk_eucl_pltp(symba_plA, symba_tpA, dt, pltp_lencounters, pltp_lvdotr, npltpenc)
! npltpenc = count(pltp_encounters > 0)
! print *,'step npltpenc: ',npltpenc
if(npltpenc>0)then

allocate(pltp_encounters_indices(npltpenc))
allocate(pltp_encounters_indices(npltpenc))

counter = 1
do k = 1,symba_tpA%helio%swiftest%num_pltp_comparisons
if(pltp_lencounters(k))then
pltp_encounters_indices(counter) = k
counter = counter + 1
endif
enddo
counter = 1
do k = 1,symba_tpA%helio%swiftest%num_pltp_comparisons
if(pltp_lencounters(k))then
pltp_encounters_indices(counter) = k
counter = counter + 1
endif
enddo

symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) = symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) + 1
symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) = symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) + 1
symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) = symba_plA%ntpenc(symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))) + 1
symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) = symba_tpA%nplenc(symba_tpA%helio%swiftest%k_pltp(2,pltp_encounters_indices(:))) + 1

pltpenc_list%status(1:npltpenc) = ACTIVE
pltpenc_list%lvdotr(1:npltpenc) = pltp_lvdotr(pltp_encounters_indices(:))
pltpenc_list%level(1:npltpenc) = 0
pltpenc_list%indexpl(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))
pltpenc_list%indextp(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))
pltpenc_list%status(1:npltpenc) = ACTIVE
pltpenc_list%lvdotr(1:npltpenc) = pltp_lvdotr(pltp_encounters_indices(:))
pltpenc_list%level(1:npltpenc) = 0
pltpenc_list%indexpl(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))
pltpenc_list%indextp(1:npltpenc) = symba_tpA%helio%swiftest%k_pltp(1,pltp_encounters_indices(:))

deallocate(pltp_encounters_indices)
endif
deallocate(pltp_encounters_indices)
endif

deallocate(pltp_lencounters, pltp_lvdotr)
endif

! END OF THINGS THAT NEED TO BE CHANGED IN THE TREE
deallocate(pltp_lencounters, pltp_lvdotr)
endif
! END OF THINGS THAT NEED TO BE CHANGED IN THE TREE

! flag to see if there was an encounter
lencounter = ((nplplenc > 0) .OR. (npltpenc > 0))
! flag to see if there was an encounter
lencounter = ((nplplenc > 0) .OR. (npltpenc > 0))

IF (lencounter) THEN ! if there was an encounter, we need to enter symba_step_interp to see if we need recursion
CALL symba_step_interp_eucl(t, npl, nplm, ntp, symba_plA, symba_tpA, dt, nplplenc, npltpenc, &
plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, param)
lfirst = .TRUE.
ELSE ! otherwise we can just advance the particles
CALL symba_step_helio(lfirst, param%lextra_force, t, npl, nplm, ntp,&
symba_plA%helio, symba_tpA%helio, &
param%j2rp2, param%j4rp4, dt)
END IF
IF (lencounter) THEN ! if there was an encounter, we need to enter symba_step_interp to see if we need recursion
CALL symba_step_interp_eucl(t, npl, nplm, ntp, symba_plA, symba_tpA, dt, nplplenc, npltpenc, &
plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list, param)
lfirst = .TRUE.
ELSE ! otherwise we can just advance the particles
CALL symba_step_helio(lfirst, param%lextra_force, t, npl, nplm, ntp,&
symba_plA%helio, symba_tpA%helio, &
param%j2rp2, param%j4rp4, dt)
END IF

RETURN
RETURN
end associate

END SUBROUTINE symba_step_eucl
!**********************************************************************************************************************************
Expand Down
31 changes: 31 additions & 0 deletions src/user/user_udio_reader.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,6 +179,37 @@ module subroutine user_udio_reader(param, unit, iotype, v_list, iostat, iomsg)
param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, i=nseeds_from_file+1, nseeds)]
end if
seed_set = .true.
case ("FIRSTKICK")
call util_toupper(param_value)
if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false.
case ("FIRSTENERGY")
call util_toupper(param_value)
if (param_value == "NO" .or. param_value == 'F') param%lfirstenergy = .false.
case("EORBIT_ORIG")
read(param_value, *) param%Eorbit_orig
case("MTOT_ORIG")
read(param_value, *) param%Mtot_orig
case("LTOT_ORIG")
read(param_value, *) param%Ltot_orig(1)
do i = 2, NDIM
ifirst = ilast + 1
param_value = user_get_token(line, ifirst, ilast, iostat)
read(param_value, *) param%Ltot_orig(i)
end do
case("LORBIT_ORIG")
read(param_value, *) param%Lorbit_orig(1)
do i = 2, NDIM
ifirst = ilast + 1
param_value = user_get_token(line, ifirst, ilast, iostat)
read(param_value, *) param%Lorbit_orig(i)
end do
case("LSPIN_ORIG")
read(param_value, *) param%Lspin_orig(1)
do i = 2, NDIM
ifirst = ilast + 1
param_value = user_get_token(line, ifirst, ilast, iostat)
read(param_value, *) param%Lspin_orig(i)
end do
case default
write(iomsg,*) "Unknown parameter -> ",param_name
iostat = -1
Expand Down
Loading

0 comments on commit 6f8359b

Please sign in to comment.