Skip to content

Commit

Permalink
Merge pull request #685 from galacticusorg/fixPositionInterpolation
Browse files Browse the repository at this point in the history
Position interpolation was broken - fix it
  • Loading branch information
abensonca authored Sep 1, 2024
2 parents fe58f3c + f30c252 commit 46273b7
Show file tree
Hide file tree
Showing 30 changed files with 866 additions and 109 deletions.
1 change: 1 addition & 0 deletions .github/workflows/cicd.yml
Original file line number Diff line number Diff line change
Expand Up @@ -1074,6 +1074,7 @@ jobs:
test-star-formation-histories-adaptive.py,
test-star-formation-histories-inSitu.py,
test-mostMassiveProgenitorIsSubhalo.py,
test-lightcone-crossing.py,
test-massProfileSatelliteBoundMass.py,
test-UniverseMachine.py,
test-stateRestore.pl,
Expand Down
95 changes: 84 additions & 11 deletions source/geometry.lightcones.square.F90
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,9 @@
(for {\normalfont \ttfamily X}$=2$ and 3). The {\normalfont \ttfamily redshift} parameters must list the redshifts of
available outputs.
</description>
<stateStore>
<stateStore variables="nodeOperator_" store="nodeOperatorStateStore_" restore="nodeOperatorStateRestore_" module="Functions_Global"/>
</stateStore>
</geometryLightcone>
!!]
type, extends(geometryLightconeClass) :: geometryLightconeSquare
Expand All @@ -73,6 +76,7 @@
class (cosmologyParametersClass), pointer :: cosmologyParameters_ => null()
class (cosmologyFunctionsClass ), pointer :: cosmologyFunctions_ => null()
class (outputTimesClass ), pointer :: outputTimes_ => null()
class (* ), pointer :: nodeOperator_ => null()
double precision , dimension(3,3) :: unitVector
double precision , dimension(3 ) :: origin , nodePositionCrossing , &
& nodeVelocityCrossing , unitVector1 , &
Expand Down Expand Up @@ -127,6 +131,9 @@
procedure :: periodicRange => squarePeriodicRange
procedure :: nodePositionReplicant => squareNodePositionReplicant
procedure :: nodeVelocityReplicant => squareNodeVelocityReplicant
procedure :: deepCopy => squareDeepCopy
procedure :: deepCopyReset => squareDeepCopyReset
procedure :: deepCopyFinalize => squareDeepCopyFinalize
end type geometryLightconeSquare

interface geometryLightconeSquare
Expand Down Expand Up @@ -156,10 +163,11 @@ function squareConstructorParameters(parameters) result(self)
Constructor for the {\normalfont \ttfamily square} lightcone geometry distribution class which takes a parameter list as
input.
!!}
use :: Cosmology_Parameters , only : cosmologyParameters , cosmologyParametersClass, hubbleUnitsLittleH
use :: Cosmology_Parameters , only : cosmologyParameters , cosmologyParametersClass, hubbleUnitsLittleH
use :: Error , only : Error_Report
use :: Input_Parameters , only : inputParameter , inputParameters
use :: Numerical_Constants_Astronomical, only : degreesToRadians , megaParsec
use :: Input_Parameters , only : inputParameter , inputParameters
use :: Numerical_Constants_Astronomical, only : degreesToRadians , megaParsec
use :: Functions_Global , only : nodeOperatorConstruct_, nodeOperatorDestruct_
implicit none
type (geometryLightconeSquare ) :: self
type (inputParameters ), intent(inout) :: parameters
Expand All @@ -169,6 +177,7 @@ function squareConstructorParameters(parameters) result(self)
class (cosmologyParametersClass), pointer :: cosmologyParameters_
class (cosmologyFunctionsClass ), pointer :: cosmologyFunctions_
class (outputTimesClass ), pointer :: outputTimes_
class (* ), pointer :: nodeOperator_
double precision :: lengthReplication , angularSize , &
& lengthUnitsInSI , unitConversionLength
integer :: lengthHubbleExponent
Expand Down Expand Up @@ -240,6 +249,7 @@ function squareConstructorParameters(parameters) result(self)
<objectBuilder class="cosmologyParameters" name="cosmologyParameters_" source="parameters"/>
<objectBuilder class="outputTimes" name="outputTimes_" source="parameters"/>
!!]
call nodeOperatorConstruct_(parameters,nodeOperator_)
! Convert angle to radians.
angularSize=angularSize*degreesToRadians
! Convert lengths units internal units.
Expand All @@ -249,19 +259,20 @@ function squareConstructorParameters(parameters) result(self)
origin =origin *unitConversionLength
lengthReplication =lengthReplication*unitConversionLength
! Construct the object.
self =geometryLightconeSquare(origin,unitVector,angularSize,lengthReplication,timeEvolvesAlongLightcone,cosmologyParameters_,cosmologyFunctions_,outputTimes_)
self =geometryLightconeSquare(origin,unitVector,angularSize,lengthReplication,timeEvolvesAlongLightcone,cosmologyParameters_,cosmologyFunctions_,outputTimes_,nodeOperator_)
self%lengthUnitsInSI =lengthUnitsInSI
self%lengthHubbleExponent=lengthHubbleExponent
!![
<inputParametersValidate source="parameters"/>
<inputParametersValidate source="parameters" extraAllowedNames="nodeOperator"/>
<objectDestructor name="cosmologyParameters_"/>
<objectDestructor name="cosmologyFunctions_" />
<objectDestructor name="outputTimes_" />
!!]
if (associated(nodeOperator_)) call nodeOperatorDestruct_(nodeOperator_)
return
end function squareConstructorParameters

function squareConstructorInternal(origin,unitVector,angularSize,lengthReplication,timeEvolvesAlongLightcone,cosmologyParameters_,cosmologyFunctions_,outputTimes_) result(self)
function squareConstructorInternal(origin,unitVector,angularSize,lengthReplication,timeEvolvesAlongLightcone,cosmologyParameters_,cosmologyFunctions_,outputTimes_,nodeOperator_) result(self)
!!{
Internal constructor for the {\normalfont \ttfamily square} lightcone geometry distribution class.
!!}
Expand All @@ -278,6 +289,7 @@ function squareConstructorInternal(origin,unitVector,angularSize,lengthReplicati
class (cosmologyParametersClass), target , intent(in ) :: cosmologyParameters_
class (cosmologyFunctionsClass ), target , intent(in ) :: cosmologyFunctions_
class (outputTimesClass ), target , intent(in ) :: outputTimes_
class (* ), target , intent(in ) :: nodeOperator_
double precision , intent(in ) :: lengthReplication , angularSize
logical , intent(in ) :: timeEvolvesAlongLightcone
double precision , parameter :: orthogonalityTolerance =1.000d-6
Expand All @@ -291,7 +303,7 @@ function squareConstructorInternal(origin,unitVector,angularSize,lengthReplicati
character (len=12 ) :: label
type (varying_string ) :: message
!![
<constructorAssign variables="origin, unitVector, angularSize, lengthReplication, timeEvolvesAlongLightcone, *cosmologyParameters_, *cosmologyFunctions_, *outputTimes_"/>
<constructorAssign variables="origin, unitVector, angularSize, lengthReplication, timeEvolvesAlongLightcone, *cosmologyParameters_, *cosmologyFunctions_, *outputTimes_, *nodeOperator_"/>
!!]

! Store unit vectors.
Expand Down Expand Up @@ -424,6 +436,7 @@ subroutine squareDestructor(self)
!!{
Destructor for the {\normalfont \ttfamily square} lightcone geometry distribution class.
!!}
use :: Functions_Global, only : nodeOperatorDestruct_
implicit none
type(geometryLightconeSquare), intent(inout) :: self

Expand All @@ -432,6 +445,7 @@ subroutine squareDestructor(self)
<objectDestructor name="self%cosmologyFunctions_" />
<objectDestructor name="self%outputTimes_" />
!!]
if (associated(self%nodeOperator_)) call nodeOperatorDestruct_(self%nodeOperator_)
return
end subroutine squareDestructor

Expand Down Expand Up @@ -1076,7 +1090,8 @@ function squareNodePositionReplicant(self,node,time,origin,replicant,setTime,pos
!!{
Compute the comoving position of the given node in the given replicant.
!!}
use :: Galacticus_Nodes, only : nodeComponentBasic, nodeComponentPosition
use :: Galacticus_Nodes, only : nodeComponentBasic , nodeComponentPosition
use :: Functions_Global, only : nodeOperatorSolveAnalytics_
implicit none
double precision , dimension(3) :: squareNodePositionReplicant
class (geometryLightconeSquare), intent(inout) :: self
Expand All @@ -1098,7 +1113,7 @@ function squareNodePositionReplicant(self,node,time,origin,replicant,setTime,pos
if (setTime_) then
basic => node %basic()
timeOriginal = basic%time ()
call basic%timeSet(time)
call nodeOperatorSolveAnalytics_(self%nodeOperator_,node,time)
end if
position => node %position ( )
expansionFactor = +self %cosmologyFunctions_%expansionFactor(time)
Expand All @@ -1114,8 +1129,8 @@ function squareNodePositionReplicant(self,node,time,origin,replicant,setTime,pos
end if
do i=1,3
squareNodePositionReplicant(i)=Dot_Product(positionComovingNode-origin+self%lengthReplication*dble(replicant),self%unitVector(:,i))
end do
if (setTime_) call basic%timeSet(timeOriginal)
end do
if (setTime_) call nodeOperatorSolveAnalytics_(self%nodeOperator_,node,timeOriginal)
return
end function squareNodePositionReplicant

Expand Down Expand Up @@ -1153,3 +1168,61 @@ function squareNodeVelocityReplicant(self,node,time,replicant,setTime)
if (setTime_) call basic%timeSet(timeOriginal)
return
end function squareNodeVelocityReplicant

subroutine squareDeepCopyReset(self)
!!{
Perform a deep copy reset of the object.
!!}
use :: Functions_Global, only : nodeOperatorDeepCopyReset_
implicit none
class(geometryLightconeSquare), intent(inout) :: self

self%copiedSelf => null()
if (associated(self%cosmologyparameters_)) call self%cosmologyparameters_%deepCopyReset()
if (associated(self%cosmologyfunctions_ )) call self%cosmologyfunctions_ %deepCopyReset()
if (associated(self%outputtimes_ )) call self%outputtimes_ %deepCopyReset()
if (associated(self%nodeOperator_ )) call nodeOperatorDeepCopyReset_(self%nodeOperator_)
return
end subroutine squareDeepCopyReset

subroutine squareDeepCopyFinalize(self)
!!{
Finalize a deep reset of the object.
!!}
use :: Functions_Global, only : nodeOperatorDeepCopyFinalize_
implicit none
class(geometryLightconeSquare), intent(inout) :: self

if (associated(self%cosmologyparameters_)) call self%cosmologyparameters_%deepCopyFinalize()
if (associated(self%cosmologyfunctions_ )) call self%cosmologyfunctions_ %deepCopyFinalize()
if (associated(self%outputtimes_ )) call self%outputtimes_ %deepCopyFinalize()
if (associated(self%nodeOperator_ )) call nodeOperatorDeepCopyFinalize_(self%nodeOperator_)
return
end subroutine squareDeepCopyFinalize

subroutine squareDeepCopy(self,destination)
!!{
Perform a deep copy of the object.
!!}
use :: Functions_Global, only : nodeOperatorDeepCopy_
implicit none
class(geometryLightconeSquare), intent(inout), target :: self
class(geometryLightconeClass ), intent(inout) :: destination

call self%deepCopy_(destination)
select type (destination)
type is (geometryLightconeSquare)
nullify(destination%nodeOperator_)
if (associated(self%nodeOperator_)) then
allocate(destination%nodeOperator_,mold=self%nodeOperator_)
call nodeOperatorDeepCopy_(self%nodeOperator_,destination%nodeOperator_)
#ifdef OBJECTDEBUG
if (debugReporting.and.mpiSelf%isMaster()) call displayMessage(var_str('functionClass[own] (class : ownerName : ownerLoc : objectLoc : sourceLoc): galacticstructure : [destination] : ')//loc(destination)//' : '//loc(destination%nodeOperator_)//' : '//{introspection:location:compact},verbosityLevelSilent)
#endif
end if
class default
call Error_Report('destination and source types do not match'//{introspection:location})
end select
return
end subroutine squareDeepCopy

116 changes: 116 additions & 0 deletions source/nodes.operators.physics.position.discrete.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018,
!! 2019, 2020, 2021, 2022, 2023, 2024
!! Andrew Benson <abenson@carnegiescience.edu>
!!
!! This file is part of Galacticus.
!!
!! Galacticus is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! Galacticus is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with Galacticus. If not, see <http://www.gnu.org/licenses/>.

!!{
Implements a node operator class that updates positions of nodes in discrete steps.
!!}

!![
<nodeOperator name="nodeOperatorPositionDiscrete">
<description>
A node operator class that interpolates positions of nodes in discrete steps.
</description>
</nodeOperator>
!!]
type, extends(nodeOperatorClass) :: nodeOperatorPositionDiscrete
!!{
A node operator class that interpolates positions of nodes in discrete steps.
!!}
private
contains
procedure :: differentialEvolutionStepFinalState => positionDiscreteDifferentialEvolutionStepFinalState
end type nodeOperatorPositionDiscrete

interface nodeOperatorPositionDiscrete
!!{
Constructors for the {\normalfont \ttfamily positionDiscrete} node operator class.
!!}
module procedure positionDiscreteConstructorParameters
end interface nodeOperatorPositionDiscrete

contains

function positionDiscreteConstructorParameters(parameters) result(self)
!!{
Constructor for the {\normalfont \ttfamily positionDiscrete} node operator class which takes a parameter set as input.
!!}
use :: Input_Parameters, only : inputParameters
implicit none
type(nodeOperatorPositionDiscrete) :: self
type(inputParameters ), intent(inout) :: parameters

self=nodeOperatorPositionDiscrete()
!![
<inputParametersValidate source="parameters"/>
!!]
return
end function positionDiscreteConstructorParameters

subroutine positionDiscreteDifferentialEvolutionAnalytics(self,node)
!!{
Mark analytically-solvable properties.
!!}
use :: Galacticus_Nodes, only : nodeComponentPosition
implicit none
class(nodeOperatorPositionDiscrete), intent(inout) :: self
type (treeNode ), intent(inout) :: node
class(nodeComponentPosition ), pointer :: position
!$GLC attributes unused :: self

position => node%position()
call position%positionAnalytic()
call position%velocityAnalytic()
return
end subroutine positionDiscreteDifferentialEvolutionAnalytics

subroutine positionDiscreteDifferentialEvolutionStepFinalState(self,node)
!!{
Compute the discrete position and velocity of the node.
!!}
use, intrinsic :: ISO_C_Binding , only : c_size_t
use :: Galacticus_Nodes , only : nodeComponentPosition, nodeComponentBasic
use :: Numerical_Interpolation, only : interpolator
use :: Histories , only : history
implicit none
class (nodeOperatorPositionDiscrete), intent(inout) :: self
type (treeNode ), intent(inout) :: node
class (nodeComponentBasic ), pointer :: basic
class (nodeComponentPosition ), pointer :: position
integer , parameter :: positionBegin=1, positionEnd=3, &
& velocityBegin=4, velocityEnd=6
type (interpolator ) :: interpolator_
type (history ) :: history_
integer (c_size_t ) :: iTime
!$GLC attributes unused :: self

! Check if this node has a position history attached to it.
position => node %position ()
history_ = position%positionHistory()
if (history_%exists()) then
basic => node%basic()
if (history_%time(1) <= basic%time()) then
interpolator_=interpolator (history_%time )
iTime =interpolator_%locate(basic %time(),closest=.true.)
call position%positionSet(history_%data(iTime,positionBegin:positionEnd))
call position%velocitySet(history_%data(iTime,velocityBegin:velocityEnd))
end if
end if
return
end subroutine positionDiscreteDifferentialEvolutionStepFinalState

8 changes: 5 additions & 3 deletions source/nodes.operators.physics.position.interpolated.F90
Original file line number Diff line number Diff line change
Expand Up @@ -253,8 +253,8 @@ subroutine positionInterpolatedDifferentialEvolutionSolveAnalytics(self,node,tim
& +coefficientsCubic(3,i)*time &
& +coefficientsCubic(4,i)
velocity_(i)=+3.0d0*coefficientsCubic(1,i)*time**2 &
& +2.0d0*coefficientsCubic(2,i)*time &
& + coefficientsCubic(3,i)
& +2.0d0*coefficientsCubic(2,i)*time &
& + coefficientsCubic(3,i)
end do
! Convert from comoving back to physical position/velocity, and to km/s.
position_=+position_ &
Expand Down Expand Up @@ -303,6 +303,9 @@ subroutine positionInterpolatedDifferentialEvolutionSolveAnalytics(self,node,tim
end do
end do
end if
! Set the position and velocity.
call position%positionSet(position_)
call position%velocitySet(velocity_)
return
end subroutine positionInterpolatedDifferentialEvolutionSolveAnalytics

Expand Down Expand Up @@ -598,7 +601,6 @@ subroutine positionInterpolatedComputeInterpolation(self,node)
positionComoving(2,:) = positionHistory %data (i+1,1:3 )
velocityComoving(1,:) = positionHistory %data (i ,4:6 )
velocityComoving(2,:) = positionHistory %data (i+1,4:6 )

end if
else
! We have no means to determine the start/end points for this interval. This node must have no future, and so we can
Expand Down
Loading

0 comments on commit 46273b7

Please sign in to comment.