Add an abstract class which automates tabulation of parametric, spherically-symmetric, collisionless mass distributions

.github/workflows/cicd.yml

@@ -839,6 +839,7 @@ jobs:
                       tests.interpolation.2D.exe,
                       tests.make_ranges.exe,
                       tests.mass_distributions.exe,
+                      tests.mass_distributions.tabulated.exe,
                       tests.math_special_functions.exe,
                       tests.math_distributions.exe,
                       tests.math.fast.exe,

scripts/build/sourceDigests.pl

@@ -223,7 +223,8 @@
 (my $outputFileName = $targetName) =~ s/\.(exe|o)$/.md5s.c/;
 open(my $outputFile,">".$ENV{'BUILDPATH'}."/".$outputFileName);
 foreach my $hashName ( sort(keys(%{$digestsPerFile->{'types'}})) ) {
-    print $outputFile "char ".$hashName."MD5[]=\"".$digestsPerFile->{'types'}->{$hashName}->{'compositeMD5'}."\";\n"
+    (my $digest = $digestsPerFile->{'types'}->{$hashName}->{'compositeMD5'}) =~ s/\//\@/g;
+    print $outputFile "char ".$hashName."MD5[]=\"".$digest."\";\n"
 }
 close($outputFile);
 

source/kinematic_distributions.collisionless.tabulated.F90

@@ -0,0 +1,139 @@
+!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018,
+!!           2019, 2020, 2021, 2022, 2023, 2024, 2025
+!!    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/>.
+
+  !!{
+  Implementation of a kinematic distribution class for collisionless mass distributions using tabulated solutions.
+  !!}
+
+  !![
+  <kinematicsDistribution name="kinematicsDistributionCollisionlessTabulated">
+   <description>A kinematic distribution class for collisionless mass distributions using tabulated solutions.</description>
+  </kinematicsDistribution>
+  !!]
+  type, public, extends(kinematicsDistributionClass) :: kinematicsDistributionCollisionlessTabulated
+     !!{
+     A kinematics distribution for collisionless distributions using tabulated solutions.
+     !!}
+   contains
+     procedure :: isCollisional        => collisionlessTabulatedIsCollisional
+     procedure :: velocityDispersion1D => collisionlessTabulatedVelocityDispersion1D
+  end type kinematicsDistributionCollisionlessTabulated
+
+  interface kinematicsDistributionCollisionlessTabulated
+     !!{
+     Constructors for the {\normalfont \ttfamily collisionless} kinematic distribution class.
+     !!}
+     module procedure collisionlessTabulatedConstructorParameters
+     module procedure collisionlessTabulatedConstructorInternal
+     module procedure collisionlessTabulatedConstructorDecorated
+  end interface kinematicsDistributionCollisionlessTabulated
+
+contains
+
+  function collisionlessTabulatedConstructorParameters(parameters) result(self)
+    !!{
+    Constructor for the {\normalfont \ttfamily collisionlessTabulated} kinematic distribution class which builds the object from a parameter
+    set.
+    !!}
+    use :: Input_Parameters, only : inputParameters
+    implicit none
+    type            (kinematicsDistributionCollisionlessTabulated)                :: self
+    type            (inputParameters                             ), intent(inout) :: parameters
+    double precision                                                              :: toleranceRelativeVelocityDispersion, toleranceRelativeVelocityDispersionMaximum
+
+    !![
+    <inputParameter>
+      <name>toleranceRelativeVelocityDispersion</name>
+      <defaultValue>1.0d-6</defaultValue>
+      <source>parameters</source>
+      <description>The relative tolerance to use in numerical solutions for the velocity dispersion in dark-matter-only density profiles.</description>
+    </inputParameter>
+    <inputParameter>
+      <name>toleranceRelativeVelocityDispersionMaximum</name>
+      <defaultValue>1.0d-3</defaultValue>
+      <source>parameters</source>
+      <description>The maximum relative tolerance to use in numerical solutions for the velocity dispersion in dark-matter-only density profiles.</description>
+    </inputParameter>
+    !!]
+    self=kinematicsDistributionCollisionlessTabulated(toleranceRelativeVelocityDispersion,toleranceRelativeVelocityDispersionMaximum)
+    !![
+    <inputParametersValidate source="parameters"/>
+    !!]
+    return
+  end function collisionlessTabulatedConstructorParameters
+
+  function collisionlessTabulatedConstructorInternal(toleranceRelativeVelocityDispersion,toleranceRelativeVelocityDispersionMaximum) result(self)
+    !!{
+    Internal constructor for the {\normalfont \ttfamily collisionlessTabulated} kinematic distribution class.
+    !!}
+    implicit none
+    type            (kinematicsDistributionCollisionlessTabulated)                          :: self
+    double precision                                              , intent(in   ), optional :: toleranceRelativeVelocityDispersion, toleranceRelativeVelocityDispersionMaximum
+    !![
+    <constructorAssign variables="toleranceRelativeVelocityDispersion, toleranceRelativeVelocityDispersionMaximum"/>
+    !!]
+
+    return
+  end function collisionlessTabulatedConstructorInternal
+
+  function collisionlessTabulatedConstructorDecorated(kinematicsDistribution_) result(self)
+    !!{
+    Internal constructor for the {\normalfont \ttfamily collisionlessTabulated} kinematic distribution class.
+    !!}
+    implicit none
+    type (kinematicsDistributionCollisionlessTabulated)                :: self
+    class(kinematicsDistributionClass                 ), intent(in   ) :: kinematicsDistribution_
+
+    self%toleranceRelativeVelocityDispersion       =kinematicsDistribution_%toleranceRelativeVelocityDispersion
+    self%toleranceRelativeVelocityDispersionMaximum=kinematicsDistribution_%toleranceRelativeVelocityDispersionMaximum
+    return
+  end function collisionlessTabulatedConstructorDecorated
+
+  logical function collisionlessTabulatedIsCollisional(self)
+    !!{
+    Return false indicating that the tabulated collisionless kinematic distribution represents collisionless particles.
+    !!}
+    implicit none
+    class(kinematicsDistributionCollisionlessTabulated), intent(inout) :: self
+
+    collisionlessTabulatedIsCollisional=.false.
+    return
+  end function collisionlessTabulatedIsCollisional
+
+  double precision function collisionlessTabulatedVelocityDispersion1D(self,coordinates,massDistributionEmbedding) result(velocityDispersion)
+    !!{
+    Return the 1D velocity dispersion at the specified {\normalfont \ttfamily coordinates} in a tabulated collisionless kinematic distribution.
+    !!}
+    implicit none
+    class(kinematicsDistributionCollisionlessTabulated), intent(inout), target :: self
+    class(coordinate                                  ), intent(in   )         :: coordinates
+    class(massDistributionClass                       ), intent(inout)         :: massDistributionEmbedding
+
+    select type (massDistributionEmbedding)
+    class is (massDistributionSphericalTabulated)
+       if (.not.massDistributionEmbedding%isTabulating()) then
+          velocityDispersion=massDistributionEmbedding%velocityDispersion1D         (coordinates                          )
+       else
+          velocityDispersion=self                     %velocityDispersion1DNumerical(coordinates,massDistributionEmbedding)
+       end if
+    class default
+       velocityDispersion   =self                     %velocityDispersion1DNumerical(coordinates,massDistributionEmbedding)
+    end select
+    return
+  end function collisionlessTabulatedVelocityDispersion1D

source/mass_distributions.F90

@@ -643,8 +643,7 @@ module Mass_Distributions
     <argument>class(coordinate           ), intent(in   ) :: coordinates              </argument>
     <argument>class(massDistributionClass), intent(inout) :: massDistributionEmbedding</argument>
     <code>
-      !$GLC attributes unused :: self, coordinates, massDistributionEmbedding
-      kinematicsDistributionVelocityDispersion1D=0.0d0
+      kinematicsDistributionVelocityDispersion1D=self%velocityDispersion1DNumerical(coordinates,massDistributionEmbedding)
     </code>
    </method>
    <method name="velocityDispersion1DNumerical" >

source/mass_distributions.spherical.F90

@@ -522,8 +522,6 @@ double precision function sphericalPotentialNumerical(self,coordinates,status) r
     double precision                                   , dimension(:) , allocatable :: potentials                  , radii               , &
          &                                                                             potentials_                 , radii_
     double precision                                   , parameter                  :: countPointsPerOctave=4.0d+0
-    double precision                                   , parameter                  :: radiusMaximumFactor =1.0d+1
-    double precision                                   , parameter                  :: toleranceRelative   =1.0d-6
     type            (integrator                       )                             :: integrator_
     integer         (c_size_t                )                                      :: countRadii                  , iMinimum            , &
          &                                                                             iMaximum                    , iPrevious           , &
@@ -707,10 +705,9 @@ double precision function sphericalPotentialDifferenceNumerical(self,coordinates
     use :: Numerical_Comparison            , only : Values_Agree
     implicit none
     class           (massDistributionSpherical        ), intent(inout), target   :: self
-    class           (coordinate                       ), intent(in   )           :: coordinates1              , coordinates2
+    class           (coordinate                       ), intent(in   )           :: coordinates1            , coordinates2
     type            (enumerationStructureErrorCodeType), intent(  out), optional :: status
-    double precision                                   , parameter               :: toleranceRelative  =1.0d-3
-    double precision                                   , parameter               :: radiusMaximumFactor=1.0d+1
+    double precision                                   , parameter               :: toleranceRelative=1.0d-3
     type            (integrator                       )                          :: integrator_
 
     if (present(status)) status=structureErrorCodeSuccess