diff --git a/source/star_formation.histories.adaptive.F90 b/source/star_formation.histories.adaptive.F90 index b058dfb610..eb4b8bf995 100644 --- a/source/star_formation.histories.adaptive.F90 +++ b/source/star_formation.histories.adaptive.F90 @@ -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 @@ -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 !![ !!] @@ -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))