Merge pull request #3 from valeriaRaffuzzi/pf-update

Conflict resolution #2

.gitignore

@@ -3,6 +3,10 @@ __pycache__/
 cream.egg-info/
 *.pyc
 
+# Ignore build folder
+Build
+build
+
 # Ignore all hidden files (except gitignore)
 .*
-!/.gitignore
+!/.gitignore

.readthedocs.yaml

@@ -1,10 +1,14 @@
 
 version: 2
 
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.11"
+
 sphinx:
   configuration: docs/conf.py
 
 python:
-  version: 3.7
   install:
-    - requirements: docs/requirements-rtd.txt
+    - requirements: docs/requirements-rtd.txt

CollisionOperator/CollisionProcessors/CMakeLists.txt

@@ -2,5 +2,6 @@
 add_sources( ./collisionProcessor_inter.f90 
              ./collisionProcessorFactory_func.f90 
              ./neutronCEstd_class.f90
-	     ./neutronCEimp_class.f90
-             ./neutronMGstd_class.f90) 
+             ./neutronCEimp_class.f90
+             ./neutronMGstd_class.f90
+             ./neutronMGimp_class.f90)

CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90

@@ -11,20 +11,21 @@ module collisionProcessorFactory_func
   use neutronCEstd_class, only : neutronCEstd
   use neutronCEimp_class, only : neutronCEimp
   use neutronMGstd_class, only : neutronMGstd
+  use neutronMGimp_class, only : neutronMGimp
 
   implicit none
   private
 
   public :: new_collisionProcessor
 
-  ! *** ADD NAME OF A NEW COLLISION PROCESSOR HERE ***!
   ! List that contains all accaptable types of collisionProcessors
   ! 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_collisionProcessors = [ 'neutronCEstd',&
                                                                                 'neutronCEimp',&
-                                                                                'neutronMGstd']
+                                                                                'neutronMGstd',&
+                                                                                'neutronMGimp']
 
 contains
 
@@ -45,31 +46,29 @@ subroutine new_collisionProcessor(new,dict)
     call dict % get(type,'type')
 
     ! Allocate approperiate subclass of collisionProcessor
-    ! *** ADD CASE STATEMENT FOR A NEW COLLISION PROCESSOR BELOW ***!
     select case(type)
       case('neutronCEstd')
         allocate(neutronCEstd :: new)
-        call new % init(dict)
 
       case('neutronCEimp')
         allocate(neutronCEimp :: new)
-        call new % init(dict)
 
       case('neutronMGstd')
         allocate(neutronMGstd :: new)
+
+      case('neutronMGimp')
+        allocate(neutronMGimp :: new)
         call new % init(dict)
 
-     !*** NEW COLLISION PROCESSOR TEMPLATE ***!
-     !case('<newcollisionProcessorName>')
-     !  allocate(<newcollisionProcessorName> :: new)
-     !  call new % init(dict)
-     !
       case default
         print *, AVALIBLE_collisionProcessors
         call fatalError(Here, 'Unrecognised type of collisionProcessor: ' // trim(type))
 
     end select
 
+    ! Initialise new processor
+    call new % init(dict)
+
   end subroutine new_collisionProcessor
 
 end module collisionProcessorFactory_func

CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90

@@ -26,6 +26,7 @@ module collisionProcessor_inter
     real(defReal)     :: muL     = ONE  !! Cosine of deflection angle in LAB-frame
     real(defReal)     :: A       = ZERO !! Target Mass [Neutron Mass]
     real(defReal)     :: kT      = ZERO !! Target temperature [MeV]
+    real(defReal)     :: E       = ZERO !! Collision energy (could be relative to target) [MeV]
   end type
 
 
@@ -86,13 +87,15 @@ module collisionProcessor_inter
     !! Procedure interface for all customisable actions associated with
     !! processing of sollision event (scatter, fission etc.)
     !!
-    subroutine collisionAction(self, p, collDat, thisCycle, nextCycle)
+    subroutine collisionAction(self, p, tally, collDat, thisCycle, nextCycle)
       import :: collisionProcessor, &
                 collisionData, &
+                tallyAdmin, &
                 particle,&
                 particleDungeon
       class(collisionProcessor), intent(inout) :: self
       class(particle), intent(inout)           :: p
+      type(tallyAdmin), intent(inout)          :: tally
       type(collisionData), intent(inout)       :: collDat
       class(particleDungeon),intent(inout)     :: thisCycle
       class(particleDungeon),intent(inout)     :: nextCycle
@@ -104,44 +107,56 @@ end subroutine collisionAction
   !!
   !! Generic flow of collision processing
   !!
-  subroutine collide(self, p, tally ,thisCycle, nextCycle)
+  subroutine collide(self, p, tally, thisCycle, nextCycle)
     class(collisionProcessor), intent(inout) :: self
     class(particle), intent(inout)           :: p
     type(tallyAdmin), intent(inout)          :: tally
     class(particleDungeon),intent(inout)     :: thisCycle
     class(particleDungeon),intent(inout)     :: nextCycle
     type(collisionData)                      :: collDat
-    character(100),parameter                 :: Here = ' collide (collisionProcessor.f90)'
+    logical(defBool)                         :: virtual
+    integer(shortInt)                        :: addCollision
+    character(100),parameter                 :: Here = 'collide (collisionProcessor.f90)'
 
     ! Load material index into data package
     collDat % matIdx = p % matIdx()
 
+    ! Choose collision nuclide and general type (Scatter, Capture or Fission)
+    call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle)
+
+    ! In case of a TMS rejection, set collision as virtual
+    if (collDat % MT == noInteraction) then
+      virtual = .true.
+      addCollision = 0
+    else
+      virtual = .false.
+      addCollision = 1
+    end if
+
     ! Report in-collision & save pre-collison state
     ! Note: the ordering must not be changed between feeding the particle to the tally
-    ! and updating the particle's preCollision state, otherwise this may cause certain 
+    ! and updating the particle's preCollision state, otherwise this may cause certain
     ! tallies (e.g., collisionProbability) to return dubious results
-    call tally % reportInColl(p)
-    call p % savePreCollision()
+    call tally % reportInColl(p, virtual)
 
-    ! Choose collision nuclide and general type (Scatter, Capture or Fission)
-    call self % sampleCollision(p, collDat, thisCycle, nextCycle)
+    call p % savePreCollision()
 
     ! Perform implicit treatment
-    call self % implicit(p, collDat, thisCycle, nextCycle)
+    if (collDat % MT /= noInteraction) call self % implicit(p, tally, collDat, thisCycle, nextCycle)
 
     ! Select physics to be processed based on MT number
     select case(collDat % MT)
       case(N_N_elastic, macroAllScatter)
-        call self % elastic(p, collDat, thisCycle, nextCycle)
+        call self % elastic(p, tally, collDat, thisCycle, nextCycle)
 
       case(N_N_inelastic, macroIEScatter)
-        call self % inelastic(p, collDat, thisCycle, nextCycle)
+        call self % inelastic(p, tally, collDat, thisCycle, nextCycle)
 
       case(N_DISAP, macroCapture)
-        call self % capture(p, collDat, thisCycle, nextCycle)
+        call self % capture(p, tally, collDat, thisCycle, nextCycle)
 
       case(N_FISSION, macroFission)
-        call self % fission(p, collDat, thisCycle, nextCycle)
+        call self % fission(p, tally, collDat, thisCycle, nextCycle)
 
       case(noInteraction)
         ! Do nothing
@@ -152,13 +167,16 @@ subroutine collide(self, p, tally ,thisCycle, nextCycle)
     end select
 
     ! Apply post collision implicit treatments
-    call self % cutoffs(p, collDat, thisCycle, nextCycle)
+    call self % cutoffs(p, tally, collDat, thisCycle, nextCycle)
+
+    ! Update particle collision counter
+    p % collisionN = p % collisionN + addCollision
 
     ! Report out-of-collision
     call tally % reportOutColl(p, collDat % MT, collDat % muL)
 
     ! Report end-of-history if particle was killed
-    if( p % isDead) then
+    if (p % isDead) then
       p % fate = ABS_FATE
       call tally % reportHist(p)
     end if