From 5e66d1c2e66ba83039ea75873440d9924931cbb1 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Mon, 23 Mar 2026 18:21:52 +0000 Subject: [PATCH 1/2] Preventing failure in delayed neutron emission --- .../endfTable/endfTable_class.f90 | 24 ++++++++++++++++++- .../fissionCE_class.f90 | 15 +++++++++--- .../releaseLawENDF/constantRelease_class.f90 | 14 +++++++++++ .../polynomialRelease_class.f90 | 14 +++++++++++ .../releaseLawENDF/releaseLawENDF_inter.f90 | 13 ++++++++++ .../releaseLawENDFslot_class.f90 | 13 ++++++++++ .../releaseLawENDF/tabularRelease_class.f90 | 13 ++++++++++ 7 files changed, 102 insertions(+), 4 deletions(-) diff --git a/NuclearData/NuclearDataStructures/endfTable/endfTable_class.f90 b/NuclearData/NuclearDataStructures/endfTable/endfTable_class.f90 index 59dce770e..537dc4eda 100644 --- a/NuclearData/NuclearDataStructures/endfTable/endfTable_class.f90 +++ b/NuclearData/NuclearDataStructures/endfTable/endfTable_class.f90 @@ -36,7 +36,8 @@ module endfTable_class !! Interface: !! at -> return y-value given x-value !! init -> initialise from components - !! intagral -> Integrate table from x_min to multiple values x + !! integral -> Integrate table from x_min to multiple values x + !! hasX -> Check on whether a given x value is within grid bounds !! reloadY -> Change values on y-axis keeping x & interpolation unchanged !! kill -> return to uninitalised state !! @@ -54,6 +55,7 @@ module endfTable_class procedure :: at procedure :: integral procedure :: reloadY + procedure :: hasX generic :: init => initSimple, initInter procedure :: kill @@ -195,6 +197,26 @@ function at(self, x) result (y) end function at + !! + !! Returns whether a given x value is within bounds. + !! Made to help deal with inconsistent nuclear data. Allows checking that a grid + !! extends up to a certain value before attempting to get a y-value in that region. + !! + function hasX(self, x) result(has) + class(endfTable), intent(in) :: self + real(defReal), intent(in) :: x + logical(defBool) :: has + integer(shortInt) :: x_idx + + x_idx = floorSearch(self % x, x) + if (x_idx < 0) then + has = .false. + else + has = .true. + end if + + end function hasX + !! !! Change values in y-axis without any othe change !! diff --git a/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 b/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 index c782f9672..26fffa57a 100644 --- a/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 +++ b/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 @@ -210,6 +210,9 @@ end function releasePrompt !! !! Returns number of particles produced on average by the reaction with delay !! + !! Has a check on energy because some nuclides in some libraries (e.g., U236 in ENDFB7) + !! have inconsistent energy grids for prompt and delayed emission!! + !! !! See uncorrelatedReactionCE for details !! function releaseDelayed(self, E) result(N) @@ -217,12 +220,18 @@ function releaseDelayed(self, E) result(N) real(defReal), intent(in) :: E real(defReal) :: N - if(allocated(self % nuBarDelayed)) then - N = self % nuBarDelayed % releaseAt(E) - else + if (.not. allocated(self % nuBarDelayed)) then + N = ZERO + return + end if + + if (.not. self % nuBarDelayed % hasEnergy(E)) then N = ZERO + return end if + N = self % nuBarDelayed % releaseAt(E) + end function releaseDelayed !! diff --git a/NuclearData/emissionENDF/releaseLawENDF/constantRelease_class.f90 b/NuclearData/emissionENDF/releaseLawENDF/constantRelease_class.f90 index 3c761ba2a..044de6cff 100644 --- a/NuclearData/emissionENDF/releaseLawENDF/constantRelease_class.f90 +++ b/NuclearData/emissionENDF/releaseLawENDF/constantRelease_class.f90 @@ -20,6 +20,7 @@ module constantRelease_class contains procedure :: init procedure :: releaseAt + procedure :: hasEnergy procedure :: kill end type constantRelease @@ -49,6 +50,19 @@ function releaseAt(self, E_in) result(release) release = self % secondaryRelease end function releaseAt + + !! + !! Return whether a given energy is available with the law. + !! Always true + !! + function hasEnergy(self, E_in) result(has) + class(constantRelease), intent(in) :: self + real(defReal), intent(in) :: E_in + logical(defBool) :: has + + has = .true. + + end function hasEnergy !! !! Return to uninitialised state diff --git a/NuclearData/emissionENDF/releaseLawENDF/polynomialRelease_class.f90 b/NuclearData/emissionENDF/releaseLawENDF/polynomialRelease_class.f90 index 3c11aa3fe..8f8fd143f 100644 --- a/NuclearData/emissionENDF/releaseLawENDF/polynomialRelease_class.f90 +++ b/NuclearData/emissionENDF/releaseLawENDF/polynomialRelease_class.f90 @@ -22,6 +22,7 @@ module polynomialrelease_class contains procedure :: init procedure :: releaseAt + procedure :: hasEnergy procedure :: kill end type polynomialRelease @@ -55,6 +56,19 @@ pure function releaseAt(self,E_in) result(release) end do end function releaseAt + + !! + !! Return whether a given energy is available with the law. + !! Always true + !! + function hasEnergy(self, E_in) result(has) + class(polynomialRelease), intent(in) :: self + real(defReal), intent(in) :: E_in + logical(defBool) :: has + + has = .true. + + end function hasEnergy !! !! Return to uninitialised state diff --git a/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDF_inter.f90 b/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDF_inter.f90 index 452b7db6c..0dfe0b745 100644 --- a/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDF_inter.f90 +++ b/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDF_inter.f90 @@ -12,6 +12,7 @@ module releaseLawENDF_inter private contains procedure(releaseAt),deferred :: releaseAt + procedure(hasEnergy),deferred :: hasEnergy procedure(kill),deferred :: kill end type releaseLawENDF @@ -28,6 +29,18 @@ function releaseAt(self,E_in) result(release) real(defReal) :: release end function releaseAt + !! + !! Return whether a given energy is available with the law + !! + function hasEnergy(self, E_in) result(has) + import :: defReal, & + defBool, & + releaseLawENDF + class(releaseLawENDF), intent(in) :: self + real(defReal), intent(in) :: E_in + logical(defBool) :: has + end function hasEnergy + !! !! Return to uninitialised state !! diff --git a/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDFslot_class.f90 b/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDFslot_class.f90 index 6f85b6047..96c406783 100644 --- a/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDFslot_class.f90 +++ b/NuclearData/emissionENDF/releaseLawENDF/releaseLawENDFslot_class.f90 @@ -13,6 +13,7 @@ module releaseLawENDFslot_class contains ! Superclass interface procedure :: releaseAt + procedure :: hasEnergy procedure :: kill ! Define assignment @@ -35,6 +36,18 @@ function releaseAt(self,E_in) result(release) release = self % slot % releaseAt(E_in) end function releaseAt + + !! + !! Returns whether an energy, E_in, is within the table + !! + function hasEnergy(self, E_in) result(has) + class(releaseLawENDFslot), intent(in) :: self + real(defReal), intent(in) :: E_in + logical(defBool) :: has + + has = self % slot % hasEnergy(E_in) + + end function hasEnergy !! !! Return to uninitialised state diff --git a/NuclearData/emissionENDF/releaseLawENDF/tabularRelease_class.f90 b/NuclearData/emissionENDF/releaseLawENDF/tabularRelease_class.f90 index 65f374df0..5cf80b5e6 100644 --- a/NuclearData/emissionENDF/releaseLawENDF/tabularRelease_class.f90 +++ b/NuclearData/emissionENDF/releaseLawENDF/tabularRelease_class.f90 @@ -28,6 +28,7 @@ module tabularRelease_class contains generic :: init => initSimple, initInter procedure :: releaseAt + procedure :: hasEnergy procedure :: kill procedure,private :: initSimple @@ -48,6 +49,18 @@ function releaseAt(self, E_in) result(release) end function releaseAt + !! + !! Returns whether an energy, E_in, is within the table + !! + function hasEnergy(self, E_in) result(has) + class(tabularRelease), intent(in) :: self + real(defReal), intent(in) :: E_in + logical(defBool) :: has + + has = self % releaseTable % hasX(E_in) + + end function hasEnergy + !! !! Return to uninitialised state !! From 250738ef204322b3e036de58588b979fccfbcb67 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Tue, 7 Apr 2026 21:51:02 +0100 Subject: [PATCH 2/2] Fixes to docs --- docs/Dictionary Input.rst | 8 +-- docs/Geometry.rst | 10 +-- docs/Input Manual.rst | 130 +++++++++++++++++++------------------- docs/Installation.rst | 13 +++- docs/index.rst | 6 +- 5 files changed, 87 insertions(+), 80 deletions(-) diff --git a/docs/Dictionary Input.rst b/docs/Dictionary Input.rst index 4df1558ee..982f3b2f2 100644 --- a/docs/Dictionary Input.rst +++ b/docs/Dictionary Input.rst @@ -18,9 +18,9 @@ of a single dictionary. The following content is available: * Subdictionary Hierarchical structure of dictionaries can be loaded from ASCII files and be -used inside the SCONE to construct different objects and set calculation +used inside SCONE to construct different objects and set calculation settings. The syntax for writing dictionaries is based on -`OpenFOAM `_. However only subset of +`OpenFOAM `_. However only a subset of OpenFOAM syntax is supported. An example of the correct input dictionary is:: // This is a line comment. C++ Style @@ -50,7 +50,7 @@ At least one space is needed before an entry name and its value :: intArray ( 1 2 3) - OK intArray(1 2 3) - WRONG -Note that this is not a case for sub-dictionaries :: +Note that this is not the case for sub-dictionaries :: sub{ entry1 1;} - OK sub { entry1 1;} - OK @@ -93,6 +93,6 @@ They have the following meaning:: ::= Real number e.g.: ^[-+]?[0-9]+[.][0-9]*([eE][-+]?[0-9]+)?$ ::= Must contain not digit character e.g.: .*[a-zA-Z]+.* -This is not yet full definition of the grammar as it does not contain limitations related to +This is not yet a full definition of the grammar as it does not contain limitations related to maximum size of integers, reals and maximum length of the characters. It may be useful if somebody will try to use a proper parser to read the SCONE dictionary files. diff --git a/docs/Geometry.rst b/docs/Geometry.rst index 33985d04a..cb48a791d 100644 --- a/docs/Geometry.rst +++ b/docs/Geometry.rst @@ -6,7 +6,7 @@ Geometry Overview '''''''' -A number of terms is going to be introduced in this section (e.g. *surface tolerance*, +A number of terms are going to be introduced in this section (e.g. *surface tolerance*, *surface transparency*). These names are specific to SCONE. However, very similar concepts can be found in other codes albeit under different names. @@ -20,14 +20,14 @@ thus only a brief overview will be given here. In principle, a surface can be represented by an equation :math:`0 = F(r)`. Thus, for any point in a geometry it is possible to evaluate :math:`F(r) = c`. Now *halfspaces* can be defined using the sign of the *remainder* (:math:`c`) of the surface expression. Positive (+ve) halfspace corresponds -to :math:`c > 0` and is considered an *outside* of the surface. Similarly, negative (-ve) halfspace -is associated with :math:`c < 0` and is considered an *inside* of the surface. Clearly each surface +to :math:`c > 0` and is considered *outside* of the surface. Similarly, negative (-ve) halfspace +is associated with :math:`c < 0` and is considered *inside* of the surface. Clearly each surface subdivides the entire space into two halfpaces. Ideally there would be no need to consider a case of :math:`c = 0`, as the probability of finding a randomly placed point exactly at the surface is zero. However, in Monte Carlo transport simulations the movement of the particles is often explicitly -resolved, which means that they need to be temporally stopped at the material interface in order +resolved, which means that they need to be temporarily stopped at the material interface in order to resample their flight distance. As a result, significant care must be taken to ensure -a correct behaviour of the geometry implementation in the vicinity of surfaces (small :math:`c`). +correct behaviour of the geometry implementation in the vicinity of surfaces (small :math:`c`). The binary subdivision of a whole space is rarely sufficient. More practical volumes can be defined with set expressions on the halfspaces. In other words, an almost arbitrary volume can be diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index e527af20a..9db1063e4 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -532,8 +532,15 @@ Example: :: collisionOperator { neutronMG { type neutronMGimp; weightWindows 1; maxSplit 50; } } +Fields +------ + +Fields allow specifying quantities which vary in space and, potentially, energy. These are +specified independently of the geometry. SCONE supports several named fields which modify +the geometry and/or particle transport. + Weight Windows --------------- +************** Weight windows can be used if, inside the collision operator ``neutronCEimp``, the keyword ``weightWindows`` is set to 1. Then, in the input file, one needs to add: :: @@ -560,7 +567,7 @@ Example: :: wUpper (2.0 1.2 1.5 1.1 2.0 4.0); Uniform Fission Sites ---------------------- +********************* Uniform Fission Sites can be used if, inside the collision operator ``neutronCEimp``, the keyword ``UFS`` is set to 1. Then, in the input file, one needs to add: :: @@ -583,6 +590,58 @@ Example: :: map { } } +Temperature and density fields +############################## + +The value of temperature and density at any point in the geometry can be overwritten by imposing +a field. In order to support surface tracking this is done with fields which are picewise constant +and endowed with a distance calculation compute the distance until the value of the field changes. +These piecewise constant fields can be initialised either with an explicit definition or with a path +to the field definition. + +A temperature field is invoked with a dictionary named ``temperature`` while a density field is invoked +with a dictionary named ``density``. The temperature field specifies local temperatures in kelvin while +the density field specifies the local dimensionless factors by which material density should be scaled. + +Currently there is only one available PieceConstantField: + +* cartesianField. This is similar to a latUniverse: the value of the field varies over a regular + Cartesian lattice with a given shape and size. The field also allows specifying different values + in different materials, or uniformly across all materials. + + - shape: (x y z) array of integers, stating the numbers of x, y and z + elements of the field. For a 2D field, the z entry has to be 0. + - pitch: (x y z) array with the x, y and z field pitches. In a 2D field, + the value entered in the third dimension is not used. [cm] + - origin (*optional*, default = (0.0 0.0 0.0)): (x y z) array with the + origin of the field. [cm] + - materials: list of material names, corresponding to materials in nuclearData. + Optionally, ``all`` can be used, applying the values of the field to all materials. + - names of each material: a map, named after every material present in the materials list. + The entries of the map are the values that the field takes in that material in that + element of the field. The order is: increasing x, increasing y and then increasing z. + - default: the value taken by the field when a point is either outside of the field or + in a material which is not included in the field. + +Example: :: + + temperature { type cartesianField; shape (3 2 2); pitch (1.0 1.0 1.5); + materials (uo2 water); + uo2 ( + 901 902 903 + 904 905 906 + 907 908 909 + 910 911 912 ); + water ( + 601 602 603 + 604 605 606 + 607 608 609 + 610 611 612); + default 302; } + + density { type cartesianField; file ./myDensityField; } + + Geometry -------- @@ -593,8 +652,6 @@ A detailed description about the geometry modelling adopted in SCONE can be foun surfaces { } cells { } universes { } - temperature { } - density { } } At the moment, the only **geometry** type available is ``geometryStd``. As for the boundary @@ -630,12 +687,6 @@ Hence, an example of a geometry input could look like: :: For more details about the graph-like structure of the nested geometry see the relevant :ref:`section `. -The geometry optionally allows the use of ``temperature`` and ``density`` fields. These -are super-imposed fields which modify the temperature and densities given in nuclear data. -The temperature field specifies local temperatures in kelvin while the density field -specifies the local dimensionless factors by which material density should be scaled. -Both of these follow the syntax of a ``PieceConstantField``. - Surfaces ######## @@ -925,52 +976,6 @@ Example: :: root { id 1000; type rootUniverse; border 10; fill u<1>; } -PieceConstantFields -################### - -These are fields which are piecewise constant and are endowed with a distance calculation to -compute the distance until the value of the field changes. These can be used for imposing -density and temperature distributions across the system in a convenient manner. Can be initialised -either with an explicit definition or with a path to the field definition. - -Currently there is only one available PieceConstantField: - -* cartesianField. This is similar to a latUniverse: the value of the field varies over a regular - Cartesian lattice with a given shape and size. The field also allows specifying different values - in different materials, or uniformly across all materials. - - - shape: (x y z) array of integers, stating the numbers of x, y and z - elements of the field. For a 2D field, the z entry has to be 0. - - pitch: (x y z) array with the x, y and z field pitches. In a 2D field, - the value entered in the third dimension is not used. [cm] - - origin (*optional*, default = (0.0 0.0 0.0)): (x y z) array with the - origin of the field. [cm] - - materials: list of material names, corresponding to materials in nuclearData. - Optionally, ``all`` can be used, applying the values of the field to all materials. - - names of each material: a map, named after every material present in the materials list. - The entries of the map are the values that the field takes in that material in that - element of the field. The order is: increasing x, increasing y and then increasing z. - - default: the value taken by the field when a point is either outside of the field or - in a material which is not included in the field. - -Example: :: - - temperature { type cartesianField; shape (3 2 2); pitch (1.0 1.0 1.5); - materials (uo2 water); - uo2 ( - 901 902 903 - 904 905 906 - 907 908 909 - 910 911 912 ); - water ( - 601 602 603 - 604 605 606 - 607 608 609 - 610 611 612); - default 302; } - - density { type cartesianField; file ./myDensityField; } - Visualiser ---------- @@ -1060,7 +1065,8 @@ aceNeutronDatabase, used for continuous energy data. In this case, the data is r from ACE files. * aceLibrary: includes the path to the *.aceXS* file, which includes the paths to - the ACE files + the ACE files. If the environmental variable ``SCONE_ACE`` is defined this can be + used instead of the path. * ures (*optional*, default = 0): 1 for true; 0 for false; activates the unresolved resonance probability tables treatment * DBRC (*optional*, default = no DBRC): list of ZAIDs of nuclides for which DBRC has @@ -1078,6 +1084,8 @@ Example: :: ceData { type aceNuclearDatabase; aceLibrary ./myFolder/ACElib/JEF311.aceXS; ures 1; DBRC (92238 94242); avgDist 32; energyPerFission 200.0;} + ceData { type aceNuclearDatabase; aceLibrary SCONE_ACE;} + .. note:: If DBRC is applied, the 0K cross section ace files of the relevant nuclides must be included in the aceLibrary file. @@ -1695,14 +1703,6 @@ Examples: :: map1 { type spaceMap; axis x; grid lin; min -50.0; max 50.0; N 100; } map2 { type spaceMap; axis z; grid unstruct; bins (0.0 0.2 0.3 0.5 0.7 0.8 1.0); } -* fieldMap (1D map), map over superimposed fields. Limited currently to pieceConstantFields. - - - field: field definition, corresponding to those in pieceConstantFields. - -Examples: :: - - map1 { type fieldMap; field {file ./myField.txt } } - * timeMap (1D map), maps particles point in time - grid: ``lin`` for linearly spaced bins diff --git a/docs/Installation.rst b/docs/Installation.rst index cd3ac5e58..9e4e6ac35 100644 --- a/docs/Installation.rst +++ b/docs/Installation.rst @@ -23,7 +23,7 @@ Requirements GNU/Linux operating system In principle any UNIX-like operating system should support SCONE. However people have experienced some significant run-time bugs when running on MacOS. Since we do not have - an access to a Mackintosh we were not able to identify the problem yet. + an access to a Macintosh we were not able to identify the problem yet. .. admonition:: Optional @@ -290,7 +290,7 @@ Obtaining Nuclear Data '''''''''''''''''''''' SCONE requires ACE-formatted nuclear data. The JEFF-3.3 evaluation can be download from the -OACD NEA `website `__. In addition SCONE requires +OECD NEA `website `__. In addition SCONE requires its own library file. An example of it is given in *IntegrationTestFiles/testLib*. Its format is:: ! This is a comment line @@ -333,3 +333,12 @@ library with thermal data we need to do the following: ./scripts/make_ace_lib.sh ./tempCE CE ace ./path_to_CE_ace_files/ ./scripts/make_ace_lib.sh ./tempSAB SAB t ./path_to_SAB_ace_files/ cat tempCE tempSAB > fullLib.xsfile + +Having made this file, one can optionally define the environmental variable ``SCONE_ACE``, e.g., +in the .basrhc file and refer to it in input files which use the aceNeutronDatabase. To do so, in +the .bashrc file include: + +.. code-block:: bash + + export SCONE_ACE=/path/to/fullLib.xsfile + diff --git a/docs/index.rst b/docs/index.rst index 9fadfa780..0bfb6771c 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,10 +3,8 @@ What is SCONE? SCONE stands for **S**\ tochastic **C**\ alculator **O**\ f **N**\ eutron Transport **E**\ quation. It is a Monte Carlo particle transport code for prototyping -of new methods. Unlike its more established cousins it does not strive for the fastest -execution speed, since it is not intended to be used for practical design calculation. Instead it -aims to offer flexibility to allow relatively inexperienced users, such as graduate students -to quickly test their ideas. +of new methods. It aims to offer flexibility to allow relatively inexperienced users, +such as graduate students to quickly test their ideas.