From 28795c71e161b5ccdce6d32138036d5ceea5b9ba Mon Sep 17 00:00:00 2001 From: Andrew Benson Date: Fri, 7 Jun 2024 16:01:20 -0700 Subject: [PATCH] feat: Add functionality to record (and output) the orbital properties of merged subhalos --- ...ysics.satellite_merging.radius_trigger.F90 | 92 +++++++- ...ty_extractor.merged_subhalo_properties.F90 | 215 ++++++++++++++++++ 2 files changed, 296 insertions(+), 11 deletions(-) create mode 100644 source/nodes.property_extractor.merged_subhalo_properties.F90 diff --git a/source/nodes.operators.physics.satellite_merging.radius_trigger.F90 b/source/nodes.operators.physics.satellite_merging.radius_trigger.F90 index 26d7d557ab..2a186276a5 100644 --- a/source/nodes.operators.physics.satellite_merging.radius_trigger.F90 +++ b/source/nodes.operators.physics.satellite_merging.radius_trigger.F90 @@ -23,6 +23,7 @@ use :: Dark_Matter_Halo_Scales, only : darkMatterHaloScaleClass use :: Galactic_Structure , only : galacticStructureClass + use :: Kepler_Orbits , only : keplerOrbitCount !![ @@ -34,9 +35,11 @@ A node operator class that triggers merging of satellites based on their orbital radius. !!} private - class (darkMatterHaloScaleClass), pointer :: darkMatterHaloScale_ => null() - class (galacticStructureClass ), pointer :: galacticStructure_ => null() + class (darkMatterHaloScaleClass), pointer :: darkMatterHaloScale_ => null() + class (galacticStructureClass ), pointer :: galacticStructure_ => null() double precision :: radiusVirialFraction + logical :: recordMergedSubhaloProperties + integer :: mergedSubhaloIDs (keplerOrbitCount) contains !![ @@ -55,6 +58,10 @@ module procedure satelliteMergingRadiusTriggerConstructorParameters module procedure satelliteMergingRadiusTriggerConstructorInternal end interface nodeOperatorSatelliteMergingRadiusTrigger + + ! Sub-module-scope pointer to self used in callback function. + class(nodeOperatorSatelliteMergingRadiusTrigger), pointer :: self_ + !$omp threadprivate(self) contains @@ -69,6 +76,7 @@ function satelliteMergingRadiusTriggerConstructorParameters(parameters) result(s class (darkMatterHaloScaleClass ), pointer :: darkMatterHaloScale_ class (galacticStructureClass ), pointer :: galacticStructure_ double precision :: radiusVirialFraction + logical :: recordMergedSubhaloProperties !![ @@ -77,10 +85,16 @@ function satelliteMergingRadiusTriggerConstructorParameters(parameters) result(s The fraction of the virial radius below which satellites are merged. parameters + + recordMergedSubhaloProperties + .false. + If true, record the orbital properties of subhalo that merge. + parameters + !!] - self=nodeOperatorSatelliteMergingRadiusTrigger(radiusVirialFraction,darkMatterHaloScale_,galacticStructure_) + self=nodeOperatorSatelliteMergingRadiusTrigger(radiusVirialFraction,recordMergedSubhaloProperties,darkMatterHaloScale_,galacticStructure_) !![ @@ -89,19 +103,31 @@ function satelliteMergingRadiusTriggerConstructorParameters(parameters) result(s return end function satelliteMergingRadiusTriggerConstructorParameters - function satelliteMergingRadiusTriggerConstructorInternal(radiusVirialFraction,darkMatterHaloScale_,galacticStructure_) result(self) + function satelliteMergingRadiusTriggerConstructorInternal(radiusVirialFraction,recordMergedSubhaloProperties,darkMatterHaloScale_,galacticStructure_) result(self) !!{ Internal constructor for the {\normalfont \ttfamily satelliteMergingRadiusTrigger} node operator class. !!} + use :: Kepler_Orbits, only : keplerOrbitTimeInitial , keplerOrbitMassSatellite, keplerOrbitMassHost, keplerOrbitRadius, & + & keplerOrbitRadiusPericenter implicit none type (nodeOperatorSatelliteMergingRadiusTrigger) :: self double precision , intent(in ) :: radiusVirialFraction + logical , intent(in ) :: recordMergedSubhaloProperties class (darkMatterHaloScaleClass ), intent(in ), target :: darkMatterHaloScale_ class (galacticStructureClass ), intent(in ), target :: galacticStructure_ !![ - + !!] + if (recordMergedSubhaloProperties) then + !![ + + + + + + !!] + end if return end function satelliteMergingRadiusTriggerConstructorInternal @@ -111,7 +137,7 @@ subroutine satelliteMergingRadiusTriggerDestructor(self) !!} implicit none type(nodeOperatorSatelliteMergingRadiusTrigger), intent(inout) :: self - + !![ @@ -149,6 +175,7 @@ subroutine satelliteMergingRadiusTriggerDifferentialEvolution(self,node,interrup ! Merging criterion met - trigger an interrupt. interrupt = .true. functionInterrupt => mergerTrigger + self_ => self end if return end subroutine satelliteMergingRadiusTriggerDifferentialEvolution @@ -157,17 +184,60 @@ subroutine mergerTrigger(node,timeEnd) !!{ Trigger a merger of the satellite by setting the time until merging to zero. !!} - use :: Galacticus_Nodes, only : nodeComponentSatellite, nodeComponentBasic, treeNode + use :: Galacticus_Nodes, only : nodeComponentSatellite , nodeComponentBasic , treeNode + use :: Kepler_Orbits , only : keplerOrbit , keplerOrbitTimeInitial, keplerOrbitMassSatellite, keplerOrbitMassHost, & + & keplerOrbitRadiusPericenter, keplerOrbitRadius implicit none - type (treeNode ), intent(inout), target :: node - double precision , intent(in ), optional :: timeEnd - class (nodeComponentBasic ) , pointer :: basic - class (nodeComponentSatellite) , pointer :: satellite + type (treeNode ), intent(inout), target :: node + double precision , intent(in ), optional :: timeEnd + type (treeNode ) , pointer :: nodeHost + class (nodeComponentBasic ) , pointer :: basic , basicHost + class (nodeComponentSatellite) , pointer :: satellite + double precision , dimension(:) , allocatable :: propertyCurrent, propertyNew + double precision :: property + type (keplerOrbit ) :: orbit + integer :: i , ID !$GLC attributes unused :: timeEnd + ! Set the time of merging to the current time. basic => node%basic () satellite => node%satellite() call satellite%timeOfMergingSet(basic%time()) + ! Record properties of the merging subhalo if necessary. + if (self_%recordMergedSubhaloProperties) then + ! Find the node to merge with. + nodeHost => node %mergesWith() + basicHost => nodeHost%basic () + ! Get the virial orbit of the halo about to merge. + orbit=satellite%virialOrbit() + ! Append the orbit data. + do i=1,5 + select case (i) + case (1) + ID =keplerOrbitTimeInitial %ID + property=basic%timeLastIsolated() + case (2) + ID =keplerOrbitMassSatellite %ID + property=orbit%massSatellite () + case (3) + ID =keplerOrbitMassHost %ID + property=orbit%massHost () + case (4) + ID =keplerOrbitRadius %ID + property=orbit%radius () + case (5) + ID =keplerOrbitRadiusPericenter%ID + property=orbit%radiusPericenter() + end select + propertyCurrent=basicHost%floatRank1MetaPropertyGet(self_%mergedSubhaloIDs(ID)) + allocate(propertyNew(size(propertyCurrent)+1_c_size_t)) + propertyNew(1_c_size_t:size(propertyCurrent))=propertyCurrent(:) + propertyNew(size(propertyNew))=property + call basicHost%floatRank1MetaPropertySet(self_%mergedSubhaloIDs(ID),propertyNew) + deallocate(propertyCurrent) + deallocate(propertyNew ) + end do + end if return end subroutine mergerTrigger diff --git a/source/nodes.property_extractor.merged_subhalo_properties.F90 b/source/nodes.property_extractor.merged_subhalo_properties.F90 new file mode 100644 index 0000000000..453da45e56 --- /dev/null +++ b/source/nodes.property_extractor.merged_subhalo_properties.F90 @@ -0,0 +1,215 @@ +!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, +!! 2019, 2020, 2021, 2022, 2023, 2024 +!! Andrew Benson +!! +!! This file is part of Galacticus. +!! +!! Galacticus 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. +!! +!! Galacticus 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 Galacticus. If not, see . + + use :: Kepler_Orbits, only : keplerOrbitCount + + !![ + + + A node property extractor which extracts properties of merged subhalo orbits. + + + !!] + type, extends(nodePropertyExtractorList) :: nodePropertyExtractorMergedSubhaloProperties + !!{ + A property extractor which extracts properties of merged subhalo orbits. + !!} + private + integer :: mergedSubhaloIDs(keplerOrbitCount) + contains + procedure :: elementCount => mergedSubhaloPropertiesElementCount + procedure :: extract => mergedSubhaloPropertiesExtract + procedure :: names => mergedSubhaloPropertiesNames + procedure :: descriptions => mergedSubhaloPropertiesDescriptions + procedure :: unitsInSI => mergedSubhaloPropertiesUnitsInSI + end type nodePropertyExtractorMergedSubhaloProperties + + interface nodePropertyExtractorMergedSubhaloProperties + !!{ + Constructors for the ``mergedSubhaloProperties'' output extractor class. + !!} + module procedure mergedSubhaloPropertiesConstructorParameters + module procedure mergedSubhaloPropertiesConstructorInternal + end interface nodePropertyExtractorMergedSubhaloProperties + +contains + + function mergedSubhaloPropertiesConstructorParameters(parameters) result(self) + !!{ + Constructor for the ``mergedSubhaloProperties'' property extractor class which takes a parameter set as input. + !!} + use :: Input_Parameters, only : inputParameter, inputParameters + implicit none + type(nodePropertyExtractorMergedSubhaloProperties) :: self + type(inputParameters ), intent(inout) :: parameters + + self=nodePropertyExtractorMergedSubhaloProperties() + !![ + + !!] + return + end function mergedSubhaloPropertiesConstructorParameters + + function mergedSubhaloPropertiesConstructorInternal() result(self) + !!{ + Internal constructor for the ``mergedSubhaloProperties'' output extractor property extractor class. + !!} + use :: Kepler_Orbits, only : keplerOrbitTimeInitial , keplerOrbitMassSatellite, keplerOrbitMassHost, keplerOrbitRadius, & + & keplerOrbitRadiusPericenter + implicit none + type(nodePropertyExtractorMergedSubhaloProperties) :: self + + !![ + + + + + + !!] + return + end function mergedSubhaloPropertiesConstructorInternal + + integer function mergedSubhaloPropertiesElementCount(self) + !!{ + Return a count of the number of properties extracted. + !!} + implicit none + class(nodePropertyExtractorMergedSubhaloProperties), intent(inout) :: self + + mergedSubhaloPropertiesElementCount=5 + return + end function mergedSubhaloPropertiesElementCount + + function mergedSubhaloPropertiesExtract(self,node,instance) result(propertiesOrbitalMergedSubhalos) + !!{ + Implement a mergedSubhaloProperties output extractor. + !!} + use :: Galacticus_Nodes, only : nodeComponentBasic + use :: Kepler_Orbits , only : keplerOrbitTimeInitial, keplerOrbitMassSatellite, keplerOrbitMassHost, keplerOrbitRadiusPericenter, & + & keplerOrbitRadius + implicit none + double precision , dimension(:,:), allocatable :: propertiesOrbitalMergedSubhalos + class (nodePropertyExtractorMergedSubhaloProperties), intent(inout) :: self + type (treeNode ), intent(inout) :: node + type (multiCounter ), intent(inout) , optional :: instance + class (nodeComponentBasic ) , pointer :: basic + double precision , dimension(: ), allocatable :: property + integer :: i , ID + !$GLC attributes unused :: instance + !$GLC attributes initialized :: timesLastIsolated + + basic => node%basic() + do i=1,5 + select case (i) + case (1) + ID=keplerOrbitTimeInitial %ID + case (2) + ID=keplerOrbitMassSatellite %ID + case (3) + ID=keplerOrbitMassHost %ID + case (4) + ID=keplerOrbitRadius %ID + case (5) + ID=keplerOrbitRadiusPericenter%ID + end select + property=basic%floatRank1MetaPropertyGet(self%mergedSubhaloIDs(ID)) + if (.not.allocated(propertiesOrbitalMergedSubhalos)) allocate(propertiesOrbitalMergedSubhalos(size(property),5)) + propertiesOrbitalMergedSubhalos(:,i)=property + deallocate(property) + end do + return + end function mergedSubhaloPropertiesExtract + + subroutine mergedSubhaloPropertiesNames(self,names) + !!{ + Return the names of the {\normalfont \ttfamily mergedSubhaloProperties} properties. + !!} + use :: Kepler_Orbits , only : enumerationKeplerOrbitDecode, keplerOrbitTimeInitial, keplerOrbitMassSatellite, keplerOrbitMassHost, & + & keplerOrbitRadiusPericenter , keplerOrbitRadius + use :: String_Handling, only : String_Upper_Case_First + implicit none + class (nodePropertyExtractorMergedSubhaloProperties), intent(inout) :: self + type (varying_string ), intent(inout), dimension(:) , allocatable :: names + integer :: i , ID + !$GLC attributes unused :: self + + allocate(names(5)) + do i=1,5 + select case (i) + case (1) + ID=keplerOrbitTimeInitial %ID + case (2) + ID=keplerOrbitMassSatellite %ID + case (3) + ID=keplerOrbitMassHost %ID + case (4) + ID=keplerOrbitRadius %ID + case (5) + ID=keplerOrbitRadiusPericenter%ID + end select + names(i)=var_str('mergedSubhalo')//String_Upper_Case_First(char(enumerationKeplerOrbitDecode(ID,includePrefix=.false.))) + end do + return + end subroutine mergedSubhaloPropertiesNames + + subroutine mergedSubhaloPropertiesDescriptions(self,descriptions) + !!{ + Return the descriptions of the {\normalfont \ttfamily mergedSubhaloProperties} properties. + !!} + use :: Kepler_Orbits, only : enumerationKeplerOrbitDescription, keplerOrbitTimeInitial, keplerOrbitMassSatellite, keplerOrbitMassHost, & + & keplerOrbitRadiusPericenter , keplerOrbitRadius + implicit none + class (nodePropertyExtractorMergedSubhaloProperties), intent(inout) :: self + type (varying_string ), intent(inout), dimension(:) , allocatable :: descriptions + integer :: i , ID + !$GLC attributes unused :: self + + allocate(descriptions(5)) + do i=1,5 + select case (i) + case (1) + ID=keplerOrbitTimeInitial %ID + case (2) + ID=keplerOrbitMassSatellite %ID + case (3) + ID=keplerOrbitMassHost %ID + case (4) + ID=keplerOrbitRadius %ID + case (5) + ID=keplerOrbitRadiusPericenter%ID + end select + descriptions(i)=var_str('Merged subhalos: ')//enumerationKeplerOrbitDescription(ID) + end do + return + end subroutine mergedSubhaloPropertiesDescriptions + + function mergedSubhaloPropertiesUnitsInSI(self) result(unitsInSI) + !!{ + Return the units of the {\normalfont \ttfamily mergedSubhaloProperties} properties in the SI system. + !!} + use :: Numerical_Constants_Astronomical, only : gigaYear, massSolar, megaParsec + implicit none + double precision , dimension(:) , allocatable :: unitsInSI + class (nodePropertyExtractorMergedSubhaloProperties), intent(inout) :: self + !$GLC attributes unused :: self + + allocate(unitsInSI(5)) + unitsInSI=[gigaYear,massSolar,massSolar,megaParsec,megaParsec] + return + end function mergedSubhaloPropertiesUnitsInSI