From 7fa96eb6b2cfcbe8f5e81046e20aa5d35bd8dee7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 20:43:30 -0400 Subject: [PATCH 01/12] Enabled conservation report --- .../1pl_1pl_encounter/init_cond.py | 1 + .../1pl_1pl_encounter/param.swiftest.in | 1 + src/io/io.f90 | 9 +++------ src/modules/swiftest_classes.f90 | 4 ++-- src/modules/symba_classes.f90 | 2 +- src/symba/symba_io.f90 | 2 +- 6 files changed, 9 insertions(+), 10 deletions(-) diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py index ece9101e0..245f5fae0 100755 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py @@ -174,6 +174,7 @@ print(f'BIG_DISCARD no') print(f'DISCARD_OUT discard.swiftest.out') print(f'ROTATION no') +print(f'ENERGY yes') print(f'GR no') print(f'MU2KG {MU2KG}') print(f'DU2M {DU2M}') diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in index d44f4df0e..3050dea4a 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in @@ -24,6 +24,7 @@ EXTRA_FORCE no BIG_DISCARD no DISCARD_OUT discard.swiftest.out ROTATION no +ENERGY yes GR no MU2KG 1.988409870698051e+30 DU2M 149597870700.0 diff --git a/src/io/io.f90 b/src/io/io.f90 index e159d019d..360088980 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -70,10 +70,6 @@ module subroutine io_conservation_report(self, param, lterminal) write(*,*) 'dpe : ',(pe_now - pe_last) / abs(Eorbit_orig) write(*,*) end if - if (Lerror > 1e-6) then - write(*,*) 'Something has gone wrong! Angular momentum is too high!' - write(*,*) 'Lerror: ', Lerror - end if end if ke_orbit_last = ke_orbit_now ke_spin_last = ke_spin_now @@ -135,7 +131,7 @@ module subroutine io_dump_swiftest(self, param, msg) implicit none ! Arguments class(swiftest_base), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation ! Internals integer(I4B) :: ierr !! Error code @@ -173,7 +169,7 @@ module subroutine io_dump_system(self, param, msg) implicit none ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation ! Internals class(swiftest_parameters), allocatable :: dump_param !! Local parameters variable used to parameters change input file names @@ -203,6 +199,7 @@ module subroutine io_dump_system(self, param, msg) ! Print the status message (format code passed in from main driver) tfrac = (param%t - param%t0) / (param%tstop - param%t0) write(*,msg) param%t, tfrac, self%pl%nbody, self%tp%nbody + if (param%lenergy) call self%conservation_report(param, lterminal=.true.) return end subroutine io_dump_system diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index dcff0f6d8..47cfc6dd1 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -516,14 +516,14 @@ end subroutine io_dump_param module subroutine io_dump_swiftest(self, param, msg) implicit none class(swiftest_base), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation end subroutine io_dump_swiftest module subroutine io_dump_system(self, param, msg) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation end subroutine io_dump_system diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index f920fff87..697356b44 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -296,7 +296,7 @@ module subroutine symba_io_dump_particle_info(self, param, msg) use swiftest_classes, only : swiftest_parameters implicit none class(symba_particle_info), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation end subroutine symba_io_dump_particle_info diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 1f8626242..2e568dd7e 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -8,7 +8,7 @@ module subroutine symba_io_dump_particle_info(self, param, msg) !! Dumps the particle information data to a file implicit none class(symba_particle_info), intent(inout) :: self !! Swiftest base object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), optional, intent(in) :: msg !! Message to display with dump operation end subroutine symba_io_dump_particle_info From aee2941a259df1eaa2bf7fa11bdc34bb11c85f04 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 22:42:46 -0400 Subject: [PATCH 02/12] Added in initial conditions values to parameter file reader/writer --- src/io/io.f90 | 529 +++++++++++++++++-------------- src/modules/swiftest_classes.f90 | 5 + 2 files changed, 301 insertions(+), 233 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 360088980..10bb911f8 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -332,270 +332,315 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) class(swiftest_parameters), intent(inout) :: self !! Collection of parameters integer, intent(in) :: unit !! File unit number character(len=*), intent(in) :: iotype !! Dummy argument passed to the input/output procedure contains the text from the char-literal-constant, prefixed with DT. - !! If you do not include a char-literal-constant, the iotype argument contains only DT. + !! If you do not include a char-literal-constant, the iotype argument contains only DT. integer, intent(in) :: v_list(:) !! The first element passes the integrator code to the reader integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals - logical :: t0_set = .false. !! Is the initial time set in the input file? - logical :: tstop_set = .false. !! Is the final time set in the input file? - logical :: dt_set = .false. !! Is the step size set in the input file? - logical :: mtiny_set = .false. !! Is the mtiny value set? - integer(I4B) :: ilength, ifirst, ilast !! Variables used to parse input file - character(STRMAX) :: line !! Line of the input file + logical :: t0_set = .false. !! Is the initial time set in the input file? + logical :: tstop_set = .false. !! Is the final time set in the input file? + logical :: dt_set = .false. !! Is the step size set in the input file? + integer(I4B) :: ilength, ifirst, ilast, i !! Variables used to parse input file + character(STRMAX) :: line !! Line of the input file character (len=:), allocatable :: line_trim,param_name, param_value !! Strings used to parse the param file - character(*),parameter :: linefmt = '(A)' !! Format code for simple text string + character(*),parameter :: linefmt = '(A)' !! Format code for simple text string ! Parse the file line by line, extracting tokens then matching them up with known parameters if possible - - do - read(unit = unit, fmt = linefmt, iostat = iostat, end = 1) line - line_trim = trim(adjustl(line)) - ilength = len(line_trim) - if ((ilength /= 0)) then - ifirst = 1 - ! Read the pair of tokens. The first one is the parameter name, the second is the value. - param_name = io_get_token(line_trim, ifirst, ilast, iostat) - if (param_name == '') cycle ! No parameter name (usually because this line is commented out) - call io_toupper(param_name) - ifirst = ilast + 1 - param_value = io_get_token(line_trim, ifirst, ilast, iostat) - select case (param_name) - case ("T0") - read(param_value, *) self%t0 - t0_set = .true. - case ("TSTOP") - read(param_value, *) self%tstop - tstop_set = .true. - case ("DT") - read(param_value, *) self%dt - dt_set = .true. - case ("CB_IN") - self%incbfile = param_value - case ("PL_IN") - self%inplfile = param_value - case ("TP_IN") - self%intpfile = param_value - case ("IN_TYPE") - call io_toupper(param_value) - self%in_type = param_value - case ("ISTEP_OUT") - read(param_value, *) self%istep_out - case ("BIN_OUT") - self%outfile = param_value - case ("OUT_TYPE") - call io_toupper(param_value) - self%out_type = param_value - case ("OUT_FORM") - call io_toupper(param_value) - self%out_form = param_value - case ("OUT_STAT") - call io_toupper(param_value) - self%out_stat = param_value - case ("ISTEP_DUMP") - read(param_value, *) self%istep_dump - case ("CHK_CLOSE") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lclose = .true. - case ("CHK_RMIN") - read(param_value, *) self%rmin - case ("CHK_RMAX") - read(param_value, *) self%rmax - case ("CHK_EJECT") - read(param_value, *) self%rmaxu - case ("CHK_QMIN") - read(param_value, *) self%qmin - case ("CHK_QMIN_COORD") - call io_toupper(param_value) - self%qmin_coord = param_value - case ("CHK_QMIN_RANGE") - read(param_value, *) self%qmin_alo + associate(param => self) + do + read(unit = unit, fmt = linefmt, iostat = iostat, end = 1) line + line_trim = trim(adjustl(line)) + ilength = len(line_trim) + if ((ilength /= 0)) then + ifirst = 1 + ! Read the pair of tokens. The first one is the parameter name, the second is the value. + param_name = io_get_token(line_trim, ifirst, ilast, iostat) + if (param_name == '') cycle ! No parameter name (usually because this line is commented out) + call io_toupper(param_name) ifirst = ilast + 1 - param_value = io_get_token(line, ifirst, ilast, iostat) - read(param_value, *) self%qmin_ahi - case ("ENC_OUT") - self%enc_out = param_value - case ("DISCARD_OUT") - self%discard_out = param_value - case ("EXTRA_FORCE") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lextra_force = .true. - case ("BIG_DISCARD") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T' ) self%lbig_discard = .true. - case ("RHILL_PRESENT") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T' ) self%lrhill_present = .true. - case ("MU2KG") - read(param_value, *) self%MU2KG - case ("TU2S") - read(param_value, *) self%TU2S - case ("DU2M") - read(param_value, *) self%DU2M - case ("ENERGY") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lenergy = .true. - case ("GR") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lgr = .true. - case ("ROTATION") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%lrotation = .true. - case ("TIDES") - call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') self%ltides = .true. - case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters - case default - write(iomsg,*) "Unknown parameter -> ",param_name - iostat = -1 - return - end select - end if - end do - 1 continue - iostat = 0 + param_value = io_get_token(line_trim, ifirst, ilast, iostat) + select case (param_name) + case ("T0") + read(param_value, *) param%t0 + t0_set = .true. + case ("TSTOP") + read(param_value, *) param%tstop + tstop_set = .true. + case ("DT") + read(param_value, *) param%dt + dt_set = .true. + case ("CB_IN") + param%incbfile = param_value + case ("PL_IN") + param%inplfile = param_value + case ("TP_IN") + param%intpfile = param_value + case ("IN_TYPE") + call io_toupper(param_value) + param%in_type = param_value + case ("ISTEP_OUT") + read(param_value, *) param%istep_out + case ("BIN_OUT") + param%outfile = param_value + case ("OUT_TYPE") + call io_toupper(param_value) + param%out_type = param_value + case ("OUT_FORM") + call io_toupper(param_value) + param%out_form = param_value + case ("OUT_STAT") + call io_toupper(param_value) + param%out_stat = param_value + case ("ISTEP_DUMP") + read(param_value, *) param%istep_dump + case ("CHK_CLOSE") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') param%lclose = .true. + case ("CHK_RMIN") + read(param_value, *) param%rmin + case ("CHK_RMAX") + read(param_value, *) param%rmax + case ("CHK_EJECT") + read(param_value, *) param%rmaxu + case ("CHK_QMIN") + read(param_value, *) param%qmin + case ("CHK_QMIN_COORD") + call io_toupper(param_value) + param%qmin_coord = param_value + case ("CHK_QMIN_RANGE") + read(param_value, *) param%qmin_alo + ifirst = ilast + 1 + param_value = io_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%qmin_ahi + case ("ENC_OUT") + param%enc_out = param_value + case ("DISCARD_OUT") + param%discard_out = param_value + case ("EXTRA_FORCE") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') param%lextra_force = .true. + case ("BIG_DISCARD") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T' ) param%lbig_discard = .true. + case ("RHILL_PRESENT") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T' ) param%lrhill_present = .true. + case ("MU2KG") + read(param_value, *) param%MU2KG + case ("TU2S") + read(param_value, *) param%TU2S + case ("DU2M") + read(param_value, *) param%DU2M + case ("ENERGY") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') param%lenergy = .true. + case ("GR") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') param%lgr = .true. + case ("ROTATION") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') param%lrotation = .true. + case ("TIDES") + call io_toupper(param_value) + if (param_value == "YES" .or. param_value == 'T') param%ltides = .true. + case ("FIRSTKICK") + call io_toupper(param_value) + if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. + case ("FIRSTENERGY") + call io_toupper(param_value) + if (param_value == "NO" .or. param_value == 'F') param%lfirstenergy = .false. + case("EORBIT_ORIG") + read(param_value, *) param%Eorbit_orig + case("MTOT_ORIG") + read(param_value, *) param%Mtot_orig + case("LTOT_ORIG") + read(param_value, *) param%Ltot_orig(1) + do i = 2, NDIM + ifirst = ilast + 1 + param_value = io_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%Ltot_orig(i) + end do + param%Lmag_orig = norm2(param%Ltot_orig(:)) + case("LORBIT_ORIG") + read(param_value, *) param%Lorbit_orig(1) + do i = 2, NDIM + ifirst = ilast + 1 + param_value = io_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%Lorbit_orig(i) + end do + case("LSPIN_ORIG") + read(param_value, *) param%Lspin_orig(1) + do i = 2, NDIM + ifirst = ilast + 1 + param_value = io_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%Lspin_orig(i) + end do + case("LESCAPE") + read(param_value, *) param%Lescape(1) + do i = 2, NDIM + ifirst = ilast + 1 + param_value = io_get_token(line, ifirst, ilast, iostat) + read(param_value, *) param%Lescape(i) + end do + case("MESCAPE") + read(param_value, *) param%Mescape + case("ECOLLISIONS") + read(param_value, *) param%Ecollisions + case("EUNTRACKED") + read(param_value, *) param%Euntracked + case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters + case default + write(iomsg,*) "Unknown parameter -> ",param_name + iostat = -1 + return + end select + end if + end do + 1 continue + iostat = 0 - !! 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 - return - end if - if (self%dt <= 0.0_DP) then - write(iomsg,*) 'Invalid timestep: ' - iostat = -1 - return - end if - if (self%inplfile == "") then - write(iomsg,*) 'No valid massive body file in input file' - iostat = -1 - return - end if - if ((self%in_type /= REAL8_TYPE) .and. (self%in_type /= "ASCII")) then - write(iomsg,*) 'Invalid input file type:',trim(adjustl(self%in_type)) - iostat = -1 - return - end if - if ((self%istep_out <= 0) .and. (self%istep_dump <= 0)) then - write(iomsg,*) 'Invalid istep' - iostat = -1 - return - end if - if ((self%istep_out > 0) .and. (self%outfile == "")) then - write(iomsg,*) 'Invalid outfile' - iostat = -1 - return - end if - if (self%outfile /= "") then - if ((self%out_type /= REAL4_TYPE) .and. (self%out_type /= REAL8_TYPE) .and. & - (self%out_type /= SWIFTER_REAL4_TYPE) .and. (self%out_type /= SWIFTER_REAL8_TYPE)) then - write(iomsg,*) 'Invalid out_type: ',trim(adjustl(self%out_type)) + !! 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 return end if - if ((self%out_form /= "EL") .and. (self%out_form /= "XV")) then - write(iomsg,*) 'Invalid out_form: ',trim(adjustl(self%out_form)) + if (param%dt <= 0.0_DP) then + write(iomsg,*) 'Invalid timestep: ' iostat = -1 return end if - if ((self%out_stat /= "NEW") .and. (self%out_stat /= "REPLACE") .and. (self%out_stat /= "APPEND") .and. (self%out_stat /= "UNKNOWN")) then - write(iomsg,*) 'Invalid out_stat: ',trim(adjustl(self%out_stat)) + if (param%inplfile == "") then + write(iomsg,*) 'No valid massive body file in input file' iostat = -1 return end if - end if - if (self%qmin > 0.0_DP) then - if ((self%qmin_coord /= "HELIO") .and. (self%qmin_coord /= "BARY")) then - write(iomsg,*) 'Invalid qmin_coord: ',trim(adjustl(self%qmin_coord)) + if ((param%in_type /= REAL8_TYPE) .and. (param%in_type /= "ASCII")) then + write(iomsg,*) 'Invalid input file type:',trim(adjustl(param%in_type)) iostat = -1 return end if - if ((self%qmin_alo <= 0.0_DP) .or. (self%qmin_ahi <= 0.0_DP)) then - write(iomsg,*) 'Invalid qmin vals' + if ((param%istep_out <= 0) .and. (param%istep_dump <= 0)) then + write(iomsg,*) 'Invalid istep' + iostat = -1 + return + end if + if ((param%istep_out > 0) .and. (param%outfile == "")) then + write(iomsg,*) 'Invalid outfile' + iostat = -1 + return + end if + if (param%outfile /= "") then + if ((param%out_type /= REAL4_TYPE) .and. (param%out_type /= REAL8_TYPE) .and. & + (param%out_type /= SWIFTER_REAL4_TYPE) .and. (param%out_type /= SWIFTER_REAL8_TYPE)) then + write(iomsg,*) 'Invalid out_type: ',trim(adjustl(param%out_type)) + iostat = -1 + return + end if + if ((param%out_form /= "EL") .and. (param%out_form /= "XV")) then + write(iomsg,*) 'Invalid out_form: ',trim(adjustl(param%out_form)) + iostat = -1 + return + end if + if ((param%out_stat /= "NEW") .and. (param%out_stat /= "REPLACE") .and. (param%out_stat /= "APPEND") .and. (param%out_stat /= "UNKNOWN")) then + write(iomsg,*) 'Invalid out_stat: ',trim(adjustl(param%out_stat)) + iostat = -1 + return + end if + end if + if (param%qmin > 0.0_DP) then + if ((param%qmin_coord /= "HELIO") .and. (param%qmin_coord /= "BARY")) then + write(iomsg,*) 'Invalid qmin_coord: ',trim(adjustl(param%qmin_coord)) + iostat = -1 + return + end if + if ((param%qmin_alo <= 0.0_DP) .or. (param%qmin_ahi <= 0.0_DP)) then + write(iomsg,*) 'Invalid qmin vals' + iostat = -1 + return + end if + end if + if (param%ltides .and. .not. param%lrotation) then + write(iomsg,*) 'Tides require rotation to be turned on' iostat = -1 return end if - end if - if (self%ltides .and. .not. self%lrotation) then - write(iomsg,*) 'Tides require rotation to be turned on' - iostat = -1 - return - end if - write(*,*) "T0 = ",self%t0 - write(*,*) "TSTOP = ",self%tstop - write(*,*) "DT = ",self%dt - write(*,*) "CB_IN = ",trim(adjustl(self%incbfile)) - write(*,*) "PL_IN = ",trim(adjustl(self%inplfile)) - write(*,*) "TP_IN = ",trim(adjustl(self%intpfile)) - write(*,*) "IN_TYPE = ",trim(adjustl(self%in_type)) - write(*,*) "ISTEP_OUT = ",self%istep_out - write(*,*) "BIN_OUT = ",trim(adjustl(self%outfile)) - write(*,*) "OUT_TYPE = ",trim(adjustl(self%out_type)) - write(*,*) "OUT_FORM = ",trim(adjustl(self%out_form)) - write(*,*) "OUT_STAT = ",trim(adjustl(self%out_stat)) - write(*,*) "ISTEP_DUMP = ",self%istep_dump - write(*,*) "CHK_CLOSE = ",self%lclose - write(*,*) "CHK_RMIN = ",self%rmin - write(*,*) "CHK_RMAX = ",self%rmax - write(*,*) "CHK_EJECT = ",self%rmaxu - write(*,*) "CHK_QMIN = ",self%qmin - write(*,*) "CHK_QMIN_COORD = ",trim(adjustl(self%qmin_coord)) - write(*,*) "CHK_QMIN_RANGE = ",self%qmin_alo, self%qmin_ahi - write(*,*) "EXTRA_FORCE = ",self%lextra_force - write(*,*) "RHILL_PRESENT = ",self%lrhill_present - write(*,*) "ROTATION = ", self%lrotation - write(*,*) "TIDES = ", self%ltides - write(*,*) "ENERGY = ",self%lenergy - write(*,*) "MU2KG = ",self%MU2KG - write(*,*) "TU2S = ",self%TU2S - write(*,*) "DU2M = ",self%DU2M - if (trim(adjustl(self%enc_out)) /= "") then - write(*,*) "ENC_OUT = ",trim(adjustl(self%enc_out)) - else - write(*,*) "! ENC_OUT not set: Encounters will not be recorded to file" - end if - if (trim(adjustl(self%discard_out)) /= "") then - write(*,*) "DISCARD_OUT = ",trim(adjustl(self%discard_out)) - write(*,*) "BIG_DISCARD = ",self%lbig_discard - else - write(*,*) "! DISCARD_OUT not set: Discards will not be recorded to file" - write(*,*) "! BIG_DISCARD = ",self%lbig_discard - 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(*,*) "OUT_STAT = ",trim(adjustl(param%out_stat)) + write(*,*) "ISTEP_DUMP = ",param%istep_dump + write(*,*) "CHK_CLOSE = ",param%lclose + write(*,*) "CHK_RMIN = ",param%rmin + write(*,*) "CHK_RMAX = ",param%rmax + write(*,*) "CHK_EJECT = ",param%rmaxu + write(*,*) "CHK_QMIN = ",param%qmin + write(*,*) "CHK_QMIN_COORD = ",trim(adjustl(param%qmin_coord)) + 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(*,*) "MU2KG = ",param%MU2KG + write(*,*) "TU2S = ",param%TU2S + write(*,*) "DU2M = ",param%DU2M + 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" + write(*,*) "! BIG_DISCARD = ",param%lbig_discard + end if - if ((self%MU2KG < 0.0_DP) .or. (self%TU2S < 0.0_DP) .or. (self%DU2M < 0.0_DP)) then - write(iomsg,*) 'Invalid unit conversion factor' - iostat = -1 - return - 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 + return + end if - ! Calculate the G for the system units - self%GU = GC / (self%DU2M**3 / (self%MU2KG * self%TU2S**2)) + ! Calculate the G for the system units + param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) - ! Calculate the inverse speed of light in the system units - self%inv_c2 = einsteinC * self%TU2S / self%DU2M - self%inv_c2 = (self%inv_c2)**(-2) + ! 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) - associate(integrator => v_list(1)) - if (integrator == RMVS) then - if (.not.self%lclose) then - write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.' - iostat = -1 - return + associate(integrator => v_list(1)) + if (integrator == RMVS) then + if (.not.param%lclose) then + write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.' + iostat = -1 + return + end if end if - end if - - ! Determine if the GR flag is set correctly for this integrator - select case(integrator) - case(WHM, RMVS, HELIO, SYMBA) - write(*,*) "GR = ", self%lgr - case default - if (self%lgr) write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.' - self%lgr = .false. - end select - end associate + + ! 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 + end associate - iostat = 0 + iostat = 0 + end associate return end subroutine io_param_reader @@ -670,6 +715,24 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, Afmt) "GR"; write(param_value, Lfmt) param%lgr; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, Afmt) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + + if (param%lenergy) then + write(param_name, Afmt) "FIRSTENERGY"; write(param_value, Lfmt) param%lfirstenergy; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "EORBIT_ORIG"; write(param_value, Rfmt) param%Eorbit_orig; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "MTOT_ORIG"; write(param_value, Rfmt) param%Mtot_orig; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(unit, '("LTOT_ORIG ",3(1X,ES25.17))') param%Ltot_orig(:) + write(unit, '("LORBIT_ORIG",3(1X,ES25.17))') param%Lorbit_orig(:) + write(unit, '("LSPIN_ORIG ",3(1X,ES25.17))') param%Lspin_orig(:) + write(unit, '("LESCAPE ",3(1X,ES25.17))') param%Lescape(:) + + write(param_name, Afmt) "MESCAPE"; write(param_value, Rfmt) param%Mescape; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "ECOLLISIONS"; write(param_value, Rfmt) param%Ecollisions; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "EUNTRACKED"; write(param_value, Rfmt) param%Euntracked; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + end if + write(param_name, Afmt) "FIRSTKICK"; write(param_value, Lfmt) param%lfirstkick; write(unit, Afmt) adjustl(param_name), adjustl(param_value) + + + iostat = 0 iomsg = "UDIO not implemented" end associate diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 47cfc6dd1..a4042eb27 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -64,6 +64,11 @@ module swiftest_classes real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector + real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector + real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) + real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) + real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions + real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies logical :: lfirstenergy = .true. !! This is the first time computing energe logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step From b4f1cb716590cfcdc94d0f2eb502403a69c11df5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 22:51:48 -0400 Subject: [PATCH 03/12] Streamlined the energy calc --- src/io/io.f90 | 7 ++++- src/modules/swiftest_classes.f90 | 13 ++++----- src/setup/setup.f90 | 1 + src/symba/symba_discard.f90 | 2 +- src/util/util_get_energy_momentum.f90 | 38 ++++++++++++--------------- 5 files changed, 30 insertions(+), 31 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 10bb911f8..d7b899475 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -39,7 +39,12 @@ module subroutine io_conservation_report(self, param, lterminal) write(EGYIU,EGYHEADER) end if end if - call system%get_energy_and_momentum(param, ke_orbit_now, ke_spin_now, pe_now, Lorbit_now, Lspin_now) + call system%get_energy_and_momentum(param) + ke_orbit_now = system%ke_orbit + ke_spin_now = system%ke_spin + pe_now = system%pe + Lorbit_now = system%Lorbit + Lspin_now = system%Lspin Eorbit_now = ke_orbit_now + ke_spin_now + pe_now Ltot_now(:) = Lorbit_now(:) + Lspin_now(:) + Lescape(:) Mtot_now = cb%mass + sum(pl%mass(1:npl)) + system%Mescape diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a4042eb27..10e60a491 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -276,10 +276,12 @@ module swiftest_classes class(swiftest_tp), allocatable :: tp !! Test particle data structure class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure real(DP) :: Gmtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion - real(DP) :: ke = 0.0_DP !! System kinetic energy + real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy + real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy real(DP) :: pe = 0.0_DP !! System potential energy real(DP) :: te = 0.0_DP !! System total energy - real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector + real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! System orbital angular momentum vector + real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! System spin angular momentum vector real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions @@ -1019,15 +1021,10 @@ module subroutine util_resize_tp(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine util_resize_tp - module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin, pe, Lorbit, Lspin) + module subroutine util_get_energy_momentum_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(out) :: ke_orbit !! Orbital kinetic energy - real(DP), intent(out) :: ke_spin !! Spin kinetic energy - real(DP), intent(out) :: pe !! Potential energy - real(DP), dimension(:), intent(out) :: Lorbit !! Orbital angular momentum - real(DP), dimension(:), intent(out) :: Lspin !! Spin angular momentum end subroutine util_get_energy_momentum_system module subroutine util_set_beg_end_pl(self, xbeg, xend, vbeg) diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index edb641907..33b5c07c4 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -133,6 +133,7 @@ module subroutine setup_initialize_system(self, param) call self%tp%set_mu(self%cb) call self%pl%eucl_index() if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) + !if (param%lfirstenergy) then return end subroutine setup_initialize_system diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 6ab835e36..5f6d3926a 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -112,7 +112,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) Ltot(:) = Ltot(:) - Lpl(:) end do Ltot(:) = Ltot(:) - cb%mass * cb%xb(:) .cross. cb%vb(:) - system%Lescape(:) = system%Lescape(:) + system%Ltot(:) + system%Lescape(:) = system%Lescape(:) + Ltot(:) if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) else diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 69936e1b2..38701229d 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -1,7 +1,7 @@ submodule (swiftest_classes) s_util_get_energy_momentum use swiftest contains - module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin, pe, Lorbit, Lspin) + module subroutine util_get_energy_momentum_system(self, param) !! author: David A. Minton !! !! Compute total system angular momentum vector and kinetic, potential and total system energy @@ -12,11 +12,6 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(out) :: ke_orbit !! Orbital kinetic energy - real(DP), intent(out) :: ke_spin !! Spin kinetic energy - real(DP), intent(out) :: pe !! Potential energy - real(DP), dimension(:), intent(out) :: Lorbit !! Orbital angular momentum - real(DP), dimension(:), intent(out) :: Lspin !! Spin angular momentum ! Internals integer(I4B) :: i, j integer(I8B) :: k @@ -28,11 +23,12 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin logical, dimension(self%pl%nplpl) :: lstatpl logical, dimension(self%pl%nbody) :: lstatus - Lorbit(:) = 0.0_DP - Lspin(:) = 0.0_DP - ke_orbit = 0.0_DP - ke_spin = 0.0_DP associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) + system%Lorbit(:) = 0.0_DP + system%Lspin(:) = 0.0_DP + system%ke_orbit = 0.0_DP + system%ke_spin = 0.0_DP + kepl(:) = 0.0_DP Lplorbitx(:) = 0.0_DP Lplorbity(:) = 0.0_DP @@ -56,6 +52,7 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin ! Kinetic energy from orbit and spin kepl(i) = pl%mass(i) * v2 end do + if (param%lrotation) then do i = 1, npl rot2 = dot_product(pl%rot(:,i), pl%rot(:,i)) @@ -90,10 +87,10 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin end associate end do - ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) - if (param%lrotation) ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) + system%ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) + if (param%lrotation) system%ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) - pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl)) + system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl)) ! Potential energy from the oblateness term if (param%loblatecb) then @@ -102,17 +99,16 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin irh(i) = 1.0_DP / norm2(pl%xh(:,i)) end do call obl_pot(npl, cb%mass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) - pe = pe + oblpot + system%pe = system%pe + oblpot end if - Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) - Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl)) - Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl)) - - Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl)) - Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl)) - Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl)) + system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) + system%Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl)) + system%Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl)) + system%Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl)) + system%Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl)) + system%Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl)) end associate return From 4fc17fd9296985550b8ad219027ee1880259f4a1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 23:05:52 -0400 Subject: [PATCH 04/12] Added collisional regime determination code --- Makefile | 6 +++ src/fragmentation/fragmentation.f90 | 32 +++++++++++++++ src/modules/swiftest_classes.f90 | 9 +++++ src/symba/symba_collision.f90 | 62 +++++++++++++++++++++++++++++ 4 files changed, 109 insertions(+) create mode 100644 src/fragmentation/fragmentation.f90 diff --git a/Makefile b/Makefile index 63cfb0ee0..162129149 100644 --- a/Makefile +++ b/Makefile @@ -96,6 +96,11 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir + cd $(SWIFTEST_HOME)/src/fragmentation; \ + 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 .; \ @@ -189,6 +194,7 @@ clean: cd $(SWIFTEST_HOME)/src/discard; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/drift; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/eucl; rm -f Makefile.Defines Makefile *.gc* + cd $(SWIFTEST_HOME)/src/fragmentation; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/gr; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/helio; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/io; rm -f Makefile.Defines Makefile *.gc* diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 new file mode 100644 index 000000000..d809bc03a --- /dev/null +++ b/src/fragmentation/fragmentation.f90 @@ -0,0 +1,32 @@ +submodule(swiftest_classes) s_fragmentation + use swiftest +contains + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Determine the collisional regime of two colliding bodies. + !! Current version requires all values to be converted to SI units prior to calling the function + !! References: + !! Kokubo, E., Genda, H., 2010. Formation of Terrestrial Planets from Protoplanets Under a Realistic Accretion + !! Condition. ApJL 714, L21. https://doi.org/10.1088/2041-8205/714/1/L21 + !! Leinhardt, Z.M., Stewart, S.T., 2012. Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling + !! Laws 745, 79. https://doi.org/10.1088/0004-637X/745/1/79 + !! Mustill, A.J., Davies, M.B., Johansen, A., 2018. The dynamical evolution of transiting planetary systems including + !! a realistic collision prescription. Mon Not R Astron Soc 478, 2896–2908. https://doi.org/10.1093/mnras/sty1273 + !! Rufu, R., Aharonson, O., 2019. Impact Dynamics of Moons Within a Planetary Potential. J. Geophys. Res. Planets 124, + !! 1008–1019. https://doi.org/10.1029/2018JE005798 + !! Stewart, S.T., Leinhardt, Z.M., 2012. Collisions between Gravity-dominated Bodies. II. The Diversity of Impact + !! Outcomes during the End Stage of Planet Formation. ApJ 751, 32. https://doi.org/10.1088/0004-637X/751/1/32 + !! + implicit none + ! Arguments + integer(I4B), intent(out) :: regime + real(DP), intent(out) :: Mlr, Mslr + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 + real(DP), intent(out) :: Qloss !! The residual energy after the collision + + return + end subroutine fragmentation_regime + +end submodule s_fragmentation \ No newline at end of file diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 10e60a491..af751fdcf 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -452,6 +452,15 @@ module subroutine eucl_dist_index_plpl(self) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object end subroutine + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) + implicit none + integer(I4B), intent(out) :: regime + real(DP), intent(out) :: Mlr, Mslr + real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny + real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 + real(DP), intent(out) :: Qloss !! The residual energy after the collision + end subroutine fragmentation_regime + module pure subroutine gr_kick_getaccb_ns_body(self, system, param) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest generic body object diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index ad0e64079..843de7dd6 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -408,6 +408,68 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions ! Internals + ! Internals + integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + real(DP), dimension(NDIM,2) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(2) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + logical :: lgoodcollision + integer(I4B) :: i, status, jtarg, jproj, regime + real(DP), dimension(2) :: radius_si, mass_si, density_si + real(DP) :: mtiny_si, Mcb_si + real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si + real(DP) :: mlr, mslr, mtot, dentot, msys, msys_new, Qloss, impact_parameter + integer(I4B), parameter :: NRES = 3 !! Number of collisional product results + real(DP), dimension(NRES) :: mass_res + + associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, cb => system%cb) + select type(pl => system%pl) + class is (symba_pl) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip) + if (.not. lgoodcollision) cycle + if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved + + ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them + ! to symba_regime + if (mass(1) > mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + mass_si(:) = (mass(:)) * param%MU2KG !! The collective mass of the parent and its children + radius_si(:) = radius(:) * param%DU2M !! The collective radius of the parent and its children + x1_si(:) = plpl_collisions%x1(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) + v1_si(:) = plpl_collisions%v1(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) + x2_si(:) = plpl_collisions%x2(:,i) * param%DU2M !! The position of the parent from inside the step (at collision) + v2_si(:) = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S !! The velocity of the parent from inside the step (at collision) + density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children + Mcb_si = cb%mass * param%MU2KG + mtiny_si = (param%MTINY / param%GU) * param%MU2KG + + mass_res(:) = 0.0_DP + + mtot = sum(mass_si(:)) + dentot = sum(mass_si(:) * density_si(:)) / mtot + + !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime + call fragmentation_regime(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& + v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), regime, mlr, mslr, mtiny_si, Qloss) + + mass_res(1) = min(max(mlr, 0.0_DP), mtot) + mass_res(2) = min(max(mslr, 0.0_DP), mtot) + mass_res(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) + mass_res(:) = (mass_res(:) / param%MU2KG) * param%GU + Qloss = Qloss * (param%GU / param%MU2KG) * (param%TU2S / param%DU2M)**2 + + !status = symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + end do + end select + end associate return end subroutine symba_collision_resolve_fragmentations From 663ad8fb6568f5e5c4cabb3ee274dd517e0f6386 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 23:08:49 -0400 Subject: [PATCH 05/12] Added the rest of fragmentation_regime (formerally symba_regime from the Fragmentation branch --- src/fragmentation/fragmentation.f90 | 251 +++++++++++++++++++++++++++- 1 file changed, 249 insertions(+), 2 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index d809bc03a..a01af206b 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -19,14 +19,261 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v !! Outcomes during the End Stage of Planet Formation. ApJ 751, 32. https://doi.org/10.1088/0004-637X/751/1/32 !! implicit none - ! Arguments + ! Arguments integer(I4B), intent(out) :: regime real(DP), intent(out) :: Mlr, Mslr real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2 real(DP), intent(out) :: Qloss !! The residual energy after the collision + ! Constants + integer(I4B), parameter :: N1 = 1 !number of objects with mass equal to the largest remnant from LS12 + integer(I4B), parameter :: N2 = 2 !number of objects with mass larger than second largest remnant from LS12 + real(DP), parameter :: DENSITY1 = 1000.0_DP !standard density parameter from LS12 [kg/m3] + real(DP), parameter :: MU_BAR = 0.37_DP !0.385#0.37#0.3333# 3.978 # 1/3 material parameter for hydrodynamic planet-size bodies (LS12) + real(DP), parameter :: BETA = 2.85_DP !slope of sfd for remnants from LS12 2.85 + real(DP), parameter :: C1 = 2.43_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: C2 = -0.0408_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: C3 = 1.86_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: C4 = 1.08_DP !! Kokubo & Genda (2010) eq. (3) + real(DP), parameter :: CRUFU = 2.0_DP - 3 * MU_BAR ! central potential variable from Rufu and Aharonson (2019) + real(DP), parameter :: SUPERCAT_QRATIO = 1.8_DP ! See Section 4.1 of LS12 + ! Internals + real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta + real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat + real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat + real(DP) :: Mint, Mtot + real(DP) :: Rp, rhill + real(DP) :: Mresidual + real(DP) :: U_binding + + Vimp = norm2(vb2(:) - vb1(:)) + b = calc_b(xh2, vb2, xh1, vb1) + l = (rad1 + rad2) * (1 - b) + egy = 0.5_DP * dot_product(vb1, vb1) - GC * Mcb / norm2(xh1) + a1 = - GC * Mcb / 2.0_DP / egy + Mtot = m1 + m2 + mu = (m1 * m2) / Mtot + if (l < 2 * rad2) then + !calculate Mint + phi = 2 * acos((l - rad2) / rad2) + aint = rad2**2 * (PI - (phi - sin(phi)) / 2.0_DP) + lint = 2 * sqrt(rad2**2 - (rad2 - l / 2.0_DP) ** 2) + Mint = aint * lint ![kg] + alpha = (l**2) * (3 * rad2 - l) / (4 * (rad2**3)) + else + alpha = 1.0_DP + Mint = m2 + end if + Rp = (3 * (m1 / den1 + alpha * m2 / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! (Mustill et al. 2018) + c_star = calc_c_star(Rp) + !calculate Vescp + Vescp = sqrt(2 * GC * Mtot / Rp) !Mustill et al. 2018 eq 6 + !calculate rhill + rhill = a1 * (m1 / 3.0_DP / (Mcb + m1))**(1.0_DP/3.0_DP) + !calculate Vhill + if ((rad2 + rad1) < rhill) then + Vhill = sqrt(2 * GC * m1 * ((rhill**2 - rhill * (rad1 + rad2)) / & + (rhill**2 - 0.5_DP * (rad1 + rad2)**2)) / (rad1 + rad2)) + else + Vhill = Vescp + end if + !calculate Qr_pstar + Qrd_pstar = calc_Qrd_pstar(m1, m2, alpha, c_star) * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3) + !calculate Verosion + Qr_erosion = 2 * (1.0_DP - m1 / Mtot) * Qrd_pstar + Verosion = (2 * Qr_erosion * Mtot / mu)** (1.0_DP / 2.0_DP) + Qr = mu*(Vimp**2) / Mtot / 2.0_DP + !calculate mass largest remnant Mlr + Mlr = (1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot ! [kg] # LS12 eq (5) + !calculate Vsupercat + Qr_supercat = SUPERCAT_QRATIO * Qrd_pstar ! See LS12 Section 4.1 + Vsupercat = sqrt(2 * Qr_supercat * Mtot / mu) + !calculate Vhr + zeta = (m1 - m2) / Mtot + theta = 1.0_DP - b + Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3) + bcrit = rad1 / (rad1 + rad2) + Qloss = 0.0_DP + U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27 + + if ((m1 < mtiny).or.(m2 < mtiny)) then + regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + write(*,*) "FORCE MERGE" + else + if( Vimp < Vescp) then + regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else if (Vimp < Verosion) then + if (b < bcrit) then + regime = COLLRESOLVE_REGIME_MERGE !partial accretion regime" + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else if ((b > bcrit) .and. (Vimp < Vhr)) then + regime = COLLRESOLVE_REGIME_MERGE ! graze and merge + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else + Mlr = m1 + Mslr = calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star) + regime = COLLRESOLVE_REGIME_HIT_AND_RUN !hit and run + Qloss = (c_star + 1.0_DP) * U_binding ! Qr + end if + else if (Vimp > Verosion .and. Vimp < Vsupercat) then + if (m2 < 0.001_DP * m1) then + regime = COLLRESOLVE_REGIME_MERGE !cratering regime" + Mlr = Mtot + Mslr = 0.0_DP + Qloss = 0.0_DP + else + Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) ! LS12 eq (37) + regime = COLLRESOLVE_REGIME_DISRUPTION !disruption + Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_erosion + end if + else if (Vimp > Vsupercat) then + Mlr = Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(-1.5_DP) !LS12 eq (44) + Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA) !LS12 eq (37) + regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic + Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_supercat + else + write(*,*) "Error no regime found in symba_regime" + end if + end if + Mresidual = Mtot - Mlr - Mslr + if (Mresidual < 0.0_DP) then ! prevents final masses from going negative + Mlr = Mlr + Mresidual + end if + + return + + ! Internal functions + contains + function calc_Qrd_pstar(Mtarg, Mp, alpha, c_star) result(Qrd_pstar) + !! author: Jennifer L.L. Pouplin and Carlisle A. Wishard + !! + !! Calculates the corrected Q* for oblique impacts. See Eq. (15) of LS12. + !! Reference: + !! Leinhardt, Z.M., Stewart, S.T., 2012. Collisions between Gravity-dominated Bodies. I. Outcome Regimes and Scaling + !! Laws 745, 79. https://doi.org/10.1088/0004-637X/745/1/79 + !! + implicit none + ! Arguments + real(DP),intent(in) :: Mtarg, Mp, alpha, c_star + ! Result + real(DP) :: Qrd_pstar + ! Internals + real(DP) :: Qrd_star1, mu_alpha, mu, Qrd_star + + ! calc mu, mu_alpha + mu = (Mtarg * Mp) / (Mtarg + Mp) ! [kg] + mu_alpha = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) ! [kg] + ! calc Qrd_star1 + Qrd_star1 = (c_star * 4 * PI * DENSITY1 * GC * Rp**2) / 5.0_DP + ! calc Qrd_star + Qrd_star = Qrd_star1 * (((Mp / Mtarg + 1.0_DP)**2) / (4 * Mp / Mtarg))**(2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 23) + ! calc Qrd_pstar, v_pstar + Qrd_pstar = ((mu / mu_alpha)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) * Qrd_star ! (eq 15) + + return + end function calc_Qrd_pstar + + function calc_Qrd_rev(Mp, Mtarg, Mint, den1, den2, Vimp, c_star) result(Mslr) + !! author: Jennifer L.L. Pouplin and Carlisle A. Wishard + !! + !! Calculates mass of second largest fragment. + !! + implicit none + ! Arguments + real(DP),intent(in) :: Mp, Mtarg, Mint, den1, den2, Vimp, c_star + ! Result + real(DP) :: Mslr + ! Internals + real(DP) :: mtot_rev, mu_rev, gamma_rev, Qrd_star1, Qrd_star, mu_alpha_rev + real(DP) :: Qrd_pstar, Rc1, Qr_rev, Qrd_pstar_rev, Qr_supercat_rev + + ! calc Mslr, Rc1, mu, gammalr + mtot_rev = Mint + Mp + Rc1 = (3 * (Mint / den1 + Mp / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! [m] Mustill et al 2018 + mu_rev = (Mint * Mp) / mtot_rev ! [kg] eq 49 LS12 + mu_alpha_rev = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) + gamma_rev = Mint / Mp ! eq 50 LS12 + !calc Qr_rev + Qr_rev = mu_rev * (Vimp**2) / (2 * mtot_rev) + ! calc Qrd_star1, v_star1 + Qrd_star1 = (c_star * 4 * PI * mtot_rev * GC ) / Rc1 / 5.0_DP + ! calc Qrd_pstar_rev + Qrd_star = Qrd_star1 * (((gamma_rev + 1.0_DP)**2) / (4 * gamma_rev)) ** (2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 52) + Qrd_pstar = Qrd_star * ((mu_rev / mu_alpha_rev)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) + Qrd_pstar_rev = Qrd_pstar * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3) + !calc Qr_supercat_rev + Qr_supercat_rev = 1.8_DP * Qrd_pstar_rev + if (Qr_rev > Qr_supercat_rev ) then + Mslr = mtot_rev * (0.1_DP * ((Qr_rev / (Qrd_pstar_rev * 1.8_DP))**(-1.5_DP))) !eq (44) + else if ( Qr_rev < Qrd_pstar_rev ) then + Mslr = Mp + else + Mslr = (1.0_DP - Qr_rev / Qrd_pstar_rev / 2.0_DP) * (mtot_rev) ! [kg] #(eq 5) + end if + + if ( Mslr > Mp ) Mslr = Mp !check conservation of mass + + return + end function calc_Qrd_rev + + function calc_b(proj_pos, proj_vel, targ_pos, targ_vel) result(sintheta) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Calculates the impact factor b = sin(theta), where theta is the angle between the relative velocity + !! and distance vectors of the target and projectile bodies. See Fig. 2 of Leinhardt and Stewart (2012) + !! + implicit none + !! Arguments + real(DP), dimension(:), intent(in) :: proj_pos, proj_vel, targ_pos, targ_vel + !! Result + real(DP) :: sintheta + !! Internals + real(DP), dimension(NDIM) :: imp_vel, distance, x_cross_v + + imp_vel(:) = proj_vel(:) - targ_vel(:) + distance(:) = proj_pos(:) - targ_pos(:) + x_cross_v(:) = distance(:) .cross. imp_vel(:) + sintheta = norm2(x_cross_v(:)) / norm2(distance(:)) / norm2(imp_vel(:)) + return + end function calc_b + + function calc_c_star(Rc1) result(c_star) + !! author: David A. Minton + !! + !! Calculates c_star as a function of impact equivalent radius. It inteRpolates between 5 for ~1 km sized bodies to + !! 1.8 for ~10000 km sized bodies. See LS12 Fig. 4 for details. + !! + implicit none + !! Arguments + real(DP), intent(in) :: Rc1 + !! Result + real(DP) :: c_star + !! Internals + real(DP), parameter :: loR = 1.0e3_DP ! Lower bound of inteRpolation size (m) + real(DP), parameter :: hiR = 1.0e7_DP ! Upper bound of inteRpolation size (m) + real(DP), parameter :: loval = 5.0_DP ! Value of C* at lower bound + real(DP), parameter :: hival = 1.9_DP ! Value of C* at upper bound + + if (Rc1 < loR) then + c_star = loval + else if (Rc1 < hiR) then + c_star = loval + (hival - loval) * log(Rc1 / loR) / log(hiR /loR) + else + c_star = hival + end if + return + end function calc_c_star - return end subroutine fragmentation_regime end submodule s_fragmentation \ No newline at end of file From b274576c252266a66312e131cd3b7734c4b7cace Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 23:12:29 -0400 Subject: [PATCH 06/12] Added in placeholder select statement to switch between fragmenation regimes --- src/symba/symba_collision.f90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 843de7dd6..97f3e415d 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -466,7 +466,19 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) mass_res(:) = (mass_res(:) / param%MU2KG) * param%GU Qloss = Qloss * (param%GU / param%MU2KG) * (param%TU2S / param%DU2M)**2 - !status = symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + select case (regime) + case (COLLRESOLVE_REGIME_DISRUPTION) + !status = symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + !status = symba_fragmentation_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + case (COLLRESOLVE_REGIME_HIT_AND_RUN) + !status = symba_fragmentation_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + status = symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) + case default + write(*,*) "Error in symba_collision, unrecognized collision regime" + call util_exit(FAILURE) + end select end do end select end associate From 11f61d6bdfd484c32b63184f1a97ad834ad4be8a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 5 Aug 2021 06:46:11 -0400 Subject: [PATCH 07/12] Added interfaces for fragmentation cases --- src/modules/symba_classes.f90 | 36 ++++++++++++++++ src/symba/symba_collision.f90 | 8 ++-- src/symba/symba_fragmentation.f90 | 72 +++++++++++++++++++++++++++++++ 3 files changed, 112 insertions(+), 4 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 697356b44..0e66ebf7c 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -275,6 +275,30 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp + module function symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + implicit none + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation + real(DP), intent(inout) :: Qloss !! Energy lost during collisionn + integer(I4B) :: status !! Status flag assigned to this outcome + end function symba_fragmentation_casedisruption + + module function symba_fragmentation_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + implicit none + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation + real(DP), intent(inout) :: Qloss !! Energy lost during collision + integer(I4B) :: status !! Status flag assigned to this outcome + end function symba_fragmentation_casehitandrun + module function symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object @@ -285,6 +309,18 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, integer(I4B) :: status !! Status flag assigned to this outcome end function symba_fragmentation_casemerge + module function symba_fragmentation_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + implicit none + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation + real(DP), intent(inout) :: Qloss !! Energy lost during collision + integer(I4B) :: status !! Status flag assigned to this outcome + end function symba_fragmentation_casesupercatastrophic + module subroutine symba_io_write_discard(self, param) use swiftest_classes, only : swiftest_parameters implicit none diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 97f3e415d..952d59709 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -420,7 +420,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si real(DP) :: mlr, mslr, mtot, dentot, msys, msys_new, Qloss, impact_parameter integer(I4B), parameter :: NRES = 3 !! Number of collisional product results - real(DP), dimension(NRES) :: mass_res + real(DP), dimension(NRES) :: mass_res associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, cb => system%cb) select type(pl => system%pl) @@ -468,11 +468,11 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) select case (regime) case (COLLRESOLVE_REGIME_DISRUPTION) - !status = symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + status = symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) case (COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - !status = symba_fragmentation_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + status = symba_fragmentation_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) case (COLLRESOLVE_REGIME_HIT_AND_RUN) - !status = symba_fragmentation_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) + status = symba_fragmentation_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) status = symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) case default diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index f8afffb85..3ad475e4d 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -2,6 +2,54 @@ use swiftest contains + module function symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a non-catastrophic disruption collision + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation + real(DP), intent(inout) :: Qloss !! Energy lost during collision + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + + status = DISRUPTION + + return + end function symba_fragmentation_casedisruption + + + module function symba_fragmentation_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a non-catastrophic hit-and-run collision + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation + real(DP), intent(inout) :: Qloss !! Energy lost during collision + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + + status = HIT_AND_RUN + + return + end function symba_fragmentation_casehitandrun + + module function symba_fragmentation_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -134,4 +182,28 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, end function symba_fragmentation_casemerge + + module function symba_fragmentation_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a supercatastrophic collision + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(:), intent(inout) :: mass_res !! The distribution of fragment mass obtained by the regime calculation + real(DP), intent(inout) :: Qloss !! Energy lost during collision + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + + status = SUPERCATASTROPHIC + + return + end function symba_fragmentation_casesupercatastrophic + end submodule s_symba_fragmentation From 8a35ffcaf20ffbd228a20e203fc3b4529f00d9d0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 5 Aug 2021 07:08:14 -0400 Subject: [PATCH 08/12] Added disruption case (without fragment initialization yet) --- src/modules/swiftest_classes.f90 | 1 + src/setup/setup.f90 | 1 + src/symba/symba_fragmentation.f90 | 115 +++++++++++++++++++++++++++++- 3 files changed, 114 insertions(+), 3 deletions(-) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index af751fdcf..1aec7dc9d 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -289,6 +289,7 @@ module swiftest_classes logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated !! separately from massive bodies. Massive body variables are saved at half steps, and passed to !! the test particles + integer(I4B) :: maxid = -1 !! The current maximum particle id number contains !> Each integrator will have its own version of the step procedure(abstract_step_system), deferred :: step diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 33b5c07c4..6cba6d27b 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -128,6 +128,7 @@ module subroutine setup_initialize_system(self, param) call self%pl%initialize(param) call self%tp%initialize(param) call util_valid(self%pl, self%tp) + self%maxid = maxval([self%pl%id(:), self%tp%id(:)]) call self%set_msys() call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index 3ad475e4d..fea7dea91 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -1,5 +1,8 @@ submodule (symba_classes) s_symba_fragmentation use swiftest + + integer(I4B), parameter :: NFRAG_DISRUPT = 12 + integer(I4B), parameter :: NFRAG_SUPERCAT = 20 contains module function symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) @@ -19,8 +22,114 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals + integer(I4B) :: i, istart, nfrag, ibiggest, nstart, nend + real(DP) :: mtot, avg_dens + real(DP), dimension(NDIM) :: xcom, vcom, Ip_new + real(DP), dimension(2) :: vol + real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag + real(DP), dimension(:), allocatable :: m_frag, rad_frag + logical :: lfailure + class(symba_pl), allocatable :: plnew + + select type(pl => system%pl) + class is (symba_pl) + associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + nfrag = NFRAG_DISRUPT + allocate(m_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Get mass weighted mean of Ip and average density + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot + vol(:) = 4._DP / 3._DP * PI * radius(:)**3 + avg_dens = mtot / sum(vol(:)) + + ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment + m_frag(1) = mass_res(1) + if (mass_res(2) > mass_res(1) / 3._DP) then + m_frag(2) = mass_res(2) + istart = 3 + else + istart = 2 + end if + ! Distribute remaining mass among the remaining bodies + do i = istart, nfrag + m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) + end do + + ! Distribute any residual mass if there is any and set the radius + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + + do i = 1, nfrag + Ip_frag(:, i) = Ip_new(:) + end do + + !call fragmentation_initialize(pl, param, family, x, v, L_spin, Ip, mass, radius, & + ! nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 + else + ! Populate the list of new bodies + write(*,'("Generating ",I2.0," fragments")') nfrag + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + + plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) + do i = 1, nfrag + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%density(:) = avg_dens + plnew%radius(:) = rad_frag(:) + plnew%info(:)%origin_type = "Disruption" + plnew%info(:)%origin_time = param%t + do i = 1, nfrag + plnew%info(i)%origin_xh(:) = plnew%xh(:,i) + plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + end do + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + + ! Append the new merged body to the list and record how many we made + nstart = mergeadd_list%nbody + 1 + nend = mergeadd_list%nbody + plnew%nbody + call mergeadd_list%append(plnew) + mergeadd_list%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) + end if - status = DISRUPTION + end associate + end select return end function symba_fragmentation_casedisruption @@ -157,8 +266,8 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, plnew%radius(1) = radius_new plnew%info(1) = pl%info(family(ibiggest)) if (param%lrotation) then - pl%Ip(:,1) = Ip_new(:) - pl%rot(:,1) = rot_new(:) + plnew%Ip(:,1) = Ip_new(:) + plnew%rot(:,1) = rot_new(:) end if if (param%ltides) then plnew%Q = pl%Q(ibiggest) From ca1d66b92d3ec8df363a10fe5af2ff5fc15cd5da Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 5 Aug 2021 07:48:38 -0400 Subject: [PATCH 09/12] Added fragment initialization code (formerly symba_frag_pos from the Fragmentation branch) --- Makefile | 2 +- Makefile.Defines | 8 +- src/fragmentation/fragmentation.f90 | 835 ++++++++++++++++++++++++++++ src/modules/swiftest_classes.f90 | 44 ++ src/util/util_minimize_bfgs.f90 | 584 +++++++++++++++++++ src/util/util_solve.f90 | 138 ++++- 6 files changed, 1605 insertions(+), 6 deletions(-) create mode 100644 src/util/util_minimize_bfgs.f90 diff --git a/Makefile b/Makefile index 162129149..94dfbceeb 100644 --- a/Makefile +++ b/Makefile @@ -45,13 +45,13 @@ #****************************************************************************** SWIFTEST_MODULES = swiftest_globals.f90 \ + lambda_function.f90\ swiftest_classes.f90 \ swiftest_operators.f90 \ whm_classes.f90 \ rmvs_classes.f90 \ helio_classes.f90 \ symba_classes.f90 \ - lambda_function.f90\ swiftest.f90 diff --git a/Makefile.Defines b/Makefile.Defines index 291f2c604..70069bb71 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -65,13 +65,13 @@ GPAR = -fopenmp -ftree-parallelize-loops=4 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries -FFLAGS = $(IDEBUG) $(HEAPARR) +#FFLAGS = $(IDEBUG) $(HEAPARR) #FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) -FORTRAN = ifort +#FORTRAN = ifort #AR = xiar -#FORTRAN = gfortran -#FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) +FORTRAN = gfortran +FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index a01af206b..d1965914f 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -1,6 +1,841 @@ submodule(swiftest_classes) s_fragmentation use swiftest contains + + subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Initialize the position and velocity of fragments to conserve energy and momentum. + use, intrinsic :: ieee_exceptions + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: system + class(swiftest_parameters), intent(in) :: param + integer(I4B), dimension(:), intent(in) :: family + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius + integer(I4B), intent(inout) :: nfrag + real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag + real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag + real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag + logical, intent(out) :: lfailure ! Answers the question: Should this have been a merger instead? + real(DP), intent(inout) :: Qloss + ! Internals + real(DP) :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system + real(DP) :: mtot + real(DP), dimension(NDIM) :: xcom, vcom + integer(I4B) :: ii + logical, dimension(:), allocatable :: lexclude + real(DP), dimension(NDIM, 2) :: rot, L_orb + real(DP), dimension(:,:), allocatable :: x_frag, v_frag, v_r_unit, v_t_unit, v_h_unit + real(DP), dimension(:), allocatable :: rmag, rotmag, v_r_mag, v_t_mag + real(DP), dimension(NDIM) :: Ltot_before + real(DP), dimension(NDIM) :: Ltot_after + real(DP) :: Etot_before, ke_orbit_before, ke_spin_before, pe_before, Lmag_before + real(DP) :: Etot_after, ke_orbit_after, ke_spin_after, pe_after, Lmag_after, dEtot, dLmag + real(DP), dimension(NDIM) :: L_frag_tot, L_frag_orb + real(DP) :: ke_frag_budget, ke_frag_orbit, ke_radial, ke_frag_spin, ke_avg_deficit, ke_avg_deficit_old + real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit + character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))" + integer(I4B) :: try, subtry + integer(I4B), parameter :: NFRAG_MIN = 7 !! The minimum allowable number of fragments (set to 6 because that's how many unknowns are needed in the tangential velocity calculation) + real(DP) :: r_max_start, r_max_start_old, r_max, f_spin + real(DP), parameter :: Ltol = 10 * epsilon(1.0_DP) + real(DP), parameter :: Etol = 1e-10_DP + integer(I4B), parameter :: MAXTRY = 3000 + integer(I4B), parameter :: TANTRY = 3 + logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes + + if (nfrag < NFRAG_MIN) then + write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." + lfailure = .true. + return + end if + + call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + fpe_quiet_modes(:) = .false. + call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes) + + f_spin = 0.05_DP + mscale = 1.0_DP + rscale = 1.0_DP + vscale = 1.0_DP + tscale = 1.0_DP + Lscale = 1.0_DP + Escale = 1.0_DP + + associate(npl => system%pl%nbody, status => system%pl%status) + allocate(lexclude(npl)) + where (status(1:npl) == INACTIVE) ! Safety check in case one of the included bodies has been previously deactivated + lexclude(1:npl) = .true. + elsewhere + lexclude(1:npl) = .false. + end where + end associate + + allocate(x_frag, source=xb_frag) + allocate(v_frag, source=vb_frag) + + call set_scale_factors() + call define_coordinate_system() + call calculate_system_energy(linclude_fragments=.false.) + + r_max_start = norm2(x(:,2) - x(:,1)) + try = 1 + lfailure = .false. + ke_avg_deficit = 0.0_DP + do while (try < MAXTRY) + lfailure = .false. + ke_avg_deficit_old = ke_avg_deficit + ke_avg_deficit = 0.0_DP + subtry = 1 + do + call set_fragment_position_vectors() + call set_fragment_tan_vel(lfailure) + ke_avg_deficit = ke_avg_deficit - ke_radial + subtry = subtry + 1 + if (.not.lfailure .or. subtry == TANTRY) exit + !write(*,*) 'Trying new arrangement' + end do + ke_avg_deficit = ke_avg_deficit / subtry + if (lfailure) write(*,*) 'Failed to find tangential velocities' + + if (.not.lfailure) then + call set_fragment_radial_velocities(lfailure) + if (lfailure) write(*,*) 'Failed to find radial velocities' + if (.not.lfailure) then + call calculate_system_energy(linclude_fragments=.true.) + !write(*,*) 'Qloss : ',Qloss + !write(*,*) '-dEtot: ',-dEtot + !write(*,*) 'delta : ',abs((dEtot + Qloss)) + if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then + write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol + lfailure = .true. + else if (abs(dLmag) / Lmag_before > Ltol) then + write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before + lfailure = .true. + end if + end if + end if + + if (.not.lfailure) exit + call restructure_failed_fragments() + try = try + 1 + end do + write(*, "(' -------------------------------------------------------------------------------------')") + write(*, "(' Final diagnostic')") + write(*, "(' -------------------------------------------------------------------------------------')") + if (lfailure) then + write(*,*) "symba_frag_pos failed after: ",try," tries" + do ii = 1, nfrag + vb_frag(:, ii) = vcom(:) + end do + else + write(*,*) "symba_frag_pos succeeded after: ",try," tries" + write(*, "(' dL_tot should be very small' )") + write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before + write(*, "(' dE_tot should be negative and equal to Qloss' )") + write(*,fmtlabel) ' dE_tot |', dEtot / abs(Etot_before) + write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) + write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) + end if + write(*, "(' -------------------------------------------------------------------------------------')") + + call restore_scale_factors() + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + + return + + contains + + ! Because of the complexity of this procedure, we have chosen to break it up into a series of nested subroutines. + + subroutine set_scale_factors() + !! author: David A. Minton + !! + !! Scales dimenional quantities to ~O(1) with respect to the collisional system. This scaling makes it easier for the non-linear minimization + !! to converge on a solution + implicit none + integer(I4B) :: i + + ! Find the center of mass of the collisional system + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Set scale factors + !! Because of the implied G, mass is actually G*mass with units of distance**3 / time**2 + Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2))) + rscale = sum(radius(:)) + mscale = sqrt(Escale * rscale) + vscale = sqrt(Escale / mscale) + tscale = rscale / vscale + Lscale = mscale * rscale * vscale + + xcom(:) = xcom(:) / rscale + vcom(:) = vcom(:) / vscale + + mtot = mtot / mscale + mass = mass / mscale + radius = radius / rscale + x = x / rscale + v = v / vscale + L_spin = L_spin / Lscale + do i = 1, 2 + rot(:,i) = L_spin(:,i) / (mass(i) * radius(i)**2 * Ip(3, i)) + end do + + m_frag = m_frag / mscale + rad_frag = rad_frag / rscale + Qloss = Qloss / Escale + + return + end subroutine set_scale_factors + + subroutine restore_scale_factors() + !! author: David A. Minton + !! + !! Restores dimenional quantities back to the system units + implicit none + integer(I4B) :: i + + call ieee_set_halting_mode(IEEE_ALL,.false.) + ! Restore scale factors + xcom(:) = xcom(:) * rscale + vcom(:) = vcom(:) * vscale + + mtot = mtot * mscale + mass = mass * mscale + radius = radius * rscale + x = x * rscale + v = v * vscale + L_spin = L_spin * Lscale + do i = 1, 2 + rot(:,i) = L_spin(:,i) * (mass(i) * radius(i)**2 * Ip(3, i)) + end do + + m_frag = m_frag * mscale + rad_frag = rad_frag * rscale + rot_frag = rot_frag / tscale + x_frag = x_frag * rscale + v_frag = v_frag * vscale + Qloss = Qloss * Escale + + do i = 1, nfrag + xb_frag(:, i) = x_frag(:, i) + xcom(:) + vb_frag(:, i) = v_frag(:, i) + vcom(:) + end do + + Etot_before = Etot_before * Escale + pe_before = pe_before * Escale + ke_spin_before = ke_spin_before * Escale + ke_orbit_before = ke_orbit_before * Escale + Ltot_before = Ltot_before * Lscale + Lmag_before = Lmag_before * Lscale + Etot_after = Etot_after * Escale + pe_after = pe_after * Escale + ke_spin_after = ke_spin_after * Escale + ke_orbit_after = ke_orbit_after * Escale + Ltot_after = Ltot_after * Lscale + Lmag_after = Lmag_after * Lscale + + mscale = 1.0_DP + rscale = 1.0_DP + vscale = 1.0_DP + tscale = 1.0_DP + Lscale = 1.0_DP + Escale = 1.0_DP + + return + end subroutine restore_scale_factors + + subroutine define_coordinate_system() + !! author: David A. Minton + !! + !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. + implicit none + integer(I4B) :: i + real(DP), dimension(NDIM) :: x_cross_v, xc, vc, delta_r, delta_v + real(DP) :: r_col_norm, v_col_norm + + allocate(rmag(nfrag)) + allocate(rotmag(nfrag)) + allocate(v_r_mag(nfrag)) + allocate(v_t_mag(nfrag)) + allocate(v_r_unit(NDIM,nfrag)) + allocate(v_t_unit(NDIM,nfrag)) + allocate(v_h_unit(NDIM,nfrag)) + + rmag(:) = 0.0_DP + rotmag(:) = 0.0_DP + v_r_mag(:) = 0.0_DP + v_t_mag(:) = 0.0_DP + v_r_unit(:,:) = 0.0_DP + v_t_unit(:,:) = 0.0_DP + v_h_unit(:,:) = 0.0_DP + + L_orb(:, :) = 0.0_DP + ! Compute orbital angular momentum of pre-impact system + do i = 1, 2 + xc(:) = x(:, i) - xcom(:) + vc(:) = v(:, i) - vcom(:) + x_cross_v(:) = xc(:) .cross. vc(:) + L_orb(:, i) = mass(i) * x_cross_v(:) + end do + + ! Compute orbital angular momentum of pre-impact system. This will be the normal vector to the collision fragment plane + L_frag_tot(:) = L_spin(:, 1) + L_spin(:, 2) + L_orb(:, 1) + L_orb(:, 2) + + delta_v(:) = v(:, 2) - v(:, 1) + v_col_norm = norm2(delta_v(:)) + delta_r(:) = x(:, 2) - x(:, 1) + r_col_norm = norm2(delta_r(:)) + + ! We will initialize fragments on a plane defined by the pre-impact system, with the z-axis aligned with the angular momentum vector + ! and the y-axis aligned with the pre-impact distance vector. + y_col_unit(:) = delta_r(:) / r_col_norm + z_col_unit(:) = L_frag_tot(:) / norm2(L_frag_tot) + ! The cross product of the y- by z-axis will give us the x-axis + x_col_unit(:) = y_col_unit(:) .cross. z_col_unit(:) + + return + end subroutine define_coordinate_system + + subroutine calculate_system_energy(linclude_fragments) + !! Author: David A. Minton + !! + !! Calculates total system energy, including all bodies in the pl list that do not have a corresponding value of the lexclude array that is true + !! and optionally including fragments. + implicit none + ! Arguments + logical, intent(in) :: linclude_fragments + ! Internals + integer(I4B) :: i, npl_new, nplm + logical, dimension(:), allocatable :: ltmp + logical :: lk_plpl + class(swiftest_nbody_system), allocatable :: tmpsys + + ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the + ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by + ! allocating two of these ginormous arrays simulteouously. This is not particularly efficient, but as this + ! subroutine should be called relatively infrequently, it shouldn't matter too much. + !if (allocated(pl%k_plpl)) deallocate(pl%k_plpl) + + ! Build the internal planet list out of the non-excluded bodies and optionally with fragments appended. This + ! will get passed to the energy calculation subroutine so that energy is computed exactly the same way is it + ! is in the main program. + associate(pl => system%pl, npl => system%pl%nbody, cb => system%cb) + lk_plpl = allocated(pl%k_plpl) + if (lk_plpl) deallocate(pl%k_plpl) + if (linclude_fragments) then ! Temporarily expand the planet list to feed it into symba_energy + lexclude(family(:)) = .true. + npl_new = npl + nfrag + else + npl_new = npl + end if + call setup_construct_system(tmpsys, param) + deallocate(tmpsys%cb) + allocate(tmpsys%cb, source=cb) + allocate(ltmp(npl)) + ltmp(:) = .true. + call tmpsys%pl%setup(npl_new, param) + call tmpsys%pl%fill(pl, ltmp) + deallocate(ltmp) + + if (linclude_fragments) then ! Append the fragments if they are included + ! Energy calculation requires the fragments to be in the system barcyentric frame, s + tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag) + tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag) + tmpsys%pl%xb(:,npl+1:npl_new) = xb_frag(:,1:nfrag) + tmpsys%pl%vb(:,npl+1:npl_new) = vb_frag(:,1:nfrag) + tmpsys%pl%status(npl+1:npl_new) = COLLISION + if (param%lrotation) then + tmpsys%pl%Ip(:,npl+1:npl_new) = Ip_frag(:,1:nfrag) + tmpsys%pl%rot(:,npl+1:npl_new) = rot_frag(:,1:nfrag) + end if + call tmpsys%pl%b2h(tmpsys%cb) + allocate(ltmp(npl_new)) + ltmp(1:npl) = lexclude(1:npl) + ltmp(npl+1:npl_new) = .false. + call move_alloc(ltmp, lexclude) + end if + + where (lexclude(1:npl_new)) + tmpsys%pl%status(1:npl_new) = INACTIVE + end where + + select type(plwksp => tmpsys%pl) + class is (symba_pl) + select type(param) + class is (symba_parameters) + plwksp%nplm = count(plwksp%Gmass > param%mtiny / mscale) + end select + end select + call tmpsys%pl%eucl_index() + call tmpsys%get_energy_and_momentum(param) + + ! Restore the big array + deallocate(tmpsys%pl%k_plpl) + select type(pl) + class is (symba_pl) + select type(param) + class is (symba_parameters) + nplm = count(pl%mass > param%mtiny) + end select + end select + if (lk_plpl) call pl%eucl_index() + + ! Calculate the current fragment energy and momentum balances + if (linclude_fragments) then + Ltot_after(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) + Lmag_after = norm2(Ltot_after(:)) + ke_orbit_after = tmpsys%ke_orbit + ke_spin_after = tmpsys%ke_spin + pe_after = tmpsys%pe + Etot_after = ke_orbit_after + ke_spin_after + pe_after + dEtot = Etot_after - Etot_before + dLmag = norm2(Ltot_after(:) - Ltot_before(:)) + else + Ltot_before(:) = tmpsys%Lorbit(:) + tmpsys%Lspin(:) + Lmag_before = norm2(Ltot_before(:)) + ke_orbit_before = tmpsys%ke_orbit + ke_spin_before = tmpsys%ke_spin + pe_before = tmpsys%pe + Etot_before = ke_orbit_before + ke_spin_before + pe_before + end if + end associate + return + end subroutine calculate_system_energy + + subroutine shift_vector_to_origin(m_frag, vec_frag) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the position or velocity of the fragments as needed to align them with the origin + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame + + ! Internals + real(DP), dimension(NDIM) :: mvec_frag, COM_offset + integer(I4B) :: i + + mvec_frag(:) = 0.0_DP + + do i = 1, nfrag + mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i) + end do + COM_offset(:) = -mvec_frag(:) / mtot + do i = 1, nfrag + vec_frag(:, i) = vec_frag(:, i) + COM_offset(:) + end do + + return + end subroutine shift_vector_to_origin + + subroutine set_fragment_position_vectors() + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the + !! pre-impact angular momentum. They are distributed on an ellipse surrounding the center of mass. + !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. + + implicit none + real(DP) :: dis, rad + real(DP), dimension(NDIM) :: L_sigma + logical, dimension(:), allocatable :: loverlap + integer(I4B) :: i, j + + allocate(loverlap(nfrag)) + + ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies + ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) + r_max = r_max_start + rad = sum(radius(:)) + + ! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector. + ! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments. + + call random_number(x_frag(:,3:nfrag)) + loverlap(:) = .true. + do while (any(loverlap(3:nfrag))) + x_frag(:, 1) = x(:, 1) - xcom(:) + x_frag(:, 2) = x(:, 2) - xcom(:) + r_max = r_max + 0.1_DP * rad + do i = 3, nfrag + if (loverlap(i)) then + call random_number(x_frag(:,i)) + x_frag(:, i) = 2 * (x_frag(:, i) - 0.5_DP) * r_max + end if + end do + loverlap(:) = .false. + do j = 1, nfrag + do i = j + 1, nfrag + dis = norm2(x_frag(:,j) - x_frag(:,i)) + loverlap(i) = loverlap(i) .or. (dis <= (rad_frag(i) + rad_frag(j))) + end do + end do + end do + call shift_vector_to_origin(m_frag, x_frag) + + do i = 1, nfrag + rmag(i) = norm2(x_frag(:, i)) + v_r_unit(:, i) = x_frag(:, i) / rmag(i) + call random_number(L_sigma(:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, + ! otherwise we can get an ill-conditioned system + v_h_unit(:, i) = z_col_unit(:) + 2e-1_DP * (L_sigma(:) - 0.5_DP) + v_h_unit(:, i) = v_h_unit(:, i) / norm2(v_h_unit(:, i)) + v_t_unit(:, i) = v_h_unit(:, i) .cross. v_r_unit(:, i) + xb_frag(:,i) = x_frag(:,i) + xcom(:) + end do + + return + end subroutine set_fragment_position_vectors + + subroutine set_fragment_tan_vel(lerr) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the tangential velocities and spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. + !! This procedure works in several stages, with a goal to solve the angular and linear momentum constraints on the fragments, while still leaving a positive balance of + !! our fragment kinetic energy (ke_frag_budget) that we can put into the radial velocity distribution. + !! + !! The first thing we'll try to do is solve for the tangential velocities of the first 6 fragments, using angular and linear momentum as constraints and an initial + !! tangential velocity distribution for the remaining bodies (if there are any) that distributes their angular momentum equally between them. + !! If that doesn't work and we blow our kinetic energy budget, we will attempt to find a tangential velocity distribution that minimizes the kinetic energy while + !! conserving momentum. + !! + !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. + implicit none + ! Arguments + logical, intent(out) :: lerr + ! Internals + integer(I4B) :: i + real(DP), parameter :: TOL = 1e-4_DP + real(DP), dimension(:), allocatable :: v_t_initial + type(lambda_obj) :: spinfunc + type(lambda_obj_err) :: objective_function + real(DP), dimension(NDIM) :: L_frag_spin, L_remainder, Li, rot_L, rot_ke + + ! Initialize the fragments with 0 velocity and spin so we can divide up the balance between the tangential, radial, and spin components while conserving momentum + lerr = .false. + vb_frag(:,:) = 0.0_DP + rot_frag(:,:) = 0.0_DP + v_t_mag(:) = 0.0_DP + v_r_mag(:) = 0.0_DP + + call calculate_system_energy(linclude_fragments=.true.) + ke_frag_budget = -dEtot - Qloss + !write(*,*) '***************************************************' + !write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) + !write(*,*) 'r_max : ',r_max + !write(*,*) 'f_spin : ',f_spin + !write(*,*) '***************************************************' + !write(*,*) 'Energy balance so far: ' + !write(*,*) 'ke_frag_budget : ',ke_frag_budget + !write(*,*) 'ke_orbit_before: ',ke_orbit_before + !write(*,*) 'ke_orbit_after : ',ke_orbit_after + !write(*,*) 'ke_spin_before : ',ke_spin_before + !write(*,*) 'ke_spin_after : ',ke_spin_after + !write(*,*) 'pe_before : ',pe_before + !write(*,*) 'pe_after : ',pe_after + !write(*,*) 'Qloss : ',Qloss + !write(*,*) '***************************************************' + if (ke_frag_budget < 0.0_DP) then + write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget + r_max_start = r_max_start / 2 + lerr = .true. + return + end if + + allocate(v_t_initial, mold=v_t_mag) + + L_frag_spin(:) = 0.0_DP + ke_frag_spin = 0.0_DP + ! Start the first two bodies with the same rotation as the original two impactors, then distribute the remaining angular momentum among the rest + do i = 1, 2 + rot_frag(:, i) = rot(:, i) + L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i) * rot_frag(:, i) + end do + L_frag_orb(:) = L_frag_tot(:) - L_frag_spin(:) + L_frag_spin(:) = 0.0_DP + do i = 1, nfrag + ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets + rot_ke(:) = sqrt(2 * f_spin * ke_frag_budget / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i))) * L_frag_orb(:) / norm2(L_frag_orb(:)) + rot_L(:) = f_spin * L_frag_orb(:) / (nfrag * m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i)) + if (norm2(rot_ke) < norm2(rot_L)) then + rot_frag(:,i) = rot_frag(:, i) + rot_ke(:) + else + rot_frag(:, i) = rot_frag(:, i) + rot_L(:) + end if + L_frag_spin(:) = L_frag_spin(:) + m_frag(i) * rad_frag(i)**2 * Ip_frag(3, i) * rot_frag(:, i) + ke_frag_spin = ke_frag_spin + m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + end do + ke_frag_spin = 0.5_DP * ke_frag_spin + ! Convert a fraction of the pre-impact angular momentum into fragment spin angular momentum + L_frag_orb(:) = L_frag_tot(:) - L_frag_spin(:) + L_remainder(:) = L_frag_orb(:) + ! Next we will solve for the tangential component of the velocities that both conserves linear momentum and uses the remaining angular momentum not used in spin. + ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, + ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution + ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. + do i = 1, nfrag + v_t_initial(i) = norm2(L_remainder(:)) / ((nfrag - i + 1) * m_frag(i) * norm2(x_frag(:,i))) + Li(:) = m_frag(i) * x_frag(:,i) .cross. v_t_initial(i) * v_t_unit(:, i) + L_remainder(:) = L_remainder(:) - Li(:) + end do + + ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum + objective_function = lambda_obj(tangential_objective_function, lerr) + v_t_mag(7:nfrag) = util_minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), TOL, lerr) + ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis + v_t_initial(7:nfrag) = v_t_mag(7:nfrag) + v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lerr=lerr) + + ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) + vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) + ! Now do a kinetic energy budget check to make sure we are still within the budget. + ke_frag_orbit = 0.0_DP + do i = 1, nfrag + v_frag(:, i) = vb_frag(:, i) - vcom(:) + ke_frag_orbit = ke_frag_orbit + m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + end do + ke_frag_orbit = 0.5_DP * ke_frag_orbit + ke_radial = ke_frag_budget - ke_frag_orbit - ke_frag_spin + + ! If we are over the energy budget, flag this as a failure so we can try again + lerr = (ke_radial < 0.0_DP) + !write(*,*) 'Tangential' + !write(*,*) 'ke_frag_budget: ',ke_frag_budget + !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit + !write(*,*) 'ke_frag_spin : ',ke_frag_spin + !write(*,*) 'ke_radial : ',ke_radial + + return + + end subroutine set_fragment_tan_vel + + function tangential_objective_function(v_t_mag_input, lerr) result(fval) + !! Author: David A. Minton + !! + !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint + logical, intent(out) :: lerr !! Error flag + ! Result + real(DP) :: fval + ! Internals + integer(I4B) :: i + real(DP), dimension(:,:), allocatable :: v_shift + real(DP), dimension(:), allocatable :: v_t_new + real(DP) :: keo + + lerr = .false. + + allocate(v_shift(NDIM, nfrag)) + allocate(v_t_new(nfrag)) + + v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lerr=lerr) + v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_new, v_t_unit, m_frag, vcom) + + keo = 0.0_DP + do i = 1, nfrag + keo = keo + m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) + end do + keo = 0.5_DP * keo + fval = keo + lerr = .false. + return + end function tangential_objective_function + + function solve_fragment_tan_vel(lerr, v_t_mag_input) result(v_t_mag_output) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum + implicit none + ! Arguments + logical, intent(out) :: lerr !! Error flag + real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag + ! Internals + integer(I4B) :: i + ! Result + real(DP), dimension(:), allocatable :: v_t_mag_output + + real(DP), dimension(2 * NDIM, 2 * NDIM) :: A ! LHS of linear equation used to solve for momentum constraint in Gauss elimination code + real(DP), dimension(2 * NDIM) :: b ! RHS of linear equation used to solve for momentum constraint in Gauss elimination code + real(DP), dimension(NDIM) :: L_lin_others, L_orb_others, L, vtmp + + v_frag(:,:) = 0.0_DP + lerr = .false. + + ! We have 6 constraint equations (2 vector constraints in 3 dimensions each) + ! The first 3 are that the linear momentum of the fragments is zero with respect to the collisional barycenter + ! The second 3 are that the sum of the angular momentum of the fragments is conserved from the pre-impact state + L_lin_others(:) = 0.0_DP + L_orb_others(:) = 0.0_DP + do i = 1, nfrag + if (i <= 2 * NDIM) then ! The tangential velocities of the first set of bodies will be the unknowns we will solve for to satisfy the constraints + A(1:3, i) = m_frag(i) * v_t_unit(:, i) + L(:) = v_r_unit(:, i) .cross. v_t_unit(:, i) + A(4:6, i) = m_frag(i) * rmag(i) * L(:) + else if (present(v_t_mag_input)) then + vtmp(:) = v_t_mag_input(i - 6) * v_t_unit(:, i) + L_lin_others(:) = L_lin_others(:) + m_frag(i) * vtmp(:) + L(:) = x_frag(:, i) .cross. vtmp(:) + L_orb_others(:) = L_orb_others(:) + m_frag(i) * L(:) + end if + end do + b(1:3) = -L_lin_others(:) + b(4:6) = L_frag_orb(:) - L_orb_others(:) + allocate(v_t_mag_output(nfrag)) + v_t_mag_output(1:6) = util_solve_linear_system(A, b, 6, lerr) + if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) + + return + end function solve_fragment_tan_vel + + subroutine set_fragment_radial_velocities(lerr) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! + !! Adjust the fragment velocities to set the fragment orbital kinetic energy. This will minimize the difference between the fragment kinetic energy and the energy budget + implicit none + ! Arguments + logical, intent(out) :: lerr + ! Internals + real(DP), parameter :: TOL = 1e-10_DP + integer(I4B) :: i, j + real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma + real(DP), dimension(:,:), allocatable :: v_r + type(lambda_obj) :: objective_function + + ! Set the "target" ke_orbit_after (the value of the orbital kinetic energy that the fragments ought to have) + + allocate(v_r_initial, source=v_r_mag) + ! Initialize radial velocity magnitudes with a random value that is approximately 10% of that found by distributing the kinetic energy equally + allocate(v_r_sigma, source=v_r_mag) + call random_number(v_r_sigma(1:nfrag)) + v_r_sigma(1:nfrag) = sqrt(1.0_DP + 2 * (v_r_sigma(1:nfrag) - 0.5_DP) * 1e-4_DP) + v_r_initial(1:nfrag) = v_r_sigma(1:nfrag) * sqrt(abs(ke_radial) / (2 * m_frag(1:nfrag) * nfrag)) + + ! Initialize the lambda function using a structure constructor that calls the init method + ! Minimize the ke objective function using the BFGS optimizer + objective_function = lambda_obj(radial_objective_function) + v_r_mag = util_minimize_bfgs(objective_function, nfrag, v_r_initial, TOL, lerr) + ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) + vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) + do i = 1, nfrag + v_frag(:, i) = vb_frag(:, i) - vcom(:) + end do + ke_frag_orbit = 0.0_DP + do i = 1, nfrag + ke_frag_orbit = ke_frag_orbit + m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + end do + ke_frag_orbit = 0.5_DP * ke_frag_orbit + !write(*,*) 'Radial' + !write(*,*) 'Failure? ',lerr + !write(*,*) 'ke_frag_budget: ',ke_frag_budget + !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit + !write(*,*) 'ke_frag_spin : ',ke_frag_spin + !write(*,*) 'ke_remainder : ',ke_frag_budget - (ke_frag_orbit + ke_frag_spin) + lerr = .false. + + return + end subroutine set_fragment_radial_velocities + + function radial_objective_function(v_r_mag_input) result(fval) + !! Author: David A. Minton + !! + !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_r_mag_input !! Unknown radial component of fragment velocity vector + ! Result + real(DP) :: fval !! The objective function result, which is the square of the difference between the calculated fragment kinetic energy and our target + !! Minimizing this brings us closer to our objective + ! Internals + integer(I4B) :: i + real(DP), dimension(:,:), allocatable :: v_shift + + allocate(v_shift, mold=vb_frag) + v_shift(:,:) = vmag_to_vb(v_r_mag_input, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) + fval = 2 * ke_frag_budget + do i = 1, nfrag + fval = fval - m_frag(i) * (Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + dot_product(v_shift(:, i), v_shift(:, i))) + end do + ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for + fval = (fval / (2 * ke_radial))**2 + + return + + end function radial_objective_function + + function vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) + !! Author: David A. Minton + !! + !! Converts radial and tangential velocity magnitudes into barycentric velocity + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector + real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint + real(DP), dimension(:,:), intent(in) :: v_r_unit, v_t_unit !! Radial and tangential unit vectors for each fragment + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass + ! Result + real(DP), dimension(:,:), allocatable :: vb + ! Internals + integer(I4B) :: i + + allocate(vb, mold=v_r_unit) + ! Make sure the velocity magnitude stays positive + do i = 1, nfrag + vb(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) + end do + ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass + call shift_vector_to_origin(m_frag, vb) + + do i = 1, nfrag + vb(:, i) = vb(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) + end do + + end function vmag_to_vb + + subroutine restructure_failed_fragments() + !! Author: David A. Minton + !! + !! We failed to find a set of positions and velocities that satisfy all the constraints, and so we will alter the fragments and try again. + implicit none + integer(I4B) :: i + real(DP), dimension(:), allocatable :: m_frag_new, rad_frag_new + real(DP), dimension(:,:), allocatable :: xb_frag_new, vb_frag_new, Ip_frag_new, rot_frag_new + real(DP) :: delta_r, delta_r_max + real(DP), parameter :: ke_avg_deficit_target = 0.0_DP + + ! Introduce a bit of noise in the radius determination so we don't just flip flop between similar failed positions + call random_number(delta_r_max) + delta_r_max = sum(radius(:)) * (1.0_DP + 2e-1_DP * (delta_r_max - 0.5_DP)) + if (try > 2) then + ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try + delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) + if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) + else + delta_r = delta_r_max + end if + r_max_start_old = r_max_start + r_max_start = r_max_start + delta_r ! The larger lever arm can help if the problem is in the angular momentum step + if (f_spin > epsilon(1.0_DP)) then + f_spin = f_spin / 2 + else + f_spin = 0.0_DP + end if + end subroutine restructure_failed_fragments + + + end subroutine fragmentation_initialize + + + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 1aec7dc9d..d00d271fb 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -966,6 +966,17 @@ end subroutine util_fill_arr_logical end interface interface + module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) + use lambda_function + implicit none + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0 + real(DP), intent(in) :: eps + logical, intent(out) :: lerr + real(DP), dimension(:), allocatable :: x1 + end function util_minimize_bfgs + module subroutine util_peri_tp(self, system, param) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object @@ -974,6 +985,39 @@ module subroutine util_peri_tp(self, system, param) end subroutine util_peri_tp end interface + interface util_solve_linear_system + module function util_solve_linear_system_d(A,b,n,lerr) result(x) + implicit none + integer(I4B), intent(in) :: n + real(DP), dimension(:,:), intent(in) :: A + real(DP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + real(DP), dimension(n) :: x + end function util_solve_linear_system_d + + module function util_solve_linear_system_q(A,b,n,lerr) result(x) + implicit none + integer(I4B), intent(in) :: n + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + real(QP), dimension(n) :: x + end function util_solve_linear_system_q + end interface + + interface + module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) + use lambda_function + implicit none + class(lambda_obj), intent(inout) :: f !! lambda function object that has been initialized to be a function of derivatives. The object will return with components lastarg and lasteval set + real(DP), dimension(:), intent(in) :: y0in !! Initial value at t=0 + real(DP), intent(in) :: t1 !! Final time + real(DP), intent(in) :: dt0 !! Initial step size guess + real(DP), intent(in) :: tol !! Tolerance on solution + real(DP), dimension(:), allocatable :: y1 !! Final result + end function util_solve_rkf45 + end interface + interface util_resize module subroutine util_resize_arr_char_string(arr, nnew) implicit none diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 new file mode 100644 index 000000000..fbd48c8c2 --- /dev/null +++ b/src/util/util_minimize_bfgs.f90 @@ -0,0 +1,584 @@ +submodule (swiftest_classes) s_util_minimize_bfgs + use swiftest +contains + module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) + !! author: David A. Minton + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !! This function implements the Broyden-Fletcher-Goldfarb-Shanno method to determine the minimum of a function of N variables. + !! It recieves as input: + !! f%eval(x) : lambda function object containing the objective function as the eval metho + !! N : Number of variables of function f + !! x0 : Initial starting value of x + !! eps : Accuracy of 1 - dimensional minimization at each step + !! The outputs include + !! lerr : Returns .true. if it could not find the minimum + !! Returns + !! x1 : Final minimum (all 0 if none found) + !! 0 = No miniumum found + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + use, intrinsic :: ieee_exceptions + implicit none + ! Arguments + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0 + real(DP), intent(in) :: eps + logical, intent(out) :: lerr + ! Result + real(DP), dimension(:), allocatable :: x1 + ! Internals + integer(I4B) :: i, j, k, l, conv, num + integer(I4B), parameter :: MAXLOOP = 1000 !! Maximum number of loops before method is determined to have failed + real(DP), parameter :: graddelta = 1e-4_DP !! Delta x for gradient calculations + real(DP), dimension(N) :: S !! Direction vectors + real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix + real(DP), dimension(N) :: grad1 !! gradient of f + real(DP), dimension(N) :: grad0 !! old value of gradient + real(DP) :: astar !! 1D minimized value + real(DP), dimension(N) :: y, P + real(DP), dimension(N,N) :: PP, PyH, HyP + real(DP) :: yHy, Py + type(ieee_status_type) :: original_fpe_status + logical, dimension(:), allocatable :: fpe_flag + + call ieee_get_status(original_fpe_status) ! Save the original floating point exception status + call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet + allocate(fpe_flag(size(ieee_usual))) + + lerr = .false. + allocate(x1, source=x0) + ! Initialize approximate Hessian with the identity matrix (i.e. begin with method of steepest descent) + ! Get initial gradient and initialize arrays for updated values of gradient and x + H(:,:) = reshape([((0._DP, i=1, j-1), 1._DP, (0._DP, i=j+1, N), j=1, N)], [N,N]) + grad0 = gradf(f, N, x0(:), graddelta, lerr) + if (lerr) then + call ieee_set_status(original_fpe_status) + return + end if + grad1(:) = grad0(:) + do i = 1, MAXLOOP + !check for convergence + conv = count(abs(grad1(:)) > eps) + if (conv == 0) then + !write(*,*) "BFGS converged on gradient after ",i," iterations" + exit + end if + S(:) = -matmul(H(:,:), grad1(:)) + astar = minimize1D(f, x1, S, N, graddelta, lerr) + if (lerr) then + !write(*,*) "Exiting BFGS with error in minimize1D step" + exit + end if + ! Get new x values + P(:) = astar * S(:) + x1(:) = x1(:) + P(:) + ! Calculate new gradient + grad0(:) = grad1(:) + grad1 = gradf(f, N, x1, graddelta, lerr) + y(:) = grad1(:) - grad0(:) + Py = sum(P(:) * y(:)) + ! set up factors for H matrix update + yHy = 0._DP + do k = 1, N + do j = 1, N + yHy = yHy + y(j) * H(j,k) * y(k) + end do + end do + ! prevent divide by zero (convergence) + if (abs(Py) < tiny(Py)) then + !write(*,*) "BFGS Converged on tiny Py after ",i," iterations" + exit + end if + ! set up update + PyH(:,:) = 0._DP + HyP(:,:) = 0._DP + do k = 1, N + do j = 1, N + PP(j, k) = P(j) * P(k) + do l = 1, N + PyH(j, k) = PyH(j, k) + P(j) * y(l) * H(l,k) + HyP(j, k) = HyP(j, k) + P(k) * y(l) * H(j,l) + end do + end do + end do + ! update H matrix + H(:,:) = H(:,:) + ((1._DP - yHy / Py) * PP(:,:) - PyH(:,:) - HyP(:,:)) / Py + ! Normalize to prevent it from blowing up if it takes many iterations to find a solution + H(:,:) = H(:,:) / norm2(H(:,:)) + ! Stop everything if there are any exceptions to allow the routine to fail gracefully + call ieee_get_flag(ieee_usual, fpe_flag) + if (any(fpe_flag)) exit + if (i == MAXLOOP) then + lerr = .true. + !write(*,*) "BFGS ran out of loops!" + end if + end do + call ieee_get_flag(ieee_usual, fpe_flag) + lerr = lerr .or. any(fpe_flag) + !if (any(fpe_flag)) write(*,*) 'BFGS did not converge due to fpe' + !if (lerr) write(*,*) "BFGS did not converge!" + call ieee_set_status(original_fpe_status) + + return + + contains + + function gradf(f, N, x1, dx, lerr) result(grad) + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !! Purpose: Estimates the gradient of a function using a central difference + !! approximation + !! Inputs: + !! f%eval(x) : lambda function object containing the objective function as the eval metho + !! N : number of variables N + !! x1 : x value array + !! dx : step size to use when calculating derivatives + !! Outputs: + !! lerr : .true. if an error occurred. Otherwise returns .false. + !! Returns + !! grad : N sized array containing estimated gradient of f at x1 + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none + ! Arguments + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x1 + real(DP), intent(in) :: dx + logical, intent(out) :: lerr + ! Result + real(DP), dimension(N) :: grad + ! Internals + integer(I4B) :: i, j + real(DP), dimension(N) :: xp, xm + real(DP) :: fp, fm + logical :: lerrp, lerrm + + do i = 1, N + do j = 1, N + if (j == i) then + xp(j) = x1(j) + dx + xm(j) = x1(j) - dx + else + xp(j) = x1(j) + xm(j) = x1(j) + end if + end do + select type (f) + class is (lambda_obj_err) + fp = f%eval(xp) + lerrp = f%lerr + fm = f%eval(xm) + lerrm = f%lerr + lerr = lerrp .or. lerrm + class is (lambda_obj) + fp = f%eval(xp) + fm = f%eval(xm) + lerr = .false. + end select + grad(i) = (fp - fm) / (2 * dx) + if (lerr) return + end do + return + end function gradf + + function minimize1D(f, x0, S, N, eps, lerr) result(astar) + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !! This program find the minimum of a function of N variables in a single direction + !! S using in sequence: + !! 1. A Bracketing method + !! 2. The golden section method + !! 3. A quadratic polynomial fit + !! Inputs + !! f%eval(x) : lambda function object containing the objective function as the eval metho + !! x0 : Array of size N of initial x values + !! S : Array of size N that determines the direction of minimization + !! N : Number of variables of function f + !! eps : Accuracy of 1 - dimensional minimization at each step + !! Output + !! lerr : .true. if an error occurred. Otherwise returns .false. + !! Returns + !! astar : Final minimum along direction S + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none + ! Arguments + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0, S + real(DP), intent(in) :: eps + logical, intent(out) :: lerr + ! Result + real(DP) :: astar + ! Internals + integer(I4B) :: num = 0 + real(DP), parameter :: step = 0.7_DP !! Bracketing method step size + real(DP), parameter :: gam = 1.2_DP !! Bracketing method expansion parameter + real(DP), parameter :: greduce = 0.2_DP !! Golden section method reduction factor + real(DP), parameter :: greduce2 = 0.1_DP ! Secondary golden section method reduction factor + real(DP) :: alo, ahi !! High and low values for 1 - D minimization routines + real(DP), parameter :: a0 = epsilon(1.0_DP) !! Initial guess of alpha + + alo = a0 + call bracket(f, x0, S, N, gam, step, alo, ahi, lerr) + if (lerr) then + !write(*,*) "BFGS bracketing step failed!" + return + end if + if (abs(alo - ahi) < eps) then + astar = alo + lerr = .false. + return + end if + call golden(f, x0, S, N, greduce, alo, ahi, lerr) + if (lerr) then + !write(*,*) "BFGS golden section step failed!" + return + end if + if (abs(alo - ahi) < eps) then + astar = alo + lerr = .false. + return + end if + call quadfit(f, x0, S, N, eps, alo, ahi, lerr) + if (lerr) then + !write(*,*) "BFGS quadfit failed!" + return + end if + if (abs(alo - ahi) < eps) then + astar = alo + lerr = .false. + return + end if + ! Quadratic fit method won't converge, so finish off with another golden section + call golden(f, x0, S, N, greduce2, alo, ahi, lerr) + if (.not. lerr) astar = (alo + ahi) / 2.0_DP + return + end function minimize1D + + function n2one(f, x0, S, N, a, lerr) result(fnew) + implicit none + ! Arguments + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0, S + real(DP), intent(in) :: a + logical, intent(out) :: lerr + + ! Return + real(DP) :: fnew + ! Internals + real(DP), dimension(N) :: xnew + integer(I4B) :: i + + xnew(:) = x0(:) + a * S(:) + fnew = f%eval(xnew(:)) + select type(f) + class is (lambda_obj_err) + lerr = f%lerr + class is (lambda_obj) + lerr = .false. + end select + return + end function n2one + + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !! This subroutine brackets the minimum. It recieves as input: + !! f%eval(x) : lambda function object containing the objective function as the eval metho + !! x0 : Array of size N of initial x values + !! S : Array of size N that determines the direction of minimization + !! gam : expansion parameter + !! step : step size + !! lo : initial guess of lo bracket value + !! The outputs include + !! lo : lo bracket + !! hi : hi bracket + !! lerr : .true. if an error occurred. Otherwise returns .false. + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none + ! Arguments + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0, S + real(DP), intent(in) :: gam, step + real(DP), intent(inout) :: lo + real(DP), intent(out) :: hi + logical, intent(out) :: lerr + ! Internals + real(DP) :: a0, a1, a2, atmp, da + real(DP) :: f0, f1, f2 + integer(I4B) :: i, j + integer(I4B), parameter :: MAXLOOP = 100 ! maximum number of loops before method is determined to have failed + real(DP), parameter :: eps = epsilon(lo) ! small number precision to test floating point equality + + ! set up initial bracket points + a0 = lo + da = step + a1 = a0 + da + a2 = a0 + 2 * da + f0 = n2one(f, x0, S, N, a0, lerr) + if (lerr) return + f1 = n2one(f, x0, S, N, a1, lerr) + if (lerr) return + f2 = n2one(f, x0, S, N, a2, lerr) + if (lerr) return + ! loop over bracket method until either min is bracketed method fails + do i = 1, MAXLOOP + if ((f0 > f1) .and. (f1 < f2)) then ! Minimum was found + lo = a0 + hi = a2 + return + else if ((f0 >= f1) .and. (f1 > f2)) then ! Function appears to decrease + da = da * gam + atmp = a2 + da + a0 = a1 + a1 = a2 + a2 = atmp + f0 = f1 + f1 = f2 + f2 = n2one(f, x0, S, N, a2, lerr) + else if ((f0 < f1) .and. (f1 <= f2)) then ! Function appears to increase + da = da * gam + atmp = a0 - da + a2 = a1 + a1 = a0 + a0 = atmp + f2 = f1 + f0 = n2one(f, x0, S, N, a0, lerr) + else if ((f0 < f1) .and. (f1 > f2)) then ! We are at a peak. Pick the direction that descends the fastest + da = da * gam + if (f2 > f0) then ! LHS is lower than RHS + atmp = a2 + da + a0 = a1 + a1 = a2 + a2 = atmp + f0 = f1 + f1 = f2 + f2 = n2one(f, x0, S, N, a2, lerr) + else ! RHS is lower than LHS + atmp = a0 - da + a2 = a1 + a1 = a0 + a0 = atmp + f2 = f1 + f1 = f2 + f0 = n2one(f, x0, S, N, a0, lerr) + end if + else if ((f0 > f1) .and. (abs(f2 - f1) <= eps)) then ! Decrasging but RHS equal + da = da * gam + atmp = a2 + da + a2 = atmp + f2 = n2one(f, x0, S, N, a2, lerr) + else if ((abs(f0 - f1) < eps) .and. (f1 < f2)) then ! Increasing but LHS equal + da = da * gam + atmp = a0 - da + a0 = atmp + f0 = n2one(f, x0, S, N, a0, lerr) + else ! all values equal. Expand in either direction and try again + a0 = a0 - da + a2 = a2 + da + f0 = n2one(f, x0, S, N, a0, lerr) + if (lerr) exit ! An error occurred while evaluating the function + f2 = n2one(f, x0, S, N, a2, lerr) + end if + if (lerr) exit ! An error occurred while evaluating the function + end do + lerr = .true. + return ! no minimum found + end subroutine bracket + + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + subroutine golden(f, x0, S, N, eps, lo, hi, lerr) + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !! This function uses the golden section method to reduce the starting interval lo, hi by some amount sigma. + !! It recieves as input: + !! f%eval(x) : lambda function object containing the objective function as the eval metho + !! x0 : Array of size N of initial x values + !! S : Array of size N that determines the direction of minimization + !! gam : expansion parameter + !! eps : reduction interval in range (0 < sigma < 1) such that: + !! hi(new) - lo(new) = eps * (hi(old) - lo(old)) + !! lo : initial guess of lo bracket value + !! The outputs include + !! lo : lo bracket + !! hi : hi bracket + !! lerr : .true. if an error occurred. Otherwise returns .false. + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none + ! Arguments + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0, S + real(DP), intent(in) :: eps + real(DP), intent(inout) :: lo + real(DP), intent(out) :: hi + logical, intent(out) :: lerr + ! Internals + real(DP), parameter :: tau = 0.5_DP * (sqrt(5.0_DP) - 1.0_DP) ! Golden section constant + integer(I4B), parameter :: MAXLOOP = 40 ! maximum number of loops before method is determined to have failed (unlikely, but could occur if no minimum exists between lo and hi) + real(DP) :: i0 ! Initial interval value + real(DP) :: a1, a2 + real(DP) :: f1, f2 + integer(I4B) :: i, j + + i0 = hi - lo + a1 = hi - tau * i0 + a2 = lo + tau * i0 + f1 = n2one(f, x0, S, N, a1, lerr) + if (lerr) return + f2 = n2one(f, x0, S, N, a2, lerr) + if (lerr) return + do i = 1, MAXLOOP + if (abs((hi - lo) / i0) <= eps) return ! interval reduced to input amount + if (f2 > f1) then + hi = a2 + a2 = a1 + f2 = f1 + a1 = hi - tau * (hi - lo) + f1 = n2one(f, x0, S, N, a1, lerr) + else + lo = a1 + a1 = a2 + f2 = f1 + a2 = hi - (1.0_DP - tau) * (hi - lo) + f2 = n2one(f, x0, S, N, a2, lerr) + end if + if (lerr) exit + end do + lerr = .true. + return ! search took too many iterations - no minimum found + end subroutine golden + + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !! This function uses a quadratic polynomial fit to locate the minimum of a function + !! to some accuracy eps. It recieves as input: + !! f%eval(x) : lambda function object containing the objective function as the eval metho + !! lo : low bracket value + !! hi : high bracket value + !! eps : desired accuracy of final minimum location + !! The outputs include + !! lo : final minimum location + !! hi : final minimum location + !! Notes: Uses the ieee_exceptions intrinsic module to allow for graceful failure due to floating point exceptions, which won't terminate the run. + !! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none + ! Arguments + integer(I4B), intent(in) :: N + class(lambda_obj), intent(inout) :: f + real(DP), dimension(:), intent(in) :: x0, S + real(DP), intent(in) :: eps + real(DP), intent(inout) :: lo + real(DP), intent(out) :: hi + logical, intent(out) :: lerr + ! Internals + integer(I4B), parameter :: MAXLOOP = 20 ! maximum number of loops before method is determined to have failed. + real(DP) :: a1, a2, a3, astar ! three points for the polynomial fit and polynomial minimum + real(DP) :: f1, f2, f3, fstar ! three function values for the polynomial and polynomial minimum + real(DP), dimension(3) :: row_1, row_2, row_3, rhs, soln ! matrix for 3 equation solver (gaussian elimination) + real(DP), dimension(3,3) :: lhs + real(DP) :: d1, d2, d3, aold, denom, errval + integer(I4B) :: i + + lerr = .false. + ! Get initial a1, a2, a3 values + a1 = lo + a2 = lo + 0.5_DP * (hi - lo) + a3 = hi + aold = a1 + astar = a2 + f1 = n2one(f, x0, S, N, a1, lerr) + if (lerr) return + f2 = n2one(f, x0, S, N, a2, lerr) + if (lerr) return + f3 = n2one(f, x0, S, N, a3, lerr) + if (lerr) return + do i = 1, MAXLOOP + ! check to see if convergence is reached and exit + errval = abs((astar - aold) / astar) + call ieee_get_flag(ieee_usual, fpe_flag) + if (any(fpe_flag)) then + !write(*,*) 'quadfit fpe' + !write(*,*) 'aold : ',aold + !write(*,*) 'astar: ',astar + lerr = .true. + exit + end if + if (errval < eps) then + lo = astar + hi = astar + exit + end if + ! Set up system for gaussian elimination equation solver + row_1 = [1.0_DP, a1, a1**2] + row_2 = [1.0_DP, a2, a2**2] + row_3 = [1.0_DP, a3, a3**2] + rhs = [f1, f2, f3] + lhs(1, :) = row_1 + lhs(2, :) = row_2 + lhs(3, :) = row_3 + ! Solve system of equations + soln(:) = util_solve_linear_system(lhs, rhs, 3, lerr) + call ieee_set_flag(ieee_all, .false.) ! Set all flags back to quiet + call ieee_set_halting_mode(ieee_divide_by_zero, .false.) + if (lerr) then + !write(*,*) 'quadfit fpe:' + !write(*,*) 'util_solve_linear_system failed' + exit + end if + aold = astar + if (soln(2) == soln(3)) then ! Handles the case where they are both 0. 0/0 is an unhandled exception + astar = -0.5_DP + else + astar = -soln(2) / (2 * soln(3)) + end if + call ieee_get_flag(ieee_usual, fpe_flag) + if (any(fpe_flag)) then + !write(*,*) 'quadfit fpe' + !write(*,*) 'soln(2:3): ',soln(2:3) + !write(*,*) 'a1, a2, a3' + !write(*,*) a1, a2, a3 + !write(*,*) 'f1, f2, f3' + !write(*,*) f1, f2, f3 + lerr = .true. + exit + end if + fstar = n2one(f, x0, S, N, astar, lerr) + if (lerr) exit + ! keep the three closest a values to astar and discard the fourth + d1 = abs(a1 - astar) + d2 = abs(a2 - astar) + d3 = abs(a3 - astar) + + if (d1 > d2) then + if (d1 > d3) then + f1 = fstar + a1 = astar + else if (d3 > d2) then + f3 = fstar + a3 = astar + end if + else + if (d2 > d3) then + f2 = fstar + a2 = astar + else if (d3 > d1) then + f3 = fstar + a3 = astar + end if + end if + end do + if (lerr) return + lo = a1 + hi = a3 + return + end subroutine quadfit + + end function util_minimize_bfgs +end submodule s_util_minimize_bfgs \ No newline at end of file diff --git a/src/util/util_solve.f90 b/src/util/util_solve.f90 index 0c3161ae2..ed9ff40ea 100644 --- a/src/util/util_solve.f90 +++ b/src/util/util_solve.f90 @@ -2,7 +2,143 @@ use swiftest contains - function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) + module function util_solve_linear_system_d(A,b,n,lerr) result(x) + !! Author: David A. Minton + !! + !! Solves the linear equation of the form A*x = b for x. + !! A is an (n,n) arrays + !! x and b are (n) arrays + !! Uses Gaussian elimination, so will have issues if system is ill-conditioned. + !! Uses quad precision intermidiate values, so works best on small arrays. + use, intrinsic :: ieee_exceptions + implicit none + ! Arguments + integer(I4B), intent(in) :: n + real(DP), dimension(:,:), intent(in) :: A + real(DP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + ! Result + real(DP), dimension(n) :: x + ! Internals + real(QP), dimension(:), allocatable :: qx + type(ieee_status_type) :: original_fpe_status + logical, dimension(:), allocatable :: fpe_flag + + call ieee_get_status(original_fpe_status) ! Save the original floating point exception status + call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet + allocate(fpe_flag(size(ieee_usual))) + + qx = solve_wbs(ge_wpp(real(A, kind=QP), real(b, kind=QP))) + + call ieee_get_flag(ieee_usual, fpe_flag) + lerr = any(fpe_flag) + if (lerr .or. (any(abs(qx) > huge(x))) .or. (any(abs(qx) < tiny(x)))) then + x = 0.0_DP + else + x = real(qx, kind=DP) + end if + call ieee_set_status(original_fpe_status) + + return + end function util_solve_linear_system_d + + module function util_solve_linear_system_q(A,b,n,lerr) result(x) + !! Author: David A. Minton + !! + !! Solves the linear equation of the form A*x = b for x. + !! A is an (n,n) arrays + !! x and b are (n) arrays + !! Uses Gaussian elimination, so will have issues if system is ill-conditioned. + !! Uses quad precision intermidiate values, so works best on small arrays. + use, intrinsic :: ieee_exceptions + use swiftest + implicit none + ! Arguments + integer(I4B), intent(in) :: n + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + ! Result + real(QP), dimension(n) :: x + ! Internals + type(ieee_status_type) :: original_fpe_status + logical, dimension(:), allocatable :: fpe_flag + + call ieee_get_status(original_fpe_status) ! Save the original floating point exception status + call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet + allocate(fpe_flag(size(ieee_usual))) + + x = solve_wbs(ge_wpp(A, b)) + + call ieee_get_flag(ieee_usual, fpe_flag) + lerr = any(fpe_flag) + if (lerr) x = 0.0_DP + call ieee_set_status(original_fpe_status) + + return + end function util_solve_linear_system_q + + function solve_wbs(u) result(x) ! solve with backward substitution + !! Based on code available on Rosetta Code: https://rosettacode.org/wiki/Gaussian_elimination#Fortran + use, intrinsic :: ieee_exceptions + use swiftest + implicit none + ! Arguments + real(QP), intent(in), dimension(:,:), allocatable :: u + ! Result + real(QP), dimension(:), allocatable :: x + ! Internals + integer(I4B) :: i,n + + n = size(u, 1) + if (allocated(x)) deallocate(x) + if (.not.allocated(x)) allocate(x(n)) + if (any(abs(u) < tiny(1._DP)) .or. any(abs(u) > huge(1._DP))) then + x(:) = 0._DP + return + end if + call ieee_set_halting_mode(ieee_divide_by_zero, .false.) + do i = n, 1, -1 + x(i) = (u(i, n + 1) - sum(u(i, i + 1:n) * x(i + 1:n))) / u(i, i) + end do + return + end function solve_wbs + + function ge_wpp(A, b) result(u) ! gaussian eliminate with partial pivoting + !! Solve Ax=b using Gaussian elimination then backwards substitution. + !! A being an n by n matrix. + !! x and b are n by 1 vectors. + !! Based on code available on Rosetta Code: https://rosettacode.org/wiki/Gaussian_elimination#Fortran + use, intrinsic :: ieee_exceptions + use swiftest + implicit none + ! Arguments + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + ! Result + real(QP), dimension(:,:), allocatable :: u + ! Internals + integer(I4B) :: i,j,n,p + real(QP) :: upi + + n = size(a, 1) + allocate(u(n, (n + 1))) + u = reshape([A, b], [n, n + 1]) + call ieee_set_halting_mode(ieee_divide_by_zero, .false.) + do j = 1, n + p = maxloc(abs(u(j:n, j)), 1) + j - 1 ! maxloc returns indices between (1, n - j + 1) + if (p /= j) u([p, j], j) = u([j, p], j) + u(j + 1:, j) = u(j + 1:, j) / u(j, j) + do i = j + 1, n + 1 + upi = u(p, i) + if (p /= j) u([p, j], i) = u([j, p], i) + u(j + 1:n, i) = u(j + 1:n, i) - upi * u(j + 1:n, j) + end do + end do + return + end function ge_wpp + + module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) !! author: David A. Minton !! !! Implements the 4th order Runge-Kutta-Fehlberg ODE solver for initial value problems of the form f=dy/dt, y0 = y(t=0), solving for y1 = y(t=t1). Uses a 5th order adaptive step size control. From 8b4d5237f4738e8dea8154ae489762ff137f7724 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 5 Aug 2021 07:50:26 -0400 Subject: [PATCH 10/12] Fixed duplicate use statement that was giving a problem in util_solve with ifort --- Makefile.Defines | 8 ++++---- src/util/util_solve.f90 | 1 - 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index 70069bb71..291f2c604 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -65,13 +65,13 @@ GPAR = -fopenmp -ftree-parallelize-loops=4 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries -#FFLAGS = $(IDEBUG) $(HEAPARR) +FFLAGS = $(IDEBUG) $(HEAPARR) #FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) -#FORTRAN = ifort +FORTRAN = ifort #AR = xiar -FORTRAN = gfortran -FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) +#FORTRAN = gfortran +#FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/util/util_solve.f90 b/src/util/util_solve.f90 index ed9ff40ea..057ed1182 100644 --- a/src/util/util_solve.f90 +++ b/src/util/util_solve.f90 @@ -51,7 +51,6 @@ module function util_solve_linear_system_q(A,b,n,lerr) result(x) !! Uses Gaussian elimination, so will have issues if system is ill-conditioned. !! Uses quad precision intermidiate values, so works best on small arrays. use, intrinsic :: ieee_exceptions - use swiftest implicit none ! Arguments integer(I4B), intent(in) :: n From d0e7fb4c150ee2ef8132c989cfe5ee6380ba4dc3 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 5 Aug 2021 07:59:37 -0400 Subject: [PATCH 11/12] Enabled the disruption and supercat fragmentation and got all of the interfaces straightened out --- src/fragmentation/fragmentation.f90 | 2 +- src/modules/swiftest_classes.f90 | 16 ++++ src/symba/symba_fragmentation.f90 | 111 +++++++++++++++++++++++++++- 3 files changed, 125 insertions(+), 4 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index d1965914f..460060183 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -2,7 +2,7 @@ use swiftest contains - subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index d00d271fb..2455e77f2 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -453,6 +453,22 @@ module subroutine eucl_dist_index_plpl(self) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object end subroutine + module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + implicit none + class(swiftest_nbody_system), intent(inout) :: system + class(swiftest_parameters), intent(in) :: param + integer(I4B), dimension(:), intent(in) :: family + real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip + real(DP), dimension(:), intent(inout) :: mass, radius + integer(I4B), intent(inout) :: nfrag + real(DP), dimension(:), allocatable, intent(inout) :: m_frag, rad_frag + real(DP), dimension(:,:), allocatable, intent(inout) :: Ip_frag + real(DP), dimension(:,:), allocatable, intent(inout) :: xb_frag, vb_frag, rot_frag + logical, intent(out) :: lfailure ! Answers the question: Should this have been a merger instead? + real(DP), intent(inout) :: Qloss + end subroutine fragmentation_initialize + module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss) implicit none integer(I4B), intent(out) :: regime diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index fea7dea91..138ab9179 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -73,8 +73,8 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, Ip_frag(:, i) = Ip_new(:) end do - !call fragmentation_initialize(pl, param, family, x, v, L_spin, Ip, mass, radius, & - ! nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) if (lfailure) then write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' @@ -83,6 +83,7 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, else ! Populate the list of new bodies write(*,'("Generating ",I2.0," fragments")') nfrag + status = DISRUPTION allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -309,8 +310,112 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals + integer(I4B) :: i, nfrag, ibiggest, nstart, nend + real(DP) :: mtot, avg_dens, min_frag_mass + real(DP), dimension(NDIM) :: xcom, vcom + real(DP), dimension(2) :: vol + real(DP), dimension(NDIM) :: Ip_new + real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag + real(DP), dimension(:), allocatable :: m_frag, rad_frag + logical :: lfailure + class(symba_pl), allocatable :: plnew + + select type(pl => system%pl) + class is (symba_pl) + associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + nfrag = NFRAG_SUPERCAT + allocate(m_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Get mass weighted mean of Ip and average density + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot + vol(:) = 4._DP / 3._DP * pi * radius(:)**3 + avg_dens = mtot / sum(vol(:)) + + ! If we are adding the first and largest fragment (lr), check to see if its mass is SMALLER than an equal distribution of + ! mass between all fragments. If so, we will just distribute the mass equally between the fragments + min_frag_mass = mtot / nfrag + if (mass_res(1) < min_frag_mass) then + m_frag(:) = min_frag_mass + else + m_frag(1) = mass_res(1) + m_frag(2:nfrag) = (mtot - mass_res(1)) / (nfrag - 1) + end if + ! Distribute any residual mass if there is any and set the radius + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + + do i = 1, nfrag + Ip_frag(:, i) = Ip_new(:) + end do - status = SUPERCATASTROPHIC + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 + else + ! Populate the list of new bodies + write(*,'("Generating ",I2.0," fragments")') nfrag + status = SUPERCATASTROPHIC + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + + plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) + do i = 1, nfrag + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%density(:) = avg_dens + plnew%radius(:) = rad_frag(:) + plnew%info(:)%origin_type = "Disruption" + plnew%info(:)%origin_time = param%t + do i = 1, nfrag + plnew%info(i)%origin_xh(:) = plnew%xh(:,i) + plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + end do + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + + ! Append the new merged body to the list and record how many we made + nstart = mergeadd_list%nbody + 1 + nend = mergeadd_list%nbody + plnew%nbody + call mergeadd_list%append(plnew) + mergeadd_list%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) + end if + + end associate + end select return end function symba_fragmentation_casesupercatastrophic From ee7bdec0544fff1e9c4af8cc1fdd22e0213cf765 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 5 Aug 2021 08:12:25 -0400 Subject: [PATCH 12/12] Added in hit and run case --- src/symba/symba_fragmentation.f90 | 166 +++++++++++++++++++++++++++++- 1 file changed, 162 insertions(+), 4 deletions(-) diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index 138ab9179..efdd8c0d7 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -22,13 +22,14 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, istart, nfrag, ibiggest, nstart, nend + integer(I4B) :: i, istart, nfrag, ibiggest, nfamily, nstart, nend real(DP) :: mtot, avg_dens real(DP), dimension(NDIM) :: xcom, vcom, Ip_new real(DP), dimension(2) :: vol real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag real(DP), dimension(:), allocatable :: m_frag, rad_frag logical :: lfailure + logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew select type(pl => system%pl) @@ -84,6 +85,18 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, ! Populate the list of new bodies write(*,'("Generating ",I2.0," fragments")') nfrag status = DISRUPTION + + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = mergesub_list%nbody + 1 + nend = mergesub_list%nbody + nfamily + call mergesub_list%append(pl, lmask) + ! Record how many bodies were subtracted in this event + mergesub_list%ncomp(nstart:nend) = nfamily + allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -153,8 +166,140 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals + integer(I4B) :: i, nfrag, jproj, jtarg, idstart, ibiggest, nfamily, nstart, nend + real(DP) :: mtot, avg_dens + real(DP), dimension(NDIM) :: xcom, vcom + real(DP), dimension(2) :: vol + real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag + real(DP), dimension(:), allocatable :: m_frag, rad_frag + integer(I4B), dimension(:), allocatable :: id_frag + logical :: lpure + logical, dimension(system%pl%nbody) :: lmask + class(symba_pl), allocatable :: plnew + + select type(pl => system%pl) + class is (symba_pl) + associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + lpure = .false. + + ! The largest body will stay untouched + if (mass(1) > mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + + if (mass_res(2) > 0.9_DP * mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + write(*,*) '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 + nfrag = NFRAG_DISRUPT - 1 + lpure = .false. + allocate(m_frag(nfrag)) + allocate(id_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + m_frag(1) = mass(jtarg) + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + id_frag(1) = pl%id(ibiggest) + rad_frag(1) = radius(jtarg) + xb_frag(:, 1) = x(:, jtarg) + vb_frag(:, 1) = v(:, jtarg) + Ip_frag(:,1) = Ip(:, jtarg) + + ! Get mass weighted mean of Ip and average density + vol(:) = 4._DP / 3._DP * pi * radius(:)**3 + avg_dens = mass(jproj) / vol(jproj) + m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) + rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + + do i = 1, nfrag + Ip_frag(:, i) = Ip(:, jproj) + end do + + ! Put the fragments on the circle surrounding the center of mass of the system + call fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & + nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lpure) + if (lpure) then + write(*,*) 'Should have been a pure hit and run instead' + nfrag = 0 + else + write(*,'("Generating ",I2.0," fragments")') nfrag + end if + end if + if (lpure) then + status = ACTIVE + else + status = HIT_AND_RUN + + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = mergesub_list%nbody + 1 + nend = mergesub_list%nbody + nfamily + call mergesub_list%append(pl, lmask) + ! Record how many bodies were subtracted in this event + mergesub_list%ncomp(nstart:nend) = nfamily + + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + + plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) + do i = 1, nfrag + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%density(:) = avg_dens + plnew%radius(:) = rad_frag(:) + plnew%info(:)%origin_type = "Hit and run fragment" + plnew%info(:)%origin_time = param%t + do i = 1, nfrag + plnew%info(i)%origin_xh(:) = plnew%xh(:,i) + plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + end do + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + + ! Append the new merged body to the list and record how many we made + nstart = mergeadd_list%nbody + 1 + nend = mergeadd_list%nbody + plnew%nbody + call mergeadd_list%append(plnew) + mergeadd_list%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) - status = HIT_AND_RUN + end if + end associate + end select return end function symba_fragmentation_casehitandrun @@ -310,7 +455,7 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, nfrag, ibiggest, nstart, nend + integer(I4B) :: i, nfrag, ibiggest, nfamily, nstart, nend real(DP) :: mtot, avg_dens, min_frag_mass real(DP), dimension(NDIM) :: xcom, vcom real(DP), dimension(2) :: vol @@ -318,6 +463,7 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag real(DP), dimension(:), allocatable :: m_frag, rad_frag logical :: lfailure + logical, dimension(system%pl%nbody) :: lmask class(symba_pl), allocatable :: plnew select type(pl => system%pl) @@ -369,6 +515,18 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, ! Populate the list of new bodies write(*,'("Generating ",I2.0," fragments")') nfrag status = SUPERCATASTROPHIC + + ! Add the family bodies to the subtraction list + nfamily = size(family(:)) + lmask(:) = .false. + lmask(family(:)) = .true. + pl%status(family(:)) = MERGED + nstart = mergesub_list%nbody + 1 + nend = mergesub_list%nbody + nfamily + call mergesub_list%append(pl, lmask) + ! Record how many bodies were subtracted in this event + mergesub_list%ncomp(nstart:nend) = nfamily + allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -387,7 +545,7 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, plnew%Gmass(:) = param%GU * m_frag(:) plnew%density(:) = avg_dens plnew%radius(:) = rad_frag(:) - plnew%info(:)%origin_type = "Disruption" + plnew%info(:)%origin_type = "Supercatastrophic" plnew%info(:)%origin_time = param%t do i = 1, nfrag plnew%info(i)%origin_xh(:) = plnew%xh(:,i)