Skip to content

Commit

Permalink
feat: Add functionality to record (and output) the orbital properties…
Browse files Browse the repository at this point in the history
… of merged subhalos
  • Loading branch information
abensonca committed Jun 7, 2024
1 parent d3a25c9 commit 28795c7
Show file tree
Hide file tree
Showing 2 changed files with 296 additions and 11 deletions.
92 changes: 81 additions & 11 deletions source/nodes.operators.physics.satellite_merging.radius_trigger.F90
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@

use :: Dark_Matter_Halo_Scales, only : darkMatterHaloScaleClass
use :: Galactic_Structure , only : galacticStructureClass
use :: Kepler_Orbits , only : keplerOrbitCount

!![
<nodeOperator name="nodeOperatorSatelliteMergingRadiusTrigger">
Expand All @@ -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
!![
<methods>
Expand All @@ -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

Expand All @@ -69,6 +76,7 @@ function satelliteMergingRadiusTriggerConstructorParameters(parameters) result(s
class (darkMatterHaloScaleClass ), pointer :: darkMatterHaloScale_
class (galacticStructureClass ), pointer :: galacticStructure_
double precision :: radiusVirialFraction
logical :: recordMergedSubhaloProperties

!![
<inputParameter>
Expand All @@ -77,10 +85,16 @@ function satelliteMergingRadiusTriggerConstructorParameters(parameters) result(s
<description>The fraction of the virial radius below which satellites are merged.</description>
<source>parameters</source>
</inputParameter>
<inputParameter>
<name>recordMergedSubhaloProperties</name>
<defaultValue>.false.</defaultValue>
<description>If true, record the orbital properties of subhalo that merge.</description>
<source>parameters</source>
</inputParameter>
<objectBuilder class="darkMatterHaloScale" name="darkMatterHaloScale_" source="parameters"/>
<objectBuilder class="galacticStructure" name="galacticStructure_" source="parameters"/>
!!]
self=nodeOperatorSatelliteMergingRadiusTrigger(radiusVirialFraction,darkMatterHaloScale_,galacticStructure_)
self=nodeOperatorSatelliteMergingRadiusTrigger(radiusVirialFraction,recordMergedSubhaloProperties,darkMatterHaloScale_,galacticStructure_)
!![
<inputParametersValidate source="parameters"/>
<objectDestructor name="darkMatterHaloScale_"/>
Expand All @@ -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_
!![
<constructorAssign variables="radiusVirialFraction, *darkMatterHaloScale_, *galacticStructure_"/>
<constructorAssign variables="radiusVirialFraction, recordMergedSubhaloProperties, *darkMatterHaloScale_, *galacticStructure_"/>
!!]

if (recordMergedSubhaloProperties) then
!![
<addMetaProperty component="basic" name="mergedSubhaloTimeInitial" id="self%mergedSubhaloIDs(keplerOrbitTimeInitial %ID)" rank="1" isCreator="yes"/>
<addMetaProperty component="basic" name="mergedSubhaloMassSatellite" id="self%mergedSubhaloIDs(keplerOrbitMassSatellite %ID)" rank="1" isCreator="yes"/>
<addMetaProperty component="basic" name="mergedSubhaloMassHost" id="self%mergedSubhaloIDs(keplerOrbitMassHost %ID)" rank="1" isCreator="yes"/>
<addMetaProperty component="basic" name="mergedSubhaloRadius" id="self%mergedSubhaloIDs(keplerOrbitRadius %ID)" rank="1" isCreator="yes"/>
<addMetaProperty component="basic" name="mergedSubhaloRadiusPericenter" id="self%mergedSubhaloIDs(keplerOrbitRadiusPericenter%ID)" rank="1" isCreator="yes"/>
!!]
end if
return
end function satelliteMergingRadiusTriggerConstructorInternal

Expand All @@ -111,7 +137,7 @@ subroutine satelliteMergingRadiusTriggerDestructor(self)
!!}
implicit none
type(nodeOperatorSatelliteMergingRadiusTrigger), intent(inout) :: self

!![
<objectDestructor name="self%darkMatterHaloScale_"/>
<objectDestructor name="self%galacticStructure_" />
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down
215 changes: 215 additions & 0 deletions source/nodes.property_extractor.merged_subhalo_properties.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,215 @@
!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018,
!! 2019, 2020, 2021, 2022, 2023, 2024
!! Andrew Benson <[email protected]>
!!
!! 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 <http://www.gnu.org/licenses/>.

use :: Kepler_Orbits, only : keplerOrbitCount

!![
<nodePropertyExtractor name="nodePropertyExtractorMergedSubhaloProperties">
<description>
A node property extractor which extracts properties of merged subhalo orbits.
</description>
</nodePropertyExtractor>
!!]
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()
!![
<inputParametersValidate source="parameters"/>
!!]
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

!![
<addMetaProperty component="basic" name="mergedSubhaloTimeInitial" id="self%mergedSubhaloIDs(keplerOrbitTimeInitial %ID)" rank="1" isCreator="no"/>
<addMetaProperty component="basic" name="mergedSubhaloMassSatellite" id="self%mergedSubhaloIDs(keplerOrbitMassSatellite %ID)" rank="1" isCreator="no"/>
<addMetaProperty component="basic" name="mergedSubhaloMassHost" id="self%mergedSubhaloIDs(keplerOrbitMassHost %ID)" rank="1" isCreator="no"/>
<addMetaProperty component="basic" name="mergedSubhaloRadius" id="self%mergedSubhaloIDs(keplerOrbitRadius %ID)" rank="1" isCreator="no"/>
<addMetaProperty component="basic" name="mergedSubhaloRadiusPericenter" id="self%mergedSubhaloIDs(keplerOrbitRadiusPericenter%ID)" rank="1" isCreator="no"/>
!!]
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

0 comments on commit 28795c7

Please sign in to comment.