Skip to content

Commit

Permalink
Merge pull request #482 from galacticusorg/bryanNorman1998GeneralCosm…
Browse files Browse the repository at this point in the history 
…ology

Allow the `virialDensityContrastBryanNorman1998` class to work for arbitrary cosmologies
  • Loading branch information
@abensonca
abensonca committed Sep 27, 2023
2 parents ceb3791 + 71469be commit 456ee65
Show file tree
Hide file tree
Showing 5 changed files with 116 additions and 61 deletions.
1 change: 1 addition & 0 deletions source/dark_matter_profiles.structure.concentration.MunozCuartas2011.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,7 @@ function munozCuartas2011ConstructorInternal(cosmologyParameters_,cosmologyFunct
<referenceConstruct owner="self" object="virialDensityContrastDefinition_">
<constructor>
virialDensityContrastBryanNorman1998 ( &amp;
&amp; allowUnsupportedCosmology = .true. , &amp;
&amp; cosmologyParameters_ =self%cosmologyParameters_ , &amp;
&amp; cosmologyFunctions_ =self%cosmologyFunctions_ &amp;
&amp; )
Expand Down
1 change: 1 addition & 0 deletions source/dark_matter_profiles.structure.concentration.Zhao2009.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,7 @@ function zhao2009ConstructorInternal(cosmologyFunctions_,cosmologyParameters_,da
<referenceConstruct owner="self" object="virialDensityContrastDefinition_">
<constructor>
virialDensityContrastBryanNorman1998 ( &amp;
&amp; allowUnsupportedCosmology = .true. , &amp;
&amp; cosmologyParameters_ =self%cosmologyParameters_ , &amp;
&amp; cosmologyFunctions_ =self%cosmologyFunctions_ &amp;
&amp; )
Expand Down
4 changes: 2 additions & 2 deletions source/output.analyses.HI_vs_halo_mass_relation.ALFALFA_Padmanabhan_2017.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -346,11 +346,11 @@ function hiVsHaloMassRelationPadmanabhan2017ConstructorInternal(likelihoodBin,sy
! lower by a factor Ωₘ^⅙ than they should be. We explicitly account for this factor when computing virial velocity below.
allocate(virialDensityContrastData )
!![
<referenceConstruct object="virialDensityContrastData" constructor="virialDensityContrastBryanNorman1998 (cosmologyParametersData ,cosmologyFunctionsData )"/>
<referenceConstruct object="virialDensityContrastData" constructor="virialDensityContrastBryanNorman1998 (.false.,cosmologyParametersData ,cosmologyFunctionsData )"/>
!!]
allocate(darkMatterHaloScaleData)
!![
<referenceConstruct object="darkMatterHaloScaleData" constructor="darkMatterHaloScaleVirialDensityContrastDefinition (cosmologyParametersData ,cosmologyFunctionsData ,virialDensityContrastData)"/>
<referenceConstruct object="darkMatterHaloScaleData" constructor="darkMatterHaloScaleVirialDensityContrastDefinition ( cosmologyParametersData ,cosmologyFunctionsData ,virialDensityContrastData)"/>
!!]
! Generate the target dataset.
allocate(massHILogarithmicTarget (massHaloCount ))
Expand Down
10 changes: 9 additions & 1 deletion source/satellites.merging.virial_orbits.Li2020.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,15 @@ function li2020ConstructorInternal(mu1,mu2,sigma1,a0,a1,a2,a3,b1,b2,c,propagateO
! Create virial density contrast definition.
allocate(self%virialDensityContrastDefinition_)
!![
<referenceConstruct isResult="yes" owner="self" object="virialDensityContrastDefinition_" constructor="virialDensityContrastBryanNorman1998(self%cosmologyParameters_,self%cosmologyFunctions_)"/>
<referenceConstruct isResult="yes" owner="self" object="virialDensityContrastDefinition_">
<constructor>
virialDensityContrastBryanNorman1998( &amp;
&amp; allowUnsupportedCosmology= .true. , &amp;
&amp; cosmologyParameters_ =self%cosmologyParameters_, &amp;
&amp; cosmologyFunctions_ =self%cosmologyFunctions_ &amp;
&amp; )
</constructor>
</referenceConstruct>
!!]
return
end function li2020ConstructorInternal
Expand Down
161 changes: 103 additions & 58 deletions source/structure_formation.virial_density_contrast.Bryan_Norman.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,20 +37,29 @@
!![
<virialDensityContrast name="virialDensityContrastBryanNorman1998">
<description>
A dark matter halo virial density contrast class using the fitting functions given by \cite{bryan_statistical_1998}. As
such, it is valid only for $\Omega_\Lambda=0$ or $\Omega_\mathrm{M}+\Omega_\Lambda=1$ cosmologies and will abort on other
cosmologies.
A dark matter halo virial density contrast class using the fitting functions given by \cite{bryan_statistical_1998}. As such,
it is valid only for $\Omega_\Lambda=0$ or $\Omega_\mathrm{M}+\Omega_\Lambda=1$ cosmologies and will either abort on other
cosmologies (if {\normalfont \ttfamily [allowUnsupportedCosmology]=false}), or revert to a numerical solution from top-hat
collapse (if {\normalfont \ttfamily [allowUnsupportedCosmology]=true}).
</description>
<deepCopy>
<functionClass variables="virialDensityContrastTopHat_"/>
</deepCopy>
<stateStorable>
<functionClass variables="virialDensityContrastTopHat_"/>
</stateStorable>
</virialDensityContrast>
!!]
type, extends(virialDensityContrastClass) :: virialDensityContrastBryanNorman1998
!!{
A dark matter halo virial density contrast class using the fitting functions of \cite{bryan_statistical_1998}.
!!}
private
class(cosmologyParametersClass ), pointer :: cosmologyParameters_ => null()
class(cosmologyFunctionsClass ), pointer :: cosmologyFunctions_ => null()
type (enumerationBryanNorman1998FitType) :: fitType
class (cosmologyParametersClass ), pointer :: cosmologyParameters_ => null()
class (cosmologyFunctionsClass ), pointer :: cosmologyFunctions_ => null()
type (virialDensityContrastSphericalCollapseClsnlssMttrCsmlgclCnstnt), pointer :: virialDensityContrastTopHat_ => null()
type (enumerationBryanNorman1998FitType ) :: fitType
logical :: useFittingFunction , allowUnsupportedCosmology
contains
final :: bryanNorman1998Destructor
procedure :: densityContrast => bryanNorman1998DensityContrast
Expand All @@ -74,16 +83,23 @@ function bryanNorman1998ConstructorParameters(parameters) result(self)
!!}
use :: Input_Parameters, only : inputParameter, inputParameters
implicit none
type (virialDensityContrastBryanNorman1998) :: self
type (inputParameters ), intent(inout) :: parameters
class(cosmologyParametersClass ), pointer :: cosmologyParameters_
class(cosmologyFunctionsClass ), pointer :: cosmologyFunctions_
type (virialDensityContrastBryanNorman1998) :: self
type (inputParameters ), intent(inout) :: parameters
class (cosmologyParametersClass ), pointer :: cosmologyParameters_
class (cosmologyFunctionsClass ), pointer :: cosmologyFunctions_
logical :: allowUnsupportedCosmology

!![
<inputParameter>
<name>allowUnsupportedCosmology</name>
<defaultValue>.false.</defaultValue>
<source>parameters</source>
<description>If true, unsupported cosmologies revert to using a numerical solution from the top-hat collapse model. Otherwise, unsupported cosmologies result in an error.</description>
</inputParameter>
<objectBuilder class="cosmologyParameters" name="cosmologyParameters_" source="parameters"/>
<objectBuilder class="cosmologyFunctions" name="cosmologyFunctions_" source="parameters"/>
!!]
self=virialDensityContrastBryanNorman1998(cosmologyParameters_,cosmologyFunctions_)
self=virialDensityContrastBryanNorman1998(allowUnsupportedCosmology,cosmologyParameters_,cosmologyFunctions_)
!![
<inputParametersValidate source="parameters"/>
<objectDestructor name="cosmologyParameters_"/>
Expand All @@ -92,25 +108,40 @@ function bryanNorman1998ConstructorParameters(parameters) result(self)
return
end function bryanNorman1998ConstructorParameters

function bryanNorman1998ConstructorInternal(cosmologyParameters_,cosmologyFunctions_) result(self)
function bryanNorman1998ConstructorInternal(allowUnsupportedCosmology,cosmologyParameters_,cosmologyFunctions_) result(self)
!!{
Internal constructor for the {\normalfont \ttfamily bryanNorman1998} dark matter halo virial density contrast class.
!!}
use :: Error , only : Error_Report
use :: Numerical_Comparison, only : Values_Differ
implicit none
type (virialDensityContrastBryanNorman1998) :: self
class(cosmologyParametersClass ), intent(in ), target :: cosmologyParameters_
class(cosmologyFunctionsClass ), intent(in ), target :: cosmologyFunctions_
type (virialDensityContrastBryanNorman1998) :: self
class (cosmologyParametersClass ), intent(in ), target :: cosmologyParameters_
class (cosmologyFunctionsClass ), intent(in ), target :: cosmologyFunctions_
logical , intent(in ) :: allowUnsupportedCosmology
!![
<constructorAssign variables="*cosmologyParameters_, *cosmologyFunctions_"/>
<constructorAssign variables="allowUnsupportedCosmology, *cosmologyParameters_, *cosmologyFunctions_"/>
!!]

! Check that fitting formulae are applicable to this cosmology.
self%useFittingFunction=.true.
if (self%cosmologyParameters_%OmegaDarkEnergy() == 0.0d0) then
self%fitType=bryanNorman1998FitZeroLambda
else if (.not.Values_Differ(self%cosmologyParameters_%OmegaMatter()+self%cosmologyParameters_%OmegaDarkEnergy(),1.0d0,absTol=1.0d-6)) then
self%fitType=bryanNorman1998FitFlatUniverse
else if (self%allowUnsupportedCosmology) then
self%useFittingFunction=.false.
allocate(self%virialDensityContrastTopHat_)
!![
<referenceConstruct isResult="yes" owner="self" object="virialDensityContrastTopHat_">
<constructor>
virialDensityContrastSphericalCollapseClsnlssMttrCsmlgclCnstnt( &amp;
&amp; tableStore =.true. , &amp;
&amp; cosmologyFunctions_=self%cosmologyFunctions_ &amp;
&amp; )
</constructor>
</referenceConstruct>
!!]
else
call Error_Report('no fitting formula available for this cosmology'//{introspection:location})
end if
Expand All @@ -128,6 +159,11 @@ subroutine bryanNorman1998Destructor(self)
<objectDestructor name="self%cosmologyParameters_" />
<objectDestructor name="self%cosmologyFunctions_" />
!!]
if (.not.self%useFittingFunction) then
!![
<objectDestructor name="self%virialDensityContrastTopHat_" />
!!]
end if
return
end subroutine bryanNorman1998Destructor

Expand All @@ -143,18 +179,21 @@ Return the virial density contrast at the given epoch, assuming the fitting func
double precision , intent(in ), optional :: time , expansionFactor
logical , intent(in ), optional :: collapsing
double precision :: x
!$GLC attributes unused :: mass

x=self%cosmologyFunctions_%omegaMatterEpochal(time,expansionFactor,collapsing)-1.0d0
select case (self%fitType%ID)
case (bryanNorman1998FitZeroLambda %ID)
bryanNorman1998DensityContrast=(18.0d0*Pi**2+60.0d0*x-32.0d0*x**2)/(x+1.0d0)
case (bryanNorman1998FitFlatUniverse%ID)
bryanNorman1998DensityContrast=(18.0d0*Pi**2+82.0d0*x-39.0d0*x**2)/(x+1.0d0)
case default
bryanNorman1998DensityContrast=0.0d0
call Error_Report('invalid fit type'//{introspection:location})
end select

if (self%useFittingFunction) then
x=self%cosmologyFunctions_%omegaMatterEpochal(time,expansionFactor,collapsing)-1.0d0
select case (self%fitType%ID)
case (bryanNorman1998FitZeroLambda %ID)
bryanNorman1998DensityContrast=(18.0d0*Pi**2+60.0d0*x-32.0d0*x**2)/(x+1.0d0)
case (bryanNorman1998FitFlatUniverse%ID)
bryanNorman1998DensityContrast=(18.0d0*Pi**2+82.0d0*x-39.0d0*x**2)/(x+1.0d0)
case default
bryanNorman1998DensityContrast=0.0d0
call Error_Report('invalid fit type'//{introspection:location})
end select
else
bryanNorman1998DensityContrast=self%virialDensityContrastTopHat_%densityContrast(mass,time,expansionFactor,collapsing)
end if
return
end function bryanNorman1998DensityContrast

Expand All @@ -170,32 +209,35 @@ Return the virial density contrast at the given epoch, assuming the fitting func
double precision , intent(in ), optional :: time , expansionFactor
logical , intent(in ), optional :: collapsing
double precision :: x
!$GLC attributes unused :: mass

x=self%cosmologyFunctions_%omegaMatterEpochal(time,expansionFactor,collapsing)-1.0d0
select case (self%fitType%ID)
case (bryanNorman1998FitZeroLambda %ID)
bryanNorman1998DensityContrastRateOfChange= &
& ( &
& +( +60.0d0 -64.0d0*x ) &
& -(18.0d0*Pi**2+60.0d0*x-32.0d0*x**2) &
& /(x+1.0d0) &
& ) &
& *self%cosmologyFunctions_%omegaMatterRateOfChange(time,expansionFactor,collapsing) &
& / (x+1.0d0)
case (bryanNorman1998FitFlatUniverse%ID)
bryanNorman1998DensityContrastRateOfChange= &
& ( &
& +( +82.0d0 -78.0d0*x ) &
& -(18.0d0*Pi**2+82.0d0*x-39.0d0*x**2) &
& /(x+1.0d0) &
& ) &
& *self%cosmologyFunctions_%omegaMatterRateOfChange(time,expansionFactor,collapsing) &
& / (x+1.0d0)
case default
bryanNorman1998DensityContrastRateOfChange=0.0d0
call Error_Report('invalid fit type'//{introspection:location})
end select

if (self%useFittingFunction) then
x=self%cosmologyFunctions_%omegaMatterEpochal(time,expansionFactor,collapsing)-1.0d0
select case (self%fitType%ID)
case (bryanNorman1998FitZeroLambda %ID)
bryanNorman1998DensityContrastRateOfChange= &
& ( &
& +( +60.0d0 -64.0d0*x ) &
& -(18.0d0*Pi**2+60.0d0*x-32.0d0*x**2) &
& /(x+1.0d0) &
& ) &
& *self%cosmologyFunctions_%omegaMatterRateOfChange(time,expansionFactor,collapsing) &
& / (x+1.0d0)
case (bryanNorman1998FitFlatUniverse%ID)
bryanNorman1998DensityContrastRateOfChange= &
& ( &
& +( +82.0d0 -78.0d0*x ) &
& -(18.0d0*Pi**2+82.0d0*x-39.0d0*x**2) &
& /(x+1.0d0) &
& ) &
& *self%cosmologyFunctions_%omegaMatterRateOfChange(time,expansionFactor,collapsing) &
& / (x+1.0d0)
case default
bryanNorman1998DensityContrastRateOfChange=0.0d0
call Error_Report('invalid fit type'//{introspection:location})
end select
else
bryanNorman1998DensityContrastRateOfChange=self%virialDensityContrastTopHat_%densityContrastRateOfChange(mass,time,expansionFactor,collapsing)
end if
return
end function bryanNorman1998DensityContrastRateOfChange

Expand All @@ -209,10 +251,13 @@ double precision function bryanNorman1998TurnAroundOverVirialRadii(self,mass,tim
double precision , intent(in ) :: mass
double precision , intent(in ), optional :: time , expansionFactor
logical , intent(in ), optional :: collapsing
!$GLC attributes unused :: self, mass, time, expansionFactor, collapsing

! In simple cosmological constant dark energy universes, this ratio is always precisely 2 (e.g. Percival 2005;
! http://adsabs.harvard.edu/abs/2005A%26A...443..819P)
bryanNorman1998TurnAroundOverVirialRadii=2.0d0
if (self%useFittingFunction) then
! In simple cosmological constant dark energy universes, this ratio is always precisely 2 (e.g. Percival 2005;
! http://adsabs.harvard.edu/abs/2005A%26A...443..819P)
bryanNorman1998TurnAroundOverVirialRadii=2.0d0
else
bryanNorman1998TurnAroundOverVirialRadii=self%turnAroundOverVirialRadii(mass,time,expansionFactor,collapsing)
end if
return
end function bryanNorman1998TurnAroundOverVirialRadii

1 comment on commit 456ee65

@github-actions
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Performance Alert ⚠️

Possible performance regression was detected for benchmark 'Milky Way model benchmarks'.
Benchmark result of this commit is worse than the previous benchmark result exceeding threshold 1.10.

Benchmark suite Current: 456ee65 Previous: ceb3791 Ratio
Milky Way model - Likelihood - localGroupMassMetallicityRelation 3.26508433888745 -logℒ 2.9594186480383 -logℒ 1.10
Milky Way model - Likelihood - localGroupMassSizeRelation 10.3101047357318 -logℒ 8.69418556042769 -logℒ 1.19

This comment was automatically generated by workflow using github-action-benchmark.

CC: @abensonca

Please sign in to comment.