Skip to content

Commit

Permalink
Merge pull request #532 from galacticusorg/fixZeroStellarHistories
Browse files Browse the repository at this point in the history 
Ensure stellar histories are zeroed when zeroing stellar mass due to ODE solver failure
  • Loading branch information
@abensonca
abensonca authored Jan 18, 2024
2 parents a34f433 + af65c7b commit 38eb56e
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 38eb56e

@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: 38eb56e Previous: a34f433 Ratio
Milky Way model - Likelihood - localGroupStellarMassFunction 190.530964778531 -logℒ 122.779370258702 -logℒ 1.55

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

CC: @abensonca

Please sign in to comment.