Skip to content

Commit

Permalink
Merge pull request #620 from galacticusorg/featSFHAdaptiveOpt
Browse files Browse the repository at this point in the history 
Store `starFormationHistoryAdaptive` timesteps and mappings to file
  • Loading branch information
@abensonca
abensonca authored Jun 6, 2024
2 parents b0a9e5c + aecf798 commit 3702136
Showing 1 changed file with 126 additions and 88 deletions.
214 changes: 126 additions & 88 deletions source/star_formation.histories.adaptive.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,11 @@ function adaptiveConstructorInternal(outputTimes_,countOutputBuffer,timeStepMini
use :: Error , only : Error_Report
use :: Galactic_Structure_Options, only : componentTypeMax, componentTypeMin
use :: Numerical_Ranges , only : Make_Range , rangeTypeLogarithmic
use :: HDF5_Access , only : hdf5Access
use :: IO_HDF5 , only : hdf5Object
use :: Input_Paths , only : inputPath , pathTypeDataDynamic
use :: File_Utilities , only : File_Exists , File_Lock , File_Unlock, lockDescriptor, &
& Directory_Make
implicit none
type (starFormationHistoryAdaptive) :: self
double precision , intent(in ), dimension(:), optional :: metallicityBoundaries
Expand All @@ -215,6 +220,10 @@ function adaptiveConstructorInternal(outputTimes_,countOutputBuffer,timeStepMini
double precision :: timeStart , timeEnd , &
& metric , metricChangeMinimum, &
& metricChange , metricMinimumGlobal
type (varying_string ) :: fileName
type (hdf5Object ) :: file
type (lockDescriptor ) :: fileLock
character (len=16 ) :: name
!![
<constructorAssign variables="countOutputBuffer, timeStepMinimum, countTimeStepsMaximum, metallicityMinimum, metallicityMaximum, countMetallicities, *outputTimes_"/>
!!]
Expand Down Expand Up @@ -255,100 +264,129 @@ function adaptiveConstructorInternal(outputTimes_,countOutputBuffer,timeStepMini
end select
end if
! Construct the time bins and rebinning strategy to be used for each output.
fileName=inputPath(pathTypeDataDynamic)//"starFormation/"//self%objectType()//"_"//self%hashedDescriptor()//".hdf5"
allocate(self%intervals(self%outputTimes_%count()))
do iOutput=1,self%outputTimes_%count()
! Our start time is either the minimum timestep size (for the first output), or the end time of the timesteps of the
! previous output, plus the minimum timestep.
if (iOutput == 1) then
timeStart=self%timeStepMinimum
else
timeStart=self%timeStepMinimum+self%intervals(iOutput-1)%time(size(self%intervals(iOutput-1)%time))
end if
timeEnd =self%outputTimes_%time(iOutput)
countTimesNew=int((timeEnd-timeStart)/self%timeStepMinimum)+1
! Construct new timesteps to span the the current output time.
!! Allocate sufficient space for timesteps and copy in the previous timesteps.
if (iOutput == 1) then
allocate(timesNew(countTimesNew))
allocate(indexMap(countTimesNew))
iStart = 1
else
allocate(timesNew(countTimesNew+size(self%intervals(iOutput-1)%time)))
allocate(indexMap(countTimesNew+size(self%intervals(iOutput-1)%time)))
timesNew(1:size(self%intervals(iOutput-1)%time))= self%intervals(iOutput-1)%time
iStart =size(self%intervals(iOutput-1)%time)+1
end if
!! Construct a mapping of indices between the current and previous interval.
indexMap=0_c_size_t
do iInterval=1,iStart-1
indexMap(iInterval)=iInterval
end do
!! Set times for the new timesteps.
do iNew=0,countTimesNew-1
timesNew(iNew+iStart)=timeStart+dble(iNew)*self%timeStepMinimum
if (File_Exists(fileName)) then
call File_Lock(char(fileName),fileLock,lockIsShared=.true.)
!$ call hdf5Access%set()
call file%openFile(char(fileName))
do iOutput=1,self%outputTimes_%count()
write (name,'(a,i4.4)') 'times' ,iOutput
call file%readDataset(name,self%intervals(iOutput)%time )
write (name,'(a,i4.4)') 'indexMap',iOutput
call file%readDataset(name,self%intervals(iOutput)%indexMap)
end do
!! Iteratively remove timesteps until we have no more than the permitted maximum.
do while (size(timesNew) > self%countTimeStepsMaximum)
! Evaluate our heuristic metric for which steps to combine. Our metric is that the logarithmic interval size, Δt/age
! should be minimized. We therefore find the global minimum of this metric across all current intervals, and then find the
! consecutive pair of intervals which, when combined, result in the smallest increase in the global minimum.
!! First, find the global minimum value of our metric.
metricMinimumGlobal=huge(0.0d0)
do iInterval=1,size(timesNew)
timeEnd=timesNew(iInterval)
if (iInterval == 1) then
timeStart=0.0d0
else
timeStart=timesNew(iInterval-1)
end if
metric =+( timeEnd -timeStart) &
& /(self%outputTimes_%time (iOutput)-timeStart)
metricMinimumGlobal=min(metric,metricMinimumGlobal)
call file%close()
!$ call hdf5Access%unset()
call File_Unlock(fileLock)
else
call Directory_Make(inputPath(pathTypeDataDynamic)//"/starFormation")
call File_Lock(char(fileName),fileLock,lockIsShared=.false.)
do iOutput=1,self%outputTimes_%count()
! Our start time is either the minimum timestep size (for the first output), or the end time of the timesteps of the
! previous output, plus the minimum timestep.
if (iOutput == 1) then
timeStart=self%timeStepMinimum
else
timeStart=self%timeStepMinimum+self%intervals(iOutput-1)%time(size(self%intervals(iOutput-1)%time))
end if
timeEnd =self%outputTimes_%time(iOutput)
countTimesNew=int((timeEnd-timeStart)/self%timeStepMinimum)+1
! Construct new timesteps to span the the current output time.
!! Allocate sufficient space for timesteps and copy in the previous timesteps.
if (iOutput == 1) then
allocate(timesNew(countTimesNew))
allocate(indexMap(countTimesNew))
iStart = 1
else
allocate(timesNew(countTimesNew+size(self%intervals(iOutput-1)%time)))
allocate(indexMap(countTimesNew+size(self%intervals(iOutput-1)%time)))
timesNew(1:size(self%intervals(iOutput-1)%time))= self%intervals(iOutput-1)%time
iStart =size(self%intervals(iOutput-1)%time)+1
end if
!! Construct a mapping of indices between the current and previous interval.
indexMap=0_c_size_t
do iInterval=1,iStart-1
indexMap(iInterval)=iInterval
end do
!! Determine which consecutive pair of intervals which, when combined, will result in the minimum increase in the global
!! minimum metric
iCombine =-1
metricChangeMinimum=huge(0.0d0)
do iInterval=2,size(timesNew)
! Metric if combined.
timeEnd=timesNew(iInterval)
if (iInterval == 2) then
timeStart=0.0d0
else
timeStart=timesNew(iInterval-2)
end if
metric=+( timeEnd -timeStart) &
& /(self%outputTimes_%time (iOutput)-timeStart)
! Change in the global metric minimum
metricChange=+metric &
& -metricMinimumGlobal
if (metricChange < metricChangeMinimum) then
metricChangeMinimum=metricChange
iCombine=iInterval
!! Set times for the new timesteps.
do iNew=0,countTimesNew-1
timesNew(iNew+iStart)=timeStart+dble(iNew)*self%timeStepMinimum
end do
!! Iteratively remove timesteps until we have no more than the permitted maximum.
do while (size(timesNew) > self%countTimeStepsMaximum)
! Evaluate our heuristic metric for which steps to combine. Our metric is that the logarithmic interval size, Δt/age
! should be minimized. We therefore find the global minimum of this metric across all current intervals, and then find the
! consecutive pair of intervals which, when combined, result in the smallest increase in the global minimum.
!! First, find the global minimum value of our metric.
metricMinimumGlobal=huge(0.0d0)
do iInterval=1,size(timesNew)
timeEnd=timesNew(iInterval)
if (iInterval == 1) then
timeStart=0.0d0
else
timeStart=timesNew(iInterval-1)
end if
metric =+( timeEnd -timeStart) &
& /(self%outputTimes_%time (iOutput)-timeStart)
metricMinimumGlobal=min(metric,metricMinimumGlobal)
end do
!! Determine which consecutive pair of intervals which, when combined, will result in the minimum increase in the global
!! minimum metric
iCombine =-1
metricChangeMinimum=huge(0.0d0)
do iInterval=2,size(timesNew)
! Metric if combined.
timeEnd=timesNew(iInterval)
if (iInterval == 2) then
timeStart=0.0d0
else
timeStart=timesNew(iInterval-2)
end if
metric=+( timeEnd -timeStart) &
& /(self%outputTimes_%time (iOutput)-timeStart)
! Change in the global metric minimum
metricChange=+metric &
& -metricMinimumGlobal
if (metricChange < metricChangeMinimum) then
metricChangeMinimum=metricChange
iCombine=iInterval
end if
end do
if (iCombine == -1) call Error_Report('no interval found - this should not happen'//{introspection:location})
! Combine the intervals.
allocate(timesNewTmp(size(timesNew)-1))
allocate(indexMapTmp(size(timesNew)-1))
if (iCombine > 2) then
timesNewTmp(1:iCombine-2)=timesNew(1:iCombine-2)
indexMapTmp(1:iCombine-2)=indexMap(1:iCombine-2)
end if
timesNewTmp(iCombine-1:size(timesNewTmp))=timesNew(iCombine:size(timesNew))
indexMapTmp(iCombine-1:size(indexMapTmp))=indexMap(iCombine:size(indexMap))
! Capture the final index if necessary.
if (iCombine > 1 .and. indexMap(iCombine) == 0) indexMapTmp(iCombine-1)=indexMap(iCombine-1)
deallocate(timesNew)
deallocate(indexMap)
call move_alloc(timesNewTmp,timesNew)
call move_alloc(indexMapTmp,indexMap)
end do
if (iCombine == -1) call Error_Report('no interval found - this should not happen'//{introspection:location})
! Combine the intervals.
allocate(timesNewTmp(size(timesNew)-1))
allocate(indexMapTmp(size(timesNew)-1))
if (iCombine > 2) then
timesNewTmp(1:iCombine-2)=timesNew(1:iCombine-2)
indexMapTmp(1:iCombine-2)=indexMap(1:iCombine-2)
end if
timesNewTmp(iCombine-1:size(timesNewTmp))=timesNew(iCombine:size(timesNew))
indexMapTmp(iCombine-1:size(indexMapTmp))=indexMap(iCombine:size(indexMap))
! Capture the final index if necessary.
if (iCombine > 1 .and. indexMap(iCombine) == 0) indexMapTmp(iCombine-1)=indexMap(iCombine-1)
deallocate(timesNew)
deallocate(indexMap)
call move_alloc(timesNewTmp,timesNew)
call move_alloc(indexMapTmp,indexMap)
! Number of intervals is now equal to or less than the maximum permitted. Store these intervals, and the map from the
! previous output's intervals.
call move_alloc(timesNew,self%intervals(iOutput)%time )
call move_alloc(indexMap,self%intervals(iOutput)%indexMap)
end do
! Number of intervals is now equal to or less than the maximum permitted. Store these intervals, and the map from the
! previous output's intervals.
call move_alloc(timesNew,self%intervals(iOutput)%time )
call move_alloc(indexMap,self%intervals(iOutput)%indexMap)
end do
!$ call hdf5Access%set()
call file%openFile(char(fileName),overWrite=.false.,readOnly=.false.)
do iOutput=1,self%outputTimes_%count()
write (name,'(a,i4.4)') 'times' ,iOutput
call file%writeDataset(self%intervals(iOutput)%time ,name)
write (name,'(a,i4.4)') 'indexMap',iOutput
call file%writeDataset(self%intervals(iOutput)%indexMap,name)
end do
call file%close()
!$ call hdf5Access%unset()
call File_Unlock(fileLock)
end if
! Construct output buffers.
allocate(self%starFormationHistoryBuffer(self%countTimeStepsMaximum,size(self%metallicityTable),self%countOutputBuffer,componentTypeMin:componentTypeMax))
allocate(self% indexOutputBuffer( componentTypeMin:componentTypeMax))
Expand Down

0 comments on commit 3702136

Please sign in to comment.