diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index e8d2b9fce..baff5c002 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -206,7 +206,7 @@ def data_stream(self, frame=0): # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(fragmentation=True, encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(collision_model="merge", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=0) print("Generating animation") diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 9666ca298..530d8b034 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -224,16 +224,18 @@ def __init__(self,read_param: bool = False, read_old_output: bool = False, simdi general_relativity : bool, default True Include the post-Newtonian correction in acceleration calculations. Parameter input file equivalent: `GR` - fragmentation : bool, default True - If set to True, this turns on the Fraggle fragment generation code and `rotation` must also be True. + collision_model: {"MERGE","BOUNCE","SIMPLE","FRAGGLE"}, default "MERGE" + This is used to set the collision/fragmentation model. [TODO: DESCRIBE THESE] This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. - Parameter input file equivalent: `FRAGMENTATION` + Parameter input file equivalent: `COLLISION_MODEL` minimum_fragment_gmass : float, optional - If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated. + If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated if a + fragmentation model is enabled. Ignored otherwise. *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass Parameter input file equivalent: None minimum_fragment_mass : float, optional - If fragmentation is turned on, this sets the mimimum mass of a collisional fragment that can be generated. + If fragmentation is turned on, this sets the mimimum mass of a collisional fragment that can be generated. if a + fragmentation model is enabled. Ignored otherwise *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass Parameter input file equivalent: `MIN_GMFRAG` rotation : bool, default False @@ -777,7 +779,7 @@ def set_parameter(self, verbose: bool = True, **kwargs): "mtiny": None, "close_encounter_check": True, "general_relativity": True, - "fragmentation": False, + "collision_model": "MERGE", "minimum_fragment_mass": None, "minimum_fragment_gmass": 0.0, "rotation": False, @@ -1013,7 +1015,7 @@ def get_integrator(self,arg_list: str | List[str] | None = None, verbose: bool | def set_feature(self, close_encounter_check: bool | None = None, general_relativity: bool | None = None, - fragmentation: bool | None = None, + collision_model: Literal["MERGE","BOUNCE","SIMPLE","FRAGGLE"] | None = None, minimum_fragment_gmass: float | None = None, minimum_fragment_mass: float | None = None, rotation: bool | None = None, @@ -1046,15 +1048,20 @@ def set_feature(self, *WARNING*: Enabling this feature could lead to very large files. general_relativity : bool, optional Include the post-Newtonian correction in acceleration calculations. - fragmentation : bool, optional - If set to True, this turns on the Fraggle fragment generation code and `rotation` must also be True. + collision_model: {"MERGE","BOUNCE","SIMPLE","FRAGGLE"}, default "MERGE" + This is used to set the collision/fragmentation model. [TODO: DESCRIBE THESE] This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. + Parameter input file equivalent: `COLLISION_MODEL` minimum_fragment_gmass : float, optional - If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated. + If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated if a + fragmentation model is enabled. Ignored otherwise. *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass + Parameter input file equivalent: None minimum_fragment_mass : float, optional - If fragmentation is turned on, this sets the mimimum mass of a collisional fragment that can be generated. + If fragmentation is turned on, this sets the mimimum mass of a collisional fragment that can be generated. if a + fragmentation model is enabled. Ignored otherwise *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass + Parameter input file equivalent: `MIN_GMFRAG` rotation : bool, optional If set to True, this turns on rotation tracking and radius, rotation vector, and moments of inertia values must be included in the initial conditions. @@ -1122,13 +1129,16 @@ def set_feature(self, self.param["GR"] = general_relativity update_list.append("general_relativity") - if fragmentation is not None: + fragmentation_models = ["FRAGGLE", "SIMPLE"] + if collision_model is not None: + collision_model = collision_model.upper() + fragmentation = collision_model in fragmentation_models if self.codename != "Swiftest" and self.integrator != "symba" and fragmentation: warnings.warn("Fragmentation is only available on Swiftest SyMBA.",stacklevel=2) - self.param['FRAGMENTATION'] = False + self.param['COLLISION_MODEL'] = "MERGE" else: - self.param['FRAGMENTATION'] = fragmentation - update_list.append("fragmentation") + self.param['COLLISION_MODEL'] = collision_model + update_list.append("collision_model") if fragmentation: if "MIN_GMFRAG" not in self.param and minimum_fragment_mass is None and minimum_fragment_gmass is None: warnings.warn("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass",stacklevel=2) @@ -1151,7 +1161,7 @@ def set_feature(self, self.param['ROTATION'] = rotation update_list.append("rotation") - if self.param['FRAGMENTATION'] and not self.param['ROTATION']: + if self.param['COLLISION_MODEL'] == "FRAGGLE" and not self.param['ROTATION']: self.param['ROTATION'] = True update_list.append("rotation") @@ -1230,7 +1240,7 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N ---------- arg_list: str | List[str], optional A single string or list of strings containing the names of the features to extract. Default is all of: - ["close_encounter_check", "general_relativity", "fragmentation", "rotation", "compute_conservation_values"] + ["close_encounter_check", "general_relativity", "collision_model", "rotation", "compute_conservation_values"] verbose: bool, optional If passed, it will override the Simulation object's verbose flag **kwargs @@ -1245,7 +1255,7 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N """ valid_var = {"close_encounter_check": "CHK_CLOSE", - "fragmentation": "FRAGMENTATION", + "collision_model": "COLLISION_MODEL", "encounter_save": "ENCOUNTER_SAVE", "minimum_fragment_gmass": "MIN_GMFRAG", "rotation": "ROTATION", diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 27d4a7bab..8cb9f6228 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -40,10 +40,11 @@ SET(FAST_MATH_FILES ${SRC}/helio/helio_module.f90 ${SRC}/symba/symba_module.f90 ${SRC}/collision/collision_check.f90 - ${SRC}/collision/collision_regime.f90 - ${SRC}/collision/collision_setup.f90 ${SRC}/collision/collision_io.f90 + ${SRC}/collision/collision_model.f90 + ${SRC}/collision/collision_regime.f90 ${SRC}/collision/collision_resolve.f90 + ${SRC}/collision/collision_setup.f90 ${SRC}/collision/collision_util.f90 ${SRC}/encounter/encounter_check.f90 ${SRC}/encounter/encounter_io.f90 diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index b85547cfd..bda2c6406 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -34,8 +34,8 @@ module base character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles character(STRMAX) :: in_netcdf = NC_INFILE !! Name of system input file for NetCDF input - character(STRMAX) :: in_type = "ASCII" !! Data representation type of input data files - character(STRMAX) :: in_form = "XV" !! Format of input data files ("EL" or "XV") + character(STRMAX) :: in_type = "NETCDF_DOUBLE" !! Data representation type of input data files + character(STRMAX) :: in_form = "XV" !! Format of input data files ("EL" or ["XV"]) integer(I4B) :: istep_out = -1 !! Number of time steps between saved outputs character(STRMAX) :: outfile = BIN_OUTFILE !! Name of output binary file character(STRMAX) :: out_type = "NETCDF_DOUBLE" !! Binary format of output file @@ -46,7 +46,7 @@ module base real(DP) :: rmax = -1.0_DP !! Maximum heliocentric radius for test particle real(DP) :: rmaxu = -1.0_DP !! Maximum unbound heliocentric radius for test particle real(DP) :: qmin = -1.0_DP !! Minimum pericenter distance for test particle - character(STRMAX) :: qmin_coord = 'HELIO' !! Coordinate frame to use for qmin + character(STRMAX) :: qmin_coord = "HELIO" !! Coordinate frame to use for qmin (["HELIO"] or "BARY") real(DP) :: qmin_alo = -1.0_DP !! Minimum semimajor axis for qmin real(DP) :: qmin_ahi = -1.0_DP !! Maximum semimajor axis for qmin real(QP) :: MU2KG = -1.0_QP !! Converts mass units to grams @@ -58,7 +58,7 @@ module base real(DP) :: min_GMfrag = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event integer(I4B), dimension(:), allocatable :: seed !! Random seeds for fragmentation modeling logical :: lmtiny_pl = .false. !! Include semi-interacting massive bodies - logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. + character(STRMAX) :: collision_model = "MERGE" !! The Coll character(STRMAX) :: encounter_save = "NONE" !! Indicate if and how encounter data should be saved logical :: lenc_save_trajectory = .false. !! Indicates that when encounters are saved, the full trajectory through recursion steps are saved logical :: lenc_save_closest = .false. !! Indicates that when encounters are saved, the closest approach distance between pairs of bodies is saved diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 new file mode 100644 index 000000000..578f3ad58 --- /dev/null +++ b/src/collision/collision_generate.f90 @@ -0,0 +1,39 @@ + +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule(collision) s_collision_model + use swiftest +contains + + module subroutine collision_generate_merge_system(self, system, param, t) + implicit none + class(collision_merge), intent(inout) :: self !! Fraggle fragment system object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + end subroutine collision_generate_merge_system + + module subroutine collision_generate_bounce_system(self, system, param, t) + implicit none + class(collision_bounce), intent(inout) :: self !! Fraggle fragment system object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + end subroutine collision_generate_bounce_system + + module subroutine collision_generate_simple_system(self, system, param, t) + implicit none + class(collision_simple), intent(inout) :: self !! Fraggle fragment system object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + end subroutine collision_generate_simple_system + +end submodule s_collision_model \ No newline at end of file diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index d04a07207..408fdfd7b 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -107,7 +107,7 @@ module collision end type collision_fragments - type :: collision_system + type, abstract :: collision_system !! This class defines a collisional system that stores impactors and fragments. This is written so that various collision models (i.e. Fraggle) could potentially be used !! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system @@ -124,7 +124,7 @@ module collision real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: Etot !! Before/after total system energy contains - procedure :: generate_fragments => abstract_generate_fragments !! Generates a system of fragments + procedure(abstract_generate_fragments), deferred :: generate_fragments !! Generates a system of fragments depending on collision model procedure :: set_mass_dist => abstract_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type procedure :: setup => collision_setup_system !! Initializer for the encounter collision system and the before/after snapshots procedure :: setup_impactors => collision_setup_impactors_system !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones @@ -134,9 +134,26 @@ module collision procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: reset => collision_util_reset_system !! Deallocates all allocatables procedure :: set_coordinate_system => collision_util_set_coordinate_system !! Sets the coordinate system of the collisional system - final :: collision_final_system !! Finalizer will deallocate all allocatables end type collision_system + type, extends(collision_system) :: collision_merge + contains + procedure :: generate_fragments => collision_generate_merge_system !! Merges the impactors to make a single final body + final :: collision_final_merge_system !! Finalizer will deallocate all allocatables + end type collision_merge + + type, extends(collision_system) :: collision_bounce + contains + procedure :: generate_fragments => collision_generate_bounce_system !! If a collision would result in a disruption, "bounce" the bodies instead. + final :: collision_final_bounce_system !! Finalizer will deallocate all allocatables + end type collision_bounce + + type, extends(collision_system) :: collision_simple + contains + procedure :: generate_fragments => collision_generate_simple_system !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012] + final :: collision_final_simple_system !! Finalizer will deallocate all allocatables + end type collision_simple + !! NetCDF dimension and variable names for the enounter save object type, extends(encounter_netcdf_parameters) :: collision_netcdf_parameters @@ -181,13 +198,13 @@ module collision abstract interface - subroutine abstract_generate_fragments(self, system, param, lfailure) - import collision_system, base_nbody_system, base_parameters + subroutine abstract_generate_fragments(self, system, param, t) + import collision_system, base_nbody_system, base_parameters, DP implicit none class(collision_system), intent(inout) :: self !! Collision system object class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters - logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? + real(DP), intent(in) :: t !! The time of the collision end subroutine abstract_generate_fragments subroutine abstract_set_mass_dist(self, param) @@ -200,6 +217,31 @@ end subroutine abstract_set_mass_dist interface + + module subroutine collision_generate_merge_system(self, system, param, t) + implicit none + class(collision_merge), intent(inout) :: self !! Fraggle fragment system object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + end subroutine collision_generate_merge_system + + module subroutine collision_generate_bounce_system(self, system, param, t) + implicit none + class(collision_bounce), intent(inout) :: self !! Fraggle fragment system object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + end subroutine collision_generate_bounce_system + + module subroutine collision_generate_simple_system(self, system, param, t) + implicit none + class(collision_simple), intent(inout) :: self !! Fraggle fragment system object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + end subroutine collision_generate_simple_system + module subroutine collision_netcdf_io_dump(self, param) implicit none class(collision_storage(*)), intent(inout) :: self !! Collision storage object @@ -354,7 +396,6 @@ module subroutine collision_util_reset_fragments(self) class(collision_fragments(*)), intent(inout) :: self end subroutine collision_util_reset_fragments - module subroutine collision_util_get_idvalues_snapshot(self, idvals) implicit none class(collision_snapshot), intent(in) :: self !! Fraggle snapshot object @@ -498,20 +539,52 @@ subroutine collision_final_storage(self) end subroutine collision_final_storage - subroutine collision_final_system(self) + subroutine collision_final_merge_system(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_merge), intent(inout) :: self !! Collision system object + + call self%reset() + if (allocated(self%impactors)) deallocate(self%impactors) + if (allocated(self%fragments)) deallocate(self%fragments) + + return + end subroutine collision_final_merge_system + + + subroutine collision_final_bounce_system(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_bounce), intent(inout) :: self !! Collision system object + + call self%reset() + if (allocated(self%impactors)) deallocate(self%impactors) + if (allocated(self%fragments)) deallocate(self%fragments) + + return + end subroutine collision_final_bounce_system + + + subroutine collision_final_simple_system(self) !! author: David A. Minton !! !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(collision_system), intent(inout) :: self !! Collision system object + type(collision_simple), intent(inout) :: self !! Collision system object call self%reset() if (allocated(self%impactors)) deallocate(self%impactors) if (allocated(self%fragments)) deallocate(self%fragments) return - end subroutine collision_final_system + end subroutine collision_final_simple_system diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index b74d76d5f..5aa170863 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -13,8 +13,116 @@ contains - subroutine collision_regime_collresolve(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & - regime, Mlr, Mslr, Qloss) + module subroutine collision_regime_impactors(self, system, param) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Determine which fragmentation regime the set of impactors will be. This subroutine is a wrapper for the non-polymorphic raggle_regime_collresolve subroutine. + !! It converts to SI units prior to calling + implicit none + ! Arguments + class(collision_impactors), intent(inout) :: self !! Collision system impactors object + class(base_nbody_system), intent(in) :: system !! Swiftest nbody system object + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + ! Internals + real (DP) :: mtot + + associate(impactors => self) + select type (system) + class is (swiftest_nbody_system) + select type(param) + class is (swiftest_parameters) + + select case(param%collision_model) + case("MERGE") + impactors%regime = COLLRESOLVE_REGIME_MERGE + mtot = sum(impactors%mass(:)) + impactors%mass_dist(1) = mtot + impactors%mass_dist(2) = 0.0_DP + impactors%mass_dist(3) = 0.0_DP + impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot + impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot + case default + call collision_regim_LS12(impactors, system, param) + end select + !call fraggle_io_log_regime(impactors, fragments) + end select + end select + end associate + + return + end subroutine collision_regime_impactors + + + subroutine collision_regime_LS12(impactors, nbody_system, param) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Determine the collisional regime of two colliding bodies based on the model by Leinhard and Stewart (2012) + !! + !! This is a wrapper subroutine that converts quantities to SI units and calls the main LS12 subroutine + implicit none + ! Arguments + class(collision_impactors), intent(inout) :: impactors !! The impactors to determine the regime for + class(swiftest_nbody_system), intent(in) :: nbody_system !! Swiftest n-body system object + class(swiftest_parameters), intent(in) :: param !! The current parameters + ! Internals + integer(I4B) :: jtarg, jproj + real(DP), dimension(2) :: radius_si, mass_si, density_si + real(DP) :: min_mfrag_si, Mcb_si + real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si, runit + real(DP) :: mlr, mslr, mtot, dentot + + ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them to the regime determination subroutine + if (impactors%mass(1) > impactors%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + mass_si(:) = impactors%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system + radius_si(:) = impactors%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system + density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The two-body equivalent density of the collider system + x1_si(:) = impactors%rb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system + v1_si(:) = impactors%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system + x2_si(:) = impactors%rb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system + v2_si(:) = impactors%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system + Mcb_si = nbody_system%cb%mass * param%MU2KG !! The central body mass of the system + min_mfrag_si = (param%min_GMfrag / param%GU) * param%MU2KG !! The minimum fragment mass to generate. Collider systems that would otherwise generate less massive fragments than this value will be forced to merge instead + + 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 collision_regime_collresolve(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), & + min_mfrag_si, impactors%regime, mlr, mslr, impactors%Qloss) + + if (allocated(impactors%mass_dist)) deallocate(impactors%mass_dist) + allocate(impactors%mass_dist(3)) + impactors%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) + impactors%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) + impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) + + ! Find the center of mass of the collisional system + mtot = sum(impactors%mass(:)) + impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot + impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot + + ! Find the point of impact between the two bodies + runit(:) = impactors%rb(:,2) - impactors%rb(:,1) + runit(:) = runit(:) / (.mag. runit(:)) + impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * runit(:) + + ! Convert quantities back to the system units and save them into the fragment system + impactors%mass_dist(:) = (impactors%mass_dist(:) / param%MU2KG) + impactors%Qloss = impactors%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG + + return + end subroutine collision_regime_LS12 + + + subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & + regime, Mlr, Mslr, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -182,8 +290,8 @@ subroutine collision_regime_collresolve(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, 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 !! @@ -305,89 +413,8 @@ function calc_c_star(Rc1) result(c_star) return end function calc_c_star - end subroutine collision_regime_collresolve - - - module subroutine collision_regime_impactors(self, system, param) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Determine which fragmentation regime the set of impactors will be. This subroutine is a wrapper for the non-polymorphic raggle_regime_collresolve subroutine. - !! It converts to SI units prior to calling - implicit none - ! Arguments - class(collision_impactors), intent(inout) :: self !! Collision system impactors object - class(base_nbody_system), intent(in) :: system !! Swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - ! Internals - integer(I4B) :: jtarg, jproj - real(DP), dimension(2) :: radius_si, mass_si, density_si - real(DP) :: min_mfrag_si, Mcb_si - real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si, runit - real(DP) :: mlr, mslr, mtot, dentot - - associate(impactors => self) - select type (system) - class is (swiftest_nbody_system) - ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them to the regime determination subroutine - if (impactors%mass(1) > impactors%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - mass_si(:) = impactors%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system - radius_si(:) = impactors%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system - density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The two-body equivalent density of the collider system - x1_si(:) = impactors%rb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system - v1_si(:) = impactors%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system - x2_si(:) = impactors%rb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system - v2_si(:) = impactors%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system - Mcb_si = system%cb%mass * param%MU2KG !! The central body mass of the system - select type(param) - class is (base_parameters) - min_mfrag_si = (param%min_GMfrag / param%GU) * param%MU2KG !! The minimum fragment mass to generate. Collider systems that would otherwise generate less massive fragments than this value will be forced to merge instead - class default - min_mfrag_si = 0.0_DP - end select - - 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 collision_regime_collresolve(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), & - min_mfrag_si, impactors%regime, mlr, mslr, impactors%Qloss) - - if (allocated(impactors%mass_dist)) deallocate(impactors%mass_dist) - allocate(impactors%mass_dist(3)) - impactors%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) - impactors%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) - impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - - ! Find the center of mass of the collisional system - mtot = sum(impactors%mass(:)) - impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot - impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot - - ! Find the point of impact between the two bodies - runit(:) = impactors%rb(:,2) - impactors%rb(:,1) - runit(:) = runit(:) / (.mag. runit(:)) - impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * runit(:) - - ! Convert quantities back to the system units and save them into the fragment system - impactors%mass_dist(:) = (impactors%mass_dist(:) / param%MU2KG) - impactors%Qloss = impactors%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG - - !call fraggle_io_log_regime(impactors, fragments) - end select - end associate - - return - end subroutine collision_regime_impactors + end subroutine collision_regime_LS12_SI - ! module subroutine collision_regime_resolve(self) - ! end subroutine collision_regime_resolve end submodule s_collision_regime \ No newline at end of file diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index ec3f8b67a..0edde0e3e 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -616,30 +616,11 @@ subroutine collision_resolve_list(plpl_collision , system, param, t) lgoodcollision = collision_resolve_consolidate_impactors(pl, cb, param, idx_parent, impactors) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle - if (param%lfragmentation) then - call impactors%get_regime(system, param) - else - impactors%regime = COLLRESOLVE_REGIME_MERGE - fragments%mtot = sum(impactors%mass(:)) - impactors%mass_dist(1) = fragments%mtot - impactors%mass_dist(2) = 0.0_DP - impactors%mass_dist(3) = 0.0_DP - impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot - impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot - end if - + call impactors%get_regime(system, param) call collision_history%take_snapshot(param,system, t, "before") - select case (impactors%regime) - case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - plpl_collision%status(i) = fraggle_resolve_disruption(system, param, t) - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - plpl_collision%status(i) = fraggle_resolve_hitandrun(system, param, t) - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - plpl_collision%status(i) = collision_resolve_merge(system, param, t) - case default - write(*,*) "Error in swiftest_collision, unrecognized collision regime" - call util_exit(FAILURE) - end select + + call collision_system%generate_fragments(system, param, t) + call collision_history%take_snapshot(param,system, t, "after") call impactors%reset() end do diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 8903c0ae1..16f93b652 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -18,7 +18,39 @@ contains - module subroutine fraggle_generate_fragments(self, system, param, lfailure) + module subroutine fraggle_generate_system(self, system, param, t) + implicit none + class(fraggle_system), intent(inout) :: self !! Fraggle fragment system object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! The time of the collision + ! Internals + integer(I4B) :: i + + select type(system) + class is (swiftest_nbody_system) + select type(param) + class is (swiftest_parameters) + associate(impactors => self%impactors, plpl_collision => system%plpl_collision) + select case (impactors%regime) + case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + plpl_collision%status(i) = fraggle_resolve_disruption(system, param, t) + case (COLLRESOLVE_REGIME_HIT_AND_RUN) + plpl_collision%status(i) = fraggle_resolve_hitandrun(system, param, t) + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + plpl_collision%status(i) = collision_resolve_merge(system, param, t) + case default + write(*,*) "Error in swiftest_collision, unrecognized collision regime" + call util_exit(FAILURE) + end select + end associate + end select + end select + + end subroutine fraggle_generate_system + + + subroutine fraggle_generate_fragments(self, system, param, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index e77dd44e2..df4f402cf 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -51,7 +51,7 @@ module fraggle real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains - procedure :: generate_fragments => fraggle_generate_fragments !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + procedure :: generate_fragments => fraggle_generate_system !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. procedure :: set_budgets => fraggle_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value procedure :: set_mass_dist => fraggle_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type procedure :: set_natural_scale => fraggle_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. @@ -65,14 +65,13 @@ module fraggle interface - module subroutine fraggle_generate_fragments(self, system, param, lfailure) - use base, only : base_nbody_system, base_parameters + module subroutine fraggle_generate_system(self, system, param, t) implicit none - class(fraggle_system), intent(inout) :: self !! Fraggle fragment system object + class(fraggle_system), intent(inout) :: self !! Fraggle fragment system object class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters - logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? - end subroutine fraggle_generate_fragments + real(DP), intent(in) :: t !! Time of collision + end subroutine fraggle_generate_system module subroutine fraggle_io_log_regime(collision_system) implicit none diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index efaf45d1f..456bef4a0 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -583,7 +583,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) implicit none ! Arguments class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to for writing a NetCDF dataset to file - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: nvar, varid, vartype real(DP) :: dfill @@ -1915,9 +1915,9 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i read(param_value, *, err = 667, iomsg = iomsg) param%maxid case ("MAXID_COLLISION") read(param_value, *, err = 667, iomsg = iomsg) param%maxid_collision - case ("FRAGMENTATION") + case ("COLLISION_MODEL") call swiftest_io_toupper(param_value) - if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true. + read(param_value, *) param%collision_model case ("GMTINY") read(param_value, *) param%GMTINY case ("MIN_GMFRAG") @@ -2076,13 +2076,24 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i param%lmtiny_pl = (integrator == INT_SYMBA) - if (param%lmtiny_pl .and. self%GMTINY < 0.0_DP) then - write(iomsg,*) "GMTINY invalid or not set: ", self%GMTINY + if (param%lmtiny_pl .and. param%GMTINY < 0.0_DP) then + write(iomsg,*) "GMTINY invalid or not set: ", param%GMTINY iostat = -1 return end if - if (param%lfragmentation) then + + if ((param%collision_model /= "MERGE") .and. & + (param%collision_model /= "SIMPLE") .and. & + (param%collision_model /= "BOUNCE") .and. & + (param%collision_model /= "FRAGGLE")) then + write(iomsg,*) 'Invalid collision_model parameter: ',trim(adjustl(param%out_type)) + write(iomsg,*) 'Valid options are NONE, TRAJECTORY, CLOSEST, or BOTH' + iostat = -1 + return + end if + + if (param%collision_model == "FRAGGLE") then if (seed_set) then call random_seed(put = param%seed) else diff --git a/src/swiftest/swiftest_setup.f90 b/src/swiftest/swiftest_setup.f90 index bbeafc6d2..ede789d7d 100644 --- a/src/swiftest/swiftest_setup.f90 +++ b/src/swiftest/swiftest_setup.f90 @@ -83,10 +83,12 @@ module subroutine swiftest_setup_construct_system(system, param) allocate(collision_list_plpl :: system%plpl_encounter) allocate(collision_list_plpl :: system%plpl_collision) - if (param%lfragmentation) then + if (param%collision_model == "FRAGGLE") then allocate(fraggle_system :: system%collision_system) - call system%collision_system%setup(system) + else + allocate(collision_system :: system%collision_system) end if + call system%collision_system%setup(system) if (param%lenc_save_trajectory .or. param%lenc_save_closest) then allocate(encounter_netcdf_parameters :: encounter_history%nc)