From d371a726eecc66218b07d296d373d983f5b24e4b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 21 Jun 2021 13:24:31 -0400 Subject: [PATCH] Fixed numerous issues that were prevening SyMBA runs from being properly restarted such that the produce identical results to non-restarted runs. --- .../original/clement_2018.in | 2 +- examples/symba_test_restart/original/param.in | 4 +- src/io/io_dump_pl_symba.f90 | 51 +++++++ src/main/swiftest_symba.f90 | 11 +- src/modules/module_interfaces.f90 | 10 ++ src/modules/module_symba.f90 | 7 + src/symba/symba_helio_getacch.f90 | 10 +- src/symba/symba_read_pl_in.f90 | 136 ++++++++++++++++++ 8 files changed, 216 insertions(+), 15 deletions(-) create mode 100644 src/io/io_dump_pl_symba.f90 create mode 100644 src/symba/symba_read_pl_in.f90 diff --git a/examples/symba_test_restart/original/clement_2018.in b/examples/symba_test_restart/original/clement_2018.in index 06e5334e2..7db3a9473 100644 --- a/examples/symba_test_restart/original/clement_2018.in +++ b/examples/symba_test_restart/original/clement_2018.in @@ -1,4 +1,4 @@ -1106 ! Solar System in unit system AU, M_sun, and years +1106 ! Solar System in unit system AU, M_sun, and years 1 3.94784176e+01 .0 .0 .0 ! x y z .0 .0 .0 ! vx vy vz diff --git a/examples/symba_test_restart/original/param.in b/examples/symba_test_restart/original/param.in index 046d19cda..68eb4df7d 100644 --- a/examples/symba_test_restart/original/param.in +++ b/examples/symba_test_restart/original/param.in @@ -7,8 +7,8 @@ DT 0.016 ! stepsize in years PL_IN clement_2018.in TP_IN tp.in IN_TYPE ASCII -ISTEP_OUT 625 ! output cadence -ISTEP_DUMP 625 ! system dump cadence +ISTEP_OUT 625 ! output cadence +ISTEP_DUMP 625 ! system dump cadence BIN_OUT bin.dat PARTICLE_FILE particle.dat OUT_TYPE REAL8 ! double precision real output diff --git a/src/io/io_dump_pl_symba.f90 b/src/io/io_dump_pl_symba.f90 new file mode 100644 index 000000000..7a16117ca --- /dev/null +++ b/src/io/io_dump_pl_symba.f90 @@ -0,0 +1,51 @@ +subroutine io_dump_pl_symba(npl, symba_plA, param) + !! Author: David A. Minton + !! + !! Dumps a SyMBA data structure to file. + !! +! modules + use swiftest + use module_symba + use module_interfaces, except_this_one => io_dump_pl_symba + implicit none + +! arguments + integer(I4B), intent(in) :: npl + type(symba_pl), intent(inout) :: symba_plA + type(user_input_parameters),intent(inout) :: param + +! internals + integer(I4B) :: ierr + integer(I4B), save :: idx = 1 + integer(I4B),parameter :: lun = 7 + + open(unit = lun, file = dump_pl_file(idx), form = "UNFORMATTED", status = 'REPLACE', iostat = ierr) + if (ierr /= 0) then + write(*, *) "Swiftest error:" + write(*, *) " Unable to open binary dump file ", trim(dump_pl_file(idx)) + call util_exit(FAILURE) + end if + write(lun) npl + write(lun) symba_plA%helio%swiftest%id(1:npl) + write(lun) symba_plA%helio%swiftest%mass(1:npl) + if (param%lrhill_present) write(lun) symba_plA%helio%swiftest%rhill(1:npl) + if (param%lclose) write(lun) symba_plA%helio%swiftest%radius(1:npl) + write(lun) symba_plA%helio%swiftest%xh(:,1:npl) + write(lun) symba_plA%helio%swiftest%vh(:,1:npl) + if (param%lrotation) then + write(lun) symba_plA%helio%swiftest%ip(:,1:npl) + write(lun) symba_plA%helio%swiftest%rot(:,1:npl) + end if + write(lun) symba_plA%helio%swiftest%vb(:,1:npl) + write(lun) symba_plA%helio%ah(:,1:npl) + write(lun) symba_plA%helio%ahi(:,1:npl) + close(lun) + if (idx == 1) then + idx = 2 + else + idx = 1 + end if + + return + +end subroutine io_dump_pl_symba \ No newline at end of file diff --git a/src/main/swiftest_symba.f90 b/src/main/swiftest_symba.f90 index 209a5f10f..848620106 100644 --- a/src/main/swiftest_symba.f90 +++ b/src/main/swiftest_symba.f90 @@ -93,7 +93,7 @@ program swiftest_symba end if ! reads in initial conditions of all massive bodies from input file - call symba_plA%helio%swiftest%read_from_file(param) + call symba_read_pl_in(symba_plA, param) call symba_tpA%helio%swiftest%read_from_file(param) !Temporary until the argument lists get fixed @@ -106,7 +106,7 @@ program swiftest_symba call symba_tp_allocate(symba_tpA,ntpfile) call symba_plA%helio%swiftest%dealloc() call symba_tpA%helio%swiftest%dealloc() - call symba_plA%helio%swiftest%read_from_file(param) + call symba_read_pl_in(symba_plA, param) call symba_tpA%helio%swiftest%read_from_file(param) !************************************************** call io_write_particle_pl(symba_plA%helio%swiftest, [(i, i=1,npl)], param) @@ -116,9 +116,9 @@ program swiftest_symba end if ! reorder by mass - if (out_stat /= "APPEND") call symba_reorder_pl(npl, symba_plA) + if (out_stat /= "APPEND") call symba_reorder_pl(npl, symba_plA) ! This is a new run, so we will sort the massive body list by mass + call util_valid(npl, ntp, symba_plA%helio%swiftest, symba_tpA%helio%swiftest) - call util_hills(npl, symba_plA%helio%swiftest) ntp0 = ntp t = t0 tbase = t0 @@ -155,6 +155,7 @@ program swiftest_symba start = clock_count / (count_rate * 1.0_DP) finish = start do while ((t < tstop) .and. ((ntp0 == 0) .or. (ntp > 0))) + call util_hills(npl, symba_plA%helio%swiftest) call symba_step_eucl(t, dt, param,npl,ntp,symba_plA, symba_tpA, nplplenc, npltpenc,& plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list) @@ -223,7 +224,7 @@ program swiftest_symba write(*,walltimefmt) finish - start, wallperstep call param%dump_to_file(t) - call io_dump_pl(npl, symba_plA%helio%swiftest, param) + call io_dump_pl_symba(npl, symba_plA, param) call io_dump_tp(ntp, symba_tpA%helio%swiftest) idump = istep_dump end if diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index de9ffac8d..7b39d79b7 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -497,6 +497,16 @@ SUBROUTINE io_dump_pl(npl, swiftest_plA, param) type(user_input_parameters),intent(inout) :: param END SUBROUTINE io_dump_pl + + subroutine io_dump_pl_symba(npl, symba_plA, param) + use swiftest_globals + use module_symba + use swiftest_data_structures + implicit none + integer(I4B), intent(in) :: npl + type(symba_pl), intent(inout):: symba_plA + type(user_input_parameters),intent(inout) :: param + end subroutine io_dump_pl_symba END INTERFACE INTERFACE diff --git a/src/modules/module_symba.f90 b/src/modules/module_symba.f90 index a93651dc3..e4e548a3c 100644 --- a/src/modules/module_symba.f90 +++ b/src/modules/module_symba.f90 @@ -103,4 +103,11 @@ MODULE module_symba end type symba_merger + interface + module subroutine symba_read_pl_in(symba_plA, param) + type(symba_pl), intent(inout) :: symba_plA !! Swiftest data structure to store massive body initial conditions + type(user_input_parameters), intent(inout) :: param !! Input collection of user-defined parameters + end subroutine symba_read_pl_in + end interface + END MODULE module_symba \ No newline at end of file diff --git a/src/symba/symba_helio_getacch.f90 b/src/symba/symba_helio_getacch.f90 index 70ed04fe5..d3cae9a1f 100644 --- a/src/symba/symba_helio_getacch.f90 +++ b/src/symba/symba_helio_getacch.f90 @@ -54,10 +54,8 @@ SUBROUTINE symba_helio_getacch(lflag, lextra_force, t, npl, nplm, helio_plA, j2r ! Executable code IF (lflag) THEN - DO i = 2, npl - helio_plA%ahi(:,i) = (/ 0.0_DP, 0.0_DP, 0.0_DP /) - END DO - CALL symba_helio_getacch_int(npl, nplm, helio_plA) + helio_plA%ahi(:,2:npl) = 0.0_DP + CALL symba_helio_getacch_int(npl, nplm, helio_plA) END IF IF (j2rp2 /= 0.0_DP) THEN DO i = 2, npl @@ -70,9 +68,7 @@ SUBROUTINE symba_helio_getacch(lflag, lextra_force, t, npl, nplm, helio_plA, j2r helio_plA%ah(:,i) = helio_plA%ahi(:,i) + aobl(:, i) - aobl(:, 1) END DO ELSE - DO i = 2, npl - helio_plA%ah(:,i) = helio_plA%ahi(:,i) - END DO + helio_plA%ah(:,2:npl) = helio_plA%ahi(:,2:npl) END IF IF (lextra_force) CALL helio_user_getacch(t, npl, helio_plA) diff --git a/src/symba/symba_read_pl_in.f90 b/src/symba/symba_read_pl_in.f90 new file mode 100644 index 000000000..e9b57bea0 --- /dev/null +++ b/src/symba/symba_read_pl_in.f90 @@ -0,0 +1,136 @@ +submodule (module_symba) s_symba_read_pl_in +contains + module subroutine symba_read_pl_in(symba_plA, param) + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Read in massive body data + !! + !! Adapted from David E. Kaufmann's Swifter routine symba_init_pl.f90 + !! Adapted from Martin Duncan's Swift routine symba_init_pl.f + use swiftest + use module_interfaces + implicit none + ! Arguments + type(symba_pl), intent(inout) :: symba_plA !! Swiftest data structure to store massive body initial conditions + type(user_input_parameters), intent(inout) :: param !! Input collection of user-defined parameters + ! Internals + integer(I4B), parameter :: LUN = 7 !! Unit number of input file + integer(I4B) :: i, ierr, npl + logical :: is_ascii + + ierr = 0 + is_ascii = (param%in_type == 'ASCII') + if (is_ascii) then + open(unit = LUN, file = param%inplfile, status = 'old', form = 'formatted', iostat = ierr) + else + open(unit = LUN, file = param%inplfile, status = 'old', form = 'unformatted', iostat = ierr) + end if + if (ierr /= 0) then + write(*,*) 'Error opening massive body initial conditions file ',trim(adjustl(param%inplfile)) + return + end if + + if (is_ascii) then + read(LUN, *, iostat = ierr) npl + else + read(LUN, iostat = ierr) npl + end if + if (npl <= 0) return + call symba_plA%helio%swiftest%alloc(npl) + + if (is_ascii) then + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%id(1), symba_plA%helio%swiftest%mass(1) + symba_plA%helio%swiftest%rhill(1) = 0.0_DP + symba_plA%helio%swiftest%radius(1) = param%rmin + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%xh(:,1) + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%vh(:,1) + if (param%lrotation) THEN + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%Ip(:,1) + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%rot(:,1) + end if + if (ierr /= 0) then + write(*,*) 'Error reading central body values in ',trim(adjustl(param%inplfile)) + return + end if + do i = 1, NDIM + if ((symba_plA%helio%swiftest%xh(i,1) /= 0.0_DP) .or. (symba_plA%helio%swiftest%vh(i,1) /= 0.0_DP)) then + write(*, *) "Swiftest error:" + write(*, *) " Input must be in heliocentric coordinates." + write(*, *) " position/velocity components of body 1 are" + write(*, *) symba_plA%helio%swiftest%xh(:,1) + write(*, *) symba_plA%helio%swiftest%vh(:,1) + end if + end do + symba_plA%helio%swiftest%status(1) = ACTIVE + do i = 2, symba_plA%helio%swiftest%nbody + if (param%lrhill_present) then + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%id(i), symba_plA%helio%swiftest%mass(i), symba_plA%helio%swiftest%rhill(i) + else + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%id(i), symba_plA%helio%swiftest%mass(i) + symba_plA%helio%swiftest%rhill(i) = 0.0_DP + end if + if (ierr /= 0 ) exit + if (param%lclose) then + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%radius(i) + if (ierr /= 0 ) exit + else + symba_plA%helio%swiftest%radius(i) = 0.0_DP + end if + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%xh(:,i) + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%vh(:,i) + if (param%lrotation) THEN + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%Ip(:,i) + read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%rot(:,i) + end if + if (ierr /= 0 ) exit + symba_plA%helio%swiftest%status(i) = ACTIVE + end do + else + read(LUN, iostat = ierr) symba_plA%helio%swiftest%id(:) + read(LUN, iostat = ierr) symba_plA%helio%swiftest%mass(:) + if (param%lrhill_present) then + read(LUN, iostat = ierr) symba_plA%helio%swiftest%rhill(:) + else + symba_plA%helio%swiftest%rhill(:) = 0.0_DP + end if + if (param%lclose) then + read(LUN, iostat = ierr) symba_plA%helio%swiftest%radius(:) + else + symba_plA%helio%swiftest%radius(:) = 0.0_DP + end if + read(LUN, iostat = ierr) symba_plA%helio%swiftest%xh(:,:) + read(LUN, iostat = ierr) symba_plA%helio%swiftest%vh(:,:) + if (param%lrotation) THEN + read(LUN, iostat = ierr) symba_plA%helio%swiftest%Ip(:,:) + read(LUN, iostat = ierr) symba_plA%helio%swiftest%rot(:,:) + end if + if (param%out_stat == 'APPEND') then + read(lun, iostat = ierr) symba_plA%helio%swiftest%vb(:,:) + read(lun, iostat = ierr) symba_plA%helio%ah(:,:) + read(lun, iostat = ierr) symba_plA%helio%ahi(:,:) + end if + if (ierr /= 0) then + write(*,*) 'An error occurred reading in ',trim(adjustl(param%inplfile)) + call util_exit(FAILURE) + end if + symba_plA%helio%swiftest%status(:) = ACTIVE + end if + close(unit = LUN) + + symba_plA%helio%swiftest%info(1)%origin_type = "Central body" + do i = 2, npl + symba_plA%helio%swiftest%info(i)%origin_xh(:) = symba_plA%helio%swiftest%xh(:,i) + symba_plA%helio%swiftest%info(i)%origin_vh(:) = symba_plA%helio%swiftest%vh(:,i) + symba_plA%helio%swiftest%info(i)%origin_time = 0.0_DP + symba_plA%helio%swiftest%info(i)%origin_type = "Initial conditions" + end do + + ! Give massive bodies a positive id + symba_plA%helio%swiftest%id(:) = abs(symba_plA%helio%swiftest%id(:)) + symba_plA%helio%swiftest%maxid = maxval(symba_plA%helio%swiftest%id(:)) + + return + end subroutine symba_read_pl_in + +end submodule s_symba_read_pl_in +