From 3b316630476c8a9c51fa0ed5150c2508cd0e9c79 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 12:10:31 -0400 Subject: [PATCH 01/15] Restructured parameter file read/write subroutines to simplify the code. --- src/io/io.f90 | 387 +++++++++++++++++++++---------- src/modules/swiftest_classes.f90 | 64 ++++- src/modules/walltime_classes.f90 | 11 +- src/symba/symba_io.f90 | 50 +--- 4 files changed, 347 insertions(+), 165 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 51d10561d..ef3f682d2 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -479,6 +479,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! !! Adapted from David E. Kaufmann's Swifter routine io_init_param.f90 !! Adapted from Martin Duncan's Swift routine io_init_param.f + use, intrinsic :: iso_fortran_env implicit none ! Arguments class(swiftest_parameters), intent(inout) :: self !! Collection of parameters @@ -500,6 +501,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ! Parse the file line by line, extracting tokens then matching them up with known parameters if possible associate(param => self) + rewind(unit) do read(unit = unit, fmt = linefmt, end = 1, err = 667, iomsg = iomsg) line line_trim = trim(adjustl(line)) @@ -601,9 +603,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("TIDES") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') param%ltides = .true. - case ("FLATTEN_INTERACTIONS") + case ("INTERACTION_LOOPS") call io_toupper(param_value) - if (param_value == "NO" .or. param_value == 'F') param%lflatten_interactions = .false. + param%interaction_loops = param_value case ("FIRSTKICK") call io_toupper(param_value) if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. @@ -654,16 +656,14 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param%particle_out = param_value case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters case default - write(*,*) "Unknown parameter -> ",param_name - iostat = -1 - return + write(*,*) "Ignoring unknown parameter -> ",param_name end select end if end do 1 continue iostat = 0 - !! Do basic sanity checks on the input values + ! Do basic sanity checks on the input values if ((.not. t0_set) .or. (.not. tstop_set) .or. (.not. dt_set)) then write(iomsg,*) 'Valid simulation time not set' iostat = -1 @@ -731,58 +731,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) return end if - write(*,*) "T0 = ",param%t0 - write(*,*) "TSTOP = ",param%tstop - write(*,*) "DT = ",param%dt - write(*,*) "CB_IN = ",trim(adjustl(param%incbfile)) - write(*,*) "PL_IN = ",trim(adjustl(param%inplfile)) - write(*,*) "TP_IN = ",trim(adjustl(param%intpfile)) - write(*,*) "IN_TYPE = ",trim(adjustl(param%in_type)) - write(*,*) "ISTEP_OUT = ",param%istep_out - write(*,*) "BIN_OUT = ",trim(adjustl(param%outfile)) - write(*,*) "OUT_TYPE = ",trim(adjustl(param%out_type)) - write(*,*) "OUT_FORM = ",trim(adjustl(param%out_form)) - write(*,*) "IN_FORM = ",trim(adjustl(param%in_form)) - write(*,*) "OUT_STAT = ",trim(adjustl(param%out_stat)) - write(*,*) "ISTEP_DUMP = ",param%istep_dump - write(*,*) "CHK_CLOSE = ",param%lclose - if (param%rmin > 0.0) then - write(*,*) "CHK_RMIN = ",param%rmin - else - write(*,*) "! CHK_RMIN value will be the central body radius" - end if - if (param%rmax > 0.0_DP) write(*,*) "CHK_RMAX = ",param%rmax - if (param%rmaxu > 0.0_DP) write(*,*) "CHK_EJECT = ",param%rmaxu - if ((param%qmin > 0.0_DP) .or. (param%qmin_alo > 0.0_DP) .or. (param%qmin_ahi > 0.0_DP)) write(*,*) "CHK_QMIN_COORD = ",trim(adjustl(param%qmin_coord)) - if (param%qmin > 0.0_DP) write(*,*) "CHK_QMIN = ",param%qmin - if ((param%qmin_alo > 0.0_DP) .or. (param%qmin_ahi > 0.0_DP)) write(*,*) "CHK_QMIN_RANGE = ",param%qmin_alo, param%qmin_ahi - write(*,*) "EXTRA_FORCE = ",param%lextra_force - write(*,*) "RHILL_PRESENT = ",param%lrhill_present - write(*,*) "ROTATION = ", param%lrotation - write(*,*) "TIDES = ", param%ltides - write(*,*) "ENERGY = ",param%lenergy - write(*,*) "FLATTEN_INTERACTIONS = ",param%lflatten_interactions - if (param%lenergy) write(*,*) "ENERGY_OUT = ",trim(adjustl(param%energy_out)) - write(*,*) "MU2KG = ",param%MU2KG - write(*,*) "TU2S = ",param%TU2S - write(*,*) "DU2M = ",param%DU2M - if (trim(adjustl(param%outfile)) == "") then - param%outfile = BIN_OUTFILE - end if - if (trim(adjustl(param%enc_out)) /= "") then - write(*,*) "ENC_OUT = ",trim(adjustl(param%enc_out)) - else - write(*,*) "! ENC_OUT not set: Encounters will not be recorded to file" - end if - if (trim(adjustl(param%discard_out)) /= "") then - write(*,*) "DISCARD_OUT = ",trim(adjustl(param%discard_out)) - write(*,*) "BIG_DISCARD = ",param%lbig_discard - else - write(*,*) "! DISCARD_OUT not set: Discards will not be recorded to file" - param%lbig_discard = .false. - write(*,*) "! BIG_DISCARD = ",param%lbig_discard - end if - if ((param%MU2KG < 0.0_DP) .or. (param%TU2S < 0.0_DP) .or. (param%DU2M < 0.0_DP)) then write(iomsg,*) 'Invalid unit conversion factor' iostat = -1 @@ -791,11 +739,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ! Calculate the G for the system units param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) - if (param%lgr) then - ! Calculate the inverse speed of light in the system units - param%inv_c2 = einsteinC * param%TU2S / param%DU2M - param%inv_c2 = (param%inv_c2)**(-2) - end if associate(integrator => v_list(1)) if ((integrator == RMVS) .or. (integrator == SYMBA)) then @@ -809,14 +752,33 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ! Determine if the GR flag is set correctly for this integrator select case(integrator) case(WHM, RMVS, HELIO, SYMBA) - write(*,*) "GR = ", param%lgr case default if (param%lgr) write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.' param%lgr = .false. end select + + if (param%lgr) then + ! Calculate the inverse speed of light in the system units + param%inv_c2 = einsteinC * param%TU2S / param%DU2M + param%inv_c2 = (param%inv_c2)**(-2) + end if + end associate + select case(trim(adjustl(param%interaction_loops))) + case("ADAPTIVE", "FLAT", "TRIANGULAR") + case default + write(*,*) "Unknown value for parameter INTERACTION_LOOPS: -> ",trim(adjustl(param%interaction_loops)) + write(*,*) "Must be one of the following: TRIANGULAR, FLAT, or ADAPTIVE" + write(*,*) "Using default value of ADAPTIVE" + param%interaction_loops = "ADAPTIVE" + end select + iostat = 0 + + ! Print the contents of the parameter file to standard output + call param%writer(unit = OUTPUT_UNIT, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg) + end associate 667 continue @@ -853,88 +815,271 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) integer(I4B) :: i associate(param => self) - write(param_name, *) "T0"; write(param_value,Rfmt) param%t0; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "TSTOP"; write(param_value, Rfmt) param%tstop; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "DT"; write(param_value, Rfmt) param%dt; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "CB_IN"; write(param_value, *) trim(adjustl(param%incbfile)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "PL_IN"; write(param_value, *) trim(adjustl(param%inplfile)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "TP_IN"; write(param_value, *) trim(adjustl(param%intpfile)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "IN_TYPE"; write(param_value, *) trim(adjustl(param%in_type)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "IN_FORM"; write(param_value, *) trim(adjustl(param%in_form)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - if (param%istep_dump > 0) write(param_name, *) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("T0", param%t0, unit) + call io_param_writer_one("TSTOP", param%tstop, unit) + call io_param_writer_one("DT", param%dt, unit) + call io_param_writer_one("CB_IN", param%incbfile, unit) + call io_param_writer_one("PL_IN", param%inplfile, unit) + call io_param_writer_one("TP_IN", param%intpfile, unit) + call io_param_writer_one("IN_TYPE", param%in_type, unit) + call io_param_writer_one("IN_FORM", param%in_form, unit) + if (param%istep_dump > 0) call io_param_writer_one("ISTEP_DUMP",param%istep_dump, unit) if (param%istep_out > 0) then - write(param_name, *) "ISTEP_OUT"; write(param_value, Ifmt) param%istep_out; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "BIN_OUT"; write(param_value, *) trim(adjustl(param%outfile)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "OUT_TYPE"; write(param_value, *) trim(adjustl(param%out_type)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "OUT_FORM"; write(param_value, *) trim(adjustl(param%out_form)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "OUT_STAT"; write(param_value, *) "APPEND"; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("ISTEP_OUT", param%istep_out, unit) + call io_param_writer_one("BIN_OUT", param%outfile, unit) + call io_param_writer_one("OUT_TYPE", param%out_type, unit) + call io_param_writer_one("OUT_FORM", param%out_form, unit) + call io_param_writer_one("OUT_STAT", "APPEND", unit) end if - write(param_name, *) "PARTICLE_OUT"; write(param_value, *) trim(adjustl(param%particle_out)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("PARTICLE_OUT", param%particle_out, unit) if (param%enc_out /= "") then - write(param_name, *) "ENC_OUT"; write(param_value, *) trim(adjustl(param%enc_out)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("ENC_OUT", param%enc_out, unit) end if - write(param_name, *) "CHK_RMIN"; write(param_value, Rfmt) param%rmin; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "CHK_RMAX"; write(param_value, Rfmt) param%rmax; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "CHK_EJECT"; write(param_value, Rfmt) param%rmaxu; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "CHK_QMIN"; write(param_value, Rfmt) param%qmin; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("CHK_RMIN", param%rmin, unit) + call io_param_writer_one("CHK_RMAX", param%rmax, unit) + call io_param_writer_one("CHK_EJECT", param%rmaxu, unit) + call io_param_writer_one("CHK_QMIN", param%qmin, unit) if (param%qmin >= 0.0_DP) then - write(param_name, *) "CHK_QMIN_COORD"; write(param_value, *) trim(adjustl(param%qmin_coord)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - allocate(param_array(2)) - write(param_array(1)%value, Rfmt) param%qmin_alo - write(param_array(2)%value, Rfmt) param%qmin_ahi - write(param_name, *) "CHK_QMIN_RANGE"; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_array(1)%value), adjustl(param_array(2)%value) + call io_param_writer_one("CHK_QMIN_COORD", param%qmin_coord, unit) + call io_param_writer_one("CHK_QMIN_RANGE", [param%qmin_alo, param%qmin_ahi], unit) end if - write(param_name, *) "MU2KG"; write(param_value, Rfmt) param%MU2KG; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "TU2S"; write(param_value, Rfmt) param%TU2S ; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "DU2M"; write(param_value, Rfmt) param%DU2M; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "RHILL_PRESENT"; write(param_value, Lfmt) param%lrhill_present; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "EXTRA_FORCE"; write(param_value, Lfmt) param%lextra_force; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("MU2KG", param%MU2KG, unit) + call io_param_writer_one("TU2S", param%TU2S , unit) + call io_param_writer_one("DU2M", param%DU2M, unit) + call io_param_writer_one("RHILL_PRESENT", param%lrhill_present, unit) + call io_param_writer_one("EXTRA_FORCE", param%lextra_force, unit) if (param%discard_out /= "") then - write(param_name, *) "DISCARD_OUT"; write(param_value, *) trim(adjustl(param%discard_out)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("DISCARD_OUT", param%discard_out, unit) end if if (param%discard_out /= "") then - write(param_name, *) "BIG_DISCARD"; write(param_value, Lfmt) param%lbig_discard; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("BIG_DISCARD", param%lbig_discard, unit) end if - write(param_name, *) "CHK_CLOSE"; write(param_value, Lfmt) param%lclose; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "ENERGY"; write(param_value, Lfmt) param%lenergy; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("CHK_CLOSE", param%lclose, unit) + call io_param_writer_one("ENERGY", param%lenergy, unit) if (param%lenergy .and. (param%energy_out /= "")) then - write(param_name, *) "ENERGY_OUT"; write(param_value, *) trim(adjustl(param%energy_out)); write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("ENERGY_OUT", param%energy_out, unit) end if - write(param_name, *) "GR"; write(param_value, Lfmt) param%lgr; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "FLATTEN_INTERACTIONS"; write(param_value, Lfmt) param%lflatten_interactions; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("GR", param%lgr, unit) + call io_param_writer_one("ROTATION", param%lrotation, unit) + call io_param_writer_one("TIDES", param%ltides, unit) + call io_param_writer_one("INTERACTION_LOOPS", param%interaction_loops, unit) if (param%lenergy) then - write(param_name, *) "FIRSTENERGY"; write(param_value, Lfmt) param%lfirstenergy; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "EORBIT_ORIG"; write(param_value, Rfmt) param%Eorbit_orig; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "GMTOT_ORIG"; write(param_value, Rfmt) param%GMtot_orig; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "LTOT_ORIG"; write(v1, Rfmt) param%Ltot_orig(1); write(v2, Rfmt) param%Ltot_orig(2); write(v3, Rfmt) param%Ltot_orig(3) - write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(v1)) // " " // trim(adjustl(v2)) // " " // trim(adjustl(v3)) - write(param_name, *) "LORBIT_ORIG"; write(v1, Rfmt) param%Lorbit_orig(1); write(v2, Rfmt) param%Lorbit_orig(2); write(v3, Rfmt) param%Lorbit_orig(3) - write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(v1)) // " " // trim(adjustl(v2)) // " " // trim(adjustl(v3)) - write(param_name, *) "LSPIN_ORIG"; write(v1, Rfmt) param%Lspin_orig(1); write(v2, Rfmt) param%Lspin_orig(2); write(v3, Rfmt) param%Lspin_orig(3) - write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(v1)) // " " // trim(adjustl(v2)) // " " // trim(adjustl(v3)) - write(param_name, *) "LESCAPE"; write(v1, Rfmt) param%Lescape(1); write(v2, Rfmt) param%Lescape(2); write(v3, Rfmt) param%Lescape(3) - write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(v1)) // " " // trim(adjustl(v2)) // " " // trim(adjustl(v3)) - - write(param_name, *) "GMESCAPE"; write(param_value, Rfmt) param%GMescape; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "ECOLLISIONS"; write(param_value, Rfmt) param%Ecollisions; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "EUNTRACKED"; write(param_value, Rfmt) param%Euntracked; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit) + call io_param_writer_one("EORBIT_ORIG", param%Eorbit_orig, unit) + call io_param_writer_one("GMTOT_ORIG", param%GMtot_orig, unit) + call io_param_writer_one("LTOT_ORIG", param%Ltot_orig(:), unit) + call io_param_writer_one("LORBIT_ORIG", param%Lorbit_orig(:), unit) + call io_param_writer_one("LSPIN_ORIG", param%Lspin_orig(:), unit) + call io_param_writer_one("LESCAPE", param%Lspin_orig(:), unit) + call io_param_writer_one("GMESCAPE",param%GMescape, unit) + call io_param_writer_one("ECOLLISIONS",param%Ecollisions, unit) + call io_param_writer_one("EUNTRACKED",param%Euntracked, unit) end if - write(param_name, *) "FIRSTKICK"; write(param_value, Lfmt) param%lfirstkick; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "MAXID"; write(param_value, Ifmt) param%maxid; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit) + call io_param_writer_one("MAXID",param%maxid, unit) iostat = 0 iomsg = "UDIO not implemented" end associate 667 continue - return end subroutine io_param_writer + module subroutine io_param_writer_one_char(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for character param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + character(len=*), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=NAMELEN) :: param_name_fixed_width !! Parameter label converted to fixed-width string + character(len=STRMAX) :: iomsg !! Message to pass if iostat /= 0 + + write(param_name_fixed_width, *) param_name + write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name_fixed_width) // " " // trim(adjustl(param_value)) + + return + 667 continue + write(*,*) 'Error writing parameter: ',trim(adjustl(iomsg)) + end subroutine io_param_writer_one_char + + + module subroutine io_param_writer_one_DP(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for real(DP) param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + real(DP), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=STRMAX) :: param_value_string !! Parameter value converted to a string + character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values + + write(param_value_string,Rfmt) param_value + call io_param_writer_one(param_name, param_value_string, unit) + + return + end subroutine io_param_writer_one_DP + + + module subroutine io_param_writer_one_DParr(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for real(DP) arrays () param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + real(DP), dimension(:), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=STRMAX) :: param_value_string !! Parameter value converted to a string + character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values + character(len=25) :: arr_val + integer(I4B) :: i, narr + + narr = size(param_value) + do i = 1, narr + write(arr_val, Rfmt) param_value(i) + if (i == 1) then + write(param_value_string, *) arr_val + else + param_value_string = trim(adjustl(param_value_string)) // " " // arr_val + end if + end do + + call io_param_writer_one(param_name, param_value_string, unit) + + return + end subroutine io_param_writer_one_DParr + + + module subroutine io_param_writer_one_I4B(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for integer(I4B) param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + integer(I4B), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=STRMAX) :: param_value_string !! Parameter value converted to a string + character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values + + write(param_value_string,Ifmt) param_value + call io_param_writer_one(param_name, param_value_string, unit) + + return + end subroutine io_param_writer_one_I4B + + + module subroutine io_param_writer_one_I8B(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for integer(I8B) param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + integer(I8B), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=STRMAX) :: param_value_string !! Parameter value converted to a string + character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values + + write(param_value_string,Ifmt) param_value + call io_param_writer_one(param_name, param_value_string, unit) + + return + end subroutine io_param_writer_one_I8B + + + module subroutine io_param_writer_one_I4Barr(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for integer(I4B) arrays param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + integer(I4B), dimension(:), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=STRMAX) :: param_value_string !! Parameter value converted to a string + character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values + character(len=25) :: arr_val + integer(I4B) :: i, narr + + narr = size(param_value) + do i = 1, narr + write(arr_val, Ifmt) param_value(i) + if (i == 1) then + write(param_value_string, *) trim(adjustl(arr_val)) + else + param_value_string = trim(adjustl(param_value_string)) // " " // trim(adjustl(arr_val)) + end if + end do + + call io_param_writer_one(param_name, param_value_string, unit) + + return + end subroutine io_param_writer_one_I4Barr + + + module subroutine io_param_writer_one_logical(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for logical param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + logical, intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=STRMAX) :: param_value_string !! Parameter value converted to a string + character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values + + write(param_value_string,Lfmt) param_value + call io_param_writer_one(param_name, param_value_string, unit) + + return + end subroutine io_param_writer_one_logical + + + module subroutine io_param_writer_one_QP(param_name, param_value, unit) + !! author: David A. Minton + !! + !! Writes a single parameter name/value pair to a file unit. + !! This version is for real(QP) param_value type + implicit none + ! Arguments + character(len=*), intent(in) :: param_name !! Name of parameter to print + real(QP), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + ! Internals + character(len=STRMAX) :: param_value_string !! Parameter value converted to a string + character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values + + write(param_value_string,Rfmt) param_value + call io_param_writer_one(param_name, param_value_string, unit) + + return + end subroutine io_param_writer_one_QP + + module 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 !! diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a1f60c0dd..af367c5ed 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -75,7 +75,6 @@ module swiftest_classes integer(I4B) :: discard_vhy_varid !! NetCDF ID for the heliocentric velocity of the body at the time of discard y variable integer(I4B) :: discard_vhz_varid !! NetCDF ID for the heliocentric velocity of the body at the time of discard z variable integer(I4B) :: discard_body_id_varid !! NetCDF ID for the id of the other body involved in the discard - contains procedure :: close => netcdf_close !! Closes an open NetCDF file procedure :: initialize => netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object @@ -122,7 +121,8 @@ module swiftest_classes real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters real(DP) :: GU = -1.0_DP !! Universal gravitational constant in the system units real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units - character(STRMAX) :: energy_out = "" !! Name of output energy and momentum report file + character(STRMAX) :: energy_out = "" !! Name of output energy and momentum report file + character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" ! Logical flags to turn on or off various features of the code logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually) @@ -728,7 +728,67 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) integer(I4B), intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine io_param_writer + end interface + + interface io_param_writer_one + module subroutine io_param_writer_one_char(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + character(len=*), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + end subroutine io_param_writer_one_char + + module subroutine io_param_writer_one_DP(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + real(DP), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + end subroutine io_param_writer_one_DP + + module subroutine io_param_writer_one_DParr(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + real(DP), dimension(:), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + end subroutine io_param_writer_one_DParr + + module subroutine io_param_writer_one_I4B(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + integer(I4B), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + end subroutine io_param_writer_one_I4B + module subroutine io_param_writer_one_I4Barr(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + integer(I4B), dimension(:), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + end subroutine io_param_writer_one_I4Barr + + module subroutine io_param_writer_one_I8B(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + integer(I8B), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + end subroutine io_param_writer_one_I8B + + module subroutine io_param_writer_one_logical(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + logical, intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + end subroutine io_param_writer_one_logical + + module subroutine io_param_writer_one_QP(param_name, param_value, unit) + implicit none + character(len=*), intent(in) :: param_name !! Name of parameter to print + real(QP), intent(in) :: param_value !! Value of parameter to print + integer(I4B), intent(in) :: unit !! Open file unit number to print parameter to + 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 diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 6425b5b5f..82fda6886 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -17,7 +17,16 @@ module walltime_classes procedure :: reset => walltime_reset !! Resets the clock ticker, settting main_start to the current ticker value procedure :: start => walltime_start !! Starts the timer, setting step_start to the current ticker value procedure :: finish => walltime_finish !! Ends the timer, setting step_finish to the current ticker value and printing the elapsed time information to the terminal - end type + end type walltimer + + type, extends(walltimer) :: interaction_timer + integer(I8B) :: max_interactions = huge(1_I8B) + integer(I4B) :: step_counter + integer(I8B) :: count_previous + character(len=NAMELEN) :: current_style + contains + + end type interaction_timer interface module subroutine walltime_finish(self, nsubsteps, message) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index c97904079..2e5e87941 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -27,7 +27,6 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms character(len=*),parameter :: linefmt = '(A)' associate(param => self) - call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) call random_seed(size = nseeds) if (allocated(param%seed)) deallocate(param%seed) @@ -85,29 +84,22 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms write(iomsg,*) "GMTINY invalid or not set: ", self%GMTINY iostat = -1 return - else - write(*,*) "GMTINY = ", self%GMTINY end if - write(*,*) "FRAGMENTATION = ", param%lfragmentation if (param%lfragmentation) then if (seed_set) then call random_seed(put = param%seed) else call random_seed(get = param%seed) end if - write(*,*) "SEED: N,VAL = ",size(param%seed), param%seed(:) if (param%min_GMfrag < 0.0_DP) param%min_GMfrag = param%GMTINY - write(*,*) "MIN_GMFRAG = ", self%min_GMfrag end if - if (.not.self%lclose) then - write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.' - iostat = -1 - return - end if ! All reporting of collision information in SyMBA (including mergers) is now recorded in the Fraggle logfile call fraggle_io_log_start(param) + + ! Call the base method (which also prints the contents to screen) + call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) end associate iostat = 0 @@ -134,43 +126,19 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals - character(*),parameter :: Ifmt = '(I0)' !! Format label for integer values - character(*),parameter :: Rfmt = '(ES25.17)' !! Format label for real values - character(*),parameter :: Rarrfmt = '(3(ES25.17,1X))' !! Format label for real values - character(*),parameter :: Lfmt = '(L1)' !! Format label for logical values - character(len=NAMELEN) :: param_name - character(len=STRMAX) :: param_value - type character_array - character(25) :: value - end type character_array - type(character_array), dimension(:), allocatable :: param_array - integer(I4B) :: i, nstr + integer(I4B) :: nseeds associate(param => self) call io_param_writer(param, unit, iotype, v_list, iostat, iomsg) ! Special handling is required for writing the random number seed array as its size is not known until runtime ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements - write(param_name, *) "GMTINY"; write(param_value, Rfmt) param%GMTINY; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "MIN_GMFRAG"; write(param_value, Rfmt) param%min_GMfrag; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) - write(param_name, *) "FRAGMENTATION"; write(param_value, Lfmt) param%lfragmentation; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + call io_param_writer_one("GMTINY",param%GMTINY, unit) + call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit) + call io_param_writer_one("FRAGMENTATION",param%lfragmentation, unit) if (param%lfragmentation) then - write(param_name, *) "SEED" - if (allocated(param_array)) deallocate(param_array) - allocate(param_array(0:size(param%seed))) - write(param_array(0)%value, Ifmt) size(param%seed) - do i = 1, size(param%seed) - write(param_array(i)%value, Ifmt) param%seed(i) - end do - write(unit, '(" ",A32)', advance='no', err = 667, iomsg = iomsg) adjustl(param_name) - do i = 0, size(param%seed) - nstr = len(trim(adjustl(param_array(i)%value))) - if (i < size(param%seed)) then - write(unit, '(A12)', advance='no', err = 667, iomsg = iomsg) trim(adjustl(param_array(i)%value)) // " " - else - write(unit, '(A12)', err = 667, iomsg = iomsg) trim(adjustl(param_array(i)%value)) - end if - end do + nseeds = size(param%seed) + call io_param_writer_one("SEED", [nseeds, param%seed(:)], unit) end if iostat = 0 From 22be9947586cfed521d230103f9fcc6eeafb1f06 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 12:25:43 -0400 Subject: [PATCH 02/15] Refactored util_index to util_flatten as it is more descriptive and consistent with the rest of the terminology. Put in safety check in case allocating k_plpl fails --- src/fraggle/fraggle_generate.f90 | 4 +- src/fraggle/fraggle_util.f90 | 4 +- src/helio/helio_setup.f90 | 2 +- src/modules/fraggle_classes.f90 | 4 +- src/modules/swiftest_classes.f90 | 22 +++++----- src/modules/symba_classes.f90 | 8 ++-- src/modules/walltime_classes.f90 | 2 + src/modules/whm_classes.f90 | 7 ---- src/setup/setup.f90 | 2 +- src/symba/symba_util.f90 | 10 ++--- src/util/{util_index.f90 => util_flatten.f90} | 40 ++++++++++--------- src/whm/whm_setup.f90 | 2 +- 12 files changed, 53 insertions(+), 54 deletions(-) rename src/util/{util_index.f90 => util_flatten.f90} (78%) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 1c189c5ca..486aed2ea 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -18,7 +18,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa class(fraggle_fragments), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? ! Internals integer(I4B) :: i @@ -140,7 +140,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call frag%set_original_scale(colliders) ! Restore the big array - if (lk_plpl) call pl%index(param) + if (lk_plpl) call pl%flatten(param) end associate call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index df8b89c41..f65bd81c5 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -135,7 +135,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa ! Internals integer(I4B) :: i, nplm @@ -172,7 +172,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) end if - call tmpsys%pl%index(param) + call tmpsys%pl%flatten(param) call tmpsys%get_energy_and_momentum(param) diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_setup.f90 index cf5f57d8a..5839628a0 100644 --- a/src/helio/helio_setup.f90 +++ b/src/helio/helio_setup.f90 @@ -16,7 +16,7 @@ module subroutine helio_setup_initialize_system(self, param) call self%pl%h2b(self%cb) call self%tp%h2b(self%cb) call self%pl%sort("mass", ascending=.false.) - call self%pl%index(param) + call self%pl%flatten(param) return end subroutine helio_setup_initialize_system diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 5a8a4a392..2f6d75a84 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -104,7 +104,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fraggle_generate_fragments @@ -249,7 +249,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa end subroutine fraggle_util_get_energy_momentum diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index af367c5ed..dc5eacc3c 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -133,7 +133,7 @@ module swiftest_classes logical :: loblatecb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input) logical :: lrotation = .false. !! Include rotation states of big bodies logical :: ltides = .false. !! Include tidal dissipation - logical :: lflatten_interactions = .true. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory) + logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory) ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values) real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy @@ -324,7 +324,7 @@ module swiftest_classes ! Massive body-specific concrete methods ! These are concrete because they are the same implemenation for all integrators procedure :: discard => discard_pl !! Placeholder method for discarding massive bodies - procedure :: index => util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix + procedure :: flatten => util_flatten_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: accel_int => kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies procedure :: accel_obl => obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: setup => setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays @@ -561,34 +561,34 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) end subroutine drift_one - module pure subroutine util_index_eucl_ij_to_k(n, i, j, k) - !$omp declare simd(util_index_eucl_ij_to_k) + module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k) + !$omp declare simd(util_flatten_eucl_ij_to_k) implicit none integer(I4B), intent(in) :: n !! Number of bodies integer(I4B), intent(in) :: i !! Index of the ith body integer(I4B), intent(in) :: j !! Index of the jth body integer(I8B), intent(out) :: k !! Index of the flattened matrix - end subroutine util_index_eucl_ij_to_k + end subroutine util_flatten_eucl_ij_to_k - module pure subroutine util_index_eucl_k_to_ij(n, k, i, j) + module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j) implicit none integer(I4B), intent(in) :: n !! Number of bodies integer(I8B), intent(in) :: k !! Index of the flattened matrix integer(I4B), intent(out) :: i !! Index of the ith body integer(I4B), intent(out) :: j !! Index of the jth body - end subroutine util_index_eucl_k_to_ij + end subroutine util_flatten_eucl_k_to_ij - module subroutine util_index_eucl_plpl(self, param) + module subroutine util_flatten_eucl_plpl(self, param) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine - module subroutine util_index_eucl_pltp(self, pl, param) + module subroutine util_flatten_eucl_pltp(self, pl, param) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine module pure subroutine gr_kick_getaccb_ns_body(self, system, param) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 0faefc3d1..52c66d293 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -67,7 +67,7 @@ module symba_classes type(symba_kinship), dimension(:), allocatable :: kin !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step contains procedure :: make_colliders => symba_collision_make_colliders_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family - procedure :: index => symba_util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix + procedure :: flatten => symba_util_flatten_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: discard => symba_discard_pl !! Process massive body discards procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other @@ -316,12 +316,12 @@ module function symba_collision_casemerge(system, param, colliders, frag) result integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casemerge - module subroutine symba_util_index_eucl_plpl(self, param) + module subroutine symba_util_flatten_eucl_plpl(self, param) use swiftest_classes, only : swiftest_parameters implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_util_index_eucl_plpl + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_util_flatten_eucl_plpl module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 82fda6886..c27a5b06a 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -6,6 +6,8 @@ module walltime_classes implicit none public + integer(I4B) :: INTERACTION_TIMER_CADENCE = 1000 !! Minimum number of steps to wait before timing an interaction loop in ADAPTIVE mode + type :: walltimer integer(I8B) :: count_rate !! Rate at wich the clock ticks integer(I8B) :: count_max !! Maximum value of the clock ticker diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index bd17f98e0..0bae1f654 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -108,13 +108,6 @@ module subroutine whm_drift_pl(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine whm_drift_pl - module subroutine whm_util_index_eucl_plpl(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine whm_util_index_eucl_plpl - !> Get heliocentric accelration of massive bodies module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_cb, swiftest_parameters diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 18d8f2189..2dc57d952 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -175,7 +175,7 @@ module subroutine setup_initialize_system(self, param) call self%pl%el2xv(self%cb) call self%tp%el2xv(self%cb) end if - call self%pl%index(param) + call self%pl%flatten(param) if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) self%pl%lfirst = param%lfirstkick self%tp%lfirst = param%lfirstkick diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 0d062894f..4b20cf0bf 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -269,7 +269,7 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) end subroutine symba_util_fill_tp - module subroutine symba_util_index_eucl_plpl(self, param) + module subroutine symba_util_flatten_eucl_plpl(self, param) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix. This also sets the lmtiny flag and computes the @@ -283,7 +283,7 @@ module subroutine symba_util_index_eucl_plpl(self, param) implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I8B) :: k, nplpl, nplplm integer(I4B) :: i, j, npl, nplm, ip, jp @@ -299,7 +299,7 @@ module subroutine symba_util_index_eucl_plpl(self, param) allocate(self%k_plpl(2, pl%nplpl)) do concurrent (i = 1:npl) do concurrent (j = i+1:npl) - call util_index_eucl_ij_to_k(npl, i, j, k) + call util_flatten_eucl_ij_to_k(npl, i, j, k) self%k_plpl(1, k) = i self%k_plpl(2, k) = j end do @@ -308,7 +308,7 @@ module subroutine symba_util_index_eucl_plpl(self, param) end associate return - end subroutine symba_util_index_eucl_plpl + end subroutine symba_util_flatten_eucl_plpl module subroutine symba_util_peri_pl(self, system, param) @@ -478,7 +478,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! Reindex the new list of bodies call pl%sort("mass", ascending=.false.) - call pl%index(param) + call pl%flatten(param) ! Reset the kinship trackers call pl%reset_kinship([(i, i=1, npl)]) diff --git a/src/util/util_index.f90 b/src/util/util_flatten.f90 similarity index 78% rename from src/util/util_index.f90 rename to src/util/util_flatten.f90 index d0dd83896..e3e97c1e7 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_flatten.f90 @@ -2,8 +2,8 @@ use swiftest contains - module pure subroutine util_index_eucl_ij_to_k(n, i, j, k) - !$omp declare simd(util_index_eucl_ij_to_k) + module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k) + !$omp declare simd(util_flatten_eucl_ij_to_k) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions. @@ -27,10 +27,10 @@ module pure subroutine util_index_eucl_ij_to_k(n, i, j, k) k = (i8 - 1_I8B) * n8 - i8 * (i8 - 1_I8B) / 2_I8B + (j8 - i8) return - end subroutine util_index_eucl_ij_to_k + end subroutine util_flatten_eucl_ij_to_k - module pure subroutine util_index_eucl_k_to_ij(n, k, i, j) + module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j) !! author: Jacob R. Elliott and David A. Minton !! !! Turns k index into i,j indices for use in the Euclidean distance matrix for pl-pl interactions. @@ -59,10 +59,10 @@ module pure subroutine util_index_eucl_k_to_ij(n, k, i, j) j = int(j8, kind=I4B) return - end subroutine util_index_eucl_k_to_ij + end subroutine util_flatten_eucl_k_to_ij - module subroutine util_index_eucl_plpl(self, param) + module subroutine util_flatten_eucl_plpl(self, param) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions for a Swiftest massive body object @@ -74,9 +74,9 @@ module subroutine util_index_eucl_plpl(self, param) implicit none ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i, j, npl + integer(I4B) :: i, j, npl, err integer(I8B) :: k npl = int(self%nbody, kind=I8B) @@ -84,20 +84,24 @@ module subroutine util_index_eucl_plpl(self, param) nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl if (param%lflatten_interactions) then if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously - allocate(self%k_plpl(2, nplpl)) - do concurrent (i=1:npl, j=1:npl, j>i) - call util_index_eucl_ij_to_k(npl, i, j, k) - self%k_plpl(1, k) = i - self%k_plpl(2, k) = j - end do + allocate(self%k_plpl(2, nplpl), stat=err) + if (err /=0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode + param%lflatten_interactions = .false. + else + do concurrent (i=1:npl, j=1:npl, j>i) + call util_flatten_eucl_ij_to_k(npl, i, j, k) + self%k_plpl(1, k) = i + self%k_plpl(2, k) = j + end do + end if end if end associate return - end subroutine util_index_eucl_plpl + end subroutine util_flatten_eucl_plpl - module subroutine util_index_eucl_pltp(self, pl, param) + module subroutine util_flatten_eucl_pltp(self, pl, param) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-tp interactions @@ -110,7 +114,7 @@ module subroutine util_index_eucl_pltp(self, pl, param) ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I8B) :: i, j, counter, npl, ntp @@ -131,6 +135,6 @@ module subroutine util_index_eucl_pltp(self, pl, param) end associate return - end subroutine util_index_eucl_pltp + end subroutine util_flatten_eucl_pltp end submodule s_util_index diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index 3c5cdeec3..60be140fd 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -82,7 +82,7 @@ module subroutine whm_setup_initialize_system(self, param) ! First we need to make sure that the massive bodies are sorted by heliocentric distance before computing jacobies call util_set_ir3h(self%pl) call self%pl%sort("ir3h", ascending=.false.) - call self%pl%index(param) + call self%pl%flatten(param) ! Make sure that the discard list gets allocated initially call self%tp_discards%setup(0, param) From 7199148a061a6ec7732628d5d2c7ccb2ca9e3232 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 13:59:56 -0400 Subject: [PATCH 03/15] Refactored to allow param to be inout in kick subroutines --- src/fraggle/fraggle_placeholder.f90 | 4 +-- src/helio/helio_kick.f90 | 8 ++--- src/kick/kick.f90 | 10 ++++-- src/main/swiftest_driver.f90 | 6 ++-- src/modules/fraggle_classes.f90 | 4 +-- src/modules/helio_classes.f90 | 10 +++--- src/modules/rmvs_classes.f90 | 2 +- src/modules/swiftest_classes.f90 | 8 ++--- src/modules/symba_classes.f90 | 8 ++--- src/modules/walltime_classes.f90 | 50 +++++++++++++++++------------ src/modules/whm_classes.f90 | 8 ++--- src/rmvs/rmvs_kick.f90 | 2 +- src/symba/symba_kick.f90 | 7 ++-- src/whm/whm_kick.f90 | 8 ++--- 14 files changed, 75 insertions(+), 60 deletions(-) diff --git a/src/fraggle/fraggle_placeholder.f90 b/src/fraggle/fraggle_placeholder.f90 index bbf08bb04..f8284b6b4 100644 --- a/src/fraggle/fraggle_placeholder.f90 +++ b/src/fraggle/fraggle_placeholder.f90 @@ -8,7 +8,7 @@ module subroutine fraggle_placeholder_accel(self, system, param, t, lbeg) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step write(*,*) "The type-bound procedure 'accel' is not defined for type fraggle_fragments" @@ -19,7 +19,7 @@ module subroutine fraggle_placeholder_kick(self, system, param, t, dt, lbeg) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index 928e7c197..2ed23654d 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -13,7 +13,7 @@ module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg) ! Arguments class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step @@ -56,7 +56,7 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) ! Arguments class(helio_tp), intent(inout) :: self !! Helio test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step @@ -89,7 +89,7 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, lbeg) ! Arguments class(helio_pl), intent(inout) :: self !! Swiftest generic body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. @@ -128,7 +128,7 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg) ! Arguments class(helio_tp), intent(inout) :: self !! Swiftest generic body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 766e829d6..0e90e73ee 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -12,8 +12,14 @@ module subroutine kick_getacch_int_pl(self, param) implicit none ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters + ! Internals + type(interaction_timer), save :: itimer + logical, save :: lfirst + if (lfirst) then + call itimer%reset(param) + end if if (param%lflatten_interactions) then call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) else @@ -34,7 +40,7 @@ module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl) implicit none ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors integer(I4B), intent(in) :: npl !! Number of active massive bodies diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index b4260f115..bd4d50855 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -70,7 +70,7 @@ program swiftest_driver !$ write(*,'(a)') ' OpenMP parameters:' !$ write(*,'(a)') ' ------------------' !$ write(*,'(a,i3,/)') ' Number of threads = ', nthreads - call timer%reset() + call timer%reset(param) write(*, *) " *************** Main Loop *************** " do iloop = 1, nloops !> Step the system forward in time @@ -86,9 +86,9 @@ program swiftest_driver iout = iout - 1 if (iout == 0) then ioutput = ioutput_t0 + iloop / istep_out - call timer%finish(nsubsteps=istep_out, message="Integration steps:") + call timer%finish(nsubsteps=istep_out, message="Integration steps:", param=param) if (t > old_t_final) call nbody_system%write_frame(param) - call timer%finish(nsubsteps=1, message="File I/O: ") + call timer%finish(nsubsteps=1, message="File I/O: ", param=param) iout = istep_out end if end if diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 2f6d75a84..63440383f 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -142,7 +142,7 @@ module subroutine fraggle_placeholder_accel(self, system, param, t, lbeg) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step end subroutine fraggle_placeholder_accel @@ -152,7 +152,7 @@ module subroutine fraggle_placeholder_kick(self, system, param, t, dt, lbeg) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 3a7de3814..a74b609e1 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -130,7 +130,7 @@ end subroutine helio_gr_p4_pl module pure subroutine helio_gr_p4_tp(self, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(helio_tp), intent(inout) :: self !! Swiftest particle object + class(helio_tp), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Step size end subroutine helio_gr_p4_tp @@ -140,7 +140,7 @@ module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg) implicit none class(helio_pl), intent(inout) :: self !! Helio massive body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine helio_kick_getacch_pl @@ -150,7 +150,7 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) implicit none class(helio_tp), intent(inout) :: self !! Helio test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine helio_kick_getacch_tp @@ -160,7 +160,7 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, lbeg) implicit none class(helio_pl), intent(inout) :: self !! Helio massive body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. @@ -171,7 +171,7 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg) implicit none class(helio_tp), intent(inout) :: self !! Helio test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 680612de5..2fa6452f9 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -144,7 +144,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) implicit none class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine rmvs_kick_getacch_tp diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index dc5eacc3c..89c5ef238 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -463,7 +463,7 @@ subroutine abstract_accel(self, system, param, t, lbeg) import swiftest_body, swiftest_nbody_system, swiftest_parameters, DP class(swiftest_body), intent(inout) :: self !! Swiftest body data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step end subroutine abstract_accel @@ -473,7 +473,7 @@ subroutine abstract_kick_body(self, system, param, t, dt, lbeg) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest generic body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. @@ -903,13 +903,13 @@ end subroutine io_write_hdr_system module subroutine kick_getacch_int_pl(self, param) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters end subroutine kick_getacch_int_pl module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors integer(I4B), intent(in) :: npl !! Number of active massive bodies diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 52c66d293..6bfef7bd9 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -355,7 +355,7 @@ end subroutine symba_io_write_discard module subroutine symba_kick_getacch_int_pl(self, param) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters end subroutine symba_kick_getacch_int_pl module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) @@ -363,7 +363,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine symba_kick_getacch_pl @@ -373,14 +373,14 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine symba_kick_getacch_tp module subroutine symba_kick_encounter(self, system, dt, irec, sgn) implicit none - class(symba_encounter), intent(in) :: self !! SyMBA pl-tp encounter list object + class(symba_encounter), intent(in) :: self !! SyMBA pl-tp encounter list object class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index c27a5b06a..76600660d 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -3,6 +3,7 @@ module walltime_classes !! !! Classes and methods used to compute elasped wall time use swiftest_globals + use swiftest_classes, only : swiftest_parameters implicit none public @@ -31,36 +32,42 @@ module walltime_classes end type interaction_timer interface - module subroutine walltime_finish(self, nsubsteps, message) + module subroutine walltime_finish(self, nsubsteps, message, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(walltimer), intent(inout) :: self !! Walltimer object - integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step - character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + class(walltimer), intent(inout) :: self !! Walltimer object + integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step + character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine walltime_finish - module subroutine walltime_reset(self) + module subroutine walltime_reset(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(walltimer), intent(inout) :: self !! Walltimer object + class(walltimer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine walltime_reset - module subroutine walltime_start(self) + module subroutine walltime_start(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(walltimer), intent(inout) :: self !! Walltimer object + class(walltimer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine walltime_start end interface contains - module subroutine walltime_finish(self, nsubsteps, message) + module subroutine walltime_finish(self, nsubsteps, message, param) !! author: David A. Minton !! !! Ends the timer, setting step_finish to the current ticker value and printing the elapsed time information to the terminal - use swiftest_globals implicit none ! Arguments - class(walltimer), intent(inout) :: self !! Walltimer object - integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step - character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + class(walltimer), intent(inout) :: self !! Walltimer object + integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step + character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals character(len=*), parameter :: walltimefmt = '" Wall time (s): ", es12.5, "; Wall time/step in this interval (s): ", es12.5' character(len=STRMAX) :: fmt @@ -85,36 +92,37 @@ module subroutine walltime_finish(self, nsubsteps, message) fmt = '("' // adjustl(message) // '",' // walltimefmt // ')' write(*,trim(adjustl(fmt))) wall_main, wall_per_substep - call self%start() + call self%start(param) return end subroutine walltime_finish - module subroutine walltime_reset(self) + + module subroutine walltime_reset(self, param) !! author: David A. Minton !! !! Resets the clock ticker, settting main_start to the current ticker value - use swiftest_globals implicit none ! Arguments - class(walltimer), intent(inout) :: self + class(walltimer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters call system_clock(self%count_start_main, self%count_rate, self%count_max) self%lmain_is_started = .true. - call self%start() + call self%start(param) return end subroutine walltime_reset - module subroutine walltime_start(self) + module subroutine walltime_start(self, param) !! author: David A. Minton !! !! Starts the timer, setting step_start to the current ticker value - use swiftest_globals implicit none ! Arguments - class(walltimer), intent(inout) :: self + class(walltimer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters if (.not.self%lmain_is_started) then write(*,*) "Wall timer error: Cannot start the step time until reset is called at least once!" diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 0bae1f654..3cad184dd 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -114,7 +114,7 @@ module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg) implicit none class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine whm_kick_getacch_pl @@ -125,7 +125,7 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) implicit none class(whm_tp), intent(inout) :: self !! WHM test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine whm_kick_getacch_tp @@ -135,7 +135,7 @@ module subroutine whm_kick_vh_pl(self, system, param, t, dt, lbeg) implicit none class(whm_pl), intent(inout) :: self !! WHM massive body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. @@ -146,7 +146,7 @@ module subroutine whm_kick_vh_tp(self, system, param, t, dt, lbeg) implicit none class(whm_tp), intent(inout) :: self !! WHM test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. diff --git a/src/rmvs/rmvs_kick.f90 b/src/rmvs/rmvs_kick.f90 index 6a8c0786e..57c04d44e 100644 --- a/src/rmvs/rmvs_kick.f90 +++ b/src/rmvs/rmvs_kick.f90 @@ -13,7 +13,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index fd0354612..e9b7d0405 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -12,7 +12,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameter + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameter if (param%lflatten_interactions) then call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) @@ -35,7 +35,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals @@ -68,6 +68,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) return end subroutine symba_kick_getacch_pl + module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) !! author: David A. Minton !! @@ -79,7 +80,7 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) ! Arguments class(symba_tp), intent(inout) :: self !! SyMBA test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index c36c1ce23..7b2b90d00 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -13,7 +13,7 @@ module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg) ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals @@ -68,7 +68,7 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) ! Arguments class(whm_tp), intent(inout) :: self !! WHM test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals @@ -201,7 +201,7 @@ module subroutine whm_kick_vh_pl(self, system, param, t, dt, lbeg) ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. @@ -243,7 +243,7 @@ module subroutine whm_kick_vh_tp(self, system, param, t, dt, lbeg) ! Arguments class(whm_tp), intent(inout) :: self !! WHM massive body object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. From eb3ce0a743fd8a94bfaab586a62d93927b0c2fbe Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 14:01:04 -0400 Subject: [PATCH 04/15] Put the walltimer_classes after swiftest_classes in the compilation sequence --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 50e5f1d89..17ef3c006 100644 --- a/Makefile +++ b/Makefile @@ -47,13 +47,13 @@ SWIFTEST_MODULES = swiftest_globals.f90 \ swiftest_operators.f90 \ lambda_function.f90\ - walltime_classes.f90 \ swiftest_classes.f90 \ fraggle_classes.f90 \ whm_classes.f90 \ rmvs_classes.f90 \ helio_classes.f90 \ symba_classes.f90 \ + walltime_classes.f90 \ swiftest.f90 From 17ba0f1879d8577a3af5f68e3e0ab0d3d759ec81 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 14:36:56 -0400 Subject: [PATCH 05/15] Moved file logging out of fraggle and into the main io submodule --- src/fraggle/fraggle_generate.f90 | 60 +++++++++++------------ src/fraggle/fraggle_io.f90 | 82 ++++++++++++++++---------------- src/fraggle/fraggle_regime.f90 | 2 +- src/io/io.f90 | 73 +++++++++++++++++++++++----- src/kick/kick.f90 | 4 +- src/modules/fraggle_classes.f90 | 11 ----- src/modules/swiftest_classes.f90 | 17 ++++++- src/modules/swiftest_globals.f90 | 7 +-- src/modules/walltime_classes.f90 | 33 ++++++++++++- src/rmvs/rmvs_io.f90 | 1 - src/symba/symba_collision.f90 | 26 +++++----- src/symba/symba_discard.f90 | 30 ++++++------ src/symba/symba_io.f90 | 3 +- 13 files changed, 214 insertions(+), 135 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 486aed2ea..af10af236 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -39,10 +39,10 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa associate(frag => self, nfrag => self%nbody, pl => system%pl) write(message,*) nfrag - call fraggle_io_log_one_message("Fraggle generating " // trim(adjustl(message)) // " fragments.") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") if (nfrag < NFRAG_MIN) then write(message,*) "Fraggle needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." - call fraggle_io_log_one_message(message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) lfailure = .true. return end if @@ -64,7 +64,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa try = 1 do while (try < MAXTRY) write(message,*) try - call fraggle_io_log_one_message("Fraggle try " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle try " // trim(adjustl(message))) if (lfailure) then call frag%restructure(colliders, try, f_spin, r_max_start) call frag%reset() @@ -87,19 +87,19 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call fraggle_generate_spins(frag, colliders, f_spin, lfailure) if (lfailure) then - call fraggle_io_log_one_message("Fraggle failed to find spins") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find spins") cycle end if call fraggle_generate_tan_vel(frag, colliders, lfailure) if (lfailure) then - call fraggle_io_log_one_message("Fraggle failed to find tangential velocities") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find tangential velocities") cycle end if call fraggle_generate_rad_vel(frag, colliders, lfailure) if (lfailure) then - call fraggle_io_log_one_message("Fraggle failed to find radial velocities") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed to find radial velocities") cycle end if @@ -110,14 +110,14 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = ((abs(dEtot + frag%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) if (lfailure) then write(message, *) dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL - call fraggle_io_log_one_message("Fraggle failed due to high energy error: " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high energy error: " // trim(adjustl(message))) cycle end if lfailure = ((abs(dLmag) / (.mag.frag%Ltot_before)) > FRAGGLE_LTOL) if (lfailure) then write(message,*) dLmag / (.mag.frag%Ltot_before(:)) - call fraggle_io_log_one_message("Fraggle failed due to high angular momentum error: " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high angular momentum error: " // trim(adjustl(message))) cycle end if @@ -126,14 +126,14 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = any(fpe_flag) if (.not.lfailure) exit write(message,*) "Fraggle failed due to a floating point exception: ", fpe_flag - call fraggle_io_log_one_message(message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) end do write(message,*) try if (lfailure) then - call fraggle_io_log_one_message("Fraggle fragment generation failed after " // trim(adjustl(message)) // " tries") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation failed after " // trim(adjustl(message)) // " tries") else - call fraggle_io_log_one_message("Fraggle fragment generation succeeded after " // trim(adjustl(message)) // " tries") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation succeeded after " // trim(adjustl(message)) // " tries") call fraggle_io_log_generate(frag) end if @@ -254,16 +254,16 @@ subroutine fraggle_generate_spins(frag, colliders, f_spin, lfailure) lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) if (lfailure) then - call fraggle_io_log_one_message(" ") - call fraggle_io_log_one_message("Spin failure diagnostics") + call io_log_one_message(FRAGGLE_LOG_OUT, " ") + call io_log_one_message(FRAGGLE_LOG_OUT, "Spin failure diagnostics") write(message, *) frag%ke_budget - call fraggle_io_log_one_message("ke_budget : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_budget : " // trim(adjustl(message))) write(message, *) frag%ke_spin - call fraggle_io_log_one_message("ke_spin : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_spin : " // trim(adjustl(message))) write(message, *) frag%ke_orbit - call fraggle_io_log_one_message("ke_orbit : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) write(message, *) frag%ke_budget - frag%ke_spin - frag%ke_orbit - call fraggle_io_log_one_message("ke_remainder : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) end if end associate @@ -355,20 +355,20 @@ subroutine fraggle_generate_tan_vel(frag, colliders, lfailure) ! If we are over the energy budget, flag this as a failure so we can try again lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP) if (lfailure) then - call fraggle_io_log_one_message(" ") - call fraggle_io_log_one_message("Tangential velocity failure diagnostics") + call io_log_one_message(FRAGGLE_LOG_OUT, " ") + call io_log_one_message(FRAGGLE_LOG_OUT, "Tangential velocity failure diagnostics") call frag%get_ang_mtm() L_frag_tot = frag%L_spin(:) + frag%L_orbit(:) write(message, *) .mag.(frag%L_budget(:) - L_frag_tot(:)) / (.mag.frag%Ltot_before(:)) - call fraggle_io_log_one_message("|L_remainder| : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "|L_remainder| : " // trim(adjustl(message))) write(message, *) frag%ke_budget - call fraggle_io_log_one_message("ke_budget : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_budget : " // trim(adjustl(message))) write(message, *) frag%ke_spin - call fraggle_io_log_one_message("ke_spin : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_spin : " // trim(adjustl(message))) write(message, *) frag%ke_orbit - call fraggle_io_log_one_message("ke_tangential : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_tangential : " // trim(adjustl(message))) write(message, *) frag%ke_budget - frag%ke_spin - frag%ke_orbit - call fraggle_io_log_one_message("ke_radial : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_radial : " // trim(adjustl(message))) end if end associate @@ -515,16 +515,16 @@ subroutine fraggle_generate_rad_vel(frag, colliders, lfailure) lfailure = abs((frag%ke_budget - (frag%ke_orbit + frag%ke_spin)) / frag%ke_budget) > FRAGGLE_ETOL if (lfailure) then - call fraggle_io_log_one_message(" ") - call fraggle_io_log_one_message("Radial velocity failure diagnostics") + call io_log_one_message(FRAGGLE_LOG_OUT, " ") + call io_log_one_message(FRAGGLE_LOG_OUT, "Radial velocity failure diagnostics") write(message, *) frag%ke_budget - call fraggle_io_log_one_message("ke_budget : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_budget : " // trim(adjustl(message))) write(message, *) frag%ke_spin - call fraggle_io_log_one_message("ke_spin : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_spin : " // trim(adjustl(message))) write(message, *) frag%ke_orbit - call fraggle_io_log_one_message("ke_orbit : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) write(message, *) frag%ke_budget - (frag%ke_orbit + frag%ke_spin) - call fraggle_io_log_one_message("ke_remainder : " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) end if end associate diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index dbd721216..9d7af73f1 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -61,24 +61,24 @@ module subroutine fraggle_io_log_generate(frag) end subroutine fraggle_io_log_generate - module subroutine fraggle_io_log_one_message(message) - !! author: David A. Minton - !! - !! Writes a single message to the fraggle log file - implicit none - ! Arguments - character(len=*), intent(in) :: message - ! Internals - character(STRMAX) :: errmsg - - open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *) trim(adjustl(message)) - close(FRAGGLE_LOG_UNIT) - - return - 667 continue - write(*,*) "Error writing Fraggle message to log file: " // trim(adjustl(errmsg)) - end subroutine fraggle_io_log_one_message + ! module subroutine io_log_one_message(FRAGGLE_LOG_OUT, message) + ! !! author: David A. Minton + ! !! + ! !! Writes a single message to the fraggle log file + ! implicit none + ! ! Arguments + ! character(len=*), intent(in) :: message + ! ! Internals + ! character(STRMAX) :: errmsg + + ! open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + ! write(FRAGGLE_LOG_UNIT, *) trim(adjustl(message)) + ! close(FRAGGLE_LOG_UNIT) + + ! return + ! 667 continue + ! write(*,*) "Error writing Fraggle message to log file: " // trim(adjustl(errmsg)) + ! end subroutine fraggle_io_log_one_message module subroutine fraggle_io_log_pl(pl, param) @@ -228,28 +228,28 @@ module subroutine fraggle_io_log_regime(colliders, frag) end subroutine fraggle_io_log_regime - module subroutine fraggle_io_log_start(param) - !! author: David A. Minton - !! - !! Checks to see if the Fraggle log file needs to be replaced if this is a new run, or appended if this is a restarted run - implicit none - ! Arguments - class(swiftest_parameters), intent(in) :: param - ! Internals - character(STRMAX) :: errmsg - logical :: fileExists - - inquire(file=FRAGGLE_LOG_OUT, exist=fileExists) - if (.not.param%lrestart .or. .not.fileExists) then - open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status="REPLACE", err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) "Fraggle logfile" - end if - close(FRAGGLE_LOG_UNIT) - - return - - 667 continue - write(*,*) "Error writing Fraggle log file: " // trim(adjustl(errmsg)) - end subroutine fraggle_io_log_start + ! module subroutine fraggle_io_log_start(param) + ! !! author: David A. Minton + ! !! + ! !! Checks to see if the Fraggle log file needs to be replaced if this is a new run, or appended if this is a restarted run + ! implicit none + ! ! Arguments + ! class(swiftest_parameters), intent(in) :: param + ! ! Internals + ! character(STRMAX) :: errmsg + ! logical :: fileExists + + ! inquire(file=FRAGGLE_LOG_OUT, exist=fileExists) + ! if (.not.param%lrestart .or. .not.fileExists) then + ! open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status="REPLACE", err = 667, iomsg = errmsg) + ! write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) "Fraggle logfile" + ! end if + ! close(FRAGGLE_LOG_UNIT) + + ! return + + ! 667 continue + ! write(*,*) "Error writing Fraggle log file: " // trim(adjustl(errmsg)) + ! end subroutine fraggle_io_log_start end submodule s_fraggle_io \ No newline at end of file diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index df9265ae7..f53888df5 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -182,7 +182,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb Mlr = Mtot Mslr = 0.0_DP Qloss = 0.0_DP - call fraggle_io_log_one_message("Fragments would have mass below the minimum. Converting this collision into a merger.") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fragments would have mass below the minimum. Converting this collision into a merger.") else if( Vimp < Vescp) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime diff --git a/src/io/io.f90 b/src/io/io.f90 index ef3f682d2..b0f30d802 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -101,7 +101,6 @@ module subroutine io_dump_param(self, param_file_name) class(swiftest_parameters),intent(in) :: self !! Output collection of parameters character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) ! Internals - integer(I4B), parameter :: LUN = 7 !! Unit number of output file character(STRMAX) :: errmsg !! Error message in UDIO procedure integer(I4B) :: ierr @@ -160,7 +159,6 @@ module subroutine io_dump_particle_info_base(self, param, idx) ! Internals logical, save :: lfirst = .true. - integer(I4B), parameter :: LUN = 22 integer(I4B) :: i character(STRMAX) :: errmsg @@ -228,7 +226,6 @@ module subroutine io_dump_base(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: ierr !! Error code - integer(I4B),parameter :: LUN = 7 !! Unit number for dump file integer(I4B) :: iu = LUN character(len=:), allocatable :: dump_file_name character(STRMAX) :: errmsg @@ -395,7 +392,6 @@ module function io_get_old_t_final_system(self, param) result(old_t_final) ! Internals class(swiftest_nbody_system), allocatable :: tmpsys class(swiftest_parameters), allocatable :: tmpparam - integer(I4B), parameter :: LUN = 76 integer(I4B) :: ierr, iu = LUN character(len=STRMAX) :: errmsg @@ -411,6 +407,7 @@ module function io_get_old_t_final_system(self, param) result(old_t_final) end do if (is_iostat_end(ierr)) then old_t_final = tmpparam%t + close(iu) return end if @@ -468,6 +465,53 @@ module function io_get_token(buffer, ifirst, ilast, ierr) result(token) return end function io_get_token + module subroutine io_log_one_message(file, message) + !! author: David A. Minton + !! + !! Writes a single message to a log file + implicit none + ! Arguments + character(len=*), intent(in) :: file !! Name of file to log + character(len=*), intent(in) :: message + ! Internals + character(STRMAX) :: errmsg + + open(unit=LUN, file=file, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(LUN, *) trim(adjustl(message)) + close(LUN) + + return + 667 continue + write(*,*) "Error writing message to log file: " // trim(adjustl(errmsg)) + end subroutine io_log_one_message + + + module subroutine io_log_start(param, file, header) + !! author: David A. Minton + !! + !! Checks to see if a log file needs to be created if this is a new run, or appended if this is a restarted run + implicit none + ! Arguments + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + character(len=*), intent(in) :: file !! Name of file to log + character(len=*), intent(in) :: header !! Header to print at top of log file + ! Internals + character(STRMAX) :: errmsg + logical :: fileExists + + inquire(file=file, exist=fileExists) + if (.not.param%lrestart .or. .not.fileExists) then + open(unit=LUN, file=file, status="REPLACE", err = 667, iomsg = errmsg) + write(LUN, *, err = 667, iomsg = errmsg) trim(adjustl(header)) + end if + close(LUN) + + return + + 667 continue + write(*,*) "Error writing log file: " // trim(adjustl(errmsg)) + end subroutine io_log_start + module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott @@ -766,12 +810,22 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) end associate select case(trim(adjustl(param%interaction_loops))) - case("ADAPTIVE", "FLAT", "TRIANGULAR") + case("ADAPTIVE") + param%ladaptive_interactions = .true. + param%lflatten_interactions = .true. + case("TRIANGULAR") + param%ladaptive_interactions = .false. + param%lflatten_interactions = .false. + case("FLAT") + param%ladaptive_interactions = .false. + param%lflatten_interactions = .true. case default write(*,*) "Unknown value for parameter INTERACTION_LOOPS: -> ",trim(adjustl(param%interaction_loops)) write(*,*) "Must be one of the following: TRIANGULAR, FLAT, or ADAPTIVE" write(*,*) "Using default value of ADAPTIVE" param%interaction_loops = "ADAPTIVE" + param%ladaptive_interactions = .true. + param%lflatten_interactions = .true. end select iostat = 0 @@ -1092,7 +1146,6 @@ module subroutine io_read_in_body(self, param) class(swiftest_body), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B), parameter :: LUN = 7 !! Unit number of input file integer(I4B) :: iu = LUN integer(I4B) :: i, nbody logical :: is_ascii, is_pl @@ -1157,7 +1210,6 @@ module subroutine io_read_in_cb(self, param) class(swiftest_cb), intent(inout) :: self class(swiftest_parameters), intent(inout) :: param ! Internals - integer(I4B), parameter :: LUN = 7 !! Unit number of input file integer(I4B) :: iu = LUN character(len=STRMAX) :: errmsg integer(I4B) :: ierr, idold @@ -1225,7 +1277,6 @@ function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & integer(I4B) :: ierr ! Internals logical , save :: lfirst = .true. - integer(I4B), parameter :: lun = 30 integer(I4B), save :: iu = lun if (lfirst) then @@ -1561,7 +1612,6 @@ module subroutine io_read_in_param(self, param_file_name) class(swiftest_parameters),intent(inout) :: self !! Current run configuration parameters character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) ! Internals - integer(I4B), parameter :: LUN = 7 !! Unit number of input file integer(I4B) :: ierr = 0 !! Input error code character(STRMAX) :: errmsg !! Error message in UDIO procedure @@ -1619,8 +1669,7 @@ module subroutine io_read_particle_info_system(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B), parameter :: LUN = 22 - integer(I4B) :: i, id, idx + integer(I4B) :: i, id, idx logical :: lmatch character(STRMAX) :: errmsg type(swiftest_particle_info), allocatable :: tmpinfo @@ -1715,7 +1764,6 @@ module subroutine io_write_discard(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B), parameter :: LUN = 40 integer(I4B) :: i logical, save :: lfirst = .true. real(DP), dimension(:,:), allocatable :: vh @@ -1800,7 +1848,6 @@ module subroutine io_write_encounter(self, pl, encbody, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals logical , save :: lfirst = .true. - integer(I4B), parameter :: LUN = 30 integer(I4B) :: k, ierr character(len=STRMAX) :: errmsg diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 0e90e73ee..2626c6d34 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -17,9 +17,7 @@ module subroutine kick_getacch_int_pl(self, param) type(interaction_timer), save :: itimer logical, save :: lfirst - if (lfirst) then - call itimer%reset(param) - end if + if (lfirst) call itimer%reset(param) if (param%lflatten_interactions) then call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) else diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 63440383f..7fefe652b 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -113,12 +113,6 @@ module subroutine fraggle_io_log_generate(frag) class(fraggle_fragments), intent(in) :: frag end subroutine fraggle_io_log_generate - module subroutine fraggle_io_log_one_message(message) - implicit none - character(len=*), intent(in) :: message - character(STRMAX) :: errmsg - end subroutine fraggle_io_log_one_message - module subroutine fraggle_io_log_pl(pl, param) implicit none class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object (only the new bodies generated in a collision) @@ -131,11 +125,6 @@ module subroutine fraggle_io_log_regime(colliders, frag) class(fraggle_fragments), intent(in) :: frag end subroutine fraggle_io_log_regime - module subroutine fraggle_io_log_start(param) - implicit none - class(swiftest_parameters), intent(in) :: param - end subroutine fraggle_io_log_start - !> The following interfaces are placeholders intended to satisfy the required abstract methods given by the parent class module subroutine fraggle_placeholder_accel(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 89c5ef238..46cf56369 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -123,6 +123,9 @@ module swiftest_classes real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units character(STRMAX) :: energy_out = "" !! Name of output energy and momentum report file character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" + ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS + logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops + logical :: ladaptive_interactions = .false. !! Adaptive interaction loop is turned on ! Logical flags to turn on or off various features of the code logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually) @@ -133,7 +136,6 @@ module swiftest_classes logical :: loblatecb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input) logical :: lrotation = .false. !! Include rotation states of big bodies logical :: ltides = .false. !! Include tidal dissipation - logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory) ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values) real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy @@ -707,6 +709,19 @@ module function io_get_token(buffer, ifirst, ilast, ierr) result(token) character(len=:), allocatable :: token !! Returned token string end function io_get_token + module subroutine io_log_one_message(file, message) + implicit none + character(len=*), intent(in) :: file !! Name of file to log + character(len=*), intent(in) :: message + end subroutine io_log_one_message + + module subroutine io_log_start(param, file, header) + implicit none + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + character(len=*), intent(in) :: file !! Name of file to log + character(len=*), intent(in) :: header !! Header to print at top of log file + end subroutine io_log_start + module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none class(swiftest_parameters), intent(inout) :: self !! Collection of parameters diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 1f9c6028c..b7fe1a0db 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -113,9 +113,9 @@ module swiftest_globals !> Standard file names integer(I4B), parameter :: NDUMPFILES = 2 - character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ] - character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ] - character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ] + character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ] + character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ] + character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ] character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.in', 'dump_param2.in'] !> Default file names that can be changed by the user in the parameters file @@ -126,6 +126,7 @@ module swiftest_globals integer(I4B), parameter :: BINUNIT = 20 !! File unit number for the binary output file character(*), parameter :: PARTICLE_OUTFILE = 'particle.dat' integer(I4B), parameter :: PARTICLEUNIT = 44 !! File unit number for the binary particle info output file + integer(I4B), parameter :: LUN = 42 !! File unit number for files that are opened and closed within a single subroutine call, and therefore should not collide !> Miscellaneous constants: integer(I4B), parameter :: NDIM = 3 !! Number of dimensions in our reality diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 76600660d..c353328d8 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -27,11 +27,25 @@ module walltime_classes integer(I4B) :: step_counter integer(I8B) :: count_previous character(len=NAMELEN) :: current_style + logical :: lflatten_interaction_old contains - + procedure :: reset => walltime_interaction_reset !! Resets the interaction loop timer, and saves the current value of the array flatten parameter end type interaction_timer interface + module subroutine walltime_interaction_reset(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(interaction_timer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine walltime_interaction_reset + + module subroutine walltime_interaction_io_log_start(param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(swiftest_parameters), intent(in) :: param + end subroutine walltime_interaction_io_log_start + module subroutine walltime_finish(self, nsubsteps, message, param) use swiftest_classes, only : swiftest_parameters implicit none @@ -58,6 +72,23 @@ end subroutine walltime_start contains + module subroutine walltime_interaction_reset(self, param) + !! author: David A. Minton + !! + !! Resets the interaction loop timer, and saves the current value of the array flatten parameter + use swiftest_classes, only : swiftest_parameters + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + self%lflatten_interaction_old = param%lflatten_interactions + call walltime_reset(self, param) + + return + end subroutine walltime_interaction_reset + + module subroutine walltime_finish(self, nsubsteps, message, param) !! author: David A. Minton !! diff --git a/src/rmvs/rmvs_io.f90 b/src/rmvs/rmvs_io.f90 index b85ce2202..51f077852 100644 --- a/src/rmvs/rmvs_io.f90 +++ b/src/rmvs/rmvs_io.f90 @@ -19,7 +19,6 @@ module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, character(*), intent(in) :: enc_out ! Internals logical , save :: lfirst = .true. - integer(I4B), parameter :: LUN = 30 integer(I4B) :: ierr if (enc_out == "") return diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 864eaa723..baa51485a 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -28,7 +28,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag) message = "Supercatastrophic disruption between" end select call symba_collision_collider_message(system%pl, colliders%idx, message) - call fraggle_io_log_one_message(message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter call frag%set_mass_dist(colliders, param) @@ -37,7 +37,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag) call frag%generate_fragments(colliders, system, param, lfailure) if (lfailure) then - call fraggle_io_log_one_message("No fragment solution found, so treat as a pure hit-and-run") + call io_log_one_message(FRAGGLE_LOG_OUT, "No fragment solution found, so treat as a pure hit-and-run") status = ACTIVE nfrag = 0 select type(pl => system%pl) @@ -50,7 +50,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag) ! Populate the list of new bodies nfrag = frag%nbody write(message, *) nfrag - call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments") + call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") select case(frag%regime) case(COLLRESOLVE_REGIME_DISRUPTION) status = DISRUPTION @@ -87,7 +87,7 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r message = "Hit and run between" call symba_collision_collider_message(system%pl, colliders%idx, message) - call fraggle_io_log_one_message(trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, trim(adjustl(message))) if (colliders%mass(1) > colliders%mass(2)) then jtarg = 1 @@ -98,7 +98,7 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r end if if (frag%mass_dist(2) > 0.9_DP * colliders%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched - call fraggle_io_log_one_message("Pure hit and run. No new fragments generated.") + call io_log_one_message(FRAGGLE_LOG_OUT, "Pure hit and run. No new fragments generated.") nfrag = 0 lpure = .true. else ! Imperfect hit and run, so we'll keep the largest body and destroy the other @@ -109,12 +109,12 @@ module function symba_collision_casehitandrun(system, param, colliders, frag) r call frag%generate_fragments(colliders, system, param, lpure) if (lpure) then - call fraggle_io_log_one_message("Should have been a pure hit and run instead") + call io_log_one_message(FRAGGLE_LOG_OUT, "Should have been a pure hit and run instead") nfrag = 0 else nfrag = frag%nbody write(message, *) nfrag - call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments") + call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") end if end if if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them @@ -164,7 +164,7 @@ module function symba_collision_casemerge(system, param, colliders, frag) resul message = "Merging" call symba_collision_collider_message(system%pl, colliders%idx, message) - call fraggle_io_log_one_message(message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) select type(pl => system%pl) class is (symba_pl) @@ -358,7 +358,7 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & // " at t = " // trim(adjustl(timestr)) - call fraggle_io_log_one_message(message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) end if end if end do @@ -984,10 +984,10 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir do write(timestr,*) t - call fraggle_io_log_one_message("") - call fraggle_io_log_one_message("***********************************************************************************************************************") - call fraggle_io_log_one_message("Collision between massive bodies detected at time t = " // trim(adjustl(timestr))) - call fraggle_io_log_one_message("***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // trim(adjustl(timestr))) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") allocate(tmp_param, source=param) tmp_param%t = t if (param%lfragmentation) then diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 704d11f85..592b146b8 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -37,11 +37,11 @@ subroutine symba_discard_cb_pl(pl, system, param) write(idstr, *) pl%id(i) write(timestr, *) param%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) - call fraggle_io_log_one_message("") - call fraggle_io_log_one_message("***********************************************************************************************************************") - call fraggle_io_log_one_message(message) - call fraggle_io_log_one_message("***********************************************************************************************************************") - call fraggle_io_log_one_message("") + call io_log_one_message(FRAGGLE_LOG_OUT, "") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, message) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "") call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. @@ -50,11 +50,11 @@ subroutine symba_discard_cb_pl(pl, system, param) write(idstr, *) pl%id(i) write(timestr, *) param%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) - call fraggle_io_log_one_message("") - call fraggle_io_log_one_message("***********************************************************************************************************************") - call fraggle_io_log_one_message(message) - call fraggle_io_log_one_message("***********************************************************************************************************************") - call fraggle_io_log_one_message("") + call io_log_one_message(FRAGGLE_LOG_OUT, "") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, message) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "") call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) @@ -67,11 +67,11 @@ subroutine symba_discard_cb_pl(pl, system, param) write(idstr, *) pl%id(i) write(timestr, *) param%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) - call fraggle_io_log_one_message("") - call fraggle_io_log_one_message("***********************************************************************************************************************") - call fraggle_io_log_one_message(message) - call fraggle_io_log_one_message("***********************************************************************************************************************") - call fraggle_io_log_one_message("") + call io_log_one_message(FRAGGLE_LOG_OUT, "") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, message) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "") call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) end if end if diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 2e5e87941..3f481d9b4 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -96,7 +96,7 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms end if ! All reporting of collision information in SyMBA (including mergers) is now recorded in the Fraggle logfile - call fraggle_io_log_start(param) + call io_log_start(param, FRAGGLE_LOG_OUT, "Fraggle logfile") ! Call the base method (which also prints the contents to screen) call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) @@ -155,7 +155,6 @@ module subroutine symba_io_write_discard(self, param) class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B), parameter :: LUN = 40 integer(I4B) :: iadd, isub, j, nsub, nadd logical, save :: lfirst = .true. real(DP), dimension(:,:), allocatable :: vh From e6c1b5ea552c75ac9f5bdc7ce19426f6e7a89283 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 14:42:39 -0400 Subject: [PATCH 06/15] Completed refactoring to make a universal temporary file unit number LUN rather than having separate random ones in a whole bunch of different places --- src/fraggle/fraggle_io.f90 | 206 +++++++++++++++---------------- src/io/io.f90 | 1 + src/modules/fraggle_classes.f90 | 1 - src/modules/walltime_classes.f90 | 1 + 4 files changed, 105 insertions(+), 104 deletions(-) diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 9d7af73f1..410592782 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -15,45 +15,45 @@ module subroutine fraggle_io_log_generate(frag) character(STRMAX) :: errmsg character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" - open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" - write(FRAGGLE_LOG_UNIT, *) " Fraggle fragment generation results" - write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" - write(FRAGGLE_LOG_UNIT, "(' dL_tot should be very small' )") - write(FRAGGLE_LOG_UNIT,fmtlabel) ' dL_tot |', (.mag.(frag%Ltot_after(:) - frag%Ltot_before(:))) / (.mag.frag%Ltot_before(:)) - write(FRAGGLE_LOG_UNIT, "(' dE_tot should be negative and equal to Qloss' )") - write(FRAGGLE_LOG_UNIT,fmtlabel) ' dE_tot |', (frag%Etot_after - frag%Etot_before) / abs(frag%Etot_before) - write(FRAGGLE_LOG_UNIT,fmtlabel) ' Qloss |', -frag%Qloss / abs(frag%Etot_before) - write(FRAGGLE_LOG_UNIT,fmtlabel) ' dE - Qloss |', (frag%Etot_after - frag%Etot_before + frag%Qloss) / abs(frag%Etot_before) - write(FRAGGLE_LOG_UNIT, "(' -------------------------------------------------------------------------------------')") - write(FRAGGLE_LOG_UNIT, *) "Individual fragment values (collisional system natural units)" - write(FRAGGLE_LOG_UNIT, *) "mass" + open(unit=LUN, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(LUN, *, err = 667, iomsg = errmsg) + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, *) " Fraggle fragment generation results" + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, "(' dL_tot should be very small' )") + write(LUN,fmtlabel) ' dL_tot |', (.mag.(frag%Ltot_after(:) - frag%Ltot_before(:))) / (.mag.frag%Ltot_before(:)) + write(LUN, "(' dE_tot should be negative and equal to Qloss' )") + write(LUN,fmtlabel) ' dE_tot |', (frag%Etot_after - frag%Etot_before) / abs(frag%Etot_before) + write(LUN,fmtlabel) ' Qloss |', -frag%Qloss / abs(frag%Etot_before) + write(LUN,fmtlabel) ' dE - Qloss |', (frag%Etot_after - frag%Etot_before + frag%Qloss) / abs(frag%Etot_before) + write(LUN, "(' -------------------------------------------------------------------------------------')") + write(LUN, *) "Individual fragment values (collisional system natural units)" + write(LUN, *) "mass" do i = 1, frag%nbody - write(FRAGGLE_LOG_UNIT, *) i, frag%mass(i) + write(LUN, *) i, frag%mass(i) end do - write(FRAGGLE_LOG_UNIT, *) "x_coll" + write(LUN, *) "x_coll" do i = 1, frag%nbody - write(FRAGGLE_LOG_UNIT, *) i, frag%x_coll(:,i) + write(LUN, *) i, frag%x_coll(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "v_coll" + write(LUN, *) "v_coll" do i = 1, frag%nbody - write(FRAGGLE_LOG_UNIT, *) i, frag%v_coll(:,i) + write(LUN, *) i, frag%v_coll(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "xb" + write(LUN, *) "xb" do i = 1, frag%nbody - write(FRAGGLE_LOG_UNIT, *) i, frag%xb(:,i) + write(LUN, *) i, frag%xb(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "vb" + write(LUN, *) "vb" do i = 1, frag%nbody - write(FRAGGLE_LOG_UNIT, *) i, frag%vb(:,i) + write(LUN, *) i, frag%vb(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "rot" + write(LUN, *) "rot" do i = 1, frag%nbody - write(FRAGGLE_LOG_UNIT, *) i, frag%rot(:,i) + write(LUN, *) i, frag%rot(:,i) end do - close(FRAGGLE_LOG_UNIT) + close(LUN) return 667 continue @@ -71,9 +71,9 @@ end subroutine fraggle_io_log_generate ! ! Internals ! character(STRMAX) :: errmsg - ! open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - ! write(FRAGGLE_LOG_UNIT, *) trim(adjustl(message)) - ! close(FRAGGLE_LOG_UNIT) + ! open(unit=LUN, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + ! write(LUN, *) trim(adjustl(message)) + ! close(LUN) ! return ! 667 continue @@ -93,68 +93,68 @@ module subroutine fraggle_io_log_pl(pl, param) integer(I4B) :: i character(STRMAX) :: errmsg - open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) + open(unit=LUN, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(LUN, *, err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" - write(FRAGGLE_LOG_UNIT, *) " Fraggle fragment final body properties" - write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" - write(FRAGGLE_LOG_UNIT, *) "id, name" + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, *) " Fraggle fragment final body properties" + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, *) "id, name" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%id(i), pl%info(i)%name + write(LUN, *) i, pl%id(i), pl%info(i)%name end do - write(FRAGGLE_LOG_UNIT, *) "mass, Gmass" + write(LUN, *) "mass, Gmass" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%mass(i), pl%Gmass(i) + write(LUN, *) i, pl%mass(i), pl%Gmass(i) end do - write(FRAGGLE_LOG_UNIT, *) "radius" + write(LUN, *) "radius" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%radius(i) + write(LUN, *) i, pl%radius(i) end do - write(FRAGGLE_LOG_UNIT, *) "xb" + write(LUN, *) "xb" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%xb(:,i) + write(LUN, *) i, pl%xb(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "vb" + write(LUN, *) "vb" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%vb(:,i) + write(LUN, *) i, pl%vb(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "xh" + write(LUN, *) "xh" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%xh(:,i) + write(LUN, *) i, pl%xh(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "vh" + write(LUN, *) "vh" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%vh(:,i) + write(LUN, *) i, pl%vh(:,i) end do if (param%lrotation) then - write(FRAGGLE_LOG_UNIT, *) "rot" + write(LUN, *) "rot" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%rot(:,i) + write(LUN, *) i, pl%rot(:,i) end do - write(FRAGGLE_LOG_UNIT, *) "Ip" + write(LUN, *) "Ip" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%Ip(:,i) + write(LUN, *) i, pl%Ip(:,i) end do end if if (param%ltides) then - write(FRAGGLE_LOG_UNIT, *) "Q" + write(LUN, *) "Q" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%Q(i) + write(LUN, *) i, pl%Q(i) end do - write(FRAGGLE_LOG_UNIT, *) "k2" + write(LUN, *) "k2" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%k2(i) + write(LUN, *) i, pl%k2(i) end do - write(FRAGGLE_LOG_UNIT, *) "tlag" + write(LUN, *) "tlag" do i = 1, pl%nbody - write(FRAGGLE_LOG_UNIT, *) i, pl%tlag(i) + write(LUN, *) i, pl%tlag(i) end do end if - close(FRAGGLE_LOG_UNIT) + close(LUN) return 667 continue @@ -173,54 +173,54 @@ module subroutine fraggle_io_log_regime(colliders, frag) ! Internals character(STRMAX) :: errmsg - open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) - write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" - write(FRAGGLE_LOG_UNIT, *) " Fraggle collisional regime determination results" - write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" - write(FRAGGLE_LOG_UNIT, *) "----------------------- Collider information -----------------------" - write(FRAGGLE_LOG_UNIT, *) "True number of colliders : ",colliders%ncoll - write(FRAGGLE_LOG_UNIT, *) "Index list of true colliders : ",colliders%idx(1:colliders%ncoll) - write(FRAGGLE_LOG_UNIT, *) "-------------------- Two-body equialent values ---------------------" - write(FRAGGLE_LOG_UNIT, *) "mass1 : ",colliders%mass(1) - write(FRAGGLE_LOG_UNIT, *) "radius1 : ",colliders%radius(1) - write(FRAGGLE_LOG_UNIT, *) "xb1 : ",colliders%xb(:,1) - write(FRAGGLE_LOG_UNIT, *) "vb1 : ",colliders%vb(:,1) - write(FRAGGLE_LOG_UNIT, *) "rot1 : ",colliders%rot(:,1) - write(FRAGGLE_LOG_UNIT, *) "Ip1 : ",colliders%Ip(:,1) - write(FRAGGLE_LOG_UNIT, *) "L_spin1 : ",colliders%L_spin(:,1) - write(FRAGGLE_LOG_UNIT, *) "L_orbit1 : ",colliders%L_orbit(:,1) - write(FRAGGLE_LOG_UNIT, *) "mass2 : ",colliders%mass(2) - write(FRAGGLE_LOG_UNIT, *) "radius2 : ",colliders%radius(2) - write(FRAGGLE_LOG_UNIT, *) "xb2 : ",colliders%xb(:,2) - write(FRAGGLE_LOG_UNIT, *) "vb2 : ",colliders%vb(:,2) - write(FRAGGLE_LOG_UNIT, *) "rot2 : ",colliders%rot(:,2) - write(FRAGGLE_LOG_UNIT, *) "Ip2 : ",colliders%Ip(:,2) - write(FRAGGLE_LOG_UNIT, *) "L_spin2 : ",colliders%L_spin(:,2) - write(FRAGGLE_LOG_UNIT, *) "L_orbit2 : ",colliders%L_orbit(:,2) - write(FRAGGLE_LOG_UNIT, *) "------------------------------ Regime -----------------------------" + open(unit=LUN, file=FRAGGLE_LOG_OUT, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + write(LUN, *, err = 667, iomsg = errmsg) + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, *) " Fraggle collisional regime determination results" + write(LUN, *) "--------------------------------------------------------------------" + write(LUN, *) "----------------------- Collider information -----------------------" + write(LUN, *) "True number of colliders : ",colliders%ncoll + write(LUN, *) "Index list of true colliders : ",colliders%idx(1:colliders%ncoll) + write(LUN, *) "-------------------- Two-body equialent values ---------------------" + write(LUN, *) "mass1 : ",colliders%mass(1) + write(LUN, *) "radius1 : ",colliders%radius(1) + write(LUN, *) "xb1 : ",colliders%xb(:,1) + write(LUN, *) "vb1 : ",colliders%vb(:,1) + write(LUN, *) "rot1 : ",colliders%rot(:,1) + write(LUN, *) "Ip1 : ",colliders%Ip(:,1) + write(LUN, *) "L_spin1 : ",colliders%L_spin(:,1) + write(LUN, *) "L_orbit1 : ",colliders%L_orbit(:,1) + write(LUN, *) "mass2 : ",colliders%mass(2) + write(LUN, *) "radius2 : ",colliders%radius(2) + write(LUN, *) "xb2 : ",colliders%xb(:,2) + write(LUN, *) "vb2 : ",colliders%vb(:,2) + write(LUN, *) "rot2 : ",colliders%rot(:,2) + write(LUN, *) "Ip2 : ",colliders%Ip(:,2) + write(LUN, *) "L_spin2 : ",colliders%L_spin(:,2) + write(LUN, *) "L_orbit2 : ",colliders%L_orbit(:,2) + write(LUN, *) "------------------------------ Regime -----------------------------" select case(frag%regime) case(COLLRESOLVE_REGIME_MERGE) - write(FRAGGLE_LOG_UNIT, *) "Merge" + write(LUN, *) "Merge" case(COLLRESOLVE_REGIME_DISRUPTION) - write(FRAGGLE_LOG_UNIT, *) "Disruption" + write(LUN, *) "Disruption" case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - write(FRAGGLE_LOG_UNIT, *) "Supercatastrophic disruption" + write(LUN, *) "Supercatastrophic disruption" case(COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - write(FRAGGLE_LOG_UNIT, *) "Graze and merge" + write(LUN, *) "Graze and merge" case(COLLRESOLVE_REGIME_HIT_AND_RUN) - write(FRAGGLE_LOG_UNIT, *) "Hit and run" + write(LUN, *) "Hit and run" end select - write(FRAGGLE_LOG_UNIT, *) "----------------------- Fragment information ----------------------" - write(FRAGGLE_LOG_UNIT, *) "Total mass of fragments : ", frag%mtot - write(FRAGGLE_LOG_UNIT, *) "Largest fragment mass : ", frag%mass_dist(1) - write(FRAGGLE_LOG_UNIT, *) "Second-largest fragment mass : ", frag%mass_dist(2) - write(FRAGGLE_LOG_UNIT, *) "Remaining fragment mass : ", frag%mass_dist(3) - write(FRAGGLE_LOG_UNIT, *) "Center of mass position : ", frag%xbcom(:) - write(FRAGGLE_LOG_UNIT, *) "Center of mass velocity : ", frag%vbcom(:) - write(FRAGGLE_LOG_UNIT, *) "Energy loss : ", frag%Qloss - write(FRAGGLE_LOG_UNIT, *) "--------------------------------------------------------------------" - close(FRAGGLE_LOG_UNIT) + write(LUN, *) "----------------------- Fragment information ----------------------" + write(LUN, *) "Total mass of fragments : ", frag%mtot + write(LUN, *) "Largest fragment mass : ", frag%mass_dist(1) + write(LUN, *) "Second-largest fragment mass : ", frag%mass_dist(2) + write(LUN, *) "Remaining fragment mass : ", frag%mass_dist(3) + write(LUN, *) "Center of mass position : ", frag%xbcom(:) + write(LUN, *) "Center of mass velocity : ", frag%vbcom(:) + write(LUN, *) "Energy loss : ", frag%Qloss + write(LUN, *) "--------------------------------------------------------------------" + close(LUN) return 667 continue @@ -241,10 +241,10 @@ end subroutine fraggle_io_log_regime ! inquire(file=FRAGGLE_LOG_OUT, exist=fileExists) ! if (.not.param%lrestart .or. .not.fileExists) then - ! open(unit=FRAGGLE_LOG_UNIT, file=FRAGGLE_LOG_OUT, status="REPLACE", err = 667, iomsg = errmsg) - ! write(FRAGGLE_LOG_UNIT, *, err = 667, iomsg = errmsg) "Fraggle logfile" + ! open(unit=LUN, file=FRAGGLE_LOG_OUT, status="REPLACE", err = 667, iomsg = errmsg) + ! write(LUN, *, err = 667, iomsg = errmsg) "Fraggle logfile" ! end if - ! close(FRAGGLE_LOG_UNIT) + ! close(LUN) ! return diff --git a/src/io/io.f90 b/src/io/io.f90 index b0f30d802..4915fba27 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -813,6 +813,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case("ADAPTIVE") param%ladaptive_interactions = .true. param%lflatten_interactions = .true. + call io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile") case("TRIANGULAR") param%ladaptive_interactions = .false. param%lflatten_interactions = .false. diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 7fefe652b..cd648c04b 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -9,7 +9,6 @@ module fraggle_classes integer(I4B), parameter :: FRAGGLE_NMASS_DIST = 3 !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder) character(len=*), parameter :: FRAGGLE_LOG_OUT = "fraggle.log" !! Name of log file for Fraggle diagnostic information - integer(I4B), parameter :: FRAGGLE_LOG_UNIT = 76 !! Unit number for Fraggle log file !******************************************************************************************************************************** ! fraggle_colliders class definitions and method interfaces diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index c353328d8..e1f05dc17 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -8,6 +8,7 @@ module walltime_classes public integer(I4B) :: INTERACTION_TIMER_CADENCE = 1000 !! Minimum number of steps to wait before timing an interaction loop in ADAPTIVE mode + character(len=*), parameter :: INTERACTION_TIMER_LOG_OUT = "interaction_timer.log" !! Name of log file for recording results of interaction loop timing type :: walltimer integer(I8B) :: count_rate !! Rate at wich the clock ticks From 473a6482e8a3d28765e66b40bcbd417feda6454c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 14:51:26 -0400 Subject: [PATCH 07/15] Added counter reset to interaction timer --- src/modules/walltime_classes.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index e1f05dc17..97f6fb520 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -27,7 +27,6 @@ module walltime_classes integer(I8B) :: max_interactions = huge(1_I8B) integer(I4B) :: step_counter integer(I8B) :: count_previous - character(len=NAMELEN) :: current_style logical :: lflatten_interaction_old contains procedure :: reset => walltime_interaction_reset !! Resets the interaction loop timer, and saves the current value of the array flatten parameter @@ -84,6 +83,7 @@ module subroutine walltime_interaction_reset(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters self%lflatten_interaction_old = param%lflatten_interactions + self%step_counter = 0 call walltime_reset(self, param) return From 78d1ca4168a4d2e59811899404fd67446b61fc0d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 18:08:28 -0400 Subject: [PATCH 08/15] Added adaptive interaction loop timing. Right now all the diagnostic logging is done in the loop. Eventually I will consolidate this inside the walltimer methods. --- Makefile.Defines | 10 +- src/io/io.f90 | 17 +-- src/kick/kick.f90 | 42 ++++++- src/modules/swiftest_classes.f90 | 3 +- src/modules/walltime_classes.f90 | 184 +++++++++++++++++++++++++------ src/symba/symba_io.f90 | 7 +- src/symba/symba_kick.f90 | 52 +++++++++ src/user/user_getacch.f90 | 2 +- 8 files changed, 268 insertions(+), 49 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index 9e06d56ba..16ce3afc3 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -67,15 +67,15 @@ GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporari GPRODUCTION = -O2 -ffree-line-length-none $(GPAR) -#FFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) -#FFASTFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) -FFLAGS = $(IPRODUCTION) $(STRICTREAL) -FFASTFLAGS = $(IPRODUCTION) -fp-model fast +#FFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) +#FFASTFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) +FFLAGS = $(IPRODUCTION) $(STRICTREAL) #$(ADVIXE_FLAGS) +FFASTFLAGS = $(IPRODUCTION) -fp-model fast #$(ADVIXE_FLAGS) FORTRAN = ifort AR = xiar #FORTRAN = gfortran -#FFLAGS = $(GDEBUG) $(GMEM) $(GPAR) +#FFLAGS = $(GDEBUG) # $(GMEM) $(GPAR) #FFLAGS = $(GPRODUCTION) -g -fbacktrace #-fcheck=all #-Wall AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only # this is done explicitly as needed in the Makefile diff --git a/src/io/io.f90 b/src/io/io.f90 index 4915fba27..cf9c6a8a9 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -465,6 +465,7 @@ module function io_get_token(buffer, ifirst, ilast, ierr) result(token) return end function io_get_token + module subroutine io_log_one_message(file, message) !! author: David A. Minton !! @@ -545,7 +546,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ! Parse the file line by line, extracting tokens then matching them up with known parameters if possible associate(param => self) - rewind(unit) + open(unit = unit, file = param%param_file_name, status = 'old', err = 667, iomsg = iomsg) do read(unit = unit, fmt = linefmt, end = 1, err = 667, iomsg = iomsg) line line_trim = trim(adjustl(line)) @@ -705,6 +706,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) end if end do 1 continue + close(unit) iostat = 0 ! Do basic sanity checks on the input values @@ -814,6 +816,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param%ladaptive_interactions = .true. param%lflatten_interactions = .true. call io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile") + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") case("TRIANGULAR") param%ladaptive_interactions = .false. param%lflatten_interactions = .false. @@ -836,8 +839,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) end associate - 667 continue return + 667 continue + write(*,*) "Error reading param file: ", trim(adjustl(iomsg)) end subroutine io_param_reader @@ -1613,19 +1617,18 @@ module subroutine io_read_in_param(self, param_file_name) class(swiftest_parameters),intent(inout) :: self !! Current run configuration parameters character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) ! Internals - integer(I4B) :: ierr = 0 !! Input error code - character(STRMAX) :: errmsg !! Error message in UDIO procedure + integer(I4B) :: ierr = 0 !! Input error code + character(STRMAX) :: errmsg !! Error message in UDIO procedure ! Read in name of parameter file write(*, *) 'Parameter input file is ', trim(adjustl(param_file_name)) write(*, *) ' ' - 100 format(A) - open(unit = LUN, file = param_file_name, status = 'old', iostat = ierr, err = 667, iomsg = errmsg) + self%param_file_name = param_file_name !! todo: Currently this procedure does not work in user-defined derived-type input mode !! as the newline characters are ignored in the input file when compiled in ifort. - !read(LUN,'(DT)', iostat= ierr, iomsg = errmsg) param + !read(LUN,'(DT)', iostat= ierr, iomsg = errmsg) self call self%reader(LUN, iotype= "none", v_list = [self%integrator], iostat = ierr, iomsg = errmsg) if (ierr == 0) return diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 2626c6d34..edc054caf 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -16,14 +16,54 @@ module subroutine kick_getacch_int_pl(self, param) ! Internals type(interaction_timer), save :: itimer logical, save :: lfirst + character(len=STRMAX) :: tstr, nstr, cstr, mstr, lstyle + character(len=1) :: schar + + if (param%ladaptive_interactions) then + if (lfirst) then + call itimer%time_this_loop(param, self, self%nplpl) + lfirst = .false. + else + if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl) + end if + end if + + if (itimer%is_on) then + write(tstr,*) param%t + write(schar,'(I1)') itimer%stage + if (itimer%stage == 1) then + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "kick_getacch_int_pl: loop timer turned on at t = " // trim(adjustl(tstr))) + end if + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "kick_getacch_int_pl: stage " // schar ) + end if - if (lfirst) call itimer%reset(param) if (param%lflatten_interactions) then call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) else call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah) end if + if (param%ladaptive_interactions) then + if (itimer%is_on) then + if (param%lflatten_interactions) then + write(lstyle,*) "FLAT" + else + write(lstyle,*) "TRIANGULAR" + end if + call itimer%adapt(param, self, self%nplpl) + write(schar,'(I1)') itimer%stage + write(nstr,*) self%nplpl + write(cstr,*) itimer%count_finish_step + select case(itimer%stage) + case(1) + write(mstr,*) itimer%stage1_metric + case(2) + write(mstr,*) itimer%stage2_metric + end select + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(lstyle)) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + end if + end if + return end subroutine kick_getacch_int_pl diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 46cf56369..2c3b4721e 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -89,6 +89,7 @@ module swiftest_classes !> Each paramter is initialized to a default values. type :: swiftest_parameters integer(I4B) :: integrator = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator used + character(STRMAX) :: param_file_name = "param.in" !! The default name of the parameter input file integer(I4B) :: maxid = -1 !! The current maximum particle id number real(DP) :: t0 = -1.0_DP !! Integration start time real(DP) :: t = -1.0_DP !! Integration current time @@ -1165,7 +1166,7 @@ module subroutine user_kick_getacch_body(self, system, param, t, lbeg) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody_system_object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step end subroutine user_kick_getacch_body diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 97f6fb520..97db74f73 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -3,7 +3,7 @@ module walltime_classes !! !! Classes and methods used to compute elasped wall time use swiftest_globals - use swiftest_classes, only : swiftest_parameters + use swiftest_classes, only : swiftest_parameters, swiftest_pl implicit none public @@ -24,28 +24,23 @@ module walltime_classes end type walltimer type, extends(walltimer) :: interaction_timer - integer(I8B) :: max_interactions = huge(1_I8B) - integer(I4B) :: step_counter - integer(I8B) :: count_previous - logical :: lflatten_interaction_old + integer(I8B) :: max_interactions = huge(1_I8B) !! Stores the number of pl-pl interactions that failed when attempting to flatten (e.g. out of memory). Adapting won't occur if ninteractions > max_interactions + integer(I8B) :: last_interactions = 0 !! Number of interactions that were computed last time. The timer is only run if there has been a change to the number of interactions + integer(I4B) :: step_counter = 0 !! Number of steps that have elapsed since the last timed loop + logical :: is_on = .false. !! The loop timer is currently active + integer(I4B) :: stage = 1 !! The stage of the loop timing (1 or 2) + logical :: stage1_is_flattened !! Logical flag indicating whether stage1 was done with a flat loop (.true.) or triangular loop (.false.) + integer(I8B) :: stage1_ninteractions !! Number of interactions computed during stage 1 + real(DP) :: stage1_metric !! Metric used to judge the performance of a timed loop (e.g. (count_finish_step - count_start_step) / ninteractions) + real(DP) :: stage2_metric !! Metric used to judge the performance of a timed loop (e.g. (count_finish_step - count_start_step) / ninteractions) contains - procedure :: reset => walltime_interaction_reset !! Resets the interaction loop timer, and saves the current value of the array flatten parameter + procedure :: adapt => walltime_interaction_adapt !! Runs the interaction loop adaptation algorithm on an interaction loop + procedure :: check => walltime_interaction_check !! Checks whether or not the loop should be timed and starts the timer if the conditions for starting are met + procedure :: flip => walltime_interaction_flip_loop_style !! Flips the interaction loop style from FLAT to TRIANGULAR or vice vers + procedure :: time_this_loop => walltime_interaction_time_this_loop !! Starts the interaction loop timer end type interaction_timer interface - module subroutine walltime_interaction_reset(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(interaction_timer), intent(inout) :: self !! Walltimer object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine walltime_interaction_reset - - module subroutine walltime_interaction_io_log_start(param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(swiftest_parameters), intent(in) :: param - end subroutine walltime_interaction_io_log_start - module subroutine walltime_finish(self, nsubsteps, message, param) use swiftest_classes, only : swiftest_parameters implicit none @@ -70,25 +65,45 @@ module subroutine walltime_start(self, param) end subroutine walltime_start end interface - contains + interface + module subroutine walltime_interaction_adapt(self, param, pl, ninteractions) + use swiftest_classes, only : swiftest_parameters + implicit none + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing + end subroutine walltime_interaction_adapt - module subroutine walltime_interaction_reset(self, param) - !! author: David A. Minton - !! - !! Resets the interaction loop timer, and saves the current value of the array flatten parameter + module function walltime_interaction_check(self, param, ninteractions) result(ltimeit) use swiftest_classes, only : swiftest_parameters implicit none - ! Arguments - class(interaction_timer), intent(inout) :: self !! Walltimer object + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing + logical :: ltimeit !! Logical flag indicating whether this loop should be timed or not + end function walltime_interaction_check + + module subroutine walltime_interaction_flip_loop_style(self, param, pl) + use swiftest_classes, only : swiftest_parameters, swiftest_pl + implicit none + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + end subroutine walltime_interaction_flip_loop_style - self%lflatten_interaction_old = param%lflatten_interactions - self%step_counter = 0 - call walltime_reset(self, param) + module subroutine walltime_interaction_time_this_loop(self, param, pl, ninteractions) + use swiftest_classes, only : swiftest_parameters, swiftest_pl + implicit none + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop) + end subroutine walltime_interaction_time_this_loop - return - end subroutine walltime_interaction_reset + end interface + contains module subroutine walltime_finish(self, nsubsteps, message, param) !! author: David A. Minton @@ -167,4 +182,111 @@ module subroutine walltime_start(self, param) end subroutine walltime_start + module subroutine walltime_interaction_adapt(self, param, pl, ninteractions) + !! author: David A. Minton + !! + !! Determines which of the two loop styles is fastest and keeps that one + implicit none + class(interaction_timer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing + + ! Record the elapsed time + call system_clock(self%count_finish_step) + + select case(self%stage) + case(1) + self%stage1_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) + case(2) + self%stage2_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) + self%is_on = .false. + self%step_counter = 0 + if (self%stage1_metric < self%stage2_metric) call self%flip(param, pl) ! Go back to the original style, otherwise keep the stage2 style + end select + + return + end subroutine walltime_interaction_adapt + + + module function walltime_interaction_check(self, param, ninteractions) result(ltimeit) + !! author: David A. Minton + !! + !! Checks whether or not the loop should be timed and starts the timer if the conditions for starting are met + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing + logical :: ltimeit !! Logical flag indicating whether this loop should be timed or not + ! Internals + character(len=STRMAX) :: tstring + + if (self%is_on) then ! Entering the second stage of the loop timing. Therefore we will swap the interaction style and time this loop + self%stage = self%stage + 1 + ltimeit = (self%stage == 2) + else + self%step_counter = max(self%step_counter + 1, INTERACTION_TIMER_CADENCE) + ltimeit = .false. + if (self%step_counter == INTERACTION_TIMER_CADENCE) then + ltimeit = (ninteractions /= self%last_interactions) + if (ltimeit) self%stage = 1 + end if + end if + self%is_on = ltimeit + + return + end function walltime_interaction_check + + + module subroutine walltime_interaction_flip_loop_style(self, param, pl) + !! author: David A. Minton + !! + !! Flips the interaction loop style from FLAT to TRIANGULAR or vice versa + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + + param%lflatten_interactions = .not. param%lflatten_interactions + if (param%lflatten_interactions) then + call pl%flatten(param) + else + if (allocated(pl%k_plpl)) deallocate(pl%k_plpl) + end if + + return + end subroutine walltime_interaction_flip_loop_style + + + module subroutine walltime_interaction_time_this_loop(self, param, pl, ninteractions) + !! author: David A. Minton + !! + !! Resets the interaction loop timer, and saves the current value of the array flatten parameter + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop) + + self%is_on = .true. + self%step_counter = 0 + select case(self%stage) + case(1) + self%stage1_ninteractions = ninteractions + self%stage1_is_flattened = param%lflatten_interactions + case(2) + param%lflatten_interactions = self%stage1_is_flattened + call self%flip(param, pl) + case default + self%stage = 1 + end select + call self%reset(param) + + return + end subroutine walltime_interaction_time_this_loop + + end module walltime_classes \ No newline at end of file diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 3f481d9b4..bea4e9de3 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -27,11 +27,10 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms character(len=*),parameter :: linefmt = '(A)' associate(param => self) - + open(unit = unit, file = param%param_file_name, status = 'old', err = 667, iomsg = iomsg) call random_seed(size = nseeds) if (allocated(param%seed)) deallocate(param%seed) allocate(param%seed(nseeds)) - rewind(unit) do read(unit = unit, fmt = linefmt, iostat = iostat, end = 1, err = 667, iomsg = iomsg) line line_trim = trim(adjustl(line)) @@ -79,6 +78,7 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms end if end do 1 continue + close(unit) if (self%GMTINY < 0.0_DP) then write(iomsg,*) "GMTINY invalid or not set: ", self%GMTINY @@ -104,8 +104,9 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms iostat = 0 - 667 continue return + 667 continue + write(*,*) "Error reading SyMBA parameters in param file: ", trim(adjustl(iomsg)) end subroutine symba_io_param_reader diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index e9b7d0405..119899ae4 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -13,6 +13,30 @@ module subroutine symba_kick_getacch_int_pl(self, param) ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameter + ! Internals + type(interaction_timer), save :: itimer + logical, save :: lfirst + character(len=STRMAX) :: tstr, nstr, cstr, mstr + character(len=10) :: lstyle + character(len=1) :: schar + + if (param%ladaptive_interactions) then + if (lfirst) then + call itimer%time_this_loop(param, self, self%nplpl) + lfirst = .false. + else + if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl) + end if + end if + + if (itimer%is_on) then + write(tstr,*) param%t + write(schar,'(I1)') itimer%stage + if (itimer%stage == 1) then + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "symba_kick_getacch_int_pl: loop timer turned on at t = " // trim(adjustl(tstr))) + end if + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "symba_kick_getacch_int_pl: stage " // schar ) + end if if (param%lflatten_interactions) then call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) @@ -20,6 +44,34 @@ module subroutine symba_kick_getacch_int_pl(self, param) call kick_getacch_int_all_triangular_pl(self%nbody, self%nplm, self%xh, self%Gmass, self%radius, self%ah) end if + if (param%ladaptive_interactions) then + if (itimer%is_on) then + if (param%lflatten_interactions) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if + call itimer%adapt(param, self, self%nplpl) + write(schar,'(I1)') itimer%stage + write(nstr,*) self%nplpl + write(cstr,*) itimer%count_finish_step - itimer%count_start_step + select case(itimer%stage) + case(1) + write(mstr,*) itimer%stage1_metric + case(2) + write(mstr,*) itimer%stage2_metric + end select + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + if (param%lflatten_interactions) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) + + end if + end if + return end subroutine symba_kick_getacch_int_pl diff --git a/src/user/user_getacch.f90 b/src/user/user_getacch.f90 index 2775de3dd..a85ec3143 100644 --- a/src/user/user_getacch.f90 +++ b/src/user/user_getacch.f90 @@ -11,7 +11,7 @@ module subroutine user_kick_getacch_body(self, system, param, t, lbeg) ! Arguments class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody_system_object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters user parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters user parameters real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the ste From 88c72d6567a801ce6ec3d58a057ff4990e5fce61 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 18:11:56 -0400 Subject: [PATCH 09/15] Fixed up some bugs in the way loop timing is reported and fixed it so that loops are not timed every time --- examples/symba_mars_disk/param.in | 2 +- src/modules/walltime_classes.f90 | 2 +- src/symba/symba_kick.f90 | 15 +++++++-------- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index b75341794..fd0cc9134 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -31,4 +31,4 @@ MU2KG 1.0 DU2M 1.0 TU2S 1.0 SEED 2 3080983 2220830 -FLATTEN_INTERACTIONS yes +INTERACTION_LOOPS ADAPTIVE diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 97db74f73..4ecc912d1 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -226,7 +226,7 @@ module function walltime_interaction_check(self, param, ninteractions) result(lt self%stage = self%stage + 1 ltimeit = (self%stage == 2) else - self%step_counter = max(self%step_counter + 1, INTERACTION_TIMER_CADENCE) + self%step_counter = min(self%step_counter + 1, INTERACTION_TIMER_CADENCE) ltimeit = .false. if (self%step_counter == INTERACTION_TIMER_CADENCE) then ltimeit = (ninteractions /= self%last_interactions) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 119899ae4..9d2e10666 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -17,7 +17,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) type(interaction_timer), save :: itimer logical, save :: lfirst character(len=STRMAX) :: tstr, nstr, cstr, mstr - character(len=10) :: lstyle + character(len=11) :: lstyle character(len=1) :: schar if (param%ladaptive_interactions) then @@ -59,16 +59,15 @@ module subroutine symba_kick_getacch_int_pl(self, param) case(1) write(mstr,*) itimer%stage1_metric case(2) + if (param%lflatten_interactions) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if write(mstr,*) itimer%stage2_metric + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) end select call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) - if (param%lflatten_interactions) then - write(lstyle,*) "FLAT " - else - write(lstyle,*) "TRIANGULAR" - end if - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) - end if end if From 3fb9f93071aaadcf3725a42dd30a49df339c1fa5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 18:19:26 -0400 Subject: [PATCH 10/15] Fixed issue in which lfirst flag was not set to .true. for the initial loop timing --- src/kick/kick.f90 | 2 +- src/modules/walltime_classes.f90 | 1 - src/symba/symba_kick.f90 | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index edc054caf..078003a2a 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -15,7 +15,7 @@ module subroutine kick_getacch_int_pl(self, param) class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters ! Internals type(interaction_timer), save :: itimer - logical, save :: lfirst + logical, save :: lfirst = .true. character(len=STRMAX) :: tstr, nstr, cstr, mstr, lstyle character(len=1) :: schar diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 4ecc912d1..a06f1103e 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -272,7 +272,6 @@ module subroutine walltime_interaction_time_this_loop(self, param, pl, ninteract integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop) self%is_on = .true. - self%step_counter = 0 select case(self%stage) case(1) self%stage1_ninteractions = ninteractions diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 9d2e10666..9d9915f7c 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -15,7 +15,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameter ! Internals type(interaction_timer), save :: itimer - logical, save :: lfirst + logical, save :: lfirst = .true. character(len=STRMAX) :: tstr, nstr, cstr, mstr character(len=11) :: lstyle character(len=1) :: schar From 53317855de7e9b21782dd3b8079d6b04c7b03c2c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 18:47:58 -0400 Subject: [PATCH 11/15] Fixed bad logging output --- src/kick/kick.f90 | 12 ++++++++++-- src/symba/symba_kick.f90 | 8 +++++--- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 078003a2a..03b6a9493 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -53,14 +53,22 @@ module subroutine kick_getacch_int_pl(self, param) call itimer%adapt(param, self, self%nplpl) write(schar,'(I1)') itimer%stage write(nstr,*) self%nplpl - write(cstr,*) itimer%count_finish_step + write(cstr,*) itimer%count_finish_step - itimer%count_start_step select case(itimer%stage) case(1) write(mstr,*) itimer%stage1_metric case(2) write(mstr,*) itimer%stage2_metric end select - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(lstyle)) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + if (itimer%stage == 2) then + if (param%lflatten_interactions) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) + end if end if end if diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 9d9915f7c..f528c7dc5 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -59,15 +59,17 @@ module subroutine symba_kick_getacch_int_pl(self, param) case(1) write(mstr,*) itimer%stage1_metric case(2) + write(mstr,*) itimer%stage2_metric + end select + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + if (itimer%stage == 2) then if (param%lflatten_interactions) then write(lstyle,*) "FLAT " else write(lstyle,*) "TRIANGULAR" end if - write(mstr,*) itimer%stage2_metric call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) - end select - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + end if end if end if From a0a179b9aca7f7cb11d55657e55d35b3c2174d2b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 21:03:54 -0400 Subject: [PATCH 12/15] Moved logging into the walltime methods. Added a walltime submodule --- Makefile | 6 + examples/symba_mars_disk/param.in | 2 +- src/kick/kick.f90 | 38 +---- src/modules/swiftest_classes.f90 | 1 + src/modules/walltime_classes.f90 | 184 +----------------------- src/symba/symba_kick.f90 | 40 +----- src/walltime/walltime.f90 | 224 ++++++++++++++++++++++++++++++ 7 files changed, 238 insertions(+), 257 deletions(-) create mode 100644 src/walltime/walltime.f90 diff --git a/Makefile b/Makefile index 17ef3c006..56df61d5b 100644 --- a/Makefile +++ b/Makefile @@ -168,6 +168,11 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir + cd $(SWIFTEST_HOME)/src/walltime; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make libdir fast: cd $(SWIFTEST_HOME)/src/fraggle; \ @@ -237,6 +242,7 @@ clean: cd $(SWIFTEST_HOME)/src/tides; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/user; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/util; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/walltime; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/whm; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/bin; rm -f swiftest_* cd $(SWIFTEST_HOME)/bin; rm -f tool_* diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index fd0cc9134..023f31647 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -1,6 +1,6 @@ !Parameter file for the SyMBA-RINGMOONS test T0 0.0 -TSTOP 12000.0 +TSTOP 1200.0 DT 600.0 CB_IN cb.in PL_IN mars.in diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 03b6a9493..b7b100e44 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -19,24 +19,17 @@ module subroutine kick_getacch_int_pl(self, param) character(len=STRMAX) :: tstr, nstr, cstr, mstr, lstyle character(len=1) :: schar + if (param%ladaptive_interactions) then if (lfirst) then call itimer%time_this_loop(param, self, self%nplpl) + write(itimer%loopname, *) "kick_getacch_int_pl" lfirst = .false. else if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl) end if end if - if (itimer%is_on) then - write(tstr,*) param%t - write(schar,'(I1)') itimer%stage - if (itimer%stage == 1) then - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "kick_getacch_int_pl: loop timer turned on at t = " // trim(adjustl(tstr))) - end if - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "kick_getacch_int_pl: stage " // schar ) - end if - if (param%lflatten_interactions) then call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) else @@ -44,32 +37,7 @@ module subroutine kick_getacch_int_pl(self, param) end if if (param%ladaptive_interactions) then - if (itimer%is_on) then - if (param%lflatten_interactions) then - write(lstyle,*) "FLAT" - else - write(lstyle,*) "TRIANGULAR" - end if - call itimer%adapt(param, self, self%nplpl) - write(schar,'(I1)') itimer%stage - write(nstr,*) self%nplpl - write(cstr,*) itimer%count_finish_step - itimer%count_start_step - select case(itimer%stage) - case(1) - write(mstr,*) itimer%stage1_metric - case(2) - write(mstr,*) itimer%stage2_metric - end select - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) - if (itimer%stage == 2) then - if (param%lflatten_interactions) then - write(lstyle,*) "FLAT " - else - write(lstyle,*) "TRIANGULAR" - end if - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) - end if - end if + if (itimer%is_on) call itimer%adapt(param, self, self%nplpl) end if return diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 2c3b4721e..055eb4a79 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -126,6 +126,7 @@ module swiftest_classes character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops + logical :: lflatten_encounters = .false. !! Use the flattened upper triangular matrix for pl-pl encounter check loops logical :: ladaptive_interactions = .false. !! Adaptive interaction loop is turned on ! Logical flags to turn on or off various features of the code diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index a06f1103e..594931a1b 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -24,6 +24,7 @@ module walltime_classes end type walltimer type, extends(walltimer) :: interaction_timer + character(len=STRMAX) :: loopname !! Stores the name of the loop being timed for logging purposes integer(I8B) :: max_interactions = huge(1_I8B) !! Stores the number of pl-pl interactions that failed when attempting to flatten (e.g. out of memory). Adapting won't occur if ninteractions > max_interactions integer(I8B) :: last_interactions = 0 !! Number of interactions that were computed last time. The timer is only run if there has been a change to the number of interactions integer(I4B) :: step_counter = 0 !! Number of steps that have elapsed since the last timed loop @@ -103,189 +104,6 @@ end subroutine walltime_interaction_time_this_loop end interface - contains - - module subroutine walltime_finish(self, nsubsteps, message, param) - !! author: David A. Minton - !! - !! Ends the timer, setting step_finish to the current ticker value and printing the elapsed time information to the terminal - implicit none - ! Arguments - class(walltimer), intent(inout) :: self !! Walltimer object - integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step - character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - character(len=*), parameter :: walltimefmt = '" Wall time (s): ", es12.5, "; Wall time/step in this interval (s): ", es12.5' - character(len=STRMAX) :: fmt - integer(I8B) :: count_delta_step, count_delta_main - real(DP) :: wall_main !! Value of total elapsed time at the end of a timed step - real(DP) :: wall_step !! Value of elapsed time since the start of a timed step - real(DP) :: wall_per_substep !! Value of time per substep - - if (.not.self%lmain_is_started) then - write(*,*) "Wall timer error: The step finish time cannot be calculated because the timer is not started!" - return - end if - - call system_clock(self%count_finish_step) - - count_delta_step = self%count_finish_step - self%count_start_step - count_delta_main = self%count_finish_step - self%count_start_main - wall_step = count_delta_step / (self%count_rate * 1.0_DP) - wall_main = count_delta_main / (self%count_rate * 1.0_DP) - wall_per_substep = wall_step / nsubsteps - - fmt = '("' // adjustl(message) // '",' // walltimefmt // ')' - write(*,trim(adjustl(fmt))) wall_main, wall_per_substep - - call self%start(param) - - return - end subroutine walltime_finish - - - module subroutine walltime_reset(self, param) - !! author: David A. Minton - !! - !! Resets the clock ticker, settting main_start to the current ticker value - implicit none - ! Arguments - class(walltimer), intent(inout) :: self !! Walltimer object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - call system_clock(self%count_start_main, self%count_rate, self%count_max) - self%lmain_is_started = .true. - call self%start(param) - - return - end subroutine walltime_reset - - - module subroutine walltime_start(self, param) - !! author: David A. Minton - !! - !! Starts the timer, setting step_start to the current ticker value - implicit none - ! Arguments - class(walltimer), intent(inout) :: self !! Walltimer object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - if (.not.self%lmain_is_started) then - write(*,*) "Wall timer error: Cannot start the step time until reset is called at least once!" - return - end if - - call system_clock(self%count_start_step) - - return - end subroutine walltime_start - - - module subroutine walltime_interaction_adapt(self, param, pl, ninteractions) - !! author: David A. Minton - !! - !! Determines which of the two loop styles is fastest and keeps that one - implicit none - class(interaction_timer), intent(inout) :: self !! Walltimer object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing - - ! Record the elapsed time - call system_clock(self%count_finish_step) - - select case(self%stage) - case(1) - self%stage1_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) - case(2) - self%stage2_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) - self%is_on = .false. - self%step_counter = 0 - if (self%stage1_metric < self%stage2_metric) call self%flip(param, pl) ! Go back to the original style, otherwise keep the stage2 style - end select - - return - end subroutine walltime_interaction_adapt - - - module function walltime_interaction_check(self, param, ninteractions) result(ltimeit) - !! author: David A. Minton - !! - !! Checks whether or not the loop should be timed and starts the timer if the conditions for starting are met - implicit none - ! Arguments - class(interaction_timer), intent(inout) :: self !! Walltimer object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing - logical :: ltimeit !! Logical flag indicating whether this loop should be timed or not - ! Internals - character(len=STRMAX) :: tstring - - if (self%is_on) then ! Entering the second stage of the loop timing. Therefore we will swap the interaction style and time this loop - self%stage = self%stage + 1 - ltimeit = (self%stage == 2) - else - self%step_counter = min(self%step_counter + 1, INTERACTION_TIMER_CADENCE) - ltimeit = .false. - if (self%step_counter == INTERACTION_TIMER_CADENCE) then - ltimeit = (ninteractions /= self%last_interactions) - if (ltimeit) self%stage = 1 - end if - end if - self%is_on = ltimeit - - return - end function walltime_interaction_check - - - module subroutine walltime_interaction_flip_loop_style(self, param, pl) - !! author: David A. Minton - !! - !! Flips the interaction loop style from FLAT to TRIANGULAR or vice versa - implicit none - ! Arguments - class(interaction_timer), intent(inout) :: self !! Interaction loop timer object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - - param%lflatten_interactions = .not. param%lflatten_interactions - if (param%lflatten_interactions) then - call pl%flatten(param) - else - if (allocated(pl%k_plpl)) deallocate(pl%k_plpl) - end if - - return - end subroutine walltime_interaction_flip_loop_style - - - module subroutine walltime_interaction_time_this_loop(self, param, pl, ninteractions) - !! author: David A. Minton - !! - !! Resets the interaction loop timer, and saves the current value of the array flatten parameter - implicit none - ! Arguments - class(interaction_timer), intent(inout) :: self !! Interaction loop timer object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop) - - self%is_on = .true. - select case(self%stage) - case(1) - self%stage1_ninteractions = ninteractions - self%stage1_is_flattened = param%lflatten_interactions - case(2) - param%lflatten_interactions = self%stage1_is_flattened - call self%flip(param, pl) - case default - self%stage = 1 - end select - call self%reset(param) - - return - end subroutine walltime_interaction_time_this_loop end module walltime_classes \ No newline at end of file diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index f528c7dc5..4a11e532a 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -16,28 +16,17 @@ module subroutine symba_kick_getacch_int_pl(self, param) ! Internals type(interaction_timer), save :: itimer logical, save :: lfirst = .true. - character(len=STRMAX) :: tstr, nstr, cstr, mstr - character(len=11) :: lstyle - character(len=1) :: schar if (param%ladaptive_interactions) then if (lfirst) then call itimer%time_this_loop(param, self, self%nplpl) + write(itimer%loopname, *) "symba_kick_getacch_int_pl" lfirst = .false. else if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl) end if end if - if (itimer%is_on) then - write(tstr,*) param%t - write(schar,'(I1)') itimer%stage - if (itimer%stage == 1) then - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "symba_kick_getacch_int_pl: loop timer turned on at t = " // trim(adjustl(tstr))) - end if - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "symba_kick_getacch_int_pl: stage " // schar ) - end if - if (param%lflatten_interactions) then call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) else @@ -45,32 +34,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) end if if (param%ladaptive_interactions) then - if (itimer%is_on) then - if (param%lflatten_interactions) then - write(lstyle,*) "FLAT " - else - write(lstyle,*) "TRIANGULAR" - end if - call itimer%adapt(param, self, self%nplpl) - write(schar,'(I1)') itimer%stage - write(nstr,*) self%nplpl - write(cstr,*) itimer%count_finish_step - itimer%count_start_step - select case(itimer%stage) - case(1) - write(mstr,*) itimer%stage1_metric - case(2) - write(mstr,*) itimer%stage2_metric - end select - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) - if (itimer%stage == 2) then - if (param%lflatten_interactions) then - write(lstyle,*) "FLAT " - else - write(lstyle,*) "TRIANGULAR" - end if - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) - end if - end if + if (itimer%is_on) call itimer%adapt(param, self, self%nplpl) end if return diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 new file mode 100644 index 000000000..325c9a2c4 --- /dev/null +++ b/src/walltime/walltime.f90 @@ -0,0 +1,224 @@ +submodule(walltime_classes) s_walltime + use swiftest +contains + + module subroutine walltime_finish(self, nsubsteps, message, param) + !! author: David A. Minton + !! + !! Ends the timer, setting step_finish to the current ticker value and printing the elapsed time information to the terminal + implicit none + ! Arguments + class(walltimer), intent(inout) :: self !! Walltimer object + integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step + character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + character(len=*), parameter :: walltimefmt = '" Wall time (s): ", es12.5, "; Wall time/step in this interval (s): ", es12.5' + character(len=STRMAX) :: fmt + integer(I8B) :: count_delta_step, count_delta_main + real(DP) :: wall_main !! Value of total elapsed time at the end of a timed step + real(DP) :: wall_step !! Value of elapsed time since the start of a timed step + real(DP) :: wall_per_substep !! Value of time per substep + + if (.not.self%lmain_is_started) then + write(*,*) "Wall timer error: The step finish time cannot be calculated because the timer is not started!" + return + end if + + call system_clock(self%count_finish_step) + + count_delta_step = self%count_finish_step - self%count_start_step + count_delta_main = self%count_finish_step - self%count_start_main + wall_step = count_delta_step / (self%count_rate * 1.0_DP) + wall_main = count_delta_main / (self%count_rate * 1.0_DP) + wall_per_substep = wall_step / nsubsteps + + fmt = '("' // adjustl(message) // '",' // walltimefmt // ')' + write(*,trim(adjustl(fmt))) wall_main, wall_per_substep + + call self%start(param) + + return + end subroutine walltime_finish + + + module subroutine walltime_reset(self, param) + !! author: David A. Minton + !! + !! Resets the clock ticker, settting main_start to the current ticker value + implicit none + ! Arguments + class(walltimer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + call system_clock(self%count_start_main, self%count_rate, self%count_max) + self%lmain_is_started = .true. + call self%start(param) + + return + end subroutine walltime_reset + + + module subroutine walltime_start(self, param) + !! author: David A. Minton + !! + !! Starts the timer, setting step_start to the current ticker value + implicit none + ! Arguments + class(walltimer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + if (.not.self%lmain_is_started) then + write(*,*) "Wall timer error: Cannot start the step time until reset is called at least once!" + return + end if + + call system_clock(self%count_start_step) + + return + end subroutine walltime_start + + + module subroutine walltime_interaction_adapt(self, param, pl, ninteractions) + !! author: David A. Minton + !! + !! Determines which of the two loop styles is fastest and keeps that one + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing + ! Internals + character(len=STRMAX) :: tstr, nstr, cstr, mstr + character(len=11) :: lstyle + character(len=1) :: schar + + ! Record the elapsed time + call system_clock(self%count_finish_step) + + if (param%lflatten_interactions) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if + write(schar,'(I1)') self%stage + write(nstr,*) ninteractions + + select case(self%stage) + case(1) + self%stage1_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) + write(mstr,*) self%stage2_metric + case(2) + self%stage2_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) + self%is_on = .false. + self%step_counter = 0 + if (self%stage1_metric < self%stage2_metric) call self%flip(param, pl) ! Go back to the original style, otherwise keep the stage2 style + write(mstr,*) self%stage1_metric + end select + + write(cstr,*) self%count_finish_step - self%count_start_step + + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + + if (self%stage == 2) then + if (param%lflatten_interactions) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) + end if + + return + end subroutine walltime_interaction_adapt + + + module function walltime_interaction_check(self, param, ninteractions) result(ltimeit) + !! author: David A. Minton + !! + !! Checks whether or not the loop should be timed and starts the timer if the conditions for starting are met + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Walltimer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop and to determine if number of interactions has changed since the last timing + logical :: ltimeit !! Logical flag indicating whether this loop should be timed or not + ! Internals + character(len=STRMAX) :: tstring + + if (self%is_on) then ! Entering the second stage of the loop timing. Therefore we will swap the interaction style and time this loop + self%stage = self%stage + 1 + ltimeit = (self%stage == 2) + else + self%step_counter = min(self%step_counter + 1, INTERACTION_TIMER_CADENCE) + ltimeit = .false. + if (self%step_counter == INTERACTION_TIMER_CADENCE) then + ltimeit = (ninteractions /= self%last_interactions) + if (ltimeit) self%stage = 1 + end if + end if + self%is_on = ltimeit + + return + end function walltime_interaction_check + + + module subroutine walltime_interaction_flip_loop_style(self, param, pl) + !! author: David A. Minton + !! + !! Flips the interaction loop style from FLAT to TRIANGULAR or vice versa + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + + param%lflatten_interactions = .not. param%lflatten_interactions + if (param%lflatten_interactions) then + call pl%flatten(param) + else + if (allocated(pl%k_plpl)) deallocate(pl%k_plpl) + end if + + return + end subroutine walltime_interaction_flip_loop_style + + + module subroutine walltime_interaction_time_this_loop(self, param, pl, ninteractions) + !! author: David A. Minton + !! + !! Resets the interaction loop timer, and saves the current value of the array flatten parameter + implicit none + ! Arguments + class(interaction_timer), intent(inout) :: self !! Interaction loop timer object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + integer(I8B), intent(in) :: ninteractions !! Current number of interactions (used to normalize the timed loop) + ! Internals + character(len=STRMAX) :: tstr + character(len=1) :: schar + + self%is_on = .true. + write(tstr,*) param%t + select case(self%stage) + case(1) + self%stage1_ninteractions = ninteractions + self%stage1_is_flattened = param%lflatten_interactions + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(self%loopname)) // ": loop timer turned on at t = " // trim(adjustl(tstr))) + case(2) + param%lflatten_interactions = self%stage1_is_flattened + call self%flip(param, pl) + case default + self%stage = 1 + end select + + write(schar,'(I1)') self%stage + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(self%loopname)) // ": stage " // schar ) + + call self%reset(param) + + return + end subroutine walltime_interaction_time_this_loop + +end submodule s_walltime \ No newline at end of file From 01b7e767c1de7615444e142f7d0173faec11e1d3 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 23:26:01 -0400 Subject: [PATCH 13/15] Restructured so that kicks and encounters can be timed and adapted separately --- src/io/io.f90 | 4 ++ src/kick/kick.f90 | 3 +- src/modules/symba_classes.f90 | 2 +- src/modules/walltime_classes.f90 | 3 +- src/symba/symba_encounter_check.f90 | 22 ++++++++++- src/symba/symba_kick.f90 | 7 ++-- src/symba/symba_util.f90 | 23 +++++++----- src/util/util_flatten.f90 | 3 +- src/walltime/walltime.f90 | 58 +++++++++++++++++++++-------- 9 files changed, 92 insertions(+), 33 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index cf9c6a8a9..226811d61 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -815,14 +815,17 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case("ADAPTIVE") param%ladaptive_interactions = .true. param%lflatten_interactions = .true. + param%lflatten_encounters = .true. call io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile") call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") case("TRIANGULAR") param%ladaptive_interactions = .false. param%lflatten_interactions = .false. + param%lflatten_encounters = .false. case("FLAT") param%ladaptive_interactions = .false. param%lflatten_interactions = .true. + param%lflatten_encounters = .true. case default write(*,*) "Unknown value for parameter INTERACTION_LOOPS: -> ",trim(adjustl(param%interaction_loops)) write(*,*) "Must be one of the following: TRIANGULAR, FLAT, or ADAPTIVE" @@ -830,6 +833,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param%interaction_loops = "ADAPTIVE" param%ladaptive_interactions = .true. param%lflatten_interactions = .true. + param%lflatten_encounters = .true. end select iostat = 0 diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index b7b100e44..637b38720 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -22,8 +22,9 @@ module subroutine kick_getacch_int_pl(self, param) if (param%ladaptive_interactions) then if (lfirst) then + write(itimer%loopname, *) "kick_getacch_int_pl" + write(itimer%looptype, *) "INTERACTION" call itimer%time_this_loop(param, self, self%nplpl) - write(itimer%loopname, *) "kick_getacch_int_pl" lfirst = .false. else if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 6bfef7bd9..a1c56a456 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -259,7 +259,7 @@ end subroutine symba_drift_tp module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter) implicit none class(symba_pl), intent(inout) :: self !! SyMBA test particle object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 594931a1b..a538baad8 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -24,7 +24,8 @@ module walltime_classes end type walltimer type, extends(walltimer) :: interaction_timer - character(len=STRMAX) :: loopname !! Stores the name of the loop being timed for logging purposes + character(len=STRMAX) :: loopname !! Stores the name of the loop being timed for logging purposes + character(len=NAMELEN) :: looptype !! Stores the type of loop (e.g. INTERACTION or ENCOUNTER) integer(I8B) :: max_interactions = huge(1_I8B) !! Stores the number of pl-pl interactions that failed when attempting to flatten (e.g. out of memory). Adapting won't occur if ninteractions > max_interactions integer(I8B) :: last_interactions = 0 !! Number of interactions that were computed last time. The timer is only run if there has been a change to the number of interactions integer(I4B) :: step_counter = 0 !! Number of steps that have elapsed since the last timed loop diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 3ba6b4adf..8215e4d85 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -130,7 +130,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA test particle object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -142,13 +142,28 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr integer(I4B), dimension(:), allocatable :: index1, index2 integer(I4B), dimension(:,:), allocatable :: k_plpl_enc + type(interaction_timer), save :: itimer + logical, save :: lfirst = .true. if (self%nbody == 0) return associate(pl => self, plplenc_list => system%plplenc_list) + + if (param%ladaptive_interactions) then + if (lfirst) then + write(itimer%loopname, *) "symba_encounter_check_pl" + write(itimer%looptype, *) "ENCOUNTERS" + call itimer%time_this_loop(param, pl, pl%nplplm) + lfirst = .false. + else + if (itimer%check(param, pl%nplplm)) call itimer%time_this_loop(param, pl, pl%nplplm) + end if + end if + npl = pl%nbody if (param%lflatten_interactions) then nplplm = pl%nplplm + allocate(lencounter(nplplm)) allocate(loc_lvdotr(nplplm)) @@ -200,6 +215,11 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l pl%nplenc(j) = pl%nplenc(j) + 1 end do end if + + if (param%ladaptive_interactions) then + if (itimer%is_on) call itimer%adapt(param, pl, pl%nplplm) + end if + end associate return diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 4a11e532a..529f00215 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -19,11 +19,12 @@ module subroutine symba_kick_getacch_int_pl(self, param) if (param%ladaptive_interactions) then if (lfirst) then - call itimer%time_this_loop(param, self, self%nplpl) write(itimer%loopname, *) "symba_kick_getacch_int_pl" + write(itimer%looptype, *) "INTERACTION" + call itimer%time_this_loop(param, self, self%nplplm) lfirst = .false. else - if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl) + if (itimer%check(param, self%nplplm)) call itimer%time_this_loop(param, self, self%nplplm) end if end if @@ -34,7 +35,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) end if if (param%ladaptive_interactions) then - if (itimer%is_on) call itimer%adapt(param, self, self%nplpl) + if (itimer%is_on) call itimer%adapt(param, self, self%nplplm) end if return diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 4b20cf0bf..8065342c6 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -285,25 +285,28 @@ module subroutine symba_util_flatten_eucl_plpl(self, param) class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I8B) :: k, nplpl, nplplm - integer(I4B) :: i, j, npl, nplm, ip, jp + integer(I8B) :: k + integer(I4B) :: i, j, npl, nplm, err - associate(pl => self) + associate(pl => self, nplpl => self%nplpl, nplplm => self%nplplm) npl = int(self%nbody, kind=I8B) nplm = count(.not. pl%lmtiny(1:npl)) pl%nplm = int(nplm, kind=I4B) - pl%nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column - pl%nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies - if (param%lflatten_interactions) then + nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column + nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies + if ((param%lflatten_interactions) .or. (param%lflatten_encounters)) then if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously - allocate(self%k_plpl(2, pl%nplpl)) - do concurrent (i = 1:npl) - do concurrent (j = i+1:npl) + allocate(self%k_plpl(2, nplpl), stat=err) + if (err /=0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode + param%lflatten_interactions = .false. + param%lflatten_encounters = .false. + else + do concurrent (i=1:npl, j=1:npl, j>i) call util_flatten_eucl_ij_to_k(npl, i, j, k) self%k_plpl(1, k) = i self%k_plpl(2, k) = j end do - end do + end if end if end associate diff --git a/src/util/util_flatten.f90 b/src/util/util_flatten.f90 index e3e97c1e7..6d2027049 100644 --- a/src/util/util_flatten.f90 +++ b/src/util/util_flatten.f90 @@ -82,11 +82,12 @@ module subroutine util_flatten_eucl_plpl(self, param) npl = int(self%nbody, kind=I8B) associate(nplpl => self%nplpl) nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl - if (param%lflatten_interactions) then + if ((param%lflatten_interactions) .or. (param%lflatten_encounters)) then if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously allocate(self%k_plpl(2, nplpl), stat=err) if (err /=0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode param%lflatten_interactions = .false. + param%lflatten_encounters = .false. else do concurrent (i=1:npl, j=1:npl, j>i) call util_flatten_eucl_ij_to_k(npl, i, j, k) diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index 325c9a2c4..2da539285 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -93,28 +93,40 @@ module subroutine walltime_interaction_adapt(self, param, pl, ninteractions) character(len=STRMAX) :: tstr, nstr, cstr, mstr character(len=11) :: lstyle character(len=1) :: schar + logical :: lflatten_final ! Record the elapsed time call system_clock(self%count_finish_step) - if (param%lflatten_interactions) then - write(lstyle,*) "FLAT " - else - write(lstyle,*) "TRIANGULAR" - end if write(schar,'(I1)') self%stage write(nstr,*) ninteractions select case(self%stage) case(1) + if (self%stage1_is_flattened) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if self%stage1_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) - write(mstr,*) self%stage2_metric + write(mstr,*) self%stage1_metric case(2) + if (.not.self%stage1_is_flattened) then + write(lstyle,*) "FLAT " + else + write(lstyle,*) "TRIANGULAR" + end if + self%stage2_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP) self%is_on = .false. self%step_counter = 0 - if (self%stage1_metric < self%stage2_metric) call self%flip(param, pl) ! Go back to the original style, otherwise keep the stage2 style - write(mstr,*) self%stage1_metric + if (self%stage1_metric < self%stage2_metric) then + lflatten_final = self%stage1_is_flattened + call self%flip(param, pl) ! Go back to the original style, otherwise keep the stage2 style + else + lflatten_final = .not.self%stage1_is_flattened + end if + write(mstr,*) self%stage2_metric end select write(cstr,*) self%count_finish_step - self%count_start_step @@ -122,12 +134,12 @@ module subroutine walltime_interaction_adapt(self, param, pl, ninteractions) call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) if (self%stage == 2) then - if (param%lflatten_interactions) then + if (lflatten_final) then write(lstyle,*) "FLAT " else write(lstyle,*) "TRIANGULAR" end if - call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle))) + call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(self%loopname)) // ": the fastest loop method tested is " // trim(adjustl(lstyle))) end if return @@ -174,8 +186,14 @@ module subroutine walltime_interaction_flip_loop_style(self, param, pl) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - param%lflatten_interactions = .not. param%lflatten_interactions - if (param%lflatten_interactions) then + select case(trim(adjustl(self%looptype))) + case("INTERACTIONS") + param%lflatten_interactions = .not. param%lflatten_interactions + case("ENCOUNTERS") + param%lflatten_encounters = .not. param%lflatten_encounters + end select + + if ((param%lflatten_interactions) .or. (param%lflatten_encounters)) then call pl%flatten(param) else if (allocated(pl%k_plpl)) deallocate(pl%k_plpl) @@ -203,11 +221,21 @@ module subroutine walltime_interaction_time_this_loop(self, param, pl, ninteract write(tstr,*) param%t select case(self%stage) case(1) - self%stage1_ninteractions = ninteractions - self%stage1_is_flattened = param%lflatten_interactions + self%stage1_ninteractions = ninteractions + select case(trim(adjustl(self%looptype))) + case("INTERACTIONS") + self%stage1_is_flattened = param%lflatten_interactions + case("ENCOUNTERS") + self%stage1_is_flattened = param%lflatten_encounters + end select call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(self%loopname)) // ": loop timer turned on at t = " // trim(adjustl(tstr))) case(2) - param%lflatten_interactions = self%stage1_is_flattened + select case(trim(adjustl(self%looptype))) + case("INTERACTIONS") + param%lflatten_interactions = self%stage1_is_flattened + case("ENCOUNTERS") + param%lflatten_encounters = self%stage1_is_flattened + end select call self%flip(param, pl) case default self%stage = 1 From 3b9999c7dbf1b65a367171c94aa51d0657422129 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 23:31:02 -0400 Subject: [PATCH 14/15] Fixed bad flag in encounter check --- src/symba/symba_encounter_check.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 8215e4d85..4896721d3 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -161,7 +161,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l end if npl = pl%nbody - if (param%lflatten_interactions) then + if (param%lflatten_encounters) then nplplm = pl%nplplm allocate(lencounter(nplplm)) From de11649f7679c4ff7c47574a61c9446a99f5b68a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 21 Sep 2021 16:56:24 -0400 Subject: [PATCH 15/15] Vectorized drift --- Makefile | 29 ++++++----- src/drift/drift.f90 | 43 +++++++++-------- src/helio/helio_drift.f90 | 39 +++++++++++++-- src/modules/swiftest_classes.f90 | 47 ++++++++---------- src/orbel/orbel.f90 | 83 ++++++++++++++++++++------------ src/rmvs/rmvs_step.f90 | 3 +- src/symba/symba_collision.f90 | 2 +- src/symba/symba_util.f90 | 5 +- src/util/util_peri.f90 | 4 +- 9 files changed, 158 insertions(+), 97 deletions(-) diff --git a/Makefile b/Makefile index 56df61d5b..c3cce41c2 100644 --- a/Makefile +++ b/Makefile @@ -93,11 +93,6 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir - cd $(SWIFTEST_HOME)/src/drift; \ - rm -f Makefile.Defines Makefile; \ - ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ - ln -s $(SWIFTEST_HOME)/Makefile .; \ - make libdir cd $(SWIFTEST_HOME)/src/gr; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ @@ -128,11 +123,6 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir - cd $(SWIFTEST_HOME)/src/orbel; \ - rm -f Makefile.Defines Makefile; \ - ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ - ln -s $(SWIFTEST_HOME)/Makefile .; \ - make libdir cd $(SWIFTEST_HOME)/src/setup; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ @@ -187,6 +177,24 @@ fast: ln -s $(SWIFTEST_HOME)/Makefile .; \ make fastdir + cd $(SWIFTEST_HOME)/src/orbel; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make fastdir + + cd $(SWIFTEST_HOME)/src/drift; \ + rm -f Makefile.Defines Makefile; \ + ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ + ln -s $(SWIFTEST_HOME)/Makefile .; \ + make fastdir + + cd $(SWIFTEST_HOME)/src/helio; \ + $(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c helio_drift.f90; \ + $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \ + $(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \ + rm -f *.o *.smod + cd $(SWIFTEST_HOME)/src/rmvs; \ $(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c rmvs_encounter_check.f90; \ $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \ @@ -211,7 +219,6 @@ fastdir: $(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \ rm -f *.o *.smod - drivers: cd $(SWIFTEST_HOME)/src/main; \ rm -f Makefile.Defines Makefile; \ diff --git a/src/drift/drift.f90 b/src/drift/drift.f90 index e929a95fa..296ceb553 100644 --- a/src/drift/drift.f90 +++ b/src/drift/drift.f90 @@ -77,18 +77,18 @@ module subroutine drift_all(mu, x, v, n, param, dt, lmask, iflag) else where(lmask(1:n)) dtp(1:n) = dt end if - !$omp parallel do default(private) & - !$omp shared(n, lmask, mu, x, v, dtp, iflag) - do i = 1,n + !$omp simd + do i = 1, n if (lmask(i)) call drift_one(mu(i), x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), dtp(i), iflag(i)) end do - !$omp end parallel do + !$omp end simd return end subroutine drift_all module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) + !$omp declare simd(drift_one) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Perform Danby drift for one body, redoing drift with smaller substeps if original accuracy is insufficient @@ -104,27 +104,22 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag ! Internals integer(I4B) :: i real(DP) :: dttmp - real(DP), dimension(NDIM) :: x, v - x = [px, py, pz] - v = [vx, vy, vz] - - call drift_dan(mu, x(:), v(:), dt, iflag) + call drift_dan(mu, px, py, pz, vx, vy, vz, dt, iflag) if (iflag /= 0) then dttmp = 0.1_DP * dt do i = 1, 10 - call drift_dan(mu, x(:), v(:), dttmp, iflag) + call drift_dan(mu, px, py, pz, vx, vy, vz, dttmp, iflag) if (iflag /= 0) exit end do end if - px = x(1); py = x(2); pz = x(3) - vx = v(1); vy = v(2); vz = v(3) return end subroutine drift_one - pure subroutine drift_dan(mu, x0, v0, dt0, iflag) + pure subroutine drift_dan(mu, px0, py0, pz0, vx0, vy0, vz0, dt0, iflag) + !$omp declare simd(drift_dan) !! author: David A. Minton !! !! Perform Kepler drift, solving Kepler's equation in appropriate variables @@ -134,14 +129,16 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) implicit none integer(I4B), intent(out) :: iflag real(DP), intent(in) :: mu, dt0 - real(DP), dimension(:), intent(inout) :: x0, v0 + real(DP), intent(inout) :: px0, py0, pz0, vx0, vy0, vz0 real(DP) :: dt, f, g, fdot, gdot, c1, c2, c3, u, alpha, fp, r0 real(DP) :: v0s, a, asq, en, dm, ec, es, esq, xkep, fchk, s, c - real(DP), dimension(NDIM) :: x, v + real(DP), dimension(NDIM) :: x, v, x0, v0 ! Executable code iflag = 0 dt = dt0 + x0 = [px0, py0, pz0] + v0 = [vx0, vy0, vz0] r0 = sqrt(dot_product(x0(:), x0(:))) v0s = dot_product(v0(:), v0(:)) u = dot_product(x0(:), v0(:)) @@ -172,8 +169,8 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) gdot = (c - 1.0_DP) / fp + 1.0_DP x(:) = x0(:) * f + v0(:) * g v(:) = x0(:) * fdot + v0(:) * gdot - x0(:) = x(:) - v0(:) = v(:) + px0 = x(1); py0 = x(2); pz0 = x(3) + vx0 = v(1); vy0 = v(2); vz0 = v(3) iflag = 0 return end if @@ -187,8 +184,8 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) gdot = 1.0_DP - mu / fp * c2 x(:) = x0(:) * f + v0(:) * g v(:) = x0(:) * fdot + v0(:) * gdot - x0(:) = x(:) - v0(:) = v(:) + px0 = x(1); py0 = x(2); pz0 = x(3) + vx0 = v(1); vy0 = v(2); vz0 = v(3) end if return @@ -196,6 +193,7 @@ end subroutine drift_dan pure subroutine drift_kepmd(dm, es, ec, x, s, c) + !$omp declare simd(drift_kepmd) !! author: David A. Minton !! !! Solve Kepler's equation in difference form for an ellipse for small input dm and eccentricity @@ -235,6 +233,7 @@ end subroutine drift_kepmd pure subroutine drift_kepu(dt,r0,mu,alpha,u,fp,c1,c2,c3,iflag) + !$omp declare simd(drift_kepu) !! author: David A. Minton !! !! Solve Kepler's equation in universal variables @@ -263,6 +262,7 @@ end subroutine drift_kepu pure subroutine drift_kepu_fchk(dt, r0, mu, alpha, u, s, f) + !$omp declare simd(drift_kepu_fchk) !! author: David A. Minton !! !! Computes the value of f, the function whose root we are trying to find in universal variables @@ -286,6 +286,7 @@ end subroutine drift_kepu_fchk pure subroutine drift_kepu_guess(dt, r0, mu, alpha, u, s) + !$omp declare simd(drift_kepu_guess) !! author: David A. Minton !! !! Compute initial guess for solving Kepler's equation using universal variables @@ -324,6 +325,7 @@ end subroutine drift_kepu_guess pure subroutine drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) + !$omp declare simd(drift_kepu_lag) !! author: David A. Minton !! !! Solve Kepler's equation in universal variables using Laguerre's method @@ -369,6 +371,7 @@ end subroutine drift_kepu_lag pure subroutine drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) + !$omp declare simd(drift_kepu_new) !! author: David A. Minton !! !! Solve Kepler's equation in universal variables using Newton's method @@ -411,6 +414,7 @@ end subroutine drift_kepu_new pure subroutine drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag) + !$omp declare simd(drift_kepu_p3solve) !! author: David A. Minton !! !! Computes real root of cubic involved in setting initial guess for solving Kepler's equation in universal variables @@ -455,6 +459,7 @@ end subroutine drift_kepu_p3solve pure subroutine drift_kepu_stumpff(x, c0, c1, c2, c3) + !$omp declare simd(drift_kepu_stumpff) !! author: David A. Minton !! !! Compute Stumpff functions needed for Kepler drift in universal variables diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index e2a55e458..03ebfbd30 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -74,6 +74,41 @@ module subroutine helio_drift_tp(self, system, param, dt) return end subroutine helio_drift_tp + + pure elemental subroutine helio_drift_linear_one(xhx, xhy, xhz, ptx, pty, ptz, dt) + !$omp declare simd(helio_drift_linear_one) + implicit none + real(DP), intent(inout) :: xhx, xhy, xhz + real(DP), intent(in) :: ptx, pty, ptz, dt + + xhx = xhx + ptx * dt + xhy = xhy + pty * dt + xhz = xhz + ptz * dt + + return + end subroutine helio_drift_linear_one + + + subroutine helio_drift_linear_all(xh, pt, dt, n, lmask) + implicit none + ! Arguments + real(DP), dimension(:,:), intent(inout) :: xh + real(DP), dimension(:), intent(in) :: pt + real(DP), intent(in) :: dt + integer(I4B), intent(in) :: n + logical, dimension(:), intent(in) :: lmask + ! Internals + integer(I4B) :: i + + !$omp parallel do simd default(shared) schedule(static) + do i = 1, n + if (lmask(i)) call helio_drift_linear_one(xh(1,i), xh(2,i), xh(3,i), pt(1), pt(2), pt(3), dt) + end do + !$omp end parallel do simd + + return + end subroutine helio_drift_linear_all + module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) !! author: David A. Minton @@ -100,9 +135,7 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl), self%lmask(1:npl)) pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl), self%lmask(1:npl)) pt(:) = pt(:) / cb%Gmass - do concurrent(i = 1:npl, self%lmask(i)) - pl%xh(:,i) = pl%xh(:,i) + pt(:) * dt - end do + call helio_drift_linear_all(pl%xh(:,:), pt(:), dt, npl, pl%lmask(:)) if (lbeg) then cb%ptbeg = pt(:) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 055eb4a79..cd68e0059 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -558,6 +558,7 @@ module subroutine drift_body(self, system, param, dt) end subroutine drift_body module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) + !$omp declare simd(drift_one) implicit none real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift real(DP), intent(inout) :: px, py, pz, vx, vy, vz !! Position and velocity of body to drift @@ -1059,45 +1060,35 @@ module subroutine orbel_el2xv_vec(self, cb) end subroutine orbel_el2xv_vec module pure subroutine orbel_scget(angle, sx, cx) + !$omp declare simd(orbel_scget) implicit none real(DP), intent(in) :: angle real(DP), intent(out) :: sx, cx end subroutine orbel_scget - module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) + module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q) + !$omp declare simd(orbel_xv2aeq) implicit none - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: e !! eccentricity - real(DP), intent(out) :: q !! periapsis + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: vx,vy,vz !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: e !! eccentricity + real(DP), intent(out) :: q !! periapsis end subroutine orbel_xv2aeq - module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) + module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi) + !$omp declare simd(orbel_xv2aqt) implicit none - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: q !! periapsis - real(DP), intent(out) :: capm !! mean anomaly - real(DP), intent(out) :: tperi !! time of pericenter passage + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: vx,vy,vz !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: q !! periapsis + real(DP), intent(out) :: capm !! mean anomaly + real(DP), intent(out) :: tperi !! time of pericenter passage end subroutine orbel_xv2aqt - module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) - implicit none - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: e !! eccentricity - real(DP), intent(out) :: inc !! inclination - real(DP), intent(out) :: capom !! longitude of ascending node - real(DP), intent(out) :: omega !! argument of periapsis - real(DP), intent(out) :: capm !! mean anomaly - end subroutine orbel_xv2el - module subroutine orbel_xv2el_vec(self, cb) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90 index 620e0f839..a815b92c7 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -129,6 +129,7 @@ end subroutine orbel_el2xv module pure subroutine orbel_scget(angle, sx, cx) + !$omp declare simd(orbel_scget) !! author: David A. Minton !! !! Efficiently compute the sine and cosine of an input angle @@ -181,7 +182,7 @@ end subroutine orbel_scget ! REVISIONS: !********************************************************************** pure subroutine orbel_schget(angle,shx,chx) - + !$omp declare simd(orbel_schget) real(DP), intent(in) :: angle real(DP), intent(out) :: shx,chx @@ -212,6 +213,7 @@ end subroutine orbel_schget ! REVISIONS: !********************************************************************** real(DP) pure function orbel_flon(e,icapn) + !$omp declare simd(orbel_flon) implicit none real(DP), intent(in) :: e, icapn integer(I4B) :: iflag,i @@ -315,6 +317,7 @@ end function orbel_flon ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_fget(e,capn) + !$omp declare simd(orbel_fget) implicit none real(DP), intent(in) :: e,capn @@ -385,6 +388,7 @@ end function orbel_fget ! series for small Q. !********************************************************************** real(DP) pure function orbel_zget(iq) + !$omp declare simd(orbel_zget) implicit none real(DP), intent(in) :: iq @@ -440,6 +444,7 @@ end function orbel_zget ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_esolmd(e,m) + !$omp declare simd(orbel_esolmd) implicit none real(DP), intent(in) :: e @@ -495,6 +500,7 @@ end function orbel_esolmd ! REVISIONS: !********************************************************************** real(DP) pure function orbel_ehie(e,im) + !$omp declare simd(orbel_ehie) implicit none real(DP), intent(in) :: e,im @@ -570,13 +576,13 @@ end function orbel_ehie ! we have an ellipse with e between 0.15 and 0.8 !********************************************************************** real(DP) pure function orbel_eget(e,m) + !$omp declare simd(orbel_eget) implicit none real(DP), intent(in) :: e,m real(DP) :: x,sm,cm,sx,cx real(DP) :: es,ec,f,fp,fpp,fppp,dx - ! function to solve kepler's eqn for e (here called ! x) for given e and m. returns value of x. ! may 21 : for e < 0.18 use esolmd for speed and sufficient accuracy @@ -644,6 +650,7 @@ end function orbel_eget ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_ehybrid(e,m) + !$omp declare simd(orbel_ehybrid) implicit none real(DP), intent(in) :: e,m @@ -683,6 +690,7 @@ end function orbel_ehybrid ! REVISIONS: 2/26/93 hfl !********************************************************************** real(DP) pure function orbel_fhybrid(e,n) + !$omp declare simd(orbel_fhybrid) implicit none real(DP), intent(in) :: e,n @@ -701,7 +709,8 @@ real(DP) pure function orbel_fhybrid(e,n) end function orbel_fhybrid - module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) + module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q) + !$omp declare simd(orbel_xv2aeq) !! author: David A. Minton !! !! Compute semimajor axis, eccentricity, and pericentric distance from relative Cartesian position and velocity @@ -710,16 +719,22 @@ module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) !! Adapted from Luke Dones' Swift routine orbel_xv2aeq.f implicit none !! Arguments - real(DP), intent(in) :: mu - real(DP), dimension(:), intent(in) :: x, v - real(DP), intent(out) :: a, e, q + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: vx,vy,vz !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: e !! eccentricity + real(DP), intent(out) :: q !! periapsis + ! Internals integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, energy, fac - real(DP), dimension(NDIM) :: hvec + real(DP), dimension(NDIM) :: hvec, x, v a = 0.0_DP e = 0.0_DP q = 0.0_DP + x = [px, py, pz] + v = [vx, vy, vz] r = sqrt(dot_product(x(:), x(:))) v2 = dot_product(v(:), v(:)) hvec(:) = x(:) .cross. v(:) @@ -760,7 +775,8 @@ module pure subroutine orbel_xv2aeq(mu, x, v, a, e, q) end subroutine orbel_xv2aeq - module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) + module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi) + !$omp declare simd(orbel_xv2aqt) !! author: David A. Minton !! !! Compute semimajor axis, pericentric distance, mean anomaly, and time to nearest pericenter passage from @@ -771,22 +787,24 @@ module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aqt.f90 implicit none ! Arguments - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: q !! periapsis - real(DP), intent(out) :: capm !! mean anomaly - real(DP), intent(out) :: tperi !! time of pericenter passage + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: vx,vy,vz !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: q !! periapsis + real(DP), intent(out) :: capm !! mean anomaly + real(DP), intent(out) :: tperi !! time of pericenter passage ! Internals integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, rdotv, energy, fac, w, face, cape, e, tmpf, capf, mm - real(DP), dimension(NDIM) :: hvec + real(DP), dimension(NDIM) :: hvec, x, v a = 0.0_DP q = 0.0_DP capm = 0.0_DP tperi = 0.0_DP + x = [px, py, pz] + v = [vx, vy, vz] r = sqrt(dot_product(x(:), x(:))) v2 = dot_product(v(:), v(:)) hvec(:) = x(:) .cross. v(:) @@ -888,13 +906,16 @@ module subroutine orbel_xv2el_vec(self, cb) if (allocated(self%omega)) deallocate(self%omega); allocate(self%omega(self%nbody)) if (allocated(self%capm)) deallocate(self%capm); allocate(self%capm(self%nbody)) do concurrent (i = 1:self%nbody) - call orbel_xv2el(self%mu(i), self%xh(:, i), self%vh(:, i), self%a(i), self%e(i), self%inc(i), & + call orbel_xv2el(self%mu(i), self%xh(1,i), self%xh(2,i), self%xh(3,i), & + self%vh(1,i), self%vh(2,i), self%vh(3,i), & + self%a(i), self%e(i), self%inc(i), & self%capom(i), self%omega(i), self%capm(i)) end do end subroutine orbel_xv2el_vec - module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) + pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm) + !$omp declare simd(orbel_xv2el) !! author: David A. Minton !! !! Compute osculating orbital elements from relative Cartesian position and velocity @@ -911,19 +932,19 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) !! Adapted from Martin Duncan's Swift routine orbel_xv2el.f implicit none ! Arguments - real(DP), intent(in) :: mu !! Gravitational constant - real(DP), dimension(:), intent(in) :: x !! Position vector - real(DP), dimension(:), intent(in) :: v !! Velocity vector - real(DP), intent(out) :: a !! semimajor axis - real(DP), intent(out) :: e !! eccentricity - real(DP), intent(out) :: inc !! inclination - real(DP), intent(out) :: capom !! longitude of ascending node - real(DP), intent(out) :: omega !! argument of periapsis - real(DP), intent(out) :: capm !! mean anomaly + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), intent(in) :: px,py,pz !! Position vector + real(DP), intent(in) :: vx,vy,vz !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: e !! eccentricity + real(DP), intent(out) :: inc !! inclination + real(DP), intent(out) :: capom !! longitude of ascending node + real(DP), intent(out) :: omega !! argument of periapsis + real(DP), intent(out) :: capm !! mean anomaly ! Internals integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, cape, tmpf, capf - real(DP), dimension(NDIM) :: hvec + real(DP), dimension(NDIM) :: hvec, x, v a = 0.0_DP e = 0.0_DP @@ -931,6 +952,8 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) capom = 0.0_DP omega = 0.0_DP capm = 0.0_DP + x = [px, py, pz] + v = [vx, vy, vz] r = .mag. x(:) v2 = dot_product(v(:), v(:)) hvec = x(:) .cross. v(:) @@ -947,11 +970,11 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) end if fac = sqrt(hvec(1)**2 + hvec(2)**2) / h if (fac**2 < VSMALL) then - u = atan2(x(2), x(1)) + u = atan2(py, px) if (hvec(3) < 0.0_DP) u = -u else capom = atan2(hvec(1), -hvec(2)) - u = atan2(x(3) / sin(inc), x(1) * cos(capom) + x(2) * sin(capom)) + u = atan2(pz / sin(inc), px * cos(capom) + py * sin(capom)) end if if (capom < 0.0_DP) capom = capom + TWOPI if (u < 0.0_DP) u = u + TWOPI diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 00b66a9f8..005b6f68e 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -542,7 +542,8 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) if (tp%isperi(i) == -1) then if (vdotr >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aqt(mu, xpc(:, i), vpc(:, i), a, peri, capm, tperi) + call orbel_xv2aqt(mu, xpc(1,i), xpc(2,i), xpc(3,i), vpc(1,i), vpc(2,i), vpc(3,i), & + a, peri, capm, tperi) r2 = dot_product(xpc(:, i), xpc(:, i)) if ((abs(tperi) > FACQDT * dt) .or. (r2 > rhill2)) peri = sqrt(r2) if (param%enc_out /= "") then diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index baa51485a..5860e4cea 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -412,7 +412,7 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt tcr2 = r2 / (vxr**2 + vyr**2 + vzr**2) dt2 = dt**2 if (tcr2 <= dt2) then - call orbel_xv2aeq(Gmtot, [xr, yr, zr], [vxr, vyr, vzr], a, e, q) + call orbel_xv2aeq(Gmtot, xr, yr, zr, vxr, vyr, vzr, a, e, q) lcollision = (q < rlim) end if end if diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 8065342c6..86250272e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -363,7 +363,8 @@ module subroutine symba_util_peri_pl(self, system, param) if (pl%isperi(i) == -1) then if (vdotr >= 0.0_DP) then pl%isperi(i) = 0 - CALL orbel_xv2aeq(pl%mu(i), pl%xh(:,i), pl%vh(:,i), pl%atp(i), e, pl%peri(i)) + CALL orbel_xv2aeq(pl%mu(i), pl%xh(1,i), pl%xh(2,i), pl%xh(3,i), pl%vh(1,i), pl%vh(2,i), pl%vh(3,i), & + pl%atp(i), e, pl%peri(i)) end if else if (vdotr > 0.0_DP) then @@ -381,7 +382,7 @@ module subroutine symba_util_peri_pl(self, system, param) if (pl%isperi(i) == -1) then if (vdotr >= 0.0_DP) then pl%isperi(i) = 0 - CALL orbel_xv2aeq(system%Gmtot, pl%xb(:,i), pl%vb(:,i), pl%atp(i), e, pl%peri(i)) + CALL orbel_xv2aeq(system%Gmtot, pl%xb(1,i), pl%xb(2,i), pl%xb(3,i), pl%vb(1,i), pl%vb(2,i), pl%vb(3,i), pl%atp(i), e, pl%peri(i)) end if else if (vdotr > 0.0_DP) then diff --git a/src/util/util_peri.f90 b/src/util/util_peri.f90 index 66f2254e1..a81748d74 100644 --- a/src/util/util_peri.f90 +++ b/src/util/util_peri.f90 @@ -29,7 +29,7 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(tp%mu(i), tp%xh(:, i), tp%vh(:, i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(tp%mu(i), tp%xh(1,i), tp%xh(2,i), tp%xh(3,i), tp%vh(1,i), tp%vh(2,i), tp%vh(3,i), tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then @@ -45,7 +45,7 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(system%Gmtot, tp%xb(:, i), tp%vb(:, i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(system%Gmtot, tp%xb(1,i), tp%xb(2,i), tp%xb(3,i), tp%vb(1,i), tp%vb(2,i), tp%vb(3,i), tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then