Skip to content

Commit

Permalink
fix: Ensure stellar histories are zeroed when zeroing stellar mass du…
Browse files Browse the repository at this point in the history 
…e to ODE solver failure

Without this, the summed star formation history will not match the stellar mass.
  • Loading branch information
@abensonca
abensonca committed Jan 18, 2024
1 parent a34f433 commit af65c7b
Show file tree
Hide file tree
Showing 3 changed files with 58 additions and 12 deletions.
19 changes: 17 additions & 2 deletions source/objects.nodes.components.disk.standard.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -490,6 +490,7 @@ subroutine Node_Component_Disk_Standard_Post_Step(node,status)
use :: Error , only : Error_Report
use :: Galacticus_Nodes , only : defaultDiskComponent, nodeComponentDisk , nodeComponentDiskStandard, nodeComponentSpin, &
& treeNode , nodeComponentBasic
use :: Histories , only : history
use :: Interface_GSL , only : GSL_Success , GSL_Continue
use :: ISO_Varying_String , only : assignment(=) , operator(//) , varying_string
use :: Stellar_Luminosities_Structure, only : abs , stellarLuminosities
Expand All @@ -509,7 +510,9 @@ subroutine Node_Component_Disk_Standard_Post_Step(node,status)
!$omp threadprivate(message)
type (stellarLuminosities), save :: luminositiesStellar
!$omp threadprivate(luminositiesStellar)

type (history ), save :: historyStellar
!$omp threadprivate(historyStellar)

! Return immediately if this class is not in use.
if (.not.defaultDiskComponent%standardIsActive()) return
! Get the disk component.
Expand Down Expand Up @@ -558,6 +561,12 @@ subroutine Node_Component_Disk_Standard_Post_Step(node,status)
specificAngularMomentum=0.0d0
call disk% massStellarSet( 0.0d0)
call disk% abundancesStellarSet( zeroAbundances)
historyStellar=disk%starFormationHistory ()
call historyStellar%reset()
call disk %starFormationHistorySet (historyStellar)
historyStellar=disk%stellarPropertiesHistory()
call historyStellar%reset()
call disk %stellarPropertiesHistorySet(historyStellar)
! We need to reset the stellar luminosities to zero. We can't simply use the "zeroStellarLuminosities" instance since
! our luminosities may have been truncated. If we were to use "zeroStellarLuminosities" then the number of stellar
! luminosities associated with the disk would change - but we are in the middle of differential evolution here and we
Expand Down Expand Up @@ -620,10 +629,16 @@ subroutine Node_Component_Disk_Standard_Post_Step(node,status)
specificAngularMomentum=disk%angularMomentum()/massDisk
if (specificAngularMomentum < 0.0d0) specificAngularMomentum=disk%radius()*disk%velocity()
end if
! Reset the stellar, abundances and angular momentum of the disk.
! Reset the stellar mass, abundances and angular momentum of the disk.
call disk% massStellarSet( 0.0d0)
call disk%abundancesStellarSet( zeroAbundances)
call disk% angularMomentumSet(specificAngularMomentum*disk%massGas())
historyStellar=disk%starFormationHistory ()
call historyStellar%reset()
call disk %starFormationHistorySet (historyStellar)
historyStellar=disk%stellarPropertiesHistory()
call historyStellar%reset()
call disk %stellarPropertiesHistorySet(historyStellar)
! Indicate that ODE evolution should continue after this state change.
if (status == GSL_Success) status=GSL_Continue
end if
Expand Down
34 changes: 25 additions & 9 deletions source/objects.nodes.components.disk.very_simple.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -763,17 +763,19 @@ subroutine satelliteMerger(self,node)
use :: Abundances_Structure , only : zeroAbundances
use :: Error , only : Error_Report
use :: Galacticus_Nodes , only : nodeComponentDisk , nodeComponentDiskVerySimple, nodeComponentSpheroid , treeNode
use :: Histories , only : history
use :: Satellite_Merging_Mass_Movements, only : destinationMergerDisk , destinationMergerSpheroid , enumerationDestinationMergerType
use :: Stellar_Luminosities_Structure , only : zeroStellarLuminosities
implicit none
class (* ), intent(inout) :: self
type (treeNode ), intent(inout) :: node
type (treeNode ), pointer :: nodeHost
class (nodeComponentDisk ), pointer :: diskHost , disk
class (nodeComponentSpheroid ), pointer :: spheroidHost
type (enumerationDestinationMergerType) :: destinationGasSatellite, destinationGasHost , &
& destinationStarsHost , destinationStarsSatellite
logical :: mergerIsMajor
class (* ), intent(inout) :: self
type (treeNode ), intent(inout) :: node
type (treeNode ), pointer :: nodeHost
class (nodeComponentDisk ), pointer :: diskHost , disk
class (nodeComponentSpheroid ), pointer :: spheroidHost
type (enumerationDestinationMergerType) :: destinationGasSatellite, destinationGasHost , &
& destinationStarsHost , destinationStarsSatellite
type (history ) :: historyHost , historyNode
logical :: mergerIsMajor
!$GLC attributes unused :: self

! Check that the disk is of the verySimple class.
Expand Down Expand Up @@ -838,7 +840,14 @@ subroutine satelliteMerger(self,node)
& diskHost %luminositiesStellar() &
& +disk %luminositiesStellar() &
& )
case (destinationMergerSpheroid%ID)
! Also add stellar properties histories.
historyNode=disk %stellarPropertiesHistory()
historyHost=diskHost%stellarPropertiesHistory()
call historyHost%interpolatedIncrement (historyNode)
call historyNode%reset ( )
call diskHost %stellarPropertiesHistorySet(historyHost)
call disk %stellarPropertiesHistorySet(historyNode)
case (destinationMergerSpheroid%ID)
call spheroidHost%massStellarSet ( &
& spheroidHost% massStellar () &
& +disk % massStellar () &
Expand All @@ -851,6 +860,13 @@ subroutine satelliteMerger(self,node)
& spheroidHost%luminositiesStellar() &
& +disk %luminositiesStellar() &
& )
! Also add stellar properties histories.
historyNode=disk %stellarPropertiesHistory()
historyHost=spheroidHost%stellarPropertiesHistory()
call historyHost %interpolatedIncrement (historyNode)
call historyNode %reset ( )
call spheroidHost%stellarPropertiesHistorySet(historyHost)
call disk %stellarPropertiesHistorySet(historyNode)
case default
call Error_Report( &
& 'unrecognized movesTo descriptor'// &
Expand Down
17 changes: 16 additions & 1 deletion source/objects.nodes.components.spheroid.standard.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -435,6 +435,7 @@ subroutine Node_Component_Spheroid_Standard_Post_Step(node,status)
use :: Abundances_Structure , only : abs , zeroAbundances
use :: Display , only : displayMessage , verbosityLevelWarn
use :: Galacticus_Nodes , only : defaultSpheroidComponent, nodeComponentSpheroid, nodeComponentSpheroidStandard, treeNode
use :: Histories , only : history
use :: Interface_GSL , only : GSL_Success , GSL_Continue
use :: ISO_Varying_String , only : assignment(=) , operator(//) , varying_string
use :: Stellar_Luminosities_Structure, only : abs , stellarLuminosities
Expand All @@ -453,6 +454,8 @@ subroutine Node_Component_Spheroid_Standard_Post_Step(node,status)
!$omp threadprivate(message)
type (stellarLuminosities ), save :: luminositiesStellar
!$omp threadprivate(luminositiesStellar)
type (history ), save :: historyStellar
!$omp threadprivate(historyStellar)

! Return immediately if this class is not in use.
if (.not.defaultSpheroidComponent%standardIsActive()) return
Expand Down Expand Up @@ -502,6 +505,12 @@ subroutine Node_Component_Spheroid_Standard_Post_Step(node,status)
specificAngularMomentum=0.0d0
call spheroid% massStellarSet( 0.0d0)
call spheroid% abundancesStellarSet( zeroAbundances)
historyStellar=spheroid%starFormationHistory ()
call historyStellar%reset()
call spheroid %starFormationHistorySet (historyStellar)
historyStellar=spheroid%stellarPropertiesHistory()
call historyStellar%reset()
call spheroid %stellarPropertiesHistorySet(historyStellar)
! We need to reset the stellar luminosities to zero. We can't simply use the "zeroStellarLuminosities" instance since
! our luminosities may have been truncated. If we were to use "zeroStellarLuminosities" then the number of stellar
! luminosities associated with the spheroid would change - but we are in the middle of differential evolution here and we
Expand Down Expand Up @@ -562,10 +571,16 @@ subroutine Node_Component_Spheroid_Standard_Post_Step(node,status)
else
specificAngularMomentum=spheroid%angularMomentum()/massSpheroid
end if
! Reset the stellar, abundances and angular momentum of the spheroid.
! Reset the stellar mass, abundances and angular momentum of the spheroid.
call spheroid% massStellarSet( 0.0d0)
call spheroid% abundancesStellarSet( zeroAbundances)
call spheroid%angularMomentumSet(specificAngularMomentum*spheroid%massGas())
historyStellar=spheroid%starFormationHistory ()
call historyStellar%reset()
call spheroid %starFormationHistorySet (historyStellar)
historyStellar=spheroid%stellarPropertiesHistory()
call historyStellar%reset()
call spheroid %stellarPropertiesHistorySet(historyStellar)
! Indicate that ODE evolution should continue after this state change.
if (status == GSL_Success) status=GSL_Continue
end if
Expand Down

1 comment on commit af65c7b

@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: af65c7b Previous: a34f433 Ratio
Milky Way model - Likelihood - localGroupStellarMassFunction 189.423032673505 -logℒ 122.779370258702 -logℒ 1.54

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

CC: @abensonca

Please sign in to comment.