Merge pull request CambridgeNuclear#126 from valeriaRaffuzzi/densityR…

…esponse

Density response

ParticleObjects/Tests/particle_test.f90

@@ -1,5 +1,6 @@
 module particle_test
   use numPrecision
+  use universalVariables
   use particle_class, only : particle, particleState, P_NEUTRON, P_PHOTON, verifyType
   use pFUnit_mod
 
@@ -247,6 +248,9 @@ subroutine testMiscAccess(this)
     @assertEqual(3, cellIdx, 'Cell Index. Level 1.')
     @assertEqual(1, uniIdx, 'Universe Index. Level 1.')
 
+    ! Verify getting velocity
+    @assertEqual(lightSpeed, this % p_CE % getVelocity())
+
   end subroutine testMiscAccess
 
   !!
@@ -306,6 +310,7 @@ subroutine testMovementProcedures(this)
 
     @assertEqual(r0_lvl1, r_lvl1, 'Global position after global teleport')
     @assertTrue(this % p_CE % coords % isAbove(), 'Particle is above geometry')
+
   end subroutine testMovementProcedures
 
   !!

ParticleObjects/particle_class.f90

@@ -5,14 +5,15 @@ module particle_class
   use genericProcedures
   use coord_class,       only : coordList
   use RNG_class,         only : RNG
+  use errors_mod,        only : fatalError
 
   implicit none
   private
 
   !!
   !! Particle types paramethers
   !!
-  integer(shortInt), parameter,public :: P_NEUTRON = 1,&
+  integer(shortInt), parameter,public :: P_NEUTRON = 1, &
                                          P_PHOTON  = 2
 
   !!
@@ -122,7 +123,7 @@ module particle_class
     generic              :: build => buildCE, buildMG
     generic              :: assignment(=) => particle_fromParticleState
 
-    ! Inquiry about coordinates
+    ! Enquiry about coordinates
     procedure                  :: rLocal
     procedure                  :: rGlobal
     procedure                  :: dirLocal
@@ -133,14 +134,17 @@ module particle_class
     procedure                  :: matIdx
     procedure, non_overridable :: getType
 
+    ! Enquiry about physical state
+    procedure :: getVelocity
+
     ! Operations on coordinates
-    procedure            :: moveGlobal
-    procedure            :: moveLocal
-    procedure            :: rotate
-    procedure            :: teleport
-    procedure            :: point
-    procedure            :: takeAboveGeom
-    procedure            :: setMatIdx
+    procedure :: moveGlobal
+    procedure :: moveLocal
+    procedure :: rotate
+    procedure :: teleport
+    procedure :: point
+    procedure :: takeAboveGeom
+    procedure :: setMatIdx
 
     ! Save particle state information
     procedure, non_overridable  :: savePreHistory
@@ -421,6 +425,48 @@ pure function getType(self) result(type)
 
   end function getType
 
+  !!
+  !! Return the particle velocity in [cm/s]
+  !! neutronMass: [MeV]
+  !! lightSpeed:  [cm/s]
+  !!
+  !! NOTE:
+  !!   The velocities are computed from non-relativistic formula for massive particles.
+  !!   A small error might appear in MeV range (e.g. for fusion applications)
+  !!
+  !! Args:
+  !!   None
+  !!
+  !! Result:
+  !!   Particle velocity
+  !!
+  !! Errors:
+  !!   fatalError if the particle type is neither P_NEUTRON nor P_PHOTON
+  !!   fatalError if the particle is MG
+  !!
+  function getVelocity(self) result(velocity)
+    class(particle), intent(in) :: self
+    real(defReal)               :: velocity
+    character(100), parameter   :: Here = 'getVelocity (particle_class.f90)'
+
+    ! Verify the particle is not MG
+    if (self % isMG) call fatalError(Here, 'Velocity cannot be calculated for MG particle')
+
+    ! Calculates the velocity for the relevant particle [cm/s]
+    if (self % type == P_NEUTRON) then
+      velocity = sqrt(TWO * self % E / neutronMass) * lightSpeed
+
+    elseif (self % type == P_PHOTON) then
+      velocity = lightSpeed
+
+    else
+      call fatalError(Here, 'Particle type requested is neither neutron (1) nor photon (2). It is: ' &
+                            & //numToChar(self % type))
+
+    end if
+
+  end function getVelocity
+
 !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
 !! Particle operations on coordinates procedures
 !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>

SharedModules/genericProcedures.f90

@@ -4,7 +4,7 @@ module genericProcedures
 
   use numPrecision
   use openmp_func, only : ompGetMaxThreads
-  use errors_mod, only: fatalError
+  use errors_mod,  only : fatalError
   use endfConstants
   use universalVariables
 

SharedModules/universalVariables.f90

@@ -75,7 +75,8 @@ module universalVariables
                                   TRACKING_XS = 3
 
   ! Physical constants
-  real(defReal), parameter :: neutronMass = 939.5654133_defReal, &   ! Neutron mass in MeV/c^2
+  ! Neutron mass and speed of light in vacuum from from https://physics.nist.gov/cuu/Constants/index.html
+  real(defReal), parameter :: neutronMass = 939.56542194_defReal,  & ! Neutron mass in MeV (m*c^2)
                               lightSpeed  = 2.99792458e10_defReal, & ! Light speed in cm/s
                               energyPerFission = 200.0_defReal       ! MeV
 

Tallies/TallyResponses/CMakeLists.txt

@@ -6,6 +6,7 @@ add_sources(./tallyResponse_inter.f90
             ./macroResponse_class.f90
             ./microResponse_class.f90
             ./weightResponse_class.f90
+            ./densityResponse_class.f90
             ./testResponse_class.f90
             )
 
@@ -15,4 +16,5 @@ add_unit_tests(./Tests/fluxResponse_test.f90
                ./Tests/macroResponse_test.f90
                ./Tests/microResponse_test.f90
                ./Tests/weightResponse_test.f90
-                )
+               ./Tests/densityResponse_test.f90
+               )

Tallies/TallyResponses/Tests/densityResponse_test.f90

@@ -0,0 +1,74 @@
+module densityResponse_test
+
+  use numPrecision
+  use universalVariables,    only : neutronMass, lightSpeed
+  use densityResponse_class, only : densityResponse
+  use particle_class,        only : particle, P_NEUTRON, P_PHOTON
+  use dictionary_class,      only : dictionary
+  use nuclearDatabase_inter, only : nuclearDatabase
+  use pFUnit_mod
+
+  implicit none
+
+@testCase
+  type, extends(TestCase) :: test_densityResponse
+    private
+    type(densityResponse) :: response
+  contains
+    procedure :: setUp
+    procedure :: tearDown
+  end type test_densityResponse
+
+
+contains
+
+  !!
+  !! Sets up test_densityResponse object we can use in a number of tests
+  !!
+  subroutine setUp(this)
+    class(test_densityResponse), intent(inout) :: this
+    type(dictionary)                           :: tempDict
+
+  end subroutine setUp
+
+  !!
+  !! Kills test_densityResponse object we can use in a number of tests
+  !!
+  subroutine tearDown(this)
+    class(test_densityResponse), intent(inout) :: this
+
+  end subroutine tearDown
+
+!!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
+!! PROPER TESTS BEGIN HERE
+!!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
+
+  !!
+  !! Test correct behaviour of the response
+  !!
+@Test
+  subroutine densityResponseing(this)
+    class(test_densityResponse), intent(inout) :: this
+    type(particle)                             :: p
+    class(nuclearDatabase),pointer             :: xsData
+    real(defReal)                              :: res
+
+    ! Test neutron density with different particle energies
+    p % type = P_NEUTRON
+    p % isMG = .false.
+    p % E = ONE
+    res   = ONE / lightSpeed / sqrt(TWO * p % E / neutronMass)
+    @assertEqual(res, this % response % get(p, xsData), res*1.0E-9_defReal)
+
+    p % E = 1.6e-06_defReal
+    res   = ONE / lightSpeed / sqrt(TWO * p % E / neutronMass)
+    @assertEqual(res, this % response % get(p, xsData), res*1.0E-9_defReal)
+
+    ! Test photon density
+    p % type = P_PHOTON
+    res   = ONE/lightSpeed
+    @assertEqual(res, this % response % get(p, xsData), res*1.0E-9_defReal)
+
+  end subroutine densityResponseing
+
+end module densityResponse_test

Tallies/TallyResponses/Tests/fluxResponse_test.f90

@@ -43,7 +43,7 @@ end subroutine tearDown
 !!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
 
   !!
-  !! Test correct behaviour of the filter
+  !! Test correct behaviour of the response
   !!
 @Test
   subroutine fluxResponseing(this)

Tallies/TallyResponses/densityResponse_class.f90

@@ -0,0 +1,77 @@
+module densityResponse_class
+
+  use numPrecision
+  use universalVariables,  only : neutronMass, lightSpeed
+  use dictionary_class,    only : dictionary
+  use particle_class,      only : particle, P_NEUTRON, P_PHOTON
+  use tallyResponse_inter, only : tallyResponse
+
+  ! Nuclear Data interface
+  use nuclearDatabase_inter, only : nuclearDatabase
+
+  implicit none
+  private
+
+  !!
+  !! tallyResponse to score particle density contribution
+  !!
+  !! Returns the inverse of the particle velocity in [cm/s]
+  !!
+  !! NOTE:
+  !!  The velocities are computed from non-relativistic formula for massive particles.
+  !!  The small error might appear in MeV range (e.g. for fusion applications)
+  !!
+  !! Interface:
+  !!   tallyResponse Interface
+  !!
+  type, public,extends(tallyResponse) :: densityResponse
+    private
+  contains
+    procedure :: init
+    procedure :: get
+    procedure :: kill
+
+  end type densityResponse
+
+contains
+
+  !!
+  !! Initialise Response from dictionary
+  !!
+  !! See tallyResponse_inter for details
+  !!
+  subroutine init(self, dict)
+    class(densityResponse), intent(inout) :: self
+    class(dictionary), intent(in)         :: dict
+
+    ! Do nothing
+
+  end subroutine init
+
+  !!
+  !! Returns the inverse of the particle velocity (response to score particle density)
+  !!
+  !! See tallyResponse_inter for details
+  !!
+  function get(self, p, xsData) result(val)
+    class(densityResponse), intent(in)    :: self
+    class(particle), intent(in)           :: p
+    class(nuclearDatabase), intent(inout) :: xsData
+    real(defReal)                         :: val
+
+    ! Gets the particle velocity from the particle
+    val = ONE / p % getVelocity()
+
+  end function get
+
+  !!
+  !! Return to uninitialised State
+  !!
+  elemental subroutine kill(self)
+    class(densityResponse), intent(inout) :: self
+
+    ! Do nothing for nothing can be done
+
+  end subroutine kill
+
+end module densityResponse_class

Tallies/TallyResponses/tallyResponseFactory_func.f90

@@ -12,6 +12,7 @@ module tallyResponseFactory_func
   use macroResponse_class,    only : macroResponse
   use microResponse_class,    only : microResponse
   use weightResponse_class,   only : weightResponse
+  use densityResponse_class,  only : densityResponse
   use testResponse_class,     only : testResponse
 
   implicit none
@@ -23,10 +24,11 @@ module tallyResponseFactory_func
   ! It is printed if type was unrecognised
   ! NOTE:
   ! For now  it is necessary to adjust trailing blanks so all enteries have the same length
-  character(nameLen),dimension(*),parameter :: AVALIBLE_tallyResponses = ['fluxResponse  ',&
-                                                                          'macroResponse ',&
-                                                                          'microResponse ',&
-                                                                          'weightResponse']
+  character(nameLen),dimension(*),parameter :: AVALIBLE_tallyResponses = ['fluxResponse   ',&
+                                                                          'macroResponse  ',&
+                                                                          'microResponse  ',&
+                                                                          'weightResponse ',&
+                                                                          'densityResponse']
 
 contains
 
@@ -61,6 +63,9 @@ subroutine new_tallyResponse(new,dict)
       case('weightResponse')
         allocate(weightResponse :: new)
 
+      case('densityResponse')
+        allocate(densityResponse :: new)
+
       case('testResponse')
         allocate(testResponse :: new)
 

docs/User Manual.rst

@@ -1062,6 +1062,15 @@ Example: ::
       collision_estimator { type collisionClerk; response (flux); flux { type fluxResponse; } }
       }
 
+* densityResponse: used to calculate the particle desnsity, i.e., the response function is 
+  the inverse of the particle velocity in [cm/s]
+
+Example: ::
+
+      tally {
+      collision_estimator { type collisionClerk; response (dens); dens { type densityResponse; } }
+      }
+
 * macroResponse: used to score macroscopic reaction rates
 
   - MT: MT number of the desired reaction. The options are: -1 total, -2 capture,