Skip to content

Commit

Permalink
Merge pull request #744 from galacticusorg/fixTreeWeights
Browse files Browse the repository at this point in the history 
Various fixes related to reading tree masses from file
  • Loading branch information
@abensonca
abensonca authored Nov 27, 2024
2 parents ddb74ce + 600e194 commit f482d6b
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 30 deletions.
52 changes: 35 additions & 17 deletions source/merger_trees.construct.build.masses.read.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@
contains
!![
<methods>
<method description="Read the halo masses, and, optionally, weights, from file.." method="read" />
<method description="Read the halo masses, and, optionally, weights, from file." method="read"/>
</methods>
!!]
procedure :: construct => readConstruct
Expand All @@ -58,13 +58,16 @@ subroutine readConstruct(self,time,mass,massMinimum,massMaximum,weight)
!!{
Construct a set of merger tree masses by reading from a file.
!!}
use :: Sorting , only : sort
use :: Sorting, only : sort
implicit none
class (mergerTreeBuildMassesRead), intent(inout) :: self
double precision , intent(in ) :: time
double precision , intent( out), allocatable, dimension(:) :: mass , weight , &
& massMinimum , massMaximum
integer :: i
double precision , intent( out), allocatable, dimension(:) :: mass , weight , &
& massMinimum , massMaximum
double precision :: massPrevious
integer :: iStart , iEnd , &
& iStartPrevious, iEndPrevious, &
& i
!$GLC attributes unused :: time

call self%read(mass,weight)
Expand All @@ -74,19 +77,34 @@ subroutine readConstruct(self,time,mass,massMinimum,massMaximum,weight)
call sort(mass )
allocate(massMinimum,mold=mass)
allocate(massMaximum,mold=mass)
do i=1,size(mass)
massMinimum(i)=+mass(i)/sqrt(1.0d0+self%massIntervalFractional)
massMaximum(i)=+mass(i)*sqrt(1.0d0+self%massIntervalFractional)
if (i > 1 .and. massMinimum(i) < massMaximum(i-1)) then
massMinimum(i )=sqrt( &
& +mass(i-1) &
& *mass(i ) &
& )
massMaximum(i-1)=sqrt( &
& +mass(i-1) &
& *mass(i ) &
& )
massPrevious =-huge(0.0d0)
iStart =-huge(0 )
iEnd =-huge(0 )
iStartPrevious=-huge(0 )
iEndPrevious =-huge(0 )
i =+ 0
do while (i <= size(mass))
i=i+1
if (i == size(mass)+1 .or. mass(i) /= massPrevious) then
! Process a block of trees of the same mass.
if (massPrevious > 0.0d0) then
massMinimum(iStart:iEnd)=+mass(iStart)/sqrt(1.0d0+self%massIntervalFractional)
massMaximum(iStart:iEnd)=+mass(iStart)*sqrt(1.0d0+self%massIntervalFractional)
if (iStart > 1 .and. massMinimum(iStart) < massMaximum(iStartPrevious)) then
massMinimum(iStart :iEnd )=+sqrt( &
& +mass (iStartPrevious) &
& *mass (iStart ) &
& )
massMaximum(iStartPrevious:iEndPrevious)=+ massMinimum(iStart )
end if
end if
! Update to the start of the next block.
iStartPrevious= iStart
iEndPrevious = iEnd
iStart = i
massPrevious =mass(i )
end if
iEnd=i
end do
end if
return
Expand Down
9 changes: 6 additions & 3 deletions source/merger_trees.construct.build.masses.read.XML.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,8 +120,8 @@ subroutine readXMLRead(self,mass,weight)
use :: File_Utilities, only : File_Name_Expand
implicit none
class (mergerTreeBuildMassesReadXML), intent(inout) :: self
double precision , intent( out), allocatable, dimension(:) :: mass, weight
type (node ), pointer :: doc , rootNode
double precision , intent( out), allocatable, dimension(:) :: mass , weight
type (node ), pointer :: doc , rootNode
integer :: ioErr

!$omp critical (FoX_DOM_Access)
Expand All @@ -131,7 +131,10 @@ subroutine readXMLRead(self,mass,weight)
! Read all tree masses.
call XML_Array_Read(doc,"treeRootMass",mass)
! Extract tree weights if available.
if (XML_Path_Exists(rootNode,"treeWeight")) call XML_Array_Read_Static(doc,"treeWeight",weight)
if (XML_Path_Exists(rootNode,"treeWeight")) then
call XML_Array_Read(doc,"treeWeight",weight)
if (size(weight) /= size(mass)) call Error_Report('`mass` and `weight` arrays must have the same size'//{introspection:location})
end if
! Finished - destroy the XML document.
call destroy(doc)
!$omp end critical (FoX_DOM_Access)
Expand Down
35 changes: 25 additions & 10 deletions source/merger_trees.construct.build.masses.replicate.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,8 +105,8 @@ subroutine replicateConstruct(self,time,mass,massMinimum,massMaximum,weight)
!!{
Construct a set of merger tree masses by sampling from a distribution.
!!}
use :: Error , only : Error_Report
use, intrinsic :: ISO_C_Binding , only : c_size_t
use :: Error , only : Error_Report
use, intrinsic :: ISO_C_Binding, only : c_size_t
implicit none
class (mergerTreeBuildMassesReplicate), intent(inout) :: self
double precision , intent(in ) :: time
Expand All @@ -117,17 +117,32 @@ subroutine replicateConstruct(self,time,mass,massMinimum,massMaximum,weight)
integer (c_size_t ) :: i

call self%mergerTreeBuildMasses_%construct(time,massTmp,massMinimumTmp,massMaximumTmp,weightTmp)
if (allocated(massTmp).and.allocated(weightTmp)) then
allocate(mass (size(massTmp )*self%replicationCount))
allocate(weight(size(weightTmp)*self%replicationCount))
if (allocated(massTmp)) then
allocate (mass (size(massTmp )*self%replicationCount))
if (allocated(massMinimumTmp)) &
& allocate(massMinimum(size(massMinimumTmp)*self%replicationCount))
if (allocated(massMaximumTmp)) &
& allocate(massMaximum(size(massMaximumTmp)*self%replicationCount))
if (allocated(weightTmp)) &
& allocate(weight (size(weightTmp )*self%replicationCount))
do i=1,size(massTmp)
mass ((i-1)*self%replicationCount+1:i*self%replicationCount)=massTmp (i)
weight((i-1)*self%replicationCount+1:i*self%replicationCount)=weightTmp(i)/dble(self%replicationCount)
mass ((i-1)*self%replicationCount+1:i*self%replicationCount)=massTmp (i)
if (allocated(massMinimumTmp)) &
& massMinimum((i-1)*self%replicationCount+1:i*self%replicationCount)=massMinimumTmp(i)
if (allocated(massMaximumTmp)) &
& massMaximum((i-1)*self%replicationCount+1:i*self%replicationCount)=massMaximumTmp(i)
if (allocated(weightTmp)) &
& weight ((i-1)*self%replicationCount+1:i*self%replicationCount)=weightTmp (i)/dble(self%replicationCount)
end do
deallocate(massTmp )
deallocate(weightTmp)
deallocate (massTmp )
if (allocated(massMinimumTmp)) &
& deallocate(massMinimumTmp)
if (allocated(massMaximumTmp)) &
& deallocate(massMaximumTmp)
if (allocated(weightTmp )) &
& deallocate(weightTmp )
else
call Error_Report('masses and weights are required'//{introspection:location})
call Error_Report('masses are required'//{introspection:location})
end if
return
end subroutine replicateConstruct

0 comments on commit f482d6b

Please sign in to comment.