DT majorant #106

valeriaRaffuzzi · 2024-01-17T18:55:07Z

This PR adds the precalculation of the majorant xss on a unionised grid for delta tracking.

Few things to notice:

I decided to have the transportOperator call the initialisation procedure in the database, so that only DT and HT call it and ST doesn't
Precalculated majorant is the new default and only option available when using DT
This way of constructing the grid is the easiest, but it's probably not the fastest. For instance, constructing the unionised energy grid (which is not really a parallelisable thing) for the CE c5g7 takes ~8 seconds (while calculating the majorant on that grid takes a fraction of a second). But then transport is x2 times faster with DT. If you have a better way in mind let me know! For now I wanted this done fast.
This was tested and it works fine, and the relevant integrationTests were adapted too to include the change.

Mikolaj-A-Kowalski

Thanks for the changes! They look great! I suppose that the optimisation of the unionised grid construction we wish to add for later (and open an issue about it)? Unless do you prefer to try to trim some fat in this PR?

The only major point is what we discussed earlier :-) My instinct would be to have a unionisation as a option in the nuclearDatabase configuration dictionary. Then initMajorant would be called (or not) in the activate procedure. Also the updateMajorantXS would decide between old and the unionised grid search based on the allocation status of the self % majorant.

Part of my fear is that with initMajorant and activate being separate functions weird things may be happening (like the example I outlined in the code attached comment) in case of some sequences of calls and we would need to add checks to protect against those situations (which is not difficult, but there is always a danger of missing some nasty, specific call sequence).

Also I might have missed something but does the unionised grid for majorant require that URR is switched off for all nuclides?

Mikolaj-A-Kowalski · 2024-01-18T13:56:17Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

@@ -3,7 +3,8 @@ module aceNeutronDatabase_class
  use numPrecision
  use endfConstants
  use universalVariables
-  use genericProcedures, only : fatalError, numToChar
+  use genericProcedures, only : fatalError, numToChar, concatenate, quickSort, &


I am guilty to keep forgetting about #37. Maybe we move the import of fatalError to use errors_mod, only : fatalError?

Mikolaj-A-Kowalski · 2024-01-18T13:58:53Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

@@ -60,6 +61,8 @@ module aceNeutronDatabase_class
  !!   nuclides   -> array of aceNeutronNuclides with data
  !!   materials  -> array of ceNeutronMaterials with data
  !!   Ebounds    -> array with bottom (1) and top (2) energy bound
+  !!   majorant   -> unionised majorant cross section
+  !!   eGrid      -> unionised energy grid


Should we change the name to make it more explicit that it is the unionised grid?
To e.g. eGridUnion? [Don't know to be honest?]

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

Mikolaj-A-Kowalski · 2024-01-18T14:09:20Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

    class(aceNeutronDatabase), intent(in) :: self
    real(defReal), intent(in)             :: E
    class(RNG), intent(inout)             :: rand
-    integer(shortInt)                     :: i, matIdx
+    integer(shortInt)                     :: idx
+    real(defReal)                         :: f
+    character(100), parameter :: Here = 'updateMajorantXS (aceNeutronDatabase_class.f90)'

    associate (maj => cache_majorantCache(1) )
      maj % E  = E
-      maj % xs = ZERO

-      do i = 1, size(self % activeMat)
-        matIdx = self % activeMat(i)
+      idx = binarySearch(self % eGrid, E)

-        ! Update if needed
-        if( cache_materialCache(matIdx) % E_tot /= E) then
-          call self % updateTotalMatXS(E, matIdx, rand)
-        end if
+      if(idx <= 0) then
+        call fatalError(Here,'Failed to find energy: '//numToChar(E)//&
+                             ' in unionised majorant grid')
+      end if
+
+      associate(E_top => self % eGrid(idx + 1), E_low  => self % eGrid(idx))
+        f = (E - E_low) / (E_top - E_low)
+      end associate
+
+      maj % xs = self % majorant(idx+1) * f + (ONE - f) * self % majorant(idx)


I have a fear that if it happens that the initMajorant is not called by the user, but for some reason a call to updateMajorantXS is made, then the code will segment on the binary search.

Similarly in the case a following calls are made:

activate with materials 1,2,3

init majorant

2nd activate (reactivate) with materials 1,2,3,4
Then it is possible that the majorant will find itself in an invalid state (not representing the expected majorant).
What we do with this will depend on what we finally decide on having initMajorant as a separate procedure in the interface, but it may involve some allocated(self % eGrid) checks to give more descriptive error message)

Mikolaj-A-Kowalski · 2024-01-18T14:22:25Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

+    do i = 1, size(self % activeMat)
+
+      ! Get current material index and number of nuclides in that material
+      matIdx = self % activeMat(i)
+      nNuc = size(self % materials(matIdx) % nuclides)
+
+      ! Loop over nuclides present in that material
+      do j = 1, nNuc
+
+        ! Get index and check if it's already been added to the grid
+        nucIdx = self % materials(matIdx) % nuclides(j)
+        isDone = nucMap % getOrDefault(nucIdx, NOT_PRESENT)
+
+        ! If it's a new nuclide, add its energy grid to the unionised one
+        if (isDone /= IN_SET) then
+
+          tmpGrid = concatenate(tmpGrid, self % nuclides(nucIdx) % eGrid)
+
+          ! Sort and remove duplicates
+          call quickSort(tmpGrid)
+          tmpGrid = removeDuplicatesSorted(tmpGrid)
+
+          ! Add nuclide to the set
+          call nucMap % add(nucIdx, IN_SET)
+
+        end if
+
+      end do
+    end do


It is perhaps not at this moment but the (relatively) low hanging fruits for optimisation here would be to create a real version of the dynArray_class and use it as the tmpGrid. I guess it could take care of a large portion of the runtime by avoiding a lot of reallocation.

Of course dynArray would than also have in-place versions of removeDuplicates and quickSort added.

It is a bit involved, so I will just open an issue. It will be needed for actual use with burnup or so.

Even though more than 70% of the time taken is spent sorting stuff, and this is unavoidable I think.

Mikolaj-A-Kowalski · 2024-01-18T14:25:51Z

NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90


-    ! Precalculate majorant xs for delta tracking
+    if (loud) print '(A)', 'Building unionised majorant cross section'


Suggested change

if (loud) print '(A)', 'Building unionised majorant cross section'

if (loud) print '(A)', 'Building MG unionised majorant cross section'

Maybe we should highlight that it is a MG case?

TransportOperator/transportOperatorHT_class.f90

ChasingNeutrons · 2024-01-18T14:58:11Z

I agree with Mikolaj's comments. Related to that, I think, until the majorant construction is optimised, it is desirable to be able to switch off precomputing a majorant in the input file - otherwise we might encounter a nasty case with many nuclides that takes impractically long to run just from constructing the grid which would run just fine at the moment. It might require switching off majorant construction and recompiling.

Add forgotten docs

…into DTmajorant

valeriaRaffuzzi · 2024-01-24T22:11:17Z

Issue #107 is addressed here. The unionised energy grid is build gradually from the nuclides grids instead of concatenating arrays and using quicksort. Now it takes a fraction of a second to precompute everything instead of ~10s.

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

ChasingNeutrons · 2024-01-31T15:34:45Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

@@ -51,19 +52,26 @@ module aceNeutronDatabase_class
  !! It's possible to use probability tables in the unresolved resonance range if
  !! ures is included in the input file
  !!
+  !! NOTE: the unionised majorant is not calculated and used if probability tables 


I think we should open an issue on this. It will be very desirable to still have a unionised majorant, e.g., if we want to accurately run fast problems. I understand it shouldn't have a big impact for thermal problems, even with burnup.

Sure, I'll open an issue. I am not sure if one can a priori find out what's the majorant when using probability tables. One approach could be to sample, say, 1000 times per each energy when calculating the majorant, and retaining the highest xs.

ChasingNeutrons · 2024-01-31T15:36:24Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

  !!   activeMat  -> array of materials present in the geometry
  !!   nucToZaid  -> map to link nuclide index to zaid index
  !!   hasUrr     -> ures probability tables flag, it's false by default
  !!   hasDBRC    -> DBRC flag, it's false by default
+  !!   has majorant -> unionised majorant cross section flag


hasMajorant, not has majorant

ChasingNeutrons · 2024-01-31T15:44:09Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

@@ -653,6 +692,9 @@ subroutine init(self, dict, ptr, silent )
      self % Ebounds(2) = min(self % Ebounds(2), self % nuclides(i) % eGrid(j))
    end do

+    ! Read unionised majorant flag
+    call dict % getOrDefault(self % hasMajorant, 'majorant', .false.)


I would be inclined to set this to true by default. We run DT often and it's one less flag to set in that case, because we would always want the precomputed majorant (as we used to do in MG). It's also a very small price to pay if we happen to be running ST.

ChasingNeutrons · 2024-01-31T15:47:38Z

NuclearData/nuclearDatabase_inter.f90

@@ -196,7 +197,7 @@ end function getMaterial
    !! Allows to retrive an access to nuclide data for all databases types


While we're here.... retrieve*

ChasingNeutrons · 2024-01-31T15:53:02Z

TransportOperator/transportOperatorHT_class.f90

@@ -98,7 +89,7 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle)
    if (abs(majorant_inv) > huge(majorant_inv)) call fatalError(Here, "Majorant is 0")

    DTLoop:do
-      distance = -log( p% pRNG % get() ) * majorant_inv
+      distance = -log( p % pRNG % get() ) * majorant_inv

      ! Move partice in the geometry


ChasingNeutrons · 2024-01-31T15:55:07Z

docs/User Manual.rst

@@ -742,7 +742,9 @@ from ACE files.
  resonance probability tables treatment
 * DBRC (*optional*, default = no DBRC): list of ZAIDs of nuclides for which DBRC has
  to be applied.
-
+* majorant (*optional*, default = 0): 1 for true; 0 for false; flag to activate the 


As above, I think this should default to true. I suppose it comes down to which would we prefer to suffer in a worst case - a slower calculation or more memory. I think the memory is more tolerable.

…into DTmajorant

valeriaRaffuzzi · 2024-02-05T17:42:05Z

Yet more changes to this PR:

now the majorant can be precalculated also with ures tables (solves issue Pre calculating majorant xs with URR #114 ): the highest entry is the ures tables at each energy is also precomputed and used to calculated the overall majorant
the majorant is nudged up by 1e-6 to avoid small floating point errors or inconsistencies that seemed to appear (e.g., majorant being smaller than the largest material xs by 1e-15)
I added to the unionised grid the energy points corresponding to the limit energies of URR and S(a,b) tables. This is because I was getting a majorant < mat_xs just below the S(a,b) upper energy boundary: two consecutive energies were one within and one outside the S(a,b) range, so interpolating gave the wrong result. Honestly nothing breaks if this happens, because it's not very frequent. Stuff like these are hard to avoid when there are discontinuities in the xss like due to Sab or ures, and they could still happen. However now for the cases I tried it behaves ok.

ChasingNeutrons

I feel bad that my only comments are tiny typos... But this looks great! We're a good bit more respectable after this!

ChasingNeutrons · 2024-02-05T18:34:07Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

+
+            needsUrr = (nuc % hasProbTab .and. E >= nuc % urrE(1) .and. E <= nuc % urrE(2))
+
+            ! Check if present nuclide uses URR tabes


ChasingNeutrons · 2024-02-05T18:35:38Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90

@@ -741,14 +742,14 @@ subroutine init(self, ACE, nucIdx, database)
    ! Create a stack of MT reactions, devide them into ones that produce 2nd-ary


ChasingNeutrons · 2024-02-05T18:35:52Z

NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90

@@ -741,14 +742,14 @@ subroutine init(self, ACE, nucIdx, database)
    ! Create a stack of MT reactions, devide them into ones that produce 2nd-ary
    ! particlues and pure absorbtion


particles*, absorption*

ChasingNeutrons · 2024-02-05T18:39:55Z

NuclearData/nuclearDatabase_inter.f90

@@ -217,15 +219,15 @@ end function getNuclide
    !!
    !! Return a pointer to a reaction
    !!
-    !! Allows to retrive an access to reaction data for all databases types
+    !! Allows to retrieve an access to reaction data for all databases types
    !! Reactions can be associated either with nuclides or materials. Thus, there is ambiguity
    !! whether material or nuclide should be asked to provide reaction data (using matIdx or nuIdx)


valeriaRaffuzzi · 2024-02-05T19:17:41Z

I feel bad that my only comments are tiny typos... But this looks great! We're a good bit more respectable after this!

Well, this is the result of several iterations already! Luckily we converged :)

Mikolaj-A-Kowalski

Looks good! Thanks for all the hard work. The URR added some complication to this one ;-)

valeriaRaffuzzi · 2024-02-06T14:51:43Z

Looks good! Thanks for all the hard work. The URR added some complication to this one ;-)

Thanks! Yes, accounting for URR was tricky, and I did spy on Serpent for that!

But interestingly the worst interpolation errors were around the S(a,b) limit energies rather than in the URR range, in the problem I tried. Hopefully now it's more robust to all of that!

V. Raffuzzi(vr339) added 2 commits January 17, 2024 18:20

Adding majorant precalculation for DT

35c9503

Adjusting intTests adding majorant precalculation

435aab7

Mikolaj-A-Kowalski requested changes Jan 18, 2024

View reviewed changes

V. Raffuzzi(vr339) and others added 4 commits January 18, 2024 16:52

Addressing comments

18fbe55

Update aceNeutronDatabase_class.f90

66efce6

Add forgotten docs

Speeding up unionised majorant cross section pre calculation

8c82ffb

Merge branch 'DTmajorant' of https://github.com/valeriaRaffuzzi/SCONE …

acaebe0

…into DTmajorant

Merge branch 'main' into DTmajorant

aa804fa

valeriaRaffuzzi commented Jan 29, 2024

View reviewed changes

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 Outdated Show resolved Hide resolved

Remove unused procedures in aceNeutronDatabase

cb9a623

ChasingNeutrons reviewed Jan 31, 2024

View reviewed changes

V. Raffuzzi(vr339) added 4 commits January 31, 2024 23:25

Typos and small fixes to DT majorant

0072e91

Merge branch 'DTmajorant' of https://github.com/valeriaRaffuzzi/SCONE …

efca20d

…into DTmajorant

Add routines to calculate majorant also with ures on

04a3774

Typos and preventing out of bound array access

b97e613

ChasingNeutrons approved these changes Feb 5, 2024

View reviewed changes

Typos

a651001

valeriaRaffuzzi requested a review from Mikolaj-A-Kowalski February 5, 2024 19:18

Mikolaj-A-Kowalski approved these changes Feb 6, 2024

View reviewed changes

valeriaRaffuzzi merged commit 4757437 into CambridgeNuclear:main Feb 6, 2024
5 checks passed

valeriaRaffuzzi deleted the DTmajorant branch February 6, 2024 14:51

This was referenced Feb 6, 2024

Pre calculating majorant xs with URR #114

Closed

Optimise unionised energy grid construction #107

Closed

Precomputing the majorant + unionised grids with S(alpha,beta) #28

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

DT majorant #106

DT majorant #106

valeriaRaffuzzi commented Jan 17, 2024 •

edited

Loading

Mikolaj-A-Kowalski left a comment

Mikolaj-A-Kowalski Jan 18, 2024

Mikolaj-A-Kowalski Jan 18, 2024

Mikolaj-A-Kowalski Jan 18, 2024

Mikolaj-A-Kowalski Jan 18, 2024

valeriaRaffuzzi Jan 18, 2024 •

edited

Loading

Mikolaj-A-Kowalski Jan 18, 2024

ChasingNeutrons commented Jan 18, 2024

valeriaRaffuzzi commented Jan 24, 2024

ChasingNeutrons Jan 31, 2024

valeriaRaffuzzi Jan 31, 2024

ChasingNeutrons Jan 31, 2024

ChasingNeutrons Jan 31, 2024 •

edited

Loading

ChasingNeutrons Jan 31, 2024

ChasingNeutrons Jan 31, 2024

ChasingNeutrons Jan 31, 2024

valeriaRaffuzzi commented Feb 5, 2024 •

edited

Loading

ChasingNeutrons left a comment

ChasingNeutrons Feb 5, 2024

ChasingNeutrons Feb 5, 2024

ChasingNeutrons Feb 5, 2024

ChasingNeutrons Feb 5, 2024

valeriaRaffuzzi commented Feb 5, 2024

Mikolaj-A-Kowalski left a comment

valeriaRaffuzzi commented Feb 6, 2024


		! Precalculate majorant xs for delta tracking
		if (loud) print '(A)', 'Building unionised majorant cross section'

	if (loud) print '(A)', 'Building unionised majorant cross section'
	if (loud) print '(A)', 'Building MG unionised majorant cross section'

		@@ -196,7 +197,7 @@ end function getMaterial
		!! Allows to retrive an access to nuclide data for all databases types


		needsUrr = (nuc % hasProbTab .and. E >= nuc % urrE(1) .and. E <= nuc % urrE(2))

		! Check if present nuclide uses URR tabes

		@@ -741,14 +742,14 @@ subroutine init(self, ACE, nucIdx, database)
		! Create a stack of MT reactions, devide them into ones that produce 2nd-ary

		@@ -741,14 +742,14 @@ subroutine init(self, ACE, nucIdx, database)
		! Create a stack of MT reactions, devide them into ones that produce 2nd-ary
		! particlues and pure absorbtion

DT majorant #106

DT majorant #106

Conversation

valeriaRaffuzzi commented Jan 17, 2024 • edited Loading

Mikolaj-A-Kowalski left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

valeriaRaffuzzi Jan 18, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChasingNeutrons commented Jan 18, 2024

valeriaRaffuzzi commented Jan 24, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChasingNeutrons Jan 31, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

valeriaRaffuzzi commented Feb 5, 2024 • edited Loading

ChasingNeutrons left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

valeriaRaffuzzi commented Feb 5, 2024

Mikolaj-A-Kowalski left a comment

Choose a reason for hiding this comment

valeriaRaffuzzi commented Feb 6, 2024

valeriaRaffuzzi commented Jan 17, 2024 •

edited

Loading

valeriaRaffuzzi Jan 18, 2024 •

edited

Loading

ChasingNeutrons Jan 31, 2024 •

edited

Loading

valeriaRaffuzzi commented Feb 5, 2024 •

edited

Loading