Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Optimize starFormationRateSurfaceDensityDisksBlitz2006 table construction #689

Merged
merged 1 commit into from
Sep 6, 2024
Merged

Optimize starFormationRateSurfaceDensityDisksBlitz2006 table construction #689

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 58 additions & 48 deletions source/star_formation.rate_surface_density.disks.Blitz2006.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -719,14 +719,15 @@ double precision function blitz2006IntegralPartiallyMolecular(self,coefficientNo
!!}
implicit none
class (starFormationRateSurfaceDensityDisksBlitz2006), intent(inout) :: self
double precision , intent(in ) :: radiusInner , radiusOuter , &
& coefficientNormalization , coefficientFactorBoost , &
double precision , intent(in ) :: radiusInner , radiusOuter , &
& coefficientNormalization , coefficientFactorBoost , &
& coefficientFactorBoostStellar
double precision , parameter :: coefficientFactorBoostStellarLarge=1.0d+4, radiusScaleFreeTiny =1.0d-3
double precision , parameter :: coefficientFactorBoostStellarLarge=1.0d+4, radiusScaleFreeTiny =1.0d-3
double precision , parameter :: coefficientFactorBoostTiny =1.0d-6, coefficientFactorBoostStellarTiny=1.0d-6
double precision , save :: coefficientFactorBoost_ , coefficientFactorBoostStellar_
!$omp threadprivate(coefficientFactorBoost_,coefficientFactorBoostStellar_)
integer :: i
double precision :: multiplier , radiusScaleFree , &
double precision :: multiplier , radiusScaleFree , &
& radiusScaleFreeInner , radiusScaleFreeOuter

! Compute scale free radii.
Expand Down Expand Up @@ -859,7 +860,6 @@ double precision function integralAnalyticPartiallyMolecularGeneric(radiusScaleF
implicit none
double precision , intent(in ) :: radiusScaleFree
double precision , parameter :: toleranceRelative =1.0d-6
double precision , parameter :: coefficientFactorBoostTiny =1.0d-6, coefficientFactorBoostStellarTiny =1.0d-6
integer , parameter :: pointsPerDecadeFactorBoost =30 , pointsPerDecadeFactorBoostStellar =30 , &
& pointsPerDecadeRadius =30
integer :: countFactorBoost , countFactorBoostStellar , &
Expand All @@ -884,9 +884,9 @@ double precision function integralAnalyticPartiallyMolecularGeneric(radiusScaleF
logical :: haveLock
!$omp threadprivate(message,file,fileLock)

! Read table if we do not have it.
! If our table is insufficient (or does not yet exist), attempt to read the table from file.
haveLock=.false.
if (.not.self%tableInitialized) then
if (tableIsInsufficient()) then
call Directory_Make(char(File_Path(char(self%filenameTable))))
call File_Lock(char(self%filenameTable),fileLock,lockIsShared=.false.)
haveLock=.true.
Expand Down Expand Up @@ -914,28 +914,8 @@ double precision function integralAnalyticPartiallyMolecularGeneric(radiusScaleF
self%tableInitialized=.true.
end if
end if
! Recompute table if it has insufficient range.
if ( &
& ( &
& coefficientFactorBoost < self%coefficientFactorBoostMinimum &
& .and. &
& coefficientFactorBoostTiny < self%coefficientFactorBoostMinimum &
& ) &
& .or. &
& coefficientFactorBoost > self%coefficientFactorBoostMaximum &
& .or. &
& ( &
& coefficientFactorBoostStellar < self%coefficientFactorBoostStellarMinimum &
& .and. &
& coefficientFactorBoostStellarTiny < self%coefficientFactorBoostStellarMinimum &
& ) &
& .or. &
& coefficientFactorBoostStellar > self%coefficientFactorBoostStellarMaximum &
& .or. &
& radiusScaleFree < self%radiusScaleFreeMinimum &
& .or. &
& radiusScaleFree > self%radiusScaleFreeMaximum &
& ) then
! Having read the table from file (if it exists), check again to see if it is sufficient. If it is not, we must retabulate.
if (tableIsInsufficient()) then
! Obtain a file lock if we don't already have one.
if (.not.haveLock) then
call File_Lock(char(self%filenameTable),fileLock,lockIsShared=.false.)
Expand Down Expand Up @@ -1064,31 +1044,61 @@ double precision function integralAnalyticPartiallyMolecularGeneric(radiusScaleF
hhi= +hi
end if
do jj=0,1
if (jj == 0) then
hhj=+1.0d0-hj
if (jj == 0) then
hhj=+1.0d0-hj
else
hhj= +hj
end if
do kk=0,1
if (kk == 0) then
hhk=+1.0d0-hk
else
hhj= +hj
hhk= +hk
end if
do kk=0,1
if (kk == 0) then
hhk=+1.0d0-hk
else
hhk= +hk
end if
integral=+integral &
& +self%integralPartiallyMolecularTable(i+ii,j+jj,k+kk) &
& * hhi &
& * hhj &
& * hhk
end do
integral=+integral &
& +self%integralPartiallyMolecularTable(i+ii,j+jj,k+kk) &
& * hhi &
& * hhj &
& * hhk
end do
end do
integralAnalyticPartiallyMolecularGeneric=+integralAnalyticPartiallyMolecularGeneric &
& +exp(integral) &
& *coefficientNormalization
end do
integralAnalyticPartiallyMolecularGeneric=+integralAnalyticPartiallyMolecularGeneric &
& +exp(integral) &
& *coefficientNormalization
return
end function integralAnalyticPartiallyMolecularGeneric


logical function tableIsInsufficient()
!!{
Determine if the current table is insufficient for our purposes.
!!}
implicit none

tableIsInsufficient= .not. self%tableInitialized &
& .or. &
& ( &
& coefficientFactorBoost < self%coefficientFactorBoostMinimum &
& .and. &
& coefficientFactorBoostTiny < self%coefficientFactorBoostMinimum &
& ) &
& .or. &
& coefficientFactorBoost > self%coefficientFactorBoostMaximum &
& .or. &
& ( &
& coefficientFactorBoostStellar < self%coefficientFactorBoostStellarMinimum &
& .and. &
& coefficientFactorBoostStellarTiny < self%coefficientFactorBoostStellarMinimum &
& ) &
& .or. &
& coefficientFactorBoostStellar > self%coefficientFactorBoostStellarMaximum &
& .or. &
& radiusScaleFree < self%radiusScaleFreeMinimum &
& .or. &
& radiusScaleFree > self%radiusScaleFreeMaximum
return
end function tableIsInsufficient

double precision function integrand(radiusScaleFree)
!!{
Integrand for the partially molecular case.
Expand Down
Loading