Moving calculating velocity to the particle

valeriaRaffuzzi · May 10, 2024 · 8fba015 · 8fba015
1 parent 9b57315
commit 8fba015
Show file tree

Hide file tree

Showing 6 changed files with 74 additions and 45 deletions.
diff --git a/ParticleObjects/Tests/particle_test.f90 b/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
 
   !!

diff --git a/ParticleObjects/particle_class.f90 b/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
 
   !!
@@ -129,7 +130,7 @@ module particle_class
     generic              :: build => buildCE, buildMG
     generic              :: assignment(=) => particle_fromParticleState
 
-    ! Inquiry about coordinates
+    ! Enquiry about coordinates
     procedure                  :: rLocal
     procedure                  :: rGlobal
     procedure                  :: dirLocal
@@ -140,14 +141,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
@@ -427,6 +431,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
 !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>

diff --git a/SharedModules/genericProcedures.f90 b/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
 

diff --git a/SharedModules/universalVariables.f90 b/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 (m*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
 

diff --git a/Tallies/TallyResponses/Tests/densityResponse_test.f90 b/Tallies/TallyResponses/Tests/densityResponse_test.f90
@@ -1,7 +1,7 @@
 module densityResponse_test
 
   use numPrecision
-  use universalVariables,    only : lightSpeed
+  use universalVariables,    only : neutronMass, lightSpeed
   use densityResponse_class, only : densityResponse
   use particle_class,        only : particle, P_NEUTRON, P_PHOTON
   use dictionary_class,      only : dictionary
@@ -55,24 +55,20 @@ subroutine densityResponseing(this)
 
     ! Test neutron density with different particle energies
     p % type = P_NEUTRON
+    p % isMG = .false.
     p % E = ONE
-    res   = ONE/1.3831592645e+09_defReal
+    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/1.7495734571e+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)
 
-    ! Test response for unknown particle type
-    p % type = 345
-    res = ZERO
-    @assertEqual(res, this % response % get(p, xsData))
-
   end subroutine densityResponseing
 
 end module densityResponse_test
diff --git a/Tallies/TallyResponses/densityResponse_class.f90 b/Tallies/TallyResponses/densityResponse_class.f90
@@ -17,8 +17,9 @@ module densityResponse_class
   !!
   !! Returns the inverse of the particle velocity in [cm/s]
   !!
-  !! The velocity is calculated from the particle energy for neutrons, and it is
-  !! the speed of light for photons
+  !! 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
@@ -29,6 +30,7 @@ module densityResponse_class
     procedure :: init
     procedure :: get
     procedure :: kill
+
   end type densityResponse
 
 contains
@@ -47,12 +49,7 @@ subroutine init(self, dict)
   end subroutine init
 
   !!
-  !! Calculates the particle velocity for neutron and photons
-  !! NOTE: neutronMass: [MeV]
-  !!       lightSpeed:  [cm/s]
-  !!
-  !! The function returns the inverse of the velocity (response to score particle density)
-  !! if the particle type is neutron or photon, and zero otherwise
+  !! Returns the inverse of the particle velocity (response to score particle density)
   !!
   !! See tallyResponse_inter for details
   !!
@@ -61,25 +58,9 @@ function get(self, p, xsData) result(val)
     class(particle), intent(in)           :: p
     class(nuclearDatabase), intent(inout) :: xsData
     real(defReal)                         :: val
-    real(defReal)                         :: velocity
-
-    ! Initialise response
-    val = ZERO
-
-    ! Calculates the velocity for the relevant particle [cm/s]
-    if (p % type == P_NEUTRON) then
-      velocity = sqrt(TWO * p % E / neutronMass) * lightSpeed
-
-    elseif (p % type == P_PHOTON) then
-      velocity = lightSpeed
-
-    else
-      return
-
-    end if
 
-    ! Returns the inverse of the velocity
-    val = ONE / velocity
+    ! Gets the particle velocity from the particle
+    val = ONE / p % getVelocity()
 
   end function get