Skip to content

Commit

Permalink
Merge pull request #734 from galacticusorg/fixTimeOfCollapseUndefined
Browse files Browse the repository at this point in the history
Catch cases where the time of collapse for a halo is undefined
  • Loading branch information
abensonca authored Nov 5, 2024
2 parents 966f99c + 2507717 commit a06d1a1
Show file tree
Hide file tree
Showing 4 changed files with 188 additions and 12 deletions.
3 changes: 2 additions & 1 deletion .github/workflows/cicd.yml
Original file line number Diff line number Diff line change
Expand Up @@ -1203,7 +1203,8 @@ jobs:
testSuite/regressions/treeWithInitialSatelliteInProgenitorlessHost.xml,
testSuite/regressions/treeWithNoPrimaryProgenitor.xml,
testSuite/regressions/outputRank2ExtendSegFault.xml,
testSuite/regressions/barInstabilityFPE.xml ]
testSuite/regressions/barInstabilityFPE.xml,
testSuite/regressions/haloMassFunctionWDMIndeterminateCollapseTime.xml ]
uses: ./.github/workflows/testModel.yml
with:
file: ${{ matrix.file }}
Expand Down
23 changes: 20 additions & 3 deletions source/dark_matter_halos.mass_accretion_history.Wechsler2002.F90
Original file line number Diff line number Diff line change
Expand Up @@ -247,6 +247,11 @@ double precision function wechsler2002MassAccretionRate(self,node,time)
case (.true.)
! Compute the expansion factor at formation.
mergerTreeFormationExpansionFactor=self%expansionFactorAtFormation(basicBase%mass(),node)
if (mergerTreeFormationExpansionFactor == huge(0.0d0)) then
! No collapse time was found - the accretion rate is indeterminate.
wechsler2002MassAccretionRate=-huge(0.0d0)
return
end if
case (.false.)
! Use the specified formation redshift.
mergerTreeFormationExpansionFactor=self%cosmologyFunctions_%expansionFactorFromRedshift(self%formationRedshift)
Expand Down Expand Up @@ -278,21 +283,33 @@ double precision function wechsler2002ExpansionFactorAtFormation(self,haloMass,n
!!{
Computes the expansion factor at formation using the simple model of \cite{bullock_profiles_2001}.
!!}
use :: Error, only : Error_Report, errorStatusSuccess, errorStatusOutOfRange
implicit none
class (darkMatterHaloMassAccretionHistoryWechsler2002), intent(inout) :: self
type (treeNode ), intent(inout) :: node
double precision , intent(in ) :: haloMass
double precision , parameter :: haloMassFraction =0.015d0 ! Wechsler et al. (2002; Astrophysical Journal, 568:52-70).
double precision :: formationTime , haloMassCharacteristic, &
& sigmaCharacteristic
integer :: status

! Compute the characteristic mass at formation time.
haloMassCharacteristic=haloMassFraction*haloMass
! Compute the corresponding rms fluctuation in the density field (i.e. σ(M)).
sigmaCharacteristic=self%cosmologicalMassVariance_%rootVariance(haloMassCharacteristic,self%cosmologyFunctions_%cosmicTime(1.0d0))
! Get the time at which this equals the critical overdensity for collapse.
formationTime=self%criticalOverdensity_%timeOfCollapse(criticalOverdensity=sigmaCharacteristic,mass=haloMass,node=node)
! Get the corresponding expansion factor.
wechsler2002ExpansionFactorAtFormation=self%cosmologyFunctions_%expansionFactor(formationTime)
formationTime=self%criticalOverdensity_%timeOfCollapse(criticalOverdensity=sigmaCharacteristic,mass=haloMass,node=node,status=status)
select case (status)
case (errorStatusSuccess)
! Get the corresponding expansion factor.
wechsler2002ExpansionFactorAtFormation=self%cosmologyFunctions_%expansionFactor(formationTime)
case (errorStatusOutOfRange)
! Collapse never occurs for this halo mass.
wechsler2002ExpansionFactorAtFormation=huge(0.0d0)
case default
! Unknown error.
wechsler2002ExpansionFactorAtFormation=huge(0.0d0)
call Error_Report('failed to find expansion factor at formation'//{introspection:location})
end select
return
end function wechsler2002ExpansionFactorAtFormation
51 changes: 43 additions & 8 deletions source/structure_formation.cosmological_density_field.F90
Original file line number Diff line number Diff line change
Expand Up @@ -323,26 +323,33 @@ module Cosmological_Density_Field

contains

double precision function criticalOverdensityTimeOfCollapse(self,criticalOverdensity,mass,node)
double precision function criticalOverdensityTimeOfCollapse(self,criticalOverdensity,mass,node,status)
!!{
Returns the time of collapse for a perturbation of linear theory overdensity {\normalfont \ttfamily criticalOverdensity}.
!!}
use :: Cosmology_Functions, only : timeToleranceRelativeBigCrunch
use :: Root_Finder , only : rangeExpandMultiplicative , rangeExpandSignExpectNegative, rangeExpandSignExpectPositive
use :: Error , only : Error_Report , errorStatusSuccess , errorStatusOutOfRange
implicit none
class (criticalOverdensityClass), intent(inout) , target :: self
double precision , intent(in ) :: criticalOverdensity
double precision , intent(in ), optional :: mass
type (treeNode ), intent(inout), optional , target :: node
double precision , parameter :: toleranceRelative =1.0d-12, toleranceAbsolute =0.0d0
integer , parameter :: countPerUnit =10000
integer , intent( out), optional :: status
double precision , parameter :: toleranceRelative =1.0d-12, toleranceAbsolute =0.0d0, &
& fractionTimeCollapseGrowthMinimum=1.0d-03
integer , parameter :: countPerUnit =10000
double precision , allocatable , dimension(:) :: threshold
double precision :: timeBigCrunch , collapseThresholdMinimum , &
& collapseThresholdMaximum , timeGuess
logical :: updateResult , remakeTable
integer :: i , countThresholds , &
& countNewLower , countNewUpper
double precision :: timeBigCrunch , timeGuess , &
& collapseThresholdMinimum , collapseThresholdMaximum , &
& collapseTimePrevious , collapseTimeUpperLimit , &
& timeUpperLimit
logical :: updateResult , remakeTable
integer :: i , countThresholds , &
& countNewLower , countNewUpper

! Assume a successful calculation by default.
if (present(status)) status=errorStatusSuccess
! Determine dependencies.
if (.not.self%dependenciesInitialized) then
self%isMassDependent_=self%isMassDependent() .or. self%cosmologicalMassVariance_%growthIsMassDependent()
Expand Down Expand Up @@ -472,6 +479,34 @@ double precision function criticalOverdensityTimeOfCollapse(self,criticalOverden
end if
self%timeOfCollapsePrevious=self%collapseThreshold%interpolate(criticalOverdensity)
else
! Check that we can find an upper bound to the collapse time. In some cosmologies (e.g. those with a cosmological
! constant) there is a limit to how much the linear growth factor can increase - and therefore some critical
! overdensities will never be reached.
timeUpperLimit =+self%cosmologyFunctions_%cosmicTime (expansionFactor=1.0d0 )
collapseTimeUpperLimit=+self %criticalOverdensityTarget &
& - collapseTimeRoot ( timeUpperLimit)
do while (collapseTimeRoot(timeUpperLimit) < 0.0d0)
timeUpperLimit =+2.0d0 &
& * timeUpperLimit
collapseTimePrevious =+ collapseTimeUpperLimit
collapseTimeUpperLimit=+self%criticalOverdensityTarget &
& - collapseTimeRoot (timeUpperLimit)
! Check if the collapsing density has increased by some minimal amount as we doubled the time of collapse.
if ( &
& collapseTimePrevious-collapseTimeUpperLimit < fractionTimeCollapseGrowthMinimum*collapseTimePrevious &
& .and. &
& collapseTimeRoot(timeUpperLimit) < 0.0d0 &
& ) then
! It did not, suggesting that there is no solution for the collapse time.
if (present(status)) then
criticalOverdensityTimeOfCollapse=-huge(0.0d0)
status =errorStatusOutOfRange
return
else
call Error_Report('unable to bracket collapse time - this can happen in cosmologies with a upper bound to the linear growth (e.g. cosmological constant models)'//{introspection:location})
end if
end if
end do
self%timeOfCollapsePrevious=self%finderTimeOfCollapse%find(rootGuess=self%timeOfCollapsePrevious)
end if
end if
Expand Down
123 changes: 123 additions & 0 deletions testSuite/regressions/haloMassFunctionWDMIndeterminateCollapseTime.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
<?xml version="1.0" encoding="UTF-8"?>
<!-- Model with a low mass WDM particle that previously failed as the collapse time for (strongly-suppressed) low mass halos becomes indeterminate. -->
<parameters>
<formatVersion>2</formatVersion>
<lastModified revision="165490f4968c699f7e5bb73ef60bcb7e626f51e6"/>

<!-- Specify tasks to perform -->
<task value="haloMassFunction">
<haloMassMinimum value="1.0e7"/>
</task>

<!-- Component selection -->
<componentBasic value="standard"/>
<componentDarkMatterProfile value="scale" />

<!-- Use a thermal WDM particle - mass is in keV -->
<darkMatterParticle value="WDMThermal">
<degreesOfFreedomEffective value="1.5" />
<mass value="0.25" />
</darkMatterParticle>

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

<!-- Power spectrum options -->
<!-- Use the Bode et al. (2001) transfer function for thermal WDM -->
<transferFunction value="bode2001">
<epsilon value="0.359" />
<eta value="3.810" />
<nu value="1.100" />
<!-- Bode2001 transfer function works by modifying a CDM transfer function - so feed it a CDM transfer function here -->
<transferFunction value="eisensteinHu1999">
<!-- Feed this transfer function a CDM particle - otherwise it will see the WDM particle defined above and complain that it
can not compute WDM transfer functions -->
<darkMatterParticle value="CDM" />
<neutrinoNumberEffective value="3.046"/>
<neutrinoMassSummed value="0.000"/>
</transferFunction>
</transferFunction>
<powerSpectrumPrimordial value="powerLaw">
<index value="0.961"/>
<wavenumberReference value="1.000"/>
<running value="0.000"/>
</powerSpectrumPrimordial>
<powerSpectrumPrimordialTransferred value="simple"/>
<!-- When computing sigma(M) for power spectra with a cut off it's better to use a filter that is sharp in k-space, instead of
the usual real-space top-hat (which introduces artificial halos below the cut-off scale -->
<cosmologicalMassVariance value="filteredPower">
<monotonicInterpolation value="true" />
<nonMonotonicIsFatal value="false" />
<powerSpectrumWindowFunction value="sharpKSpace">
<normalization value="2.5" />
</powerSpectrumWindowFunction>
<sigma_8 value="0.807" />
<tolerance value="3.0e-4" />
<toleranceTopHat value="3.0e-4" />
</cosmologicalMassVariance>

<!-- Structure formation options -->
<linearGrowth value="collisionlessMatter"/>
<haloMassFunction value="tinker2008"/>
<virialDensityContrast value="sphericalCollapseClsnlssMttrCsmlgclCnstnt"/>
<!-- Use the Barkana et al. (2001) method for the critical overdensity for collapse for WDM -->
<criticalOverdensity value="barkana2001WDM">
<!-- Barkana2001 critical overdensity works by modifying a CDM critical overdensity - so feed it a CDM critical overdensity
here -->
<criticalOverdensity value="sphericalCollapseClsnlssMttrCsmlgclCnstnt">
<!-- Feed this critical overdensity a CDM particle - otherwise it will see the WDM particle defined above and complain that
it can not compute WDM critical overdensities-->
<darkMatterParticle value="CDM" />
</criticalOverdensity>
</criticalOverdensity>

<!-- Dark matter halo structure options -->
<darkMatterProfileDMO value="NFW"/>
<darkMatterProfileConcentration value="schneider2015">
<!-- Define a reference CDM universe - the Schneider algorithm works by finding halos with the same formation epoch in this
reference universe -->
<reference>
<darkMatterParticle value="CDM" />
<darkMatterProfileConcentration value="gao2008"/>
<criticalOverdensity value="sphericalCollapseClsnlssMttrCsmlgclCnstnt"/>
<cosmologicalMassVariance value="filteredPower">
<sigma_8 value="0.807"/>
<monotonicInterpolation value="true" />
<nonMonotonicIsFatal value="false" />
<powerSpectrumPrimordial value="powerLaw">
<index value="0.961"/>
<wavenumberReference value="1.000"/>
<running value="0.000"/>
</powerSpectrumPrimordial>
<powerSpectrumPrimordialTransferred value="simple"/>
<powerSpectrumWindowFunction value="sharpKSpace">
<normalization value="2.5" />
</powerSpectrumWindowFunction>
<transferFunction value="eisensteinHu1999">
<darkMatterParticle value="CDM" />
<neutrinoNumberEffective value="3.046"/>
<neutrinoMassSummed value="0.000"/>
</transferFunction>
</cosmologicalMassVariance>
</reference>
</darkMatterProfileConcentration>
<darkMatterProfileScaleRadius value="concentrationLimiter">
<concentrationMinimum value=" 1.0"/>
<concentrationMaximum value="100.0"/>
<darkMatterProfileScaleRadius value="concentration"/>
</darkMatterProfileScaleRadius>

<!-- Output options -->
<outputFileName value="testSuite/outputs/regressions/haloMassFunctionWDMIndeterminateCollapseTime.hdf5"/>
<outputTimes value="list">
<redshifts value="0.0 1.0"/>
</outputTimes>

</parameters>

0 comments on commit a06d1a1

Please sign in to comment.