Skip to content

Commit

Permalink
Merge pull request #686 from cgannonucm/UniverseMachine-fix
Browse files Browse the repository at this point in the history
fix spheroidRadiusPowerLaw nodeOperator
  • Loading branch information
abensonca authored Sep 1, 2024
2 parents 46273b7 + 6815728 commit bcf66cb
Show file tree
Hide file tree
Showing 3 changed files with 244 additions and 15 deletions.
33 changes: 18 additions & 15 deletions source/nodes.operators.empirical.spheroid_radius_power_law.F90
Original file line number Diff line number Diff line change
Expand Up @@ -25,13 +25,13 @@
!!}

!![
<nodeOperator name="nodeOperatorspheroidRadiusPowerLaw">
<nodeOperator name="nodeOperatorSpheroidRadiusPowerLaw">
<description>
A node operator that implements an an empirical power law relationship between spheroid stellar radius and stellar mass.
</description>
</description>
</nodeOperator>
!!]
type, extends(nodeOperatorClass) :: nodeOperatorspheroidRadiusPowerLaw
type, extends(nodeOperatorClass) :: nodeOperatorSpheroidRadiusPowerLaw
!!{
Implements a power law prescription for the stellar mass--stellar radius relation of spheroids. Specificially:
\begin{equation}
Expand All @@ -52,15 +52,15 @@
procedure :: nodeInitialize => spheroidRadiusPowerLawNodeInitialize
procedure :: differentialEvolutionSolveAnalytics => spheroidRadiusPowerLawSolveAnalytics
procedure :: nodesMerge => spheroidRadiusPowerLawNodesMerge
end type nodeOperatorspheroidRadiusPowerLaw
end type nodeOperatorSpheroidRadiusPowerLaw

interface nodeOperatorspheroidRadiusPowerLaw
interface nodeOperatorSpheroidRadiusPowerLaw
!!{
Constructors for the {\normalfont \ttfamily spheroidRadiusPowerLaw} node operator class.
!!}
module procedure spheroidRadiusPowerLawConstructorParameters
module procedure spheroidRadiusPowerLawConstructorInternal
end interface nodeOperatorspheroidRadiusPowerLaw
end interface nodeOperatorSpheroidRadiusPowerLaw

contains

Expand All @@ -70,7 +70,7 @@ function spheroidRadiusPowerLawConstructorParameters(parameters) result(self)
!!}
use :: Input_Parameters, only : inputParameters
implicit none
type (nodeOperatorspheroidRadiusPowerLaw) :: self
type (nodeOperatorSpheroidRadiusPowerLaw) :: self
type (inputParameters ), intent(inout) :: parameters
double precision :: alpha , beta

Expand All @@ -86,23 +86,26 @@ function spheroidRadiusPowerLawConstructorParameters(parameters) result(self)
<name>beta</name>
<source>parameters</source>
<description>Coefficient $\beta$ in the power law fit.</description>
<defaultValue>1.19d-6</defaultValue>
<defaultValue>1.19d-9</defaultValue>
<defaultSource>\cite[][table J1: Parameter $b$, for early type galaxies and re-scaled from half light radius to Hernquist radius \protect\citep{hernquist_analytical_1990}---Note: there was a typo in the originally provided value, see \protect\cite{shen_erratum_2007} for the corrected value]{shen_size_2003}</defaultSource>
</inputParameter>
!!]
self=spheroidRadiusPowerLawConstructorInternal(alpha, beta)
!![
<inputParametersValidate source="parameters"/>
!!]
end function spheroidRadiusPowerLawConstructorParameters

function spheroidRadiusPowerLawConstructorInternal(alpha, beta) result(self)
!!{
Internal constructor for the {\normalfont \ttfamily spheroidRadiusPowerLaw} node operator class.
!!}
implicit none
type (nodeOperatorspheroidRadiusPowerLaw) :: self
type (nodeOperatorSpheroidRadiusPowerLaw) :: self
double precision , intent(in) :: alpha, beta
!![
<constructorAssign variables="alpha, beta"/>
!!]

return
end function spheroidRadiusPowerLawConstructorInternal

Expand All @@ -112,15 +115,15 @@ subroutine spheroidRadiusPowerLawUpdate(self, node)
!!}
use :: Galacticus_Nodes, only : nodeComponentSpheroid
implicit none
class (nodeOperatorspheroidRadiusPowerLaw), intent(inout) :: self
class (nodeOperatorSpheroidRadiusPowerLaw), intent(inout) :: self
type (treeNode ), intent(inout) :: node
class (nodeComponentSpheroid ), pointer :: spheroid
double precision :: radiusStellar

if (.not.node%isOnMainBranch()) return
spheroid => node %spheroid ()
radiusStellar = +self %beta &
& +spheroid%massStellar()**self%alpha
& *spheroid%massStellar()**self%alpha
call spheroid%radiusSet(radiusStellar)
return
end subroutine spheroidRadiusPowerLawUpdate
Expand All @@ -131,7 +134,7 @@ subroutine spheroidRadiusPowerLawNodeInitialize(self,node)
!!}
use :: Galacticus_Nodes, only : nodeComponentSpheroid
implicit none
class(nodeOperatorspheroidRadiusPowerLaw), intent(inout), target :: self
class(nodeOperatorSpheroidRadiusPowerLaw), intent(inout), target :: self
type (treeNode ), intent(inout), target :: node
class (nodeComponentSpheroid ), pointer :: spheroid

Expand All @@ -147,7 +150,7 @@ subroutine spheroidRadiusPowerLawSolveAnalytics(self,node,time)
Set the radius of the spheroid.
!!}
implicit none
class (nodeOperatorspheroidRadiusPowerLaw), intent(inout) :: self
class (nodeOperatorSpheroidRadiusPowerLaw), intent(inout) :: self
type (treeNode ), intent(inout) :: node
double precision , intent(in ) :: time
!$GLC attributes unused :: time
Expand All @@ -161,7 +164,7 @@ subroutine spheroidRadiusPowerLawNodesMerge(self,node)
Update the radius of the spheroid after a merger.
!!}
implicit none
class(nodeOperatorspheroidRadiusPowerLaw), intent(inout) :: self
class(nodeOperatorSpheroidRadiusPowerLaw), intent(inout) :: self
type (treeNode ), intent(inout) :: node

call self%update(node)
Expand Down
168 changes: 168 additions & 0 deletions testSuite/parameters/test-spheroid-power-law.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,168 @@
<?xml version="1.0" encoding="UTF-8"?>
<!--Parameters for testing univerMachine nodeOperator -->
<parameters>
<formatVersion>2</formatVersion>
<version>0.9.4</version>
<randomNumberGenerator value="GSL">
<seed value="153" />
</randomNumberGenerator>

<!-- Specify tasks to perform -->
<task value="evolveForests"/>

<!-- Component selection -->
<componentBasic value="standard"/>
<componentBlackHole value="null"/>
<componentDarkMatterProfile value="scale"/>
<componentDisk value="standard"/>
<componentHotHalo value="null"/>
<componentSatellite value="standard"/>
<componentSpheroid value="standard"/>
<componentSpin value="null"/>

<!-- Cosmological parameters and options -->
<cosmologyFunctions value="matterLambda"/>
<cosmologyParameters value="simple">
<HubbleConstant value="70.20000"/>
<OmegaMatter value=" 0.27250"/>
<OmegaDarkEnergy value=" 0.72750"/>
<OmegaBaryon value=" 0.00000"/>
<temperatureCMB value=" 2.72548"/>
</cosmologyParameters>

<!-- Power spectrum options -->
<transferFunction value="eisensteinHu1999">
<neutrinoNumberEffective value="3.046"/>
<neutrinoMassSummed value="0.000"/>
<cosmologyParameters value="simple">
<HubbleConstant value="70.20000"/>
<OmegaMatter value=" 0.27250"/>
<OmegaDarkEnergy value=" 0.72750"/>
<OmegaBaryon value=" 0.04550"/>
<temperatureCMB value=" 2.72548"/>
</cosmologyParameters>
</transferFunction>
<powerSpectrumPrimordial value="powerLaw">
<index value="0.961"/>
<wavenumberReference value="1.000"/>
<running value="0.000"/>
</powerSpectrumPrimordial>
<powerSpectrumPrimordialTransferred value="simple"/>
<cosmologicalMassVariance value="filteredPower">
<sigma_8 value="0.807"/>
</cosmologicalMassVariance>

<!-- Structure formation options -->
<linearGrowth value="collisionlessMatter"/>
<haloMassFunction value="tinker2008"/>
<criticalOverdensity value="sphericalCollapseClsnlssMttrCsmlgclCnstnt"/>
<virialDensityContrast value="sphericalCollapseClsnlssMttrCsmlgclCnstnt"/>

<!-- Merger tree building options -->
<mergerTreeConstructor value="build">
<redshiftBase value="0.0" />
</mergerTreeConstructor>
<mergerTreeBuilder value="cole2000">
<accretionLimit value="0.1"/>
<mergeProbability value="0.1"/>
</mergerTreeBuilder>
<mergerTreeBranchingProbability value="parkinsonColeHelly">
<G0 value="+0.57"/>
<gamma1 value="+0.38"/>
<gamma2 value="-0.01"/>
<accuracyFirstOrder value="+0.10"/>
</mergerTreeBranchingProbability>

<mergerTreeBuildMasses value="readXML">
<fileName value="testSuite/parameters/UniverseMachine-treemass.xml" />
</mergerTreeBuildMasses>

<mergerTreeMassResolution value="scaled">
<massResolutionFractional value="1.0e-1"/>
</mergerTreeMassResolution>


<!-- Substructure hierarchy options -->
<mergerTreeNodeMerger value="singleLevelHierarchy"/>

<!-- Dark matter halo structure options -->
<darkMatterProfileDMO value="NFW"/>
<darkMatterProfileConcentration value="gao2008"/>
<darkMatterProfileScaleRadius value="concentrationLimiter">
<concentrationMinimum value=" 4.0"/>
<concentrationMaximum value="100.0"/>
<darkMatterProfileScaleRadius value="concentration"/>
</darkMatterProfileScaleRadius>

<!-- Switch off baryonic physics -->
<hotHaloMassDistribution value="null"/>

<!-- Tree evolution -->
<mergerTreeNodeEvolver value="standard">
<odeToleranceAbsolute value="0.001"/>
<odeToleranceRelative value="0.001"/>
</mergerTreeNodeEvolver>
<mergerTreeEvolver value="standard">
<timestepHostAbsolute value="0.01"/>
<timestepHostRelative value="0.01"/>
</mergerTreeEvolver>

<!-- Output options -->
<outputFileName value="testSuite/outputs/test-spheroid-power-law/test-spheroid-power-law.hdf5"/>
<outputTimes value="list">
<redshifts value="0.0"/>
</outputTimes>
<mergerTreeOutputter value="standard">
<outputReferences value="false"/>
</mergerTreeOutputter>
<nodePropertyExtractor value="multi">
<nodePropertyExtractor value="massHalo">
<virialDensityContrastDefinition value="bryanNorman1998" />
</nodePropertyExtractor>
<nodePropertyExtractor value="nodeIndices" />
<nodePropertyExtractor value="redshiftLastIsolated" />
</nodePropertyExtractor>

<!-- Node physics -->
<nodeOperator value="multi">
<!-- Cosmological epoch -->
<nodeOperator value="cosmicTime"/>
<!-- DMO evolution -->
<nodeOperator value="DMOInterpolate"/>
<!-- Halo concentrations -->
<nodeOperator value="darkMatterProfileScaleSet" />
<nodeOperator value="darkMatterProfileScaleInterpolate"/>
<!-- Satellite evolution -->
<nodeOperator value="satelliteMergingTime"/>
<nodeOperator value="satelliteMassLoss" />
<nodeOperator value="empiricalGalaxyUniverseMachine">
<massStellarFinal value="-1.0d0" />
<fractionMassSpheroid value="1.0d0" />
<fractionMassDisk value="0.0d0" />
<epsilon_0 value="-1.435" />
<epsilon_a value= "+1.831" />
<epsilon_lna value= "+1.368"/>
<epsilon_z value= "-0.217"/>
<M_0 value="+12.035"/>
<M_a value="+4.556"/>
<M_lna value="+4.417"/>
<M_z value="-0.731"/>
<alpha_0 value="+1.963"/>
<alpha_a value="-2.316"/>
<alpha_lna value="-1.733"/>
<alpha_z value="+0.178"/>
<beta_0 value= "+0.482"/>
<beta_a value = "-0.841"/>
<beta_z value= "-0.471"/>
<delta_0 value="+0.411"/>
<gamma_0 value="-1.034"/>
<gamma_a value="-3.10"/>
<gamma_z value="-1.055"/>
</nodeOperator>
<nodeOperator value="spheroidRadiusPowerLaw">
<alpha value="0.56" />
<beta value="1.19E-9"/>
</nodeOperator>
</nodeOperator>

</parameters>
58 changes: 58 additions & 0 deletions testSuite/test-spheroid-power-law.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
#!/usr/bin/env python3
import subprocess
import os
import sys
import h5py
import numpy as np
import xml.etree.ElementTree as ET

# Test spheroidRadiusPowerLaw nodeOperator.
# Charles Gannon (31-August-2024)
path_in_param = os.path.abspath(f"parameters/test-spheroid-power-law.xml" )
path_out_dir = os.path.abspath(f"outputs/test-spheroid-power-law" )
path_out_log = os.path.abspath(f"outputs/test-spheroid-power-law/test-spheroid-power-law.log" )
path_out_hdf = os.path.abspath(f"outputs/test-spheroid-power-law/test-spheroid-power-law.hdf5" )

status = subprocess.run(f"mkdir -p {path_out_dir}",shell=True)

# Run the model and check for completion
print("Running model...")

log = open(path_out_log,"w")
status = subprocess.run(f"cd ..; ./Galacticus.exe {path_in_param}",stdout=log,stderr=log,shell=True)
log.close()

print("...done ("+str(status)+")")

if status.returncode != 0:
print("FAILED: model run:")
subprocess.run(f"cat {path_out_log}",shell=True)
sys.exit()
print("Checking for errors...")
status = subprocess.run(f"grep -q -i -e fatal -e aborted -e \"Galacticus experienced an error in the GSL library\" {path_out_log}",shell=True)
print("...done ("+str(status)+")")
if status.returncode == 0:
print("FAILED: model run (errors):")
subprocess.run(f"cat {path_out_log}",shell=True)
sys.exit()
print("SUCCESS: model run")

# Open the model and extract the recycled fraction.
model = h5py.File(path_out_hdf,'r')
massStellar = model["Outputs/Output1/nodeData/spheroidMassStellar" ][:]
radiusStellar = model["Outputs/Output1/nodeData/spheroidRadius" ]
isIsolated = model["Outputs/Output1/nodeData/nodeIsIsolated" ][:].astype(bool)

massStellarHost, radiusStellarHost = massStellar[isIsolated], radiusStellar[isIsolated]

alpha = 0.56
beta = 1.19E-9

radiusStellarHostPython = (
+beta
*massStellarHost**alpha
)
if np.allclose(np.log10(radiusStellarHost), np.log10(radiusStellarHostPython), rtol=1e-6):
print(f"SUCCESS: results do agree" )
else:
print(f"FAILED: results do not agree")

0 comments on commit bcf66cb

Please sign in to comment.