From e0f0e44c7a98bcf9e90836096682c97ef4f62ce0 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Tue, 24 Feb 2026 19:54:38 +0000 Subject: [PATCH 01/15] Added fission * energy release as a cross section --- .../Reactions/reactionMG/fissionMG_class.f90 | 21 ++++++++++++++++- .../fissionCE_class.f90 | 17 ++++++++++++++ .../aceDatabase/aceNeutronDatabase_class.f90 | 3 +++ .../aceDatabase/aceNeutronNuclide_class.f90 | 23 +++++++++++++++++-- .../ceNeutronData/ceNeutronDatabase_inter.f90 | 6 +++++ .../baseMgNeutronMaterial_class.f90 | 2 ++ .../testNeutronDatabase_class.f90 | 11 ++++++++- .../xsPackages/neutronXsPackages_class.f90 | 12 ++++++++-- SharedModules/endfConstants.f90 | 8 ++++--- SharedModules/universalVariables.f90 | 3 +-- .../Tests/macroResponse_test.f90 | 13 +++++++++-- .../Tests/microResponse_test.f90 | 14 +++++++++-- .../TallyResponses/macroResponse_class.f90 | 3 +++ .../TallyResponses/microResponse_class.f90 | 3 +++ docs/Input Manual.rst | 7 +++--- 15 files changed, 128 insertions(+), 18 deletions(-) diff --git a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 index 3d51ee384..47d41d89f 100644 --- a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 +++ b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 @@ -26,14 +26,17 @@ module fissionMG_class !! Fission spectrum is independent of the energy group !! !! Public members: - !! data -> data for fissionMG [energyGroup, reaction (nu or chi)] + !! data -> data for fissionMG [energyGroup, reaction (nu or chi)] + !! kappa -> Energy per fission. 200 MeV by default. !! !! Interface: !! reactionMG interface !! buildFromDict -> builds fissionMG from a SCONE dictionary + !! getKappa -> obtains the energy per fission in MeV !! type, public, extends(reactionMG) :: fissionMG real(defReal),dimension(:,:),allocatable :: data + real(defReal) :: kappa = 200.0_defReal contains ! Superclass procedures procedure :: init @@ -46,6 +49,7 @@ module fissionMG_class ! Local procedures procedure :: buildFromDict + procedure :: getKappa end type fissionMG @@ -216,6 +220,9 @@ subroutine buildFromDict(self, dict) ! Get number of groups call dict % get(nG, 'numberOfGroups') + + ! Get energy per fission + call dict % getOrDefault(self % kappa, 'kappa', 200.0_defReal) ! Allocate space allocate(self % data(nG, 2)) @@ -246,6 +253,18 @@ subroutine buildFromDict(self, dict) end subroutine buildFromDict + !! + !! Get the energy release from fission + !! + pure function getKappa(self) result(kappa) + class(fissionMG), intent(in) :: self + real(defReal) :: kappa + + kappa = self % kappa + + end function getKappa + + !! !! Cast reactionHandle pointer to fissionMG pointer !! diff --git a/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 b/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 index bb64e01a2..c782f9672 100644 --- a/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 +++ b/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 @@ -54,6 +54,7 @@ module fissionCE_class !! eLawPrompt -> Energy Law for the prompt fission neutrons !! nuBarDelayed -> Delayed release table for incindent energy [MeV] !! delayed -> Information about Delayed emission precursors + !! Q -> Q value for the reaction [MeV] !! !! Interface: !! uncorrelatedReactionCE interface @@ -65,6 +66,7 @@ module fissionCE_class class(energyLawENDF), allocatable :: eLawPrompt class(releaseLawENDF),allocatable :: nuBarDelayed type(precursor),dimension(:),allocatable :: delayed + real(defReal) :: Q = ZERO contains ! Superclass procedures @@ -79,6 +81,7 @@ module fissionCE_class ! Type specific procedures procedure :: buildFromACE + procedure :: getQ end type fissionCE contains @@ -148,8 +151,21 @@ elemental subroutine kill(self) deallocate(self % delayed) end if + self % Q = ZERO + end subroutine kill + !! + !! Returns the Q-value + !! + pure function getQ(self) result(Q) + class(fissionCE), intent(in) :: self + real(defReal) :: Q + + Q = self % Q + + end function getQ + !! !! Returns true if reaction is in Centre-Of-Mass frame !! @@ -349,6 +365,7 @@ subroutine buildFromACE(self, ACE, MT) ! Read basic data call new_totalNU(self % nuBarTotal, ACE) call new_energyLawENDF(self % eLawPrompt, ACE, MT) + self % Q = ACE % QforMT(MT) ! Read Delayed Data if (withDelayed) then diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 6899dd076..f6d7c8528 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -854,6 +854,9 @@ subroutine init(self, dict, ptr, silent ) end if + ! Obtain the energy per fission with which to scale fission heating + call dict % getOrDefault(self % H235, 'energyPerFission', 202.27_defReal) + ! Get path to ACE library call dict % get(aceLibPath,'aceLibrary') diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index 2f1979edf..c92cbee1f 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -47,6 +47,7 @@ module aceNeutronNuclide_class integer(shortInt), parameter :: CAPTURE_XS = 4 integer(shortInt), parameter :: FISSION_XS = 5 integer(shortInt), parameter :: NU_FISSION = 6 + integer(shortInt), parameter :: KAPPA_XS = 7 !! @@ -436,9 +437,11 @@ elemental subroutine microXSs(self, xss, idx, f) if (self % isFissile()) then xss % fission = data(FISSION_XS, 2) * f + (ONE-f) * data(FISSION_XS, 1) xss % nuFission = data(NU_FISSION, 2) * f + (ONE-f) * data(NU_FISSION, 1) + xss % kappaXS = data(KAPPA_XS, 2) * f + (ONE-f) * data(KAPPA_XS, 1) else xss % fission = ZERO xss % nuFission = ZERO + xss % kappaXS = ZERO end if end associate @@ -483,9 +486,11 @@ subroutine getThXSs(self, xss, idx, f, E, kT, rand) if (self % isFissile()) then xss % fission = data(FISSION_XS, 2) * f + (ONE-f) * data(FISSION_XS, 1) xss % nuFission = data(NU_FISSION, 2) * f + (ONE-f) * data(NU_FISSION, 1) + xss % kappaXS = data(KAPPA_XS, 2) * f + (ONE-f) * data(KAPPA_XS, 1) else xss % fission = ZERO xss % nuFission = ZERO + xss % kappaXS = ZERO end if ! Read S(a,b) tables for elastic scatter: return zero if elastic scatter is off. @@ -555,9 +560,11 @@ subroutine getUrrXSs(self, xss, idx, f, E, xi) if (self % isFissile()) then xss % fission = data(FISSION_XS, 2) * f + (ONE-f) * data(FISSION_XS, 1) xss % nuFission = data(NU_FISSION, 2) * f + (ONE-f) * data(NU_FISSION, 1) + xss % kappaXS = data(KAPPA_XS, 2) * f + (ONE-f) * data(KAPPA_XS, 1) else xss % fission = ZERO xss % nuFission = ZERO + xss % kappaXS = ZERO end if ! Check if flag for multiplication factor (IFF) is true, and apply it to elastic scattering, @@ -587,6 +594,7 @@ subroutine getUrrXSs(self, xss, idx, f, E, xi) if (self % isFissile()) then xss % nuFission = xss % nuFission/xss % fission * val(3) + xss % kappaXS = xss % kappaXS / xss % fission * val(3) xss % fission = val(3) end if @@ -725,6 +733,7 @@ subroutine init(self, ACE, nucIdx, database) top, firstIdxMT4 real(defReal), dimension(:), allocatable :: xsMT4 type(stackInt) :: scatterMT, absMT + real(defReal) :: H_Q character(100), parameter :: Here = "init (aceNeutronNuclide_class.f90)" ! Reset nuclide just in case @@ -745,7 +754,7 @@ subroutine init(self, ACE, nucIdx, database) ! Allocate space for main XSs if (self % isFissile()) then - N = 6 + N = 7 else N = 4 end if @@ -807,10 +816,20 @@ subroutine init(self, ACE, nucIdx, database) call self % fission % init(ACE, N_f) end if - ! Calculate nuFission + ! Obtain Heating/Q scaling ratio + ! Check if database is associated in order to satisfy tests where it might not be! + if (associated(database)) then + H_Q = database % H235 / database % Q235 + else + H_Q = ONE + end if + + ! Calculate nuFission and kappaXS do i = bottom, Ngrid self % mainData(NU_FISSION,i) = self % mainData(FISSION_XS,i) * & self % fission % release(self % eGrid(i)) + self % mainData(KAPPA_XS,i) = self % mainData(FISSION_XS,i) * & + self % fission % getQ() * H_Q end do end if diff --git a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 index c044c707e..ef368c55b 100644 --- a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 @@ -39,6 +39,8 @@ module ceNeutronDatabase_inter !! !! Public Members: !! mapDBRCnuc -> map to link indexes of DBRC nuclides with their corresponding 0K + !! H235 -> heating fission of U235, used for scaling fission energy + !! Q235 -> Q-value of U235 fission, used for scaling fission energy !! !! Interface: !! nuclearDatabase Interface @@ -53,6 +55,10 @@ module ceNeutronDatabase_inter !! type, public, abstract, extends(nuclearDatabase) :: ceNeutronDatabase type(intMap) :: mapDBRCnuc + + ! Default fission scaling [MeV] + real(defReal) :: H235 = 202.27_defReal + real(defReal) :: Q235 = 193.406_defReal contains diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 index dc98bdc35..068a16549 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 @@ -129,9 +129,11 @@ subroutine getMacroXSs_byG(self, xss, G, rand) if(self % isFissile()) then xss % fission = self % data(FISSION_XS, G) xss % nuFission = self % data(NU_FISSION, G) + xss % kappaXS = xss % fission * self % fission % getKappa() else xss % fission = ZERO xss % nuFission = ZERO + xss % kappaXS = ZERO end if end subroutine getMacroXSs_byG diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index e6032e73b..dbe12118f 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -70,11 +70,12 @@ module testNeutronDatabase_class !! captureXS [in] -> Optional. Value of Capture XS !! fissionXS [in] -> Optional. Value of Fission XS !! nuFissionXS [in] -> Optional Value of nuFission + !! kappaXS [in] -> Optional Value of kappa XS !! !! Errors: !! None !! - subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuFissionXS) + subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuFissionXS, kappaXS) class(testNeutroNDatabase), intent(inout) :: self real(defReal), intent(in) :: xsVal real(defReal), intent(in),optional :: eScatterXS @@ -82,6 +83,7 @@ subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuF real(defReal), intent(in),optional :: captureXS real(defReal), intent(in),optional :: fissionXS real(defReal), intent(in),optional :: nuFissionXS + real(defReal), intent(in),optional :: kappaXS self % xsVal = xsVal @@ -125,6 +127,13 @@ subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuF else self % mat % xss % nuFission = xsVal end if + + ! kappa * fission + if(present(kappaXS)) then + self % mat % xss % kappaXS = kappaXS + else + self % mat % xss % kappaXS = xsVal + end if end subroutine build diff --git a/NuclearData/xsPackages/neutronXsPackages_class.f90 b/NuclearData/xsPackages/neutronXsPackages_class.f90 index efb3cf159..6ad3e9771 100644 --- a/NuclearData/xsPackages/neutronXsPackages_class.f90 +++ b/NuclearData/xsPackages/neutronXsPackages_class.f90 @@ -2,7 +2,6 @@ !! This module brakes standard rules !! It contains a library of XS Packages for Neutron particle type !! -!! module neutronXsPackages_class use numPrecision @@ -22,6 +21,7 @@ module neutronXsPackages_class !! capture -> sum of all reactions without secendary neutrons excluding fission [1/cm] !! fission -> total Fission MT=18 Cross-section [1/cm] !! nuFission -> total average neutron production Cross-section [1/cm] + !! kappaXS -> heating cross-section [MeV/cm] !! !! Interface: !! clean -> Set all XSs to 0.0 @@ -35,6 +35,7 @@ module neutronXsPackages_class real(defReal) :: capture = ZERO real(defReal) :: fission = ZERO real(defReal) :: nuFission = ZERO + real(defReal) :: kappaXS = ZERO contains procedure :: clean => clean_neutronMacroXSs procedure :: add => add_neutronMacroXSs @@ -54,6 +55,7 @@ module neutronXsPackages_class !! capture -> all reactions without secendary neutrons excluding fission [barn] !! fission -> total Fission MT=18 Cross-section [barn] !! nuFission -> total average neutron production Cross-section [barn] + !! kappaXS -> heating Cross-section [MeV*barn] !! type, public :: neutronMicroXSs real(defReal) :: total = ZERO @@ -62,6 +64,7 @@ module neutronXsPackages_class real(defReal) :: capture = ZERO real(defReal) :: fission = ZERO real(defReal) :: nuFission = ZERO + real(defReal) :: kappaXS = ZERO contains procedure :: invert => invert_microXSs end type neutronMicroXSs @@ -88,6 +91,7 @@ elemental subroutine clean_neutronMacroXSs(self) self % capture = ZERO self % fission = ZERO self % nuFission = ZERO + self % kappaXS = ZERO end subroutine clean_neutronMacroXSs @@ -114,6 +118,7 @@ elemental subroutine add_neutronMacroXSs(self, micro, dens) self % capture = self % capture + dens * micro % capture self % fission = self % fission + dens * micro % fission self % nuFission = self % nuFission + dens * micro % nuFission + self % kappaXS = self % kappaXS + dens * micro % kappaXS end subroutine add_neutronMacroXSs @@ -159,6 +164,9 @@ elemental function get(self, MT) result(xs) case(macroNuFission) xs = self % nuFission + case(macroKappaFission) + xs = self % kappaXS + case(macroAbsorbtion) xs = self % fission + self % capture @@ -233,7 +241,7 @@ end function invert_macroXSs !! !! Use a real r in <0;1> to sample reaction from Microscopic XSs !! - !! This function involves a bit of code so is written for conviniance + !! This function involves a bit of code so is written for convenience !! !! Args: !! r [in] -> Real number in <0;1> diff --git a/SharedModules/endfConstants.f90 b/SharedModules/endfConstants.f90 index e5ed867a2..97d68b000 100644 --- a/SharedModules/endfConstants.f90 +++ b/SharedModules/endfConstants.f90 @@ -87,6 +87,7 @@ module endfConstants N_Xt = 205 ,& N_X3He = 206 ,& N_Xa = 207 ,& + N_KAPPA = 301 ,& nubar_tot = 452 ,& nubar_del = 455 ,& nubar_prompt = 456 ,& @@ -114,7 +115,8 @@ module endfConstants macroEscatter = -3 ,& macroIEscatter = -4 ,& macroFission = -6 ,& - macroNuFission = -7 + macroNuFission = -7 ,& + macroKappaFission = -80 ! List of Macro MT numbers for macroscopic XSs. Unique to SCONE (not from Serpent) integer(shortInt), parameter :: macroAllScatter = -20 ,& @@ -129,10 +131,10 @@ module endfConstants 36, 37, 38, 41, 42, 44, 45, 91, 101, 102, 103, & 104, 105, 106, 107, 108, 109, 111, 112, 113, & 114, 115, 116, 117, 203, 204, 205, 206, 207, & - 891, N_Nl, N_2Nl] + 301, 891, N_Nl, N_2Nl] integer(shortInt), dimension(*), parameter :: availableMacroMTs = & - [-1, -2, -3, -4, -6, -7, -20, -21, -22] + [-1, -2, -3, -4, -6, -7, -20, -21, -22, -80] diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 790fff3e5..14a6d15b8 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -80,8 +80,7 @@ module universalVariables ! 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 - kBoltzmann = 1.380649e-23_defReal, & ! Bolztmann constant in J/K - energyPerFission = 200.0_defReal ! MeV + kBoltzmann = 1.380649e-23_defReal ! Bolztmann constant in J/K ! Unit conversion real(defReal), parameter :: joulesPerMeV = 1.60218e-13_defReal ,& ! Convert MeV to J diff --git a/Tallies/TallyResponses/Tests/macroResponse_test.f90 b/Tallies/TallyResponses/Tests/macroResponse_test.f90 index bfa87292e..103c61007 100644 --- a/Tallies/TallyResponses/Tests/macroResponse_test.f90 +++ b/Tallies/TallyResponses/Tests/macroResponse_test.f90 @@ -18,6 +18,7 @@ module macroResponse_test type(macroResponse) :: response_fission type(macroResponse) :: response_nuFission type(macroResponse) :: response_absorbtion + type(macroResponse) :: response_kappa type(macroResponse) :: response_MT type(testNeutronDatabase) :: xsData contains @@ -36,8 +37,8 @@ subroutine setUp(this) type(dictionary) :: tempDict ! Allocate and initialise test nuclearData - ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission - call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal) + ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappFission + call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal, 9.0_defReal) ! Set up responses ! Total @@ -75,6 +76,13 @@ subroutine setUp(this) call this % response_absorbtion % init(tempDict) call tempDict % kill() + ! Energy deposition + call tempDict % init(2) + call tempDict % store('type','macroResponse') + call tempDict % store('MT', macroKappaFission) + call this % response_kappa % init(tempDict) + call tempDict % kill() + ! MT number call tempDict % init(2) call tempDict % store('type','macroResponse') @@ -117,6 +125,7 @@ subroutine testGettingResponse(this) @assertEqual(1.0_defReal, this % response_fission % get(p, this % xsData), TOL) @assertEqual(1.5_defReal, this % response_nuFission % get(p, this % xsData), TOL) @assertEqual(3.0_defReal, this % response_absorbtion % get(p, this % xsData), TOL) + @assertEqual(9.0_defReal, this % response_kappa % get(p, this % xsData), TOL) @assertEqual(105.0_defReal, this % response_MT % get(p, this % xsData), TOL) p % isMG = .true. diff --git a/Tallies/TallyResponses/Tests/microResponse_test.f90 b/Tallies/TallyResponses/Tests/microResponse_test.f90 index ec1e810f9..8f96e8adf 100644 --- a/Tallies/TallyResponses/Tests/microResponse_test.f90 +++ b/Tallies/TallyResponses/Tests/microResponse_test.f90 @@ -19,6 +19,7 @@ module microResponse_test type(microResponse) :: response_capture type(microResponse) :: response_fission type(microResponse) :: response_absorbtion + type(microResponse) :: response_kappa type(microResponse) :: response_MT type(testNeutronDatabase) :: xsData contains @@ -38,8 +39,8 @@ subroutine setUp(this) ! Allocate and initialise test nuclearData - ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission - call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal) + ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappFission + call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal, 9.0_defReal) ! Set dictionaries to initialise material call dictMat1 % init(1) @@ -97,6 +98,14 @@ subroutine setUp(this) call this % response_absorbtion % init(tempDict) call tempDict % kill() + ! Energy deposition + call tempDict % init(3) + call tempDict % store('type', 'microResponse') + call tempDict % store('MT', N_KAPPA) + call tempDict % store('material', 'Xenon') + call this % response_kappa % init(tempDict) + call tempDict % kill() + ! MT number call tempDict % init(2) call tempDict % store('type','microResponse') @@ -140,6 +149,7 @@ subroutine testGettingResponse(this) @assertEqual(0.5_defReal, this % response_fission % get(p, this % xsData), TOL) @assertEqual(1.5_defReal, this % response_eScatter % get(p, this % xsData), TOL) @assertEqual(1.5_defReal, this % response_absorbtion % get(p, this % xsData), TOL) + @assertEqual(4.5_defReal, this % response_kappa % get(p, this % xsData), TOL) @assertEqual(8.0_defReal, this % response_MT % get(p, this % xsData), TOL) p % isMG = .true. diff --git a/Tallies/TallyResponses/macroResponse_class.f90 b/Tallies/TallyResponses/macroResponse_class.f90 index 44ba0c9f8..d5d488d84 100644 --- a/Tallies/TallyResponses/macroResponse_class.f90 +++ b/Tallies/TallyResponses/macroResponse_class.f90 @@ -115,6 +115,9 @@ subroutine build(self, MT) case(N_ABSORPTION) self % MT = macroAbsorbtion + case(N_KAPPA) + self % MT = macroKappaFission + case default self % mainData = .false. self % MT = MT diff --git a/Tallies/TallyResponses/microResponse_class.f90 b/Tallies/TallyResponses/microResponse_class.f90 index fd7c63808..9bd5bb72f 100644 --- a/Tallies/TallyResponses/microResponse_class.f90 +++ b/Tallies/TallyResponses/microResponse_class.f90 @@ -143,6 +143,9 @@ subroutine build(self, MT) case(N_ABSORPTION) self % MT = macroAbsorbtion + case(N_KAPPA) + self % MT = macroKappaFission + case default self % MT = MT self % mainData = .false. diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index cc4ceb1bd..874d54cb2 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -1370,8 +1370,9 @@ Example: :: * macroResponse: used to score macroscopic reaction rates - MT: MT number of the desired reaction. The options are: -1 (total), -2 (disappearance), - -3 (elastic scattering), -4 (total inelastic scattering), -6 (fission), -7 nu*fission), - -20 (total scattering), -21 (absorption), -22 (total non elastic, i.e., absorption + inelastic). + -3 (elastic scattering), -4 (total inelastic scattering), -6 (fission), -7 (nu*fission), + -20 (total scattering), -21 (absorption), -22 (total non elastic, i.e., absorption + inelastic), + -80 (kappa*fission). Additionally, all the MT numbers allowed by microResponse can be used here. Example: :: @@ -1385,7 +1386,7 @@ Example: :: * microResponse: used to score microscopic reaction rates - MT: MT number of the desired reaction. The options are: 1, 2, 3, 4, 5, 11, 16-25, 27-30, - 32-38, 41, 42, 44, 45, 51-90, 91, 101-109, 111-117, 203-207, 875-890. These MT numbers + 32-38, 41, 42, 44, 45, 51-90, 91, 101-109, 111-117, 203-207, 301, 875-890. These MT numbers are defined in the conventional way, i.e., following the ENDF standard - material: material name where to score the reaction. The material must be defined to include only one nuclide; its density could be anything, it doesn't From de6b78eb5b7581f7e833deacf386e4c7ba3203b5 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Sun, 1 Mar 2026 23:02:54 +0000 Subject: [PATCH 02/15] Most of the coupling admin --- Coupling/CMakeLists.txt | 6 + Coupling/couplingAdmin_class.f90 | 365 ++++++++++++++++++ PhysicsPackages/eigenPhysicsPackage_class.f90 | 11 + Tallies/scoreMemory_class.f90 | 16 + Tallies/tallyAdmin_class.f90 | 35 +- 5 files changed, 422 insertions(+), 11 deletions(-) create mode 100644 Coupling/CMakeLists.txt create mode 100644 Coupling/couplingAdmin_class.f90 diff --git a/Coupling/CMakeLists.txt b/Coupling/CMakeLists.txt new file mode 100644 index 000000000..fa0bfa5bd --- /dev/null +++ b/Coupling/CMakeLists.txt @@ -0,0 +1,6 @@ +# Add Source Files to a global List +add_sources( ./couplingAdmin_class.f90) + +# Add Unit Tests to a global List +add_unit_tests( ./Tests/couplingAdmin_test.f90) + diff --git a/Coupling/couplingAdmin_class.f90 b/Coupling/couplingAdmin_class.f90 new file mode 100644 index 000000000..1d3e3988b --- /dev/null +++ b/Coupling/couplingAdmin_class.f90 @@ -0,0 +1,365 @@ +module couplingAdmin_class + + use numPrecision + use universalVariables only : nameDensity, nameTemperature + use dictionary_class, only : dictionary + use dictParser_func, only : fileToDict + use outputFile_class, only : outputFile + use errors_mod, only : fatalError + + use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + + implicit none + private + + character(nameLen), parameter, dimension(2) :: ALLOWABLE_FIELDS = ['temperature', 'density'] + character(nameLen), parameter :: END_SIGNAL = "SIGUSR2", & + CONTINUE_SIGNAL = "SIGUSR1" + + !! + !! Object to manage coupling between radiation and external physics solvers/driver scripts + !! + !! This object is responsible for managing the status of coupled calculations. It should do the following: + !! - say whether coupling is happening (default: no) + !! - know the means of communication by which coupling occurs (currently files are used as flags) + !! - send a signal to the driver that it has finished a number of iterations + !! - wait for the driver/other physics to finish their calculations + !! - restart on receiving a signal from the driver + !! - create any tallies necessary for coupling (such as power tallies) + !! - output the results of those tallies to an appropriate location + !! - know the location of fields provided by other physics and update them in SCONE + !! + !! Private members: + !! isCoupled -> flag for whether coupling is being performed + !! commFileIn -> path to the file which triggers a SCONE calculation by the driver + !! commFileOut -> path to the file which SCONE outputs to the driver on completing a calculation + !! outputFile -> path to the file which SCONE outputs containing results to be used by the driver + !! updateFreq -> frequency with which SCONE pauses to exchange data with the driver + !! maxIt -> maximum number of iterations before finishing coupling + !! tally -> Tally admin containing necessary tallies for the driver + !! fieldPaths -> paths to files which SCONE reads to update the given fields + !! + !! Public interface: + !! init -> initialises the coupling setup from a dictionary + !! kill -> clean up + !! couple -> performs the main coupling logic + !! endCoupling -> deactivates coupling + !! + !! doCoupling -> returns a logical if coupling is happening + !! doUpdate -> returns a logical if a physics update should be performed on a given iteration. + !! + !! makeTallies -> returns a tally admin used for coupling + !! outputTallies -> outputs tally data to be read by other physics + !! updateFields -> updates the relevant physics fields before resuming a calculation + !! + !! waitForSignal -> pauses the calculation until a 'resume' signal is received + !! signalIterationOver -> signals that an iteration is complete + !! signalCalculationOver -> signals that the coupled calculation is complete + !! + type, public :: couplingAdmin + private + logical(defBool) :: isCoupled = .false. + character(pathLen) :: commFileIn = '' + character(pathLen) :: commFileOut = '' + character(pathLen) :: outputFile = '' + character(nameLen) :: outputFormat = '' + integer(shortInt) :: updateFreq = 0 + integer(shortInt) :: maxIt = huge(shortInt) + type(tallyAdmin) :: tally + character(pathLen), dimension(:), allocatable :: fieldPaths + contains + procedure :: init + procedure :: kill + procedure :: couple + procedure :: endCoupling + + ! Status checks + procedure :: doCoupling + procedure, private :: doUpdate + + ! Results handling and manipulation + procedure :: attachTally + procedure :: outputTallies + procedure, private :: updateFields + + ! Communication + procedure, private :: waitForSignal + procedure, private :: signalIterationOver + procedure, private :: signalCalculationOver + + end type couplingAdmin + +contains + + !! + !! Initialises coupling from a dictionary + !! + subroutine init(self, dict) + class(couplingAdmin), intent(inout) :: self + class(dictionary), intent(in) :: dict + integer(shortInt) :: i, nFields + logical(defBool) :: found, duplicate + class(dictionary), pointer :: tempDict + type(outputFile) :: test_out + character(100), parameter :: Here ='init (couplingAdmin_class.f90)' + + ! Read communication details + call dict % get(self % commFileIn, 'receiveFile') + call dict % get(self % commFileOut, 'sendFile') + + ! Read how often to exchange data + call dict % get(self % updateFreq, 'updateFreq') + if (self % updateFreq < 1) call fatalError(Here,& + 'Physics update frequency is silly: '//numToChar(self % updateFreq)) + + ! Read maximum number of iterations over which to perform coupling + call dict % getOrDefault(self % maxIt, 'maxIt', huge(shortInt)) + if (self % maxIt < 1) call fatalError(Here, & + 'Maximum iteration number is silly: '//numToChar(self % maxIt)) + + ! Read tally info needed by other solvers + tallyDict => dict % getDictPtr('tally') + call self % tally % init(tallyDict) + + ! Read the output file name + call dict % get(self % outputFile, 'outputFile') + + ! Initialise output file before calculation (so mistake in format will be caught early) + call dict % getOrDefault(self % outputFormat, 'outputFormat', 'asciiMATLAB') + call test_out % init(self % outputFormat) + + ! Read the field update location(s) + call dict % get(self % fieldPaths, 'fieldPaths') + + ! Read the field name(s) + call dict % get(self % fieldNames, 'fieldNames') + + ! Ensure the paths correspond to the name + if (size(self % fieldPaths) /= size(self % fieldNames)) call fatalError(Here,& + 'The number of field paths must match the number of field names.') + + ! Ensure the field names are unique and allowable + nFields = size(self % fieldNames) + do i = 1, nFields + + found = any(trim(self % fieldNames(i)) == trim(ALLOWABLE_FIELDS)) + if (.not. found) call fatalError(Here,'Field names must be one of: '//ALLOWABLE_FIELDS) + + if (i < nFields) then + duplicate = any(trim(self % fieldNames(i)) == trim(self % fieldNames(i+1:nFields))) + if (duplicate) call fatalError(Here,'Field names must be unique: '//trim(self % fieldNames(i))) + end if + + end do + + end subroutine init + + !! + !! Clean-up + !! + subroutine kill(self) + class(couplingAdmin), intent(inout) :: self + + self % isCoupled = .false. + self % updateFreq = 0 + self % commFileIn = '' + self % commFileOut = '' + self % outputFile = '' + self % outputFormat = '' + self % maxIt = huge(shortInt) + call self % tally % kill() + if (allocated(fieldPaths)) deallocate(fieldPaths) + + end subroutine kill + + !! + !! Subroutine which performs coupling logic during transport. + !! Evaluates whether a physics update should be performed. + !! If so, outputs tally information, cleans up tallies, signals and waits, + !! and updates fields. Deactivates coupling if later than the maximum + !! iteration. + !! + subroutine couple(self, it) + class(couplingAdmin), intent(inout) :: self + integer(shortInt), intent(in) :: it + + if (.not. self % doCoupling()) return + if (.not. self % doUpdate(it)) return + + call self % outputTallies() + call self % tally % resetMemory() + + if (it < self % maxIt) then + call self % signalIterationOver() + call self % waitForSignal() + call self % updateFields() + else + call self % signalCalculationOver() + self % isCoupled = .false. + end if + + end subroutine couple + + !! + !! Returns whether coupling should be occurring or not + !! + function doCoupling(self) result(isCoupled) + class(couplingAdmin), intent(in) :: self + logical(defBool) :: isCoupled + + isCoupled = self % isCoupled + + end function doCoupling + + !! + !! End coupling + !! + subroutine endCoupling(self) + class(couplingAdmin), intent(inout) :: self + + if (self % isCoupled) then + self % doCoupling = .false. + call self % signalCalculationOver() + end if + + end subroutine endCoupling + + !! + !! Returns whether to perform a physics update based on iteration number + !! and update frequency + !! + function doUpdate(self, it) result(update) + class(couplingAdmin), intent(in) :: self + integer(shortInt), intent(in) :: it + logical(defBool) :: update + + update = mod(it, self % updateFreq) == 0 + + end function doUpdate + + !! + !! Update the fields from the paths to the new field info + !! + subroutine updateFields(self) + class(couplingAdmin), intent(in) :: self + integer(shortInt) :: i + class(dictionary) :: dict + character(nameLen) :: fieldName + class(field), pointer :: ptr + character(100), parameter :: Here = 'updateFields (couplingAdmin_class.f90)' + + do i = 1, size(self % fieldPaths) + call fileToDict(dict, self % fieldPaths(i)) + + ! Obtain the pointer to the corresponding field. + ! Reinitialise the field from the dictionary. + select case(trim(self % fieldNames(i))) + case('temperature') + fieldName = nameTemperature + + case('density') + fieldName = nameDensity + + case default + call fatalError(Here, 'Unrecognised field type requested') + end select + + ptr => gr_fieldPtr(gr_fieldIdx(fieldName)) + + call ptr % kill() + call ptr % init(dict) + + end do + + end subroutine updateFields + + !! + !! Writes a file to signal that the current iteration is finished + !! + subroutine signalIterationOver(self) + class(couplingAdmin), intent(in) :: self + integer(shortInt) :: unit + + open(newunit=unit, file=self % commFileOut, status="replace", action="write") + write(unit, '(A)') CONTINUE_SIGNAL + close(unit) + + end subroutine signalIterationOver + + !! + !! Writes a file to signal that the entire calculation is finished + !! + subroutine signalCalculationOver(self) + class(couplingAdmin), intent(in) :: self + integer(shortInt) :: unit + + open(newunit=unit, file=self % commFileOut, status="replace", action="write") + write(unit, '(A)') END_SIGNAL + close(unit) + + end subroutine signalCalculationOver + + !! + !! Pauses the calculation while waiting for the existence + !! of a signal file. + !! + subroutine waitForSignal(self) + class(couplingAdmin), intent(in) :: self + logical(defBool) :: exists + integer(shortInt) :: unit + character(*) :: line + + do + + ! Check for the presence of a signal file + inquire(file=self % commFileIn, exist=exists) + + if (exists) then + ! Read and delete the communication file + open(newunit=unit, file=self % commFileIn, status="old") + + read(unit, '(A)') line + + ! Check whether to continue iteration or end coupling + if (trim(adjustl(line)) == END_SIGNAL) then + self % isCoupled = .false. + elseif (trim(adjustl(line)) /= CONTINUE_SIGNAL) then + call fatalError(Here, 'Unrecognised signal: '//trim(adjustl(line))) + end if + + close(unit, status="delete") + exit + end if + + call sleep(1) + + end do + + end subroutine waitForSignal + + !! + !! Attach coupling tallies to another tally + !! + subroutine attachTally(self, tallyPtr) + class(couplingAdmin), intent(in) :: self + class(tallyAdmin), pointer, intent(inout) :: tallyPtr + + if (self % doCoupling()) + call tallyPtr % push(self % tally) + end if + + end subroutine attachTally + + !! + !! Output tally results + !! + subroutine outputTallies(self) + class(couplingAdmin), intent(in) :: self + type(outputFile) :: out + + call out % init(self % outputFormat, filename = self % outputFile) + call self % tally % print(out) + + end subroutine outputTallies + +end module couplingAdmin_class diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index 878780e4f..a479466dd 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -445,6 +445,12 @@ subroutine init(self, dict) ! Parallel buffer size call dict % getOrDefault(self % bufferSize, 'buffer', 1000) + ! Check if the calculation is coupled + if (dict % isPresent('coupling')) then + tempDict => dict % getDictPtr('coupling') + call self % couplingInfo % init(tempDict) + end if + ! Process type of data select case(energy) case('mg') @@ -556,6 +562,7 @@ subroutine init(self, dict) allocate(self % inactiveTally) call self % inactiveTally % init(tempDict) + tempDict => dict % getDictPtr('activeTally') allocate(self % activeTally) call self % activeTally % init(tempDict) @@ -609,7 +616,11 @@ subroutine init(self, dict) ! Attach attachments to result tallies call self % inactiveTally % push(self % inactiveAtch) call self % activeTally % push(self % activeAtch) + + ! Attach a tally admin for coupling + if (self % couplingInfo % isCoupled()) then + end if call self % printSettings() diff --git a/Tallies/scoreMemory_class.f90 b/Tallies/scoreMemory_class.f90 index bf67940c5..c125ed1fc 100644 --- a/Tallies/scoreMemory_class.f90 +++ b/Tallies/scoreMemory_class.f90 @@ -66,6 +66,8 @@ module scoreMemory_class !! !! reduceBins(): Move the scores from parallelBins and different processes to bins. !! + !! reset(): Resets the memory: zeros all scores and reverts cycles/batches to 0. + !! !! Example use case: !! !! do batches=1,20 @@ -112,6 +114,7 @@ module scoreMemory_class procedure :: getBatchSize procedure :: reduceBins procedure :: collectDistributed + procedure :: reset ! Private procedures procedure, private :: score_defReal @@ -619,4 +622,17 @@ elemental function getScore(self, idx) result (score) end function getScore + !! + !! Resets score memory to having no scores, cycles, or batches + !! + subroutine reset(self) + class(scoreMemory), intent(inout) :: self + + self % bins = ZERO + self % parallelBins = ZERO + self % batchN = 0 + self % cycles = 0 + + end subroutine reset + end module scoreMemory_class diff --git a/Tallies/tallyAdmin_class.f90 b/Tallies/tallyAdmin_class.f90 index feba84a6d..839c5d111 100644 --- a/Tallies/tallyAdmin_class.f90 +++ b/Tallies/tallyAdmin_class.f90 @@ -67,17 +67,19 @@ module tallyAdmin_class !! push -> Push new attachment to the end of the list !! pop -> Remove attachment from the end of the list !! getEnd -> Returns pointer to the ned of the list (does not remove it like pop) - !! reportInColl -> Process pre-collision reports in all clerks - !! reportOutColl -> Process post-collision reports in all clerks - !! reportPath -> Process pathlength reports in all clerks - !! reportTrans -> Process transition reports in all clerks - !! reportSpawn -> Process particle generation reports in all clerks - !! reportHist -> Process history reports in all clerks - !! reportCycleStart -> Process start of cycle reports in all clerks - !! reportCycleEnd -> Process end of cycle reports in all clerks - !! getResult -> Return tallyResult object from a named Clerk + !! reportInColl -> Process pre-collision reports in all clerks + !! reportOutColl -> Process post-collision reports in all clerks + !! reportPath -> Process pathlength reports in all clerks + !! reportTrans -> Process transition reports in all clerks + !! reportSpawn -> Process particle generation reports in all clerks + !! reportHist -> Process history reports in all clerks + !! reportCycleStart -> Process start of cycle reports in all clerks + !! reportCycleEnd -> Process end of cycle reports in all clerks + !! getResult -> Return tallyResult object from a named Clerk + !! collectDistributed -> Sums tally results from all ranks into the master + !! resetMemory -> Resets score memory for all tallies !! display -> Call "display" on all clerks registered to display - !! isConverged -> Return .true. if all convergance targets have been reached + !! isConverged -> Return .true. if all convergence targets have been reached !! print -> Prints results to an output file object !! !! SAMPLE DICTIOANRY INPUT: @@ -151,6 +153,7 @@ module tallyAdmin_class ! Interaction procedures procedure :: getResult procedure :: collectDistributed + procedure :: resetMemory ! Display procedures procedure :: display @@ -410,7 +413,7 @@ end subroutine display !! None !! !! Result: - !! True is all convergance targets have been reached. False otherwise. + !! True if all convergence targets have been reached. False otherwise. !! !! Errors: !! None @@ -893,4 +896,14 @@ subroutine addToReports(self, reportCode, idx) end subroutine addToReports + !! + !! Resets memory to a fresh state + !! + subroutine resetMemory(self) + class(tallyAdmin),intent(inout) :: self + + call self % mem % reset() + + end subroutine resetMemory + end module tallyAdmin_class From 4e99450252734d100154d563e50c790a68b55c92 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Mon, 2 Mar 2026 17:38:44 +0000 Subject: [PATCH 03/15] Mostly there --- CMakeLists.txt | 1 + Coupling/CMakeLists.txt | 2 +- Coupling/couplingAdmin_class.f90 | 74 ++++++++++++------- PhysicsPackages/eigenPhysicsPackage_class.f90 | 18 ++++- SharedModules/timer_mod.f90 | 14 ++++ 5 files changed, 78 insertions(+), 31 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2861bb10d..b140c3b3f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -135,6 +135,7 @@ add_subdirectory(UserInterface) add_subdirectory(PhysicsPackages) add_subdirectory(DataStructures) +add_subdirectory(Coupling) #################################################################################################### # Compile SCONE static library diff --git a/Coupling/CMakeLists.txt b/Coupling/CMakeLists.txt index fa0bfa5bd..bf93af224 100644 --- a/Coupling/CMakeLists.txt +++ b/Coupling/CMakeLists.txt @@ -2,5 +2,5 @@ add_sources( ./couplingAdmin_class.f90) # Add Unit Tests to a global List -add_unit_tests( ./Tests/couplingAdmin_test.f90) +#add_unit_tests( ./Tests/couplingAdmin_test.f90) diff --git a/Coupling/couplingAdmin_class.f90 b/Coupling/couplingAdmin_class.f90 index 1d3e3988b..6a4f85457 100644 --- a/Coupling/couplingAdmin_class.f90 +++ b/Coupling/couplingAdmin_class.f90 @@ -1,18 +1,23 @@ module couplingAdmin_class use numPrecision - use universalVariables only : nameDensity, nameTemperature - use dictionary_class, only : dictionary - use dictParser_func, only : fileToDict - use outputFile_class, only : outputFile - use errors_mod, only : fatalError + use universalVariables, only : nameDensity, nameTemperature + use dictionary_class, only : dictionary + use dictParser_func, only : fileToDict + use outputFile_class, only : outputFile + use errors_mod, only : fatalError + use genericProcedures, only : numToChar + use timer_mod, only : c_sleep - use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + use field_inter, only : field + use tallyAdmin_class, only : tallyAdmin + use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr implicit none private - character(nameLen), parameter, dimension(2) :: ALLOWABLE_FIELDS = ['temperature', 'density'] + character(nameLen), parameter, dimension(2) :: ALLOWABLE_FIELDS = ['temperature',& + 'density '] character(nameLen), parameter :: END_SIGNAL = "SIGUSR2", & CONTINUE_SIGNAL = "SIGUSR1" @@ -38,6 +43,7 @@ module couplingAdmin_class !! maxIt -> maximum number of iterations before finishing coupling !! tally -> Tally admin containing necessary tallies for the driver !! fieldPaths -> paths to files which SCONE reads to update the given fields + !! fieldNames -> names of fields which are updated by other physics !! !! Public interface: !! init -> initialises the coupling setup from a dictionary @@ -46,9 +52,11 @@ module couplingAdmin_class !! endCoupling -> deactivates coupling !! !! doCoupling -> returns a logical if coupling is happening + !! attachTally -> attach coupling tallies to another tallyAdmin used by a physics package + !! + !! Private procedures: !! doUpdate -> returns a logical if a physics update should be performed on a given iteration. !! - !! makeTallies -> returns a tally admin used for coupling !! outputTallies -> outputs tally data to be read by other physics !! updateFields -> updates the relevant physics fields before resuming a calculation !! @@ -65,8 +73,9 @@ module couplingAdmin_class character(nameLen) :: outputFormat = '' integer(shortInt) :: updateFreq = 0 integer(shortInt) :: maxIt = huge(shortInt) - type(tallyAdmin) :: tally + type(tallyAdmin), pointer :: tally character(pathLen), dimension(:), allocatable :: fieldPaths + character(nameLen), dimension(:), allocatable :: fieldNames contains procedure :: init procedure :: kill @@ -79,7 +88,7 @@ module couplingAdmin_class ! Results handling and manipulation procedure :: attachTally - procedure :: outputTallies + procedure, private :: outputTallies procedure, private :: updateFields ! Communication @@ -97,9 +106,9 @@ module couplingAdmin_class subroutine init(self, dict) class(couplingAdmin), intent(inout) :: self class(dictionary), intent(in) :: dict - integer(shortInt) :: i, nFields + integer(shortInt) :: i, nFields, j logical(defBool) :: found, duplicate - class(dictionary), pointer :: tempDict + class(dictionary), pointer :: tallyDict type(outputFile) :: test_out character(100), parameter :: Here ='init (couplingAdmin_class.f90)' @@ -119,6 +128,7 @@ subroutine init(self, dict) ! Read tally info needed by other solvers tallyDict => dict % getDictPtr('tally') + allocate(self % tally) call self % tally % init(tallyDict) ! Read the output file name @@ -142,11 +152,19 @@ subroutine init(self, dict) nFields = size(self % fieldNames) do i = 1, nFields - found = any(trim(self % fieldNames(i)) == trim(ALLOWABLE_FIELDS)) - if (.not. found) call fatalError(Here,'Field names must be one of: '//ALLOWABLE_FIELDS) + found = .false. + do j = 1, size(ALLOWABLE_FIELDS) + found = found .or. (trim(self % fieldNames(i)) == trim(ALLOWABLE_FIELDS(j))) + end do + + if (.not. found) then + print '(A)', "ALLOWABLE FIELDS:" + print '(A)', ALLOWABLE_FIELDS + call fatalError(Here,'Field names is invalid. See list above.') + end if if (i < nFields) then - duplicate = any(trim(self % fieldNames(i)) == trim(self % fieldNames(i+1:nFields))) + duplicate = any(trim(self % fieldNames(i)) == [ (trim(self % fieldNames(j)), j=i+1,nFields) ]) if (duplicate) call fatalError(Here,'Field names must be unique: '//trim(self % fieldNames(i))) end if @@ -168,7 +186,8 @@ subroutine kill(self) self % outputFormat = '' self % maxIt = huge(shortInt) call self % tally % kill() - if (allocated(fieldPaths)) deallocate(fieldPaths) + if (allocated(self % fieldPaths)) deallocate(self % fieldPaths) + if (allocated(self % fieldNames)) deallocate(self % fieldNames) end subroutine kill @@ -218,7 +237,7 @@ subroutine endCoupling(self) class(couplingAdmin), intent(inout) :: self if (self % isCoupled) then - self % doCoupling = .false. + self % isCoupled = .false. call self % signalCalculationOver() end if @@ -243,7 +262,7 @@ end function doUpdate subroutine updateFields(self) class(couplingAdmin), intent(in) :: self integer(shortInt) :: i - class(dictionary) :: dict + class(dictionary), allocatable :: dict character(nameLen) :: fieldName class(field), pointer :: ptr character(100), parameter :: Here = 'updateFields (couplingAdmin_class.f90)' @@ -304,10 +323,11 @@ end subroutine signalCalculationOver !! of a signal file. !! subroutine waitForSignal(self) - class(couplingAdmin), intent(in) :: self - logical(defBool) :: exists - integer(shortInt) :: unit - character(*) :: line + class(couplingAdmin), intent(inout) :: self + logical(defBool) :: exists + integer(shortInt) :: unit + character(nameLen) :: line + character(100), parameter :: Here ='waitForSignal (couplingAdmin_class.f90)' do @@ -328,10 +348,10 @@ subroutine waitForSignal(self) end if close(unit, status="delete") - exit + return end if - call sleep(1) + call c_sleep(1) end do @@ -341,10 +361,10 @@ end subroutine waitForSignal !! Attach coupling tallies to another tally !! subroutine attachTally(self, tallyPtr) - class(couplingAdmin), intent(in) :: self - class(tallyAdmin), pointer, intent(inout) :: tallyPtr + class(couplingAdmin), intent(in) :: self + type(tallyAdmin), pointer, intent(inout) :: tallyPtr - if (self % doCoupling()) + if (self % doCoupling()) then call tallyPtr % push(self % tally) end if diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index bd88e96a8..f5aab91fc 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -65,6 +65,9 @@ module eigenPhysicsPackage_class use tallyResult_class, only : tallyResult use keffAnalogClerk_class, only : keffResult + ! Coupling + use couplingAdmin_class, only : couplingAdmin + ! Factories use transportOperatorFactory_func, only : new_transportOperator @@ -107,6 +110,9 @@ module eigenPhysicsPackage_class integer(shortInt) :: bufferSize logical(defBool) :: UFS = .false. logical(defBool) :: reproducible = .true. + + ! Coupling handler + type(couplingAdmin) :: couplingInfo ! Calculation components type(particleDungeon), pointer :: thisCycle => null() @@ -144,6 +150,10 @@ subroutine run(self) call self % generateInitialState() call self % cycles(self % inactiveTally, self % inactiveAtch, self % N_inactive) + + ! Deactivate coupling for active cycles + call self % couplingInfo % endCoupling() + call self % cycles(self % activeTally, self % activeAtch, self % N_active) ! Collect results from other processes @@ -334,6 +344,9 @@ subroutine cycles(self, tally, tallyAtch, N_cycles) call statusMsg("End time: " // trim(secToChar(end_T))) call statusMsg("Time to end: " // trim(secToChar(T_toEnd))) call tally % display() + + ! Perform coupling operations + call self % couplingInfo % couple(i) end do @@ -591,7 +604,6 @@ subroutine init(self, dict) allocate(self % inactiveTally) call self % inactiveTally % init(tempDict) - tempDict => dict % getDictPtr('activeTally') allocate(self % activeTally) call self % activeTally % init(tempDict) @@ -647,8 +659,8 @@ subroutine init(self, dict) call self % activeTally % push(self % activeAtch) ! Attach a tally admin for coupling - if (self % couplingInfo % isCoupled()) then - + if (self % couplingInfo % doCoupling()) then + call self % couplingInfo % attachTally(self % inactiveTally) end if call self % printSettings() diff --git a/SharedModules/timer_mod.f90 b/SharedModules/timer_mod.f90 index c87e57ba8..eac4f16c7 100644 --- a/SharedModules/timer_mod.f90 +++ b/SharedModules/timer_mod.f90 @@ -22,11 +22,14 @@ !! timerReset -> Reset a defined timer !! timerTime -> Get total elapsed time in Seconds !! +!! c_sleep -> pauses SCONE for several seconds +!! !! module timer_mod use numPrecision use genericProcedures, only : fatalError + use iso_c_binding, only : c_int implicit none private @@ -38,6 +41,7 @@ module timer_mod public :: timerStop public :: timerReset public :: timerTime + public :: c_sleep !! !! This derived type allows to measure time in a program @@ -77,6 +81,15 @@ module timer_mod real(defReal),parameter :: GROWTH_RATIO = 1.6_defReal integer(shortInt),parameter :: MIN_SIZE = 5 + ! Interface to C sleep procedure + interface + subroutine c_sleep(seconds) bind(C, name="sleep") + import :: c_int + integer(c_int), value :: seconds + end subroutine + end interface + + contains !! @@ -338,4 +351,5 @@ function toSec_stopWatch(self) result(sec) end function toSec_stopWatch + end module timer_mod From c3683540a978082e3098b3a62c09d3201a0ea2b1 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Tue, 3 Mar 2026 20:05:31 +0000 Subject: [PATCH 04/15] Fixes --- .../Reactions/reactionMG/fissionMG_class.f90 | 6 ++- SharedModules/endfConstants.f90 | 2 +- Tallies/TallyClerks/Tests/mgXsClerk_test.f90 | 18 ++++--- Tallies/TallyClerks/mgXsClerk_class.f90 | 49 +++++++++++++------ .../Tests/macroResponse_test.f90 | 2 +- .../Tests/microResponse_test.f90 | 2 +- 6 files changed, 50 insertions(+), 29 deletions(-) diff --git a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 index 47d41d89f..ed7bd9d9c 100644 --- a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 +++ b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 @@ -18,6 +18,8 @@ module fissionMG_class !! public :: fissionMG_TptrCast + real(defReal), private, parameter :: KAPPA_DEFAULT = 202.27 ! [MeV] + !! !! A special type of MG reaction that contains data related to fission !! @@ -36,7 +38,7 @@ module fissionMG_class !! type, public, extends(reactionMG) :: fissionMG real(defReal),dimension(:,:),allocatable :: data - real(defReal) :: kappa = 200.0_defReal + real(defReal) :: kappa = KAPPA_DEFAULT contains ! Superclass procedures procedure :: init @@ -222,7 +224,7 @@ subroutine buildFromDict(self, dict) call dict % get(nG, 'numberOfGroups') ! Get energy per fission - call dict % getOrDefault(self % kappa, 'kappa', 200.0_defReal) + call dict % getOrDefault(self % kappa, 'kappa', KAPPA_DEFAULT) ! Allocate space allocate(self % data(nG, 2)) diff --git a/SharedModules/endfConstants.f90 b/SharedModules/endfConstants.f90 index ff3b1ad34..de8ee3dba 100644 --- a/SharedModules/endfConstants.f90 +++ b/SharedModules/endfConstants.f90 @@ -120,7 +120,7 @@ module endfConstants macroFission = -6 ,& macroNuFission = -7 ,& macroPromptNuFission = -8 ,& - macroDelayedNuFission = -9 + macroDelayedNuFission = -9 ,& macroKappaFission = -80 ! List of Macro MT numbers for macroscopic XSs. Unique to SCONE (not from Serpent) diff --git a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 index d14e284c4..d193cdb13 100644 --- a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 +++ b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 @@ -74,7 +74,7 @@ subroutine setUp(this) ! Build test neutronDatabase call this % nucData % build(ONE, captureXS = 2.0_defReal, & - fissionXS = 1.5_defReal, nuFissionXS = 3.0_defReal) + fissionXS = 1.5_defReal, nuFissionXS = 3.0_defReal, kappaXS = 300.0_defReal) end subroutine setUp @@ -105,8 +105,8 @@ subroutine testScoring_clerk1(this) type(particleState) :: pFiss type(outputFile) :: out real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & - nu, chi, P0, P1, P2, P3, P4, & - P5, P6, P7, prod + nu, chi, kappa, P0, P1, P2, P3, & + P4, P5, P6, P7, prod real(defReal), parameter :: TOL = 1.0E-9 ! Configure memory @@ -143,7 +143,7 @@ subroutine testScoring_clerk1(this) call mem % closeCycle(ONE) ! Process and get results - call this % clerk_test1 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + call this % clerk_test1 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) call this % clerk_test1 % processPN(mem, P2, P3, P4, P5, P6, P7) ! Verify results of scoring @@ -151,6 +151,7 @@ subroutine testScoring_clerk1(this) @assertEqual([ZERO, ZERO, 1.5_defReal, ZERO], fiss(1,:), TOL, 'Fission XS' ) @assertEqual([ZERO, ZERO, TWO, ZERO], nu(1,:), TOL, 'NuFission XS' ) @assertEqual([ZERO, ZERO, HALF, HALF], chi(1,:), TOL, 'Chi' ) + @assertEqual([ZERO, ZERO, 300.0_defReal, ZERO], kappa(1,:), TOL, 'Kappa XS' ) @assertEqual([ZERO, ZERO, 4.0_defReal, ZERO], transOS(1,:), TOL, 'Transport XS O.S.' ) @assertEqual([ZERO, ZERO, 5.5_defReal, ZERO], transFL(1,:), TOL, 'Transport XS F.L.' ) @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, TWO, ZERO, ZERO], P0(1,:), TOL, 'P0' ) @@ -165,7 +166,7 @@ subroutine testScoring_clerk1(this) @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, -0.0683670044_defReal, ZERO, ZERO], P7(1,:), TOL, 'P7' ) ! Test getting size - @assertEqual(100, this % clerk_test1 % getSize(),'Test getSize():') + @assertEqual(104, this % clerk_test1 % getSize(),'Test getSize():') ! Test correctness of output calls call out % init('dummyPrinter', fatalErrors = .false.) @@ -185,7 +186,7 @@ subroutine testScoring_clerk2(this) type(particleState) :: pFiss type(outputFile) :: out real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & - nu, chi, P0, P1, prod + nu, chi, kappa, P0, P1, prod real(defReal), parameter :: TOL = 1.0E-9 ! Configure memory @@ -221,13 +222,14 @@ subroutine testScoring_clerk2(this) call mem % closeCycle(ONE) ! Process and get results - call this % clerk_test2 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + call this % clerk_test2 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) ! Verify results of scoring @assertEqual([ZERO, ZERO, TWO], capt(1,:), TOL, 'Capture XS' ) @assertEqual([ZERO, ZERO, 1.5_defReal], fiss(1,:), TOL, 'Fission XS' ) @assertEqual([ZERO, ZERO, TWO], nu(1,:), TOL, 'NuFission XS' ) @assertEqual([HALF, ZERO, HALF], chi(1,:), TOL, 'Chi' ) + @assertEqual([ZERO, ZERO, 300.0_defReal], kappa(1,:), TOL, 'Kappa XS' ) @assertEqual([ZERO, ZERO, 4.0_defReal], transOS(1,:), TOL, 'Transport XS O.S.' ) @assertEqual([ZERO, ZERO, 5.5_defReal], transFL(1,:), TOL, 'Transport XS F.L.' ) @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, TWO, ZERO], P0(1,:), TOL, 'P0' ) @@ -235,7 +237,7 @@ subroutine testScoring_clerk2(this) @assertEqual([ONE, ONE, ONE, ONE, ONE, ONE, ONE, TWO, ONE], prod(1,:), TOL, 'prod' ) ! Test getting size - @assertEqual(48, this % clerk_test2 % getSize(),'Test getSize():') + @assertEqual(51, this % clerk_test2 % getSize(),'Test getSize():') ! Test correctness of output calls call out % init('dummyPrinter', fatalErrors = .false.) diff --git a/Tallies/TallyClerks/mgXsClerk_class.f90 b/Tallies/TallyClerks/mgXsClerk_class.f90 index a65d4aaf0..567cbc1e6 100644 --- a/Tallies/TallyClerks/mgXsClerk_class.f90 +++ b/Tallies/TallyClerks/mgXsClerk_class.f90 @@ -29,7 +29,7 @@ module mgXsClerk_class private !! Size of clerk memory - integer(shortInt), parameter :: ARRAY_SCORE_SIZE = 7 ,& ! Size of data to store as 1D arrays + integer(shortInt), parameter :: ARRAY_SCORE_SIZE = 8 ,& ! Size of data to store as 1D arrays MATRIX_SCORE_SMALL = 3 ,& ! Size of data to store as 2D arrays when scoring up to P1 MATRIX_SCORE_FULL = 9 ! Size of data to store as 2D arrays when scoring up to P7 @@ -40,14 +40,16 @@ module mgXsClerk_class FISS_idx = 4 ,& ! Fission macroscopic reaction rate NUBAR_idx = 5 ,& ! NuBar CHI_idx = 6 ,& ! Fission neutron spectrum - SCATT_EV_idx = 7 ! Analog: number of scattering events + KAPPAFISS_idx = 7 ,& ! Fission heating macroscopic reaction rate + SCATT_EV_idx = 8 ! Analog: number of scattering events !! !! Multi-group macroscopic cross section calculation !! !! It prints out: - !! capture xs, fission xs, transport xs, nu, chi, the P0 and P1 scattering matrices, - !! the P0 scattering production matrix. On request, also the P2 -> P7 scattering matrices. + !! capture xs, fission xs, transport xs, nu, chi, kappa * fission xs, the P0 and P1 scattering + !! matrices, the P0 scattering production matrix. On request, also the P2 -> P7 scattering + !! matrices. !! !! NOTE: !! - the cross sections are tallied with a collision estimator; @@ -230,7 +232,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) type(particleState) :: state type(neutronMacroXSs) :: xss class(neutronMaterial), pointer :: mat - real(defReal) :: nuFissXS, captXS, fissXS, scattXS, flux + real(defReal) :: nuFissXS, captXS, fissXS, scattXS, kappaXS, flux integer(shortInt) :: enIdx, locIdx, binIdx integer(longInt) :: addr character(100), parameter :: Here =' reportInColl (mgXsClerk_class.f90)' @@ -288,6 +290,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) captXS = xss % capture * flux fissXS = xss % fission * flux scattXS = (xss % elasticScatter + xss % inelasticScatter) * flux + kappaXS = xss % kappaXS * flux else @@ -296,6 +299,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) captXS = ZERO fissXS = ZERO scattXS = ZERO + kappaXS = ZERO end if @@ -305,6 +309,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) call mem % score(captXS, addr + CAPT_idx) call mem % score(fissXS, addr + FISS_idx) call mem % score(scattXS, addr + SCATT_idx) + call mem % score(kappaXS, addr + KAPPAFISS_idx) end subroutine reportInColl @@ -500,7 +505,7 @@ end subroutine reportSpawn !! none !! pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_res, & - nu_res, chi_res, P0_res, P1_res, prod_res) + nu_res, chi_res, kappa_res, P0_res, P1_res, prod_res) class(mgXsClerk), intent(in) :: self type(scoreMemory), intent(in) :: mem real(defReal), dimension(:,:), allocatable, intent(out) :: capt_res @@ -509,6 +514,7 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r real(defReal), dimension(:,:), allocatable, intent(out) :: transOS_res real(defReal), dimension(:,:), allocatable, intent(out) :: nu_res real(defReal), dimension(:,:), allocatable, intent(out) :: chi_res + real(defReal), dimension(:,:), allocatable, intent(out) :: kappa_res real(defReal), dimension(:,:), allocatable, intent(out) :: P0_res real(defReal), dimension(:,:), allocatable, intent(out) :: P1_res real(defReal), dimension(:,:), allocatable, intent(out) :: prod_res @@ -518,7 +524,7 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r integer(shortInt) :: N, M, i, j, k, g1, gEnd, idx real(defReal) :: capt, fiss, scatt, nu, chi, P0, P1, prod, sumChi, flux, scattProb, & captStd, fissStd, scattStd, nuStd, chiStd, P0std, P1std, prodStd, & - fluxStd, scattProbStd, scattXS, scattXSstd + fluxStd, scattProbStd, scattXS, scattXSstd, kappa, kappaStd ! Get number of bins N = self % energyN @@ -526,9 +532,9 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r ! Allocate arrays for MG xss allocate( capt_res(2,N*M), fiss_res(2,N*M), transFL_res(2,N*M), transOS_res(2,N*M), & - nu_res(2,N*M), chi_res(2,N*M), P0_res(2,N*N*M), P1_res(2,N*N*M), & - prod_res(2,N*N*M), tot(N*M), fluxG(N*M), delta(M, N), totStd(N*M), & - fluxGstd(N*M), deltaStd(M, N) ) + nu_res(2,N*M), chi_res(2,N*M), kappa_res(2,N*M), P0_res(2,N*N*M), & + P1_res(2,N*N*M), prod_res(2,N*N*M), tot(N*M), fluxG(N*M), delta(M, N), & + totStd(N*M), fluxGstd(N*M), deltaStd(M, N) ) ! Initialise values sumChi = 0 ! to normalise chi @@ -548,6 +554,7 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r call mem % getResult(nu, nuStd, addr + NUBAR_idx) call mem % getResult(chi, chiStd, addr + CHI_idx) call mem % getResult(scattProb, scattProbStd, addr + SCATT_EV_idx) + call mem % getResult(kappa, kappaStd, addr + KAPPAFISS_idx) ! Calculate MG constants, being careful to avoid division by zero ! If flux is zero all cross sections must be set to zero @@ -564,13 +571,16 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r ! Calculate fission production term and uncertainties if (fiss == ZERO) then - fiss_res(1:2,i) = ZERO - nu_res(1:2,i) = ZERO + fiss_res(1:2,i) = ZERO + nu_res(1:2,i) = ZERO + kappa_res(1:2,i) = ZERO else fiss_res(1,i) = fiss/flux fiss_res(2,i) = fiss_res(1,i) * sqrt((fissStd/fiss)**2 + (fluxStd/flux)**2) nu_res(1,i) = nu/fiss nu_res(2,i) = nu_res(1,i) * sqrt((nuStd/nu)**2 + (fissStd/fiss)**2) + kappa_res(1,i) = kappa/flux + kappa_res(2,i) = kappa_res(1,i) * sqrt((kappaStd/kappa)**2 + (fluxStd/flux)**2) end if ! Store fission spectrum @@ -758,8 +768,8 @@ subroutine print(self, outFile, mem) type(scoreMemory), intent(in) :: mem integer(shortInt),dimension(:),allocatable :: resArrayShape real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & - nu, chi, P0, P1, P2, P3, P4, & - P5, P6, P7, prod + nu, chi, kappa, P0, P1, P2, P3, & + P4, P5, P6, P7, prod character(nameLen) :: name integer(shortInt) :: i @@ -784,7 +794,7 @@ subroutine print(self, outFile, mem) end if ! Process and get results - call self % processRes(mem, capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + call self % processRes(mem, capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) ! Print results name = 'capture' @@ -801,6 +811,13 @@ subroutine print(self, outFile, mem) end do call outFile % endArray() + name = 'kappaFission' + call outFile % startArray(name, resArrayShape) + do i=1,product(resArrayShape) + call outFile % addResult(kappa(1,i),kappa(2,i)) + end do + call outFile % endArray() + name = 'transportFluxLimited' call outFile % startArray(name, resArrayShape) do i=1,product(resArrayShape) @@ -853,7 +870,7 @@ subroutine print(self, outFile, mem) call outFile % endArray() ! Deallocate to limit memory consumption when writing to the output file - deallocate(capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + deallocate(capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) ! If high order scattering is requested, print the other matrices if (self % PN) then diff --git a/Tallies/TallyResponses/Tests/macroResponse_test.f90 b/Tallies/TallyResponses/Tests/macroResponse_test.f90 index 103c61007..69c3815f0 100644 --- a/Tallies/TallyResponses/Tests/macroResponse_test.f90 +++ b/Tallies/TallyResponses/Tests/macroResponse_test.f90 @@ -37,7 +37,7 @@ subroutine setUp(this) type(dictionary) :: tempDict ! Allocate and initialise test nuclearData - ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappFission + ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappaFission call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal, 9.0_defReal) ! Set up responses diff --git a/Tallies/TallyResponses/Tests/microResponse_test.f90 b/Tallies/TallyResponses/Tests/microResponse_test.f90 index 8f96e8adf..462fac181 100644 --- a/Tallies/TallyResponses/Tests/microResponse_test.f90 +++ b/Tallies/TallyResponses/Tests/microResponse_test.f90 @@ -39,7 +39,7 @@ subroutine setUp(this) ! Allocate and initialise test nuclearData - ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappFission + ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappaFission call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal, 9.0_defReal) ! Set dictionaries to initialise material From 221eee7fe1bb10e23646e30f112a04de0a956f11 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Sat, 14 Mar 2026 22:03:53 +0000 Subject: [PATCH 05/15] Finishing off coupling --- Coupling/CMakeLists.txt | 5 +- Coupling/Tests/couplingAdmin_iTest.f90 | 180 ++++++++++++++++++ Coupling/couplingAdmin_class.f90 | 78 ++++++-- .../cartesianField_class.f90 | 3 +- Geometry/Tests/geometryStd_iTest.f90 | 56 ++++++ Geometry/geometryReg_mod.f90 | 39 ++-- Geometry/geometryStd_class.f90 | 4 +- .../Benchmarks/Transients/LatticeBrunner | 2 +- IntegrationTestFiles/Coupling/test_coupling | 19 ++ IntegrationTestFiles/Geometry/test_lat | 11 +- IntegrationTestFiles/dens | 7 + IntegrationTestFiles/temp | 7 + .../aceDatabase/aceNeutronDatabase_class.f90 | 75 ++++++-- .../ceNeutronData/ceNeutronCache_mod.f90 | 11 +- .../ceNeutronData/ceNeutronDatabase_inter.f90 | 7 +- .../ceNeutronData/ceNeutronMaterial_class.f90 | 21 +- PhysicsPackages/alphaPhysicsPackage_class.f90 | 46 ++++- PhysicsPackages/eigenPhysicsPackage_class.f90 | 9 +- SharedModules/universalVariables.f90 | 8 +- .../transportOperatorHT_class.f90 | 4 +- .../transportOperatorST_class.f90 | 4 +- docs/Input Manual.rst | 65 ++++++- 22 files changed, 566 insertions(+), 95 deletions(-) create mode 100644 Coupling/Tests/couplingAdmin_iTest.f90 create mode 100644 IntegrationTestFiles/Coupling/test_coupling create mode 100644 IntegrationTestFiles/dens create mode 100644 IntegrationTestFiles/temp diff --git a/Coupling/CMakeLists.txt b/Coupling/CMakeLists.txt index bf93af224..044952249 100644 --- a/Coupling/CMakeLists.txt +++ b/Coupling/CMakeLists.txt @@ -1,6 +1,3 @@ -# Add Source Files to a global List add_sources( ./couplingAdmin_class.f90) -# Add Unit Tests to a global List -#add_unit_tests( ./Tests/couplingAdmin_test.f90) - +add_integration_tests( ./Tests/couplingAdmin_iTest.f90) diff --git a/Coupling/Tests/couplingAdmin_iTest.f90 b/Coupling/Tests/couplingAdmin_iTest.f90 new file mode 100644 index 000000000..d0a2f632b --- /dev/null +++ b/Coupling/Tests/couplingAdmin_iTest.f90 @@ -0,0 +1,180 @@ +module couplingAdmin_iTest + + use numPrecision + use universalVariables + use dictionary_class, only : dictionary + use dictParser_func, only : fileToDict + use fieldFactory_func, only : new_field + use couplingAdmin_class, only : couplingAdmin + use funit + + implicit none + +contains + + !! + !! Coupling admin integration test + !! +@Test + subroutine test_coupling() + type(couplingAdmin) :: coupler + character(pathLen), parameter :: path = './IntegrationTestFiles/Coupling/test_coupling' + character(pathLen), parameter :: sendPath = './IntegrationTestFiles/Coupling/toDriver.dat' + character(pathLen), parameter :: recvPath = './IntegrationTestFiles/Coupling/toSCONE.dat' + character(pathLen), parameter :: tallyPath = './IntegrationTestFiles/Coupling/sconeOut.m' + type(dictionary) :: dict + integer(shortInt) :: it, unit + logical(defBool) :: exists, correctFile + character(nameLen) :: line + + ! Create a recv file to prevent SCONE from hanging + open(newunit=unit, file=recvPath, status="replace", action="write") + write(unit, '(A)') 'SIGUSR1' + close(unit) + + ! Attempt a coupling update before initialisation: no files should be produced + ! or deleted + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! Initialise + call fileToDict(dict, path) + call coupler % init(dict) + + ! Feed in different iteration numbers and see whether files are produced + ! This one should not result in any information exchange + it = 1 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! This one should be coupled + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertTrue(exists) + inquire(file=sendPath, exist=exists) + @assertTrue(exists) + + ! Send file should contain SIGUSR1 + open(newunit=unit, file=sendPath, status='old', action='read') + read(unit,'(A)') line + close(unit) + correctFile = (trim(line) == 'SIGUSR1') + @assertTrue(correctFile) + + ! Reset the files + open(newunit=unit, file=tallyPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=sendPath, status='old') + close(unit, status='delete') + + ! This one should end coupling due to hitting the maximum iteration + ! No recv file is necessary + it = 10 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertTrue(exists) + inquire(file=sendPath, exist=exists) + @assertTrue(exists) + + ! Send file should contain SIGTERM + open(newunit=unit, file=sendPath, status='old', action='read') + read(unit,'(A)') line + close(unit) + correctFile = (trim(line) == 'SIGTERM') + @assertTrue(correctFile) + + ! Reset the files + open(newunit=unit, file=tallyPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=sendPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=recvPath, status="replace", action="write") + write(unit, '(A)') 'SIGTERM' + close(unit) + + ! Coupling should not occur + it = 15 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! Reinitialise the coupler, ensuring that no coupling occurs successfully + ! until reinitialisation + call coupler % kill() + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + call coupler % init(dict) + + ! Given the recv file, coupling should occur on this iteration, but not on the 10th + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertTrue(exists) + inquire(file=sendPath, exist=exists) + @assertTrue(exists) + + ! Send file should contain SIGUSR1 + open(newunit=unit, file=sendPath, status='old', action='read') + read(unit,'(A)') line + close(unit) + correctFile = (trim(line) == 'SIGUSR1') + @assertTrue(correctFile) + + ! Reset the files - no need for the recv file + open(newunit=unit, file=tallyPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=sendPath, status='old') + close(unit, status='delete') + + it = 10 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! Kill coupling admin + call coupler % kill() + + end subroutine test_coupling + +end module couplingAdmin_iTest diff --git a/Coupling/couplingAdmin_class.f90 b/Coupling/couplingAdmin_class.f90 index 6a4f85457..acede20d2 100644 --- a/Coupling/couplingAdmin_class.f90 +++ b/Coupling/couplingAdmin_class.f90 @@ -1,24 +1,28 @@ module couplingAdmin_class use numPrecision - use universalVariables, only : nameDensity, nameTemperature + use universalVariables, only : nameDensity, nameTemperature, NO_TEMPERATURE use dictionary_class, only : dictionary use dictParser_func, only : fileToDict use outputFile_class, only : outputFile use errors_mod, only : fatalError + use display_func, only : statusMsg use genericProcedures, only : numToChar use timer_mod, only : c_sleep - use field_inter, only : field - use tallyAdmin_class, only : tallyAdmin - use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + use field_inter, only : field + use fieldFactory_func, only : new_field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast + use tallyAdmin_class, only : tallyAdmin + use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr, & + gr_hasField => hasField implicit none private character(nameLen), parameter, dimension(2) :: ALLOWABLE_FIELDS = ['temperature',& 'density '] - character(nameLen), parameter :: END_SIGNAL = "SIGUSR2", & + character(nameLen), parameter :: END_SIGNAL = "SIGTERM", & CONTINUE_SIGNAL = "SIGUSR1" !! @@ -33,6 +37,7 @@ module couplingAdmin_class !! - create any tallies necessary for coupling (such as power tallies) !! - output the results of those tallies to an appropriate location !! - know the location of fields provided by other physics and update them in SCONE + !! - possibly provide checks on validity of information fed back to SCONE !! !! Private members: !! isCoupled -> flag for whether coupling is being performed @@ -41,6 +46,7 @@ module couplingAdmin_class !! outputFile -> path to the file which SCONE outputs containing results to be used by the driver !! updateFreq -> frequency with which SCONE pauses to exchange data with the driver !! maxIt -> maximum number of iterations before finishing coupling + !! maxTemperature -> Assumed maximum temperature present in the problem !! tally -> Tally admin containing necessary tallies for the driver !! fieldPaths -> paths to files which SCONE reads to update the given fields !! fieldNames -> names of fields which are updated by other physics @@ -51,11 +57,12 @@ module couplingAdmin_class !! couple -> performs the main coupling logic !! endCoupling -> deactivates coupling !! - !! doCoupling -> returns a logical if coupling is happening - !! attachTally -> attach coupling tallies to another tallyAdmin used by a physics package + !! doCoupling -> returns a logical if coupling is happening + !! attachTally -> attach coupling tallies to another tallyAdmin used by a physics package !! !! Private procedures: - !! doUpdate -> returns a logical if a physics update should be performed on a given iteration. + !! doUpdate -> returns a logical if a physics update should be performed on a given iteration. + !! maxTemperatureCheck -> checks whether the maximum field temperature is lower than that assumed initially !! !! outputTallies -> outputs tally data to be read by other physics !! updateFields -> updates the relevant physics fields before resuming a calculation @@ -73,6 +80,7 @@ module couplingAdmin_class character(nameLen) :: outputFormat = '' integer(shortInt) :: updateFreq = 0 integer(shortInt) :: maxIt = huge(shortInt) + real(defReal) :: maxTemperature = NO_TEMPERATURE type(tallyAdmin), pointer :: tally character(pathLen), dimension(:), allocatable :: fieldPaths character(nameLen), dimension(:), allocatable :: fieldNames @@ -84,6 +92,7 @@ module couplingAdmin_class ! Status checks procedure :: doCoupling + procedure, private :: checkMaxTemperature procedure, private :: doUpdate ! Results handling and manipulation @@ -112,6 +121,8 @@ subroutine init(self, dict) type(outputFile) :: test_out character(100), parameter :: Here ='init (couplingAdmin_class.f90)' + self % isCoupled = .true. + ! Read communication details call dict % get(self % commFileIn, 'receiveFile') call dict % get(self % commFileOut, 'sendFile') @@ -127,10 +138,14 @@ subroutine init(self, dict) 'Maximum iteration number is silly: '//numToChar(self % maxIt)) ! Read tally info needed by other solvers + ! Should this be modified to ensure MPI sync? tallyDict => dict % getDictPtr('tally') allocate(self % tally) call self % tally % init(tallyDict) + ! Read maximum temperature that may occur + call dict % getOrDefault(self % maxTemperature, 'maxTemp', NO_TEMPERATURE) + ! Read the output file name call dict % get(self % outputFile, 'outputFile') @@ -160,7 +175,7 @@ subroutine init(self, dict) if (.not. found) then print '(A)', "ALLOWABLE FIELDS:" print '(A)', ALLOWABLE_FIELDS - call fatalError(Here,'Field names is invalid. See list above.') + call fatalError(Here,'Field name is invalid. See list above.') end if if (i < nFields) then @@ -185,6 +200,7 @@ subroutine kill(self) self % outputFile = '' self % outputFormat = '' self % maxIt = huge(shortInt) + self % maxTemperature = NO_TEMPERATURE call self % tally % kill() if (allocated(self % fieldPaths)) deallocate(self % fieldPaths) if (allocated(self % fieldNames)) deallocate(self % fieldNames) @@ -219,6 +235,25 @@ subroutine couple(self, it) end subroutine couple + !! + !! Compares the maximum temperature field value with the assumed maximum + !! temperature given in the input. If this is exceeded, emits a warning + !! concerning the validity of the majorant and the resulting simulation + !! outputs. + !! + subroutine checkMaxTemperature(self, temp) + class(couplingAdmin), intent(in) :: self + real(defReal), intent(in) :: temp + + if (temp > self % maxTemperature) then + call statusMsg('WARNING: Coupled field temperature exceeds the assumed maximum temperature '& + //'for the simulation. If using delta tracking, results may be incorrect! '& + //'Assumed maximum: '//numToChar(self % maxTemperature)//'K. '& + //'Actual maximum: '//numToChar(temp)//'K.') + end if + + end subroutine checkMaxTemperature + !! !! Returns whether coupling should be occurring or not !! @@ -260,11 +295,13 @@ end function doUpdate !! Update the fields from the paths to the new field info !! subroutine updateFields(self) - class(couplingAdmin), intent(in) :: self - integer(shortInt) :: i - class(dictionary), allocatable :: dict - character(nameLen) :: fieldName - class(field), pointer :: ptr + class(couplingAdmin), intent(in) :: self + integer(shortInt) :: i + type(dictionary) :: dict + character(nameLen) :: fieldName + class(field), pointer :: ptr + class(pieceConstantField), pointer :: piecePtr + real(defReal) :: maxTemp character(100), parameter :: Here = 'updateFields (couplingAdmin_class.f90)' do i = 1, size(self % fieldPaths) @@ -283,11 +320,24 @@ subroutine updateFields(self) call fatalError(Here, 'Unrecognised field type requested') end select + if (.not. gr_hasField(fieldName)) then + call new_field(dict, fieldName) + end if + ptr => gr_fieldPtr(gr_fieldIdx(fieldName)) call ptr % kill() call ptr % init(dict) + ! Ensure assumed maximum temperature is not violated + ! For now: down cast to pieceConstantField + ! TODO: implement getMaxValue for all fields + if (fieldName == nameTemperature) then + piecePtr => pieceConstantField_CptrCast(ptr) + maxTemp = piecePtr % getMaxValue() + call self % checkMaxTemperature(maxTemp) + end if + end do end subroutine updateFields diff --git a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 index e1b109513..74a54fe58 100644 --- a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 +++ b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 @@ -130,7 +130,7 @@ subroutine init_dict(self, dict) ! Check for invalid pitch if (any(self % pitch < 10 * SURF_TOL)) then call fatalError(Here, 'Pitch size must be larger than: '//numToChar( 10 * SURF_TOL)) - end if + end if ! Calculate halfwidth and corner self % a_bar = self % pitch * HALF - SURF_TOL @@ -347,6 +347,7 @@ elemental subroutine kill(self) self % nMat = 0 self % corner = ZERO self % a_bar = ZERO + if (allocated(self % matIdxs)) deallocate(self % matIdxs) call self % outline % kill() end subroutine kill diff --git a/Geometry/Tests/geometryStd_iTest.f90 b/Geometry/Tests/geometryStd_iTest.f90 index 976e821ac..5a9e91eaf 100644 --- a/Geometry/Tests/geometryStd_iTest.f90 +++ b/Geometry/Tests/geometryStd_iTest.f90 @@ -8,6 +8,8 @@ module geometryStd_iTest use coord_class, only : coordList use geometryStd_class, only : geometryStd use materialMenu_mod, only : mm_init => init + use geometryReg_mod, only : gr_kill => kill + use fieldFactory_func, only : new_field use funit implicit none @@ -213,8 +215,62 @@ subroutine test_lattice_geom() @assertEqual(idx, coords % matIdx) @assertEqual(0.08_defReal, maxDist, TOL) + ! Add superimposed density field to complicate movement + tempDict => dict % getDictPtr('density') + call new_field(tempDict, nameDensity) + + ! Place a particle and move it, ensuring the event is as expected + r = [-0.64_defReal, -1.18_defReal, 0.0_defReal] + u = [ONE, ZERO, ZERO] + call coords % init(r, u) + call geom % placeCoord(coords) + maxDist = 10.0_defReal + + ! Cross the field + r_ref = [-0.63_defReal, -1.18_defReal, ZERO] + u_ref = [ONE, ZERO, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(FIELD_EV, event) + @assertEqual(0.01_defReal, maxDist, TOL) + + ! Cross a cell rather than hit the field + maxDist = 1.0_defReal + r_ref = [ZERO, -1.18_defReal, ZERO] + u_ref = [ONE, ZERO, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(CROSS_EV, event) + @assertEqual(0.63_defReal, maxDist, TOL) + + ! Collision close to the field crossing + maxDist = 0.629_defReal + r_ref = [0.629_defReal, -1.18_defReal, ZERO] + u_ref = [ONE, ZERO, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(COLL_EV, event) + @assertEqual(0.629_defReal, maxDist, TOL) + + ! Boundary crossing + u = [ZERO, -ONE, ZERO] + call coords % assignDirection(u) + maxDist = 2.0_defReal + + r_ref = [0.629_defReal, 1.26_defReal, ZERO] + u_ref = [ZERO, -ONE, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(BOUNDARY_EV, event) + @assertEqual(0.08_defReal, maxDist, TOL) + ! Kill geometry call geom % kill() + call gr_kill() end subroutine test_lattice_geom diff --git a/Geometry/geometryReg_mod.f90 b/Geometry/geometryReg_mod.f90 index a4ceb805b..f9107cdc8 100644 --- a/Geometry/geometryReg_mod.f90 +++ b/Geometry/geometryReg_mod.f90 @@ -43,6 +43,8 @@ module geometryReg_mod implicit none private + + integer(shortInt), parameter :: NO_GEOM = -7, NO_FIELD = -8 !! !! Small type to hold a single polymorphic geometry @@ -86,8 +88,8 @@ module geometryReg_mod type(charMap) :: fieldNameMap ! Save indices for speed. Avoids repeated calls to the character map - integer(shortInt) :: temperatureIdx = -8 - integer(shortInt) :: densityIdx = -8 + integer(shortInt) :: temperatureIdx = NO_FIELD + integer(shortInt) :: densityIdx = NO_FIELD contains @@ -107,13 +109,12 @@ subroutine addGeom(geom, name) class(geometry), allocatable, intent(inout) :: geom character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -7 character(100), parameter :: Here = 'addGeom (geometryReg_mod.f90)' ! Get free index - idx = geometryNameMap % getOrDefault(name, NOT_PRESENT) + idx = geometryNameMap % getOrDefault(name, NO_GEOM) - if (idx /= NOT_PRESENT) then + if (idx /= NO_GEOM) then call fatalError(Here, 'Geometry with name: '//trim(name)//' has already been defined with & &idx: '//numToChar(idx)) else @@ -145,11 +146,10 @@ end subroutine addGeom function geomIdx(name) result(idx) character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -8 character(*), parameter :: Here = 'geomIdx (geometryReg_mod.f90)' - idx = geometryNameMap % getOrDefault(name, NOT_PRESENT) - if (idx == NOT_PRESENT) then + idx = geometryNameMap % getOrDefault(name, NO_GEOM) + if (idx == NO_GEOM) then call fatalError(Here, 'Geometry: '//trim(name)//' was not found.') end if @@ -217,13 +217,12 @@ subroutine addField(kentta, name) class(field), allocatable, intent(inout) :: kentta character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -7 character(100), parameter :: Here = 'addField (geometryReg_mod.f90)' ! Get free index - idx = fieldNameMap % getOrDefault(name, NOT_PRESENT) + idx = fieldNameMap % getOrDefault(name, NO_FIELD) - if (idx /= NOT_PRESENT) then + if (idx /= NO_FIELD) then call fatalError(Here, 'Field with name: '//trim(name)//' has already been defined with & &idx: '//numToChar(idx)) else @@ -241,7 +240,7 @@ subroutine addField(kentta, name) elseif (name == nameDensity) then densityIdx = idx end if - + ! Point field call move_alloc(kentta, fields(idx) % kentta) @@ -259,7 +258,6 @@ end subroutine addField pure function hasField(name) result(exists) character(nameLen), intent(in) :: name logical(defBool) :: exists - integer(shortInt), parameter :: NOT_PRESENT = -8 integer(shortInt) :: idx select case(name) @@ -268,10 +266,10 @@ pure function hasField(name) result(exists) case(nameDensity) idx = densityIdx case default - idx = fieldNameMap % getOrDefault(name, NOT_PRESENT) + idx = fieldNameMap % getOrDefault(name, NO_FIELD) end select - exists = (idx /= NOT_PRESENT) + exists = (idx /= NO_FIELD) end function hasField @@ -290,11 +288,10 @@ end function hasField function fieldIdx(name) result(idx) character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -8 character(*), parameter :: Here = 'fieldIdx (geometryReg_mod.f90)' - idx = fieldNameMap % getOrDefault(name, NOT_PRESENT) - if (idx == NOT_PRESENT) then + idx = fieldNameMap % getOrDefault(name, NO_FIELD) + if (idx == NO_FIELD) then call fatalError(Here, 'Field: '//trim(name)//' was not found.') end if @@ -319,7 +316,6 @@ function fieldPtrName(name) result(ptr) character(nameLen), intent(in) :: name class(field), pointer :: ptr integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -8 character(*), parameter :: Here = 'fieldPtrName (geometryReg_mod.f90)' select case(name) @@ -328,7 +324,7 @@ function fieldPtrName(name) result(ptr) case(nameDensity) idx = densityIdx case default - idx = NOT_PRESENT + idx = NO_FIELD end select if (idx < 1 .or. idx > fieldTop) then @@ -394,6 +390,9 @@ subroutine kill() deallocate(fields) end if fieldTop = 0 + + temperatureIdx = NO_FIELD + densityIdx = NO_FIELD end subroutine kill diff --git a/Geometry/geometryStd_class.f90 b/Geometry/geometryStd_class.f90 index d3a0e03f4..f7fc3c1be 100644 --- a/Geometry/geometryStd_class.f90 +++ b/Geometry/geometryStd_class.f90 @@ -238,7 +238,7 @@ subroutine move_noCache(self, coords, maxDist, event) event = COLL_EV maxDist = maxDist ! Left for explicitness. Compiler will not stand it anyway - else if (maxDist < dist .and. maxDist >= fieldDist) then ! Stays within the same cell, but crosses field boundary + else if (fieldDist < dist - NUDGE) then ! Stays within the same cell, but crosses field boundary call coords % moveLocal(fieldDist, level0) event = FIELD_EV maxDist = fieldDist @@ -311,7 +311,7 @@ subroutine move_withCache(self, coords, maxDist, event, cache) maxDist = maxDist ! Left for explicitness. Compiler will not stand it anyway cache % lvl = 0 - else if (maxDist < dist .and. maxDist >= fieldDist) then ! Stays within the same cell, but crosses field boundary + else if (fieldDist < dist - NUDGE) then ! Stays within the same cell, but crosses field boundary call coords % moveLocal(fieldDist, level0) event = FIELD_EV maxDist = fieldDist diff --git a/InputFiles/Benchmarks/Transients/LatticeBrunner b/InputFiles/Benchmarks/Transients/LatticeBrunner index bcdb236ab..0e82a9d08 100644 --- a/InputFiles/Benchmarks/Transients/LatticeBrunner +++ b/InputFiles/Benchmarks/Transients/LatticeBrunner @@ -109,7 +109,7 @@ nuclearData { temp 12345; composition { } - xsFile ./XS/abs; + xsFile ./XS/air; } diff --git a/IntegrationTestFiles/Coupling/test_coupling b/IntegrationTestFiles/Coupling/test_coupling new file mode 100644 index 000000000..703382c07 --- /dev/null +++ b/IntegrationTestFiles/Coupling/test_coupling @@ -0,0 +1,19 @@ +tally { + norm totalPow; + normVal 70000; + totalPow { type collisionClerk; response (fission); fission { type macroResponse; MT -80;}} + axialPow {type collisionClerk; + map {type spaceMap; axis z; grid lin; min -200; max 200.0; N 40;} + response (fission); fission { type macroResponse; MT -80;} + } +} + +maxTemp 1100.0; +receiveFile ./IntegrationTestFiles/Coupling/toSCONE.dat; +sendFile ./IntegrationTestFiles/Coupling/toDriver.dat; +outputFile ./IntegrationTestFiles/Coupling/sconeOut; +updateFreq 5; +maxIt 10; +outputFormat asciiMATLAB; +fieldPaths (./IntegrationTestFiles/temp ./IntegrationTestFiles/dens); +fieldNames (temperature density); diff --git a/IntegrationTestFiles/Geometry/test_lat b/IntegrationTestFiles/Geometry/test_lat index df07f69ff..96a93bd9e 100644 --- a/IntegrationTestFiles/Geometry/test_lat +++ b/IntegrationTestFiles/Geometry/test_lat @@ -45,5 +45,12 @@ nuclearData { } } - - +density { +type cartesianField; +materials (all); +origin (0 0 0); +shape (4 4 0); +pitch (0.63 0.63 0); +default -1; +all ( 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1); +} diff --git a/IntegrationTestFiles/dens b/IntegrationTestFiles/dens new file mode 100644 index 000000000..8670240f7 --- /dev/null +++ b/IntegrationTestFiles/dens @@ -0,0 +1,7 @@ +type cartesianField; +materials (cool); +origin (0 0 0); +shape (1 1 0); +pitch (1.26 1.26 0); +default -1; +cool ( 0.9 ); diff --git a/IntegrationTestFiles/temp b/IntegrationTestFiles/temp new file mode 100644 index 000000000..d53eaac1d --- /dev/null +++ b/IntegrationTestFiles/temp @@ -0,0 +1,7 @@ +type cartesianField; +materials (fuel); +origin (0 0 0); +shape (1 1 0); +pitch (1.26 1.26 0); +default -1; +fuel ( 900 ); diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 8492876d5..a5a99e86c 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -123,6 +123,7 @@ module aceNeutronDatabase_class procedure :: updateTotalTempMajXS procedure :: updateRelEnMacroXSs procedure :: makeNuclideName + procedure :: computeMaxTrackXS end type aceNeutronDatabase @@ -416,14 +417,14 @@ subroutine updateTrackMatXS(self, E, matIdx, temp, rho, rand) if (present(temp)) then T = temp else - T = -ONE + T = NO_TEMPERATURE end if matCache % T_track = T if (present(rho)) then densityFactor = rho else - densityFactor = -ONE + densityFactor = NO_DENSITY end if matCache % rho_track = densityFactor @@ -452,6 +453,8 @@ end subroutine updateTrackMatXS !! Args: !! E [in] -> Incident neutron energy for which temperature majorant is found !! matIdx [in] -> Index of material for which the material temperature majorant is found + !! temp [in] -> Temperature, in Kelvin, locally + !! rho [in] -> Density scaling factor (dimensionless) locally !! subroutine updateTotalTempMajXS(self, E, matIdx, temp, rho) class(aceNeutronDatabase), intent(in) :: self @@ -469,7 +472,7 @@ subroutine updateTotalTempMajXS(self, E, matIdx, temp, rho) matCache % trackXS = ZERO ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = mat % kT else kT = temp * kBoltzmannMeV @@ -499,6 +502,52 @@ subroutine updateTotalTempMajXS(self, E, matIdx, temp, rho) end associate end subroutine updateTotalTempMajXS + + !! + !! This subroutine exists to simplify the construction of the majorant in the presence + !! of a temperature field. It updates the materialCache to have the largest cross section + !! of two if evaluated at two different temperatures: an imposed temperature and the input + !! material temperature. + !! + !! The majorant cross section at a given energy point may be given by either the higher + !! or lower of the possible temperatures. If a temperature field is imposed (where the + !! temperature will vary spatially across a given material) then the majorant should + !! be chosen appropriately/conservatively. + !! + subroutine computeMaxTrackXS(self, E, matIdx, maxTemp, rho) + class(aceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: maxTemp + real(defReal), intent(in) :: rho + real(defReal) :: xs1, xs2 + + associate (matCache => cache_materialCache(matIdx)) + + ! Clean current total XS + matCache % trackXS = ZERO + + ! If no temperature is given then the usual majorant from the material + ! temperature is fine and we needn't go further + if (maxTemp < ZERO) then + call self % updateTotalTempMajXS(E, matIdx, maxTemp, rho) + + ! Otherwise, compute the xs for both the material temperature and the given + ! max temperature + else + call self % updateTotalTempMajXS(E, matIdx, NO_TEMPERATURE, rho) + xs1 = matCache % trackXS + + call self % updateTotalTempMajXS(E, matIdx, maxTemp, rho) + xs2 = matCache % trackXS + + matCache % trackXS = max(xs1, xs2) + + end if + + end associate + + end subroutine computeMaxTrackXS !! !! Make sure that totalXS of material with matIdx is at energy E @@ -525,7 +574,7 @@ subroutine updateTotalMatXS(self, E, matIdx, temp, rho, rand) matCache % rho_tot = rho ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = mat % kT else kT = temp * kBoltzmannMeV @@ -563,7 +612,6 @@ subroutine updateTotalMatXS(self, E, matIdx, temp, rho, rand) matCache % xss % total = matCache % xss % total * densityFactor end if - matCache % T_tot = temp end associate @@ -594,7 +642,7 @@ subroutine updateMacroXSs(self, E, matIdx, temp, rho, rand) matCache % rho_tot = rho ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = mat % kT else kT = temp * kBoltzmannMeV @@ -832,7 +880,7 @@ subroutine updateTotalTempNucXS(self, E, kT, nucIdx) ! Catch negative dkT - perhaps due to input field values if (deltakT < ZERO) call fatalError(Here,& - 'Invalid input temperatureL '//numToChar(kT/kBoltzmannMeV)) + 'Invalid input temperature: '//numToChar(kT/kBoltzmannMeV)) ! Check if an update is required if (nucCache % E_maj /= E .or. nucCache % deltakT /= deltakT) then @@ -1353,12 +1401,12 @@ subroutine initMajorant(self, loud, maxTemp, scaleDensity) ! Check maxTemp is present and sensible if (present(maxTemp)) then if (maxTemp < ZERO) then - broadenTemp = -ONE + broadenTemp = NO_TEMPERATURE else broadenTemp = maxTemp end if else - broadenTemp = -ONE + broadenTemp = NO_TEMPERATURE end if ! Check scaleDensity is present and sensible @@ -1568,8 +1616,9 @@ subroutine initMajorant(self, loud, maxTemp, scaleDensity) matIdx = self % activeMat(j) ! Get material tracking cross section - if (self % materials(matIdx) % hasTMS .and. broadenTemp > ZERO) then - call self % updateTrackMatXS(E, matIdx, temp = broadenTemp, rho = densityFactor, rand = rand) + if (self % materials(matIdx) % useTMS(E) .and. broadenTemp > ZERO) then + call self % computeMaxTrackXS(E, matIdx, maxTemp = broadenTemp, rho = densityFactor) + !call self % updateTrackMatXS(E, matIdx, temp = broadenTemp, rho = densityFactor, rand = rand) else call self % updateTrackMatXS(E, matIdx, rho = densityFactor, rand = rand) end if @@ -1597,8 +1646,8 @@ subroutine initMajorant(self, loud, maxTemp, scaleDensity) nucXS = urrMaj end if - ! Update total material cross section - trackXS = trackXS + dens * (nucXS - cache_nuclideCache(nucIdx) % xss % total) + ! Update total material cross section + trackXS = trackXS + dens * (nucXS - cache_nuclideCache(nucIdx) % xss % total) end if diff --git a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 index 60c5402c8..b3b8e7154 100644 --- a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 +++ b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 @@ -19,6 +19,7 @@ module ceNeutronCache_mod use numPrecision + use universalVariables, only : NO_TEMPERATURE, NO_DENSITY use genericProcedures, only : fatalError, numToChar use neutronXsPackages_class, only : neutronMicroXSs, neutronMacroXSs @@ -37,16 +38,16 @@ module ceNeutronCache_mod !! idx -> Index on a nuclide grid for energy E_tot !! xss -> Cached cross-section values !! E_track -> Energy of the tracking xs + !! trackXS -> Cached tracking xs; this can be different to xss % total when using TMS !! T_track -> Temperature at which tracking XS was evaluated !! rho_track -> Density scaling at which tracking XS was evaluated - !! trackXS -> Cached tracking xs; this can be different to xss % total when using TMS !! E_rel -> Base energy for which relative energy cross sections are found (for TMS) !! xssRel -> Cached effective cross-section values at energy relative to E_rel (for TMS) !! type, public :: cacheMatDat real(defReal) :: E_tot = ZERO - real(defReal) :: T_tot = -ONE - real(defReal) :: rho_tot = -ONE + real(defReal) :: T_tot = NO_TEMPERATURE + real(defReal) :: rho_tot = NO_DENSITY real(defReal) :: E_tail = ZERO real(defReal) :: f = ZERO integer(shortInt) :: idx = 0 @@ -55,8 +56,8 @@ module ceNeutronCache_mod ! Tracking data real(defReal) :: E_track = ZERO real(defReal) :: trackXS = ZERO - real(defReal) :: T_track = -ONE - real(defReal) :: rho_track = -ONE + real(defReal) :: T_track = NO_TEMPERATURE + real(defReal) :: rho_track = NO_DENSITY ! TMS data real(defReal) :: E_rel = ZERO diff --git a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 index 96d2c41fd..20ef6f473 100644 --- a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 @@ -86,8 +86,8 @@ module ceNeutronDatabase_inter !! !! Return energy bounds for data in the database !! - !! eMin and eMax are minimun and maximumum energy such that data - !! for ALL nuclides if avalible + !! eMin and eMax are minimum and maximumum energy such that data + !! for ALL nuclides is available !! !! Args: !! eMin [out] -> minimum value of energy [MeV] @@ -477,7 +477,7 @@ function getMajorantXS(self, p) result(xs) ! Check dynamic type of the particle if (p % isMG .or. p % type /= P_NEUTRON) then - call fatalError(Here, 'Dynamic type of the partcle is not CE Neutron but:'//p % typeToChar()) + call fatalError(Here, 'Dynamic type of the particle is not CE Neutron but:'//p % typeToChar()) end if ! Check Cache and update if needed @@ -515,5 +515,4 @@ pure function ceNeutronDatabase_CptrCast(source) result(ptr) end function ceNeutronDatabase_CptrCast - end module ceNeutronDatabase_inter diff --git a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 index 05395c74f..4175c10c3 100644 --- a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 +++ b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 @@ -360,7 +360,7 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) trackMatXS = materialCache(self % matIdx) % trackXS * rand % get() ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV @@ -373,8 +373,8 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) densityFactor = rho end if - ! Loop over nuclides - do i = 1,size(self % nuclides) + ! Loop over nuclides until one is sampled + nucLoop: do i = 1,size(self % nuclides) nucIdx = self % nuclides(i) dens = self % dens(i) * densityFactor @@ -393,9 +393,9 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) else ! Update nuclide cache if needed - if (E /= nucCache % E_tot .or. temp /= materialCache(self % matIdx) % T_tot) then + !!if (E /= nucCache % E_tot .or. temp /= materialCache(self % matIdx) % T_tot) then call self % data % updateTotalNucXS(E, nucIdx, kT, rand) - end if + !!end if totNucXS = nucCache % xss % total end if @@ -440,17 +440,18 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) end if - ! Exit function, return the sampled nucIdx + materialCache(self % matIdx) % T_tot = temp return end if end associate - end do + end do nucLoop ! Print error message as the inversion failed call fatalError(Here,'Nuclide sampling loop failed to terminate') + end subroutine sampleNuclide @@ -503,7 +504,7 @@ function sampleFission(self, E, temp, rho, rand) result(nucIdx) end if ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV @@ -599,7 +600,7 @@ function sampleScatter(self, E, temp, rho, rand) result(nucIdx) end if ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV @@ -691,7 +692,7 @@ function sampleScatterWithFission(self, E, temp, rho, rand) result(nucIdx) end if ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV diff --git a/PhysicsPackages/alphaPhysicsPackage_class.f90 b/PhysicsPackages/alphaPhysicsPackage_class.f90 index 449e4d561..f2aa219cc 100644 --- a/PhysicsPackages/alphaPhysicsPackage_class.f90 +++ b/PhysicsPackages/alphaPhysicsPackage_class.f90 @@ -30,13 +30,15 @@ module alphaPhysicsPackage_class ! Geometry use geometry_inter, only : geometry use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx, & - gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr, & + gr_fieldPtrName => fieldPtrName use geometryFactory_func, only : new_geometry ! Fields use field_inter, only : field use uniFissSitesField_class, only : uniFissSitesField, uniFissSitesField_TptrCast use fieldFactory_func, only : new_field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast ! Nuclear Data use materialMenu_mod, only : mm_nMat => nMat @@ -62,6 +64,9 @@ module alphaPhysicsPackage_class use tallyResult_class, only : tallyResult use kAlphaAnalogClerk_class, only : kAlphaResult + ! Coupling + use couplingAdmin_class, only : couplingAdmin + ! Factories use transportOperatorFactory_func, only : new_transportOperator @@ -105,6 +110,9 @@ module alphaPhysicsPackage_class integer(shortInt) :: bufferSize logical(defBool) :: UFS = .false. + ! Coupling handler + type(couplingAdmin) :: couplingInfo + ! Calculation components type(particleDungeon), pointer :: thisCycle => null() type(particleDungeon), pointer :: nextCycle => null() @@ -142,6 +150,7 @@ subroutine run(self) call self % generateInitialState() call self % cycles(self % inactiveTally, self % inactiveAtch, self % N_inactive) call self % buildActiveAtch() + call self % couplingInfo % endCoupling() call self % cycles(self % activeTally, self % activeAtch, self % N_active) ! Collect results from other processes @@ -328,6 +337,9 @@ subroutine cycles(self, tally, tallyAtch, N_cycles) call statusMsg("Time to end: " // trim(secToChar(T_toEnd))) call tally % display() + ! Perform coupling operations + call self % couplingInfo % couple(i) + end do ! Load elapsed time @@ -422,6 +434,8 @@ subroutine init(self, dict) type(outputFile) :: test_out type(visualiser) :: viz class(field), pointer :: field + class(pieceConstantField), pointer :: pcField + real(defReal) :: maxTemperature, maxDensityScale character(100), parameter :: Here ='init (alphaPhysicsPackage_class.f90)' call cpu_time(self % CPU_time_start) @@ -503,6 +517,31 @@ subroutine init(self, dict) ! Activate Nuclear Data *** All materials are active call ndReg_activate(self % particleType, nucData, self % geom % activeMats()) self % nucData => ndReg_get(self % particleType) + + ! If present, build temperature field + if (dict % isPresent('temperature')) then + tempDict => dict % getDictPtr('temperature') + call new_field(tempDict, nameTemperature) + field => gr_fieldPtrName(nameTemperature) + pcField => pieceConstantField_CptrCast(field) + maxTemperature = pcField % getMaxValue() + else + maxTemperature = NO_TEMPERATURE + end if + + ! If present, build density field + if (dict % isPresent('density')) then + tempDict => dict % getDictPtr('density') + call new_field(tempDict, nameDensity) + field => gr_fieldPtrName(nameDensity) + pcField => pieceConstantField_CptrCast(field) + maxDensityScale = pcField % getMaxValue() + else + maxDensityScale = NO_DENSITY + end if + + ! Update majorant in case of density and temperature fields + call self % nucData % initMajorant(.false., maxTemp = maxTemperature, scaleDensity = maxDensityScale) ! Call visualisation if (dict % isPresent('viz')) then @@ -590,6 +629,11 @@ subroutine init(self, dict) ! Don't create the active tally yet as the guess value ! strongly affects the convergence rate and may bias ! the active cycle estimate + + ! Attach a tally admin for coupling + if (self % couplingInfo % doCoupling()) then + call self % couplingInfo % attachTally(self % inactiveTally) + end if call self % printSettings() diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index f5aab91fc..993424c17 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -442,7 +442,7 @@ subroutine init(self, dict) type(visualiser) :: viz class(field), pointer :: field class(pieceConstantField), pointer :: pcField - real(defReal) :: maxDensityScale, maxTemperature + real(defReal) :: maxDensityScale, maxTemperature, maxCouplingTemp character(100), parameter :: Here ='init (eigenPhysicsPackage_class.f90)' call cpu_time(self % CPU_time_start) @@ -463,9 +463,15 @@ subroutine init(self, dict) call dict % getOrDefault(self % bufferSize, 'buffer', 1000) ! Check if the calculation is coupled + maxCouplingTemp = NO_TEMPERATURE if (dict % isPresent('coupling')) then tempDict => dict % getDictPtr('coupling') call self % couplingInfo % init(tempDict) + call tempDict % getOrDefault(maxCouplingTemp, 'maxTemp', NO_TEMPERATURE) + if (.not. tempDict % isPresent('maxTemp')) then + call statusMsg('If coupling with temperature fields, it is advisable to include a maximum '& + //'possible temperature (maxTemp) in order to correctly precompute the majorant.') + end if end if ! Process type of data @@ -556,6 +562,7 @@ subroutine init(self, dict) else maxTemperature = NO_TEMPERATURE end if + maxTemperature = max(maxTemperature, maxCouplingTemp) ! If present, build density field if (dict % isPresent('density')) then diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 5a60419ab..c9e93b82c 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -18,11 +18,9 @@ module universalVariables integer(shortInt), parameter, public :: MAX_COL = 70 ! Maximum number of columns in console display ! Define variables which are important for tracking neutrons in the geometry - real(defReal), parameter, public :: INFINITY = 2.0_defReal**63, & - surface_tol = 1.0e-12_defReal, & ! Tol. on closeness to surface - SURF_TOL = 1.0E-12_defReal, & - INF = 2.0_defReal**63, & - NUDGE = 1.0e-8_defReal ! Distance to poke neutrons across boundaries for surface tracking + real(defReal), parameter, public :: INF = 2.0_defReal**63, & + SURF_TOL = 1.0E-12_defReal, & + NUDGE = 1.0e-8_defReal ! Distance to poke neutrons across boundaries for surface tracking ! Flags for different possible events in movement in geometry integer(shortInt), parameter, public :: COLL_EV = 1, & diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index c1461b4a0..1bc9c0015 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -195,8 +195,8 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! This branch is called in the case of voids with no imposed XS if (sigmaTrack < tol) then - dist = INFINITY - invSigmaTrack = INFINITY + dist = INF + invSigmaTrack = INF sigmaT = ZERO else diff --git a/TransportOperator/transportOperatorST_class.f90 b/TransportOperator/transportOperatorST_class.f90 index f73d27563..80a35b7ef 100644 --- a/TransportOperator/transportOperatorST_class.f90 +++ b/TransportOperator/transportOperatorST_class.f90 @@ -70,8 +70,8 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! This branch is called in the case of voids with no imposed XS if (sigmaTrack < tol) then - dist = INFINITY - invSigmaTrack = INFINITY + dist = INF + invSigmaTrack = INF sigmaT = ZERO else diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index 4cc8998c1..d7e274fc0 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -62,6 +62,9 @@ Example: :: uniformFissionSites { } varianceReduction { } source { } + coupling { } + temperature { } + density { } .. note:: Although a ``source`` definition is not required, it can be included to replace @@ -121,6 +124,8 @@ Example: :: *Optional entries* :: varianceReduction { } + temperature { } + density { } kineticPhysicsPackage ##################### @@ -190,6 +195,8 @@ Example: :: *Optional entries* :: varianceReduction { } + temperature { } + density { } alphaPhysicsPackage ################### @@ -249,6 +256,8 @@ Example: :: uniformFissionSites { } varianceReduction { } source { } + temperature { } + density { } .. note:: Although a ``source`` definition is not required, it can be included to replace @@ -971,6 +980,54 @@ Example: :: density { type cartesianField; file ./myDensityField; } +Coupling +-------- + +To performed coupled physics simulations (currently limited to k and alpha eigenvalue calculations), the keyword +``coupling`` must be present in the input file with an appropriately defined dictionary. The entries +in this dictionary are: +* tally: a tallyAdmin definition, used to produce tallies and normalisation required for coupling, e.g., power. +* receiveFile: tex file that SCONE can expect to receive from external codes. This will be read to signal either that + coupling is continuing (containing 'SIGUSR1') or that coupling is finished (containing 'SIGTERM'). + The latter could be based on some convergence criteria being met by the external solver. +* sendFile: text file that SCONE should write to send signals to other codes. This will signal whether coupling should + continue (containing 'SIGUSR1') or that coupling has finished (containing 'SIGTERM'). +* updateFreq: sets how often SCONE exchanges data with coupled codes. More frequently may mean faster convergence + of the coupled system, but results sent to other codes will be noisier due to poorer statistics. +* outputFile: name of the file containing tallies that SCONE will produce. +* outputFormat (*optional*, default = ``asciiMATLAB``): type of output file. + Choices are ``asciiMATLAB`` and ``asciiJSON``. +* fieldPaths: array of paths to files that SCONE will read which contain field definitions, produced by external solvers. + There must be as many paths as fields, and their ordering must match. +* fields: array of names of fields that SCONE will update during coupling. Can be ``temperature`` and/or ``density``. +* maxTemp: (*optional*, defaults to the highest material temperature initially present in the geometry): + the maximum temperature expected to be produced during coupling. This is important to correctly initialise the + majorant if using delta tracking. If this temperature is exceeded during coupling, a warning will be produced + regarding the validity of the results. +* maxIt: (*optional*, defaults to infinity) the maximum number of MC cycles that will take place before SCONE + terminates the coupling. SCONE will continue to run, but it will signal to the external solver that coupling + has finished and fields will no longer be updated. + +An example of a coupling dictionary is: :: + + coupling { + tally {norm power; normVal 100; + power {type collisionClerk; response (kappa) kappa {type macroResponse; MT -80;}} + powerLocal {type collisionClerk; response (kappa) kappa {type macroResponse; MT -80;} + map {type spaceMap; axis z; grid lin; min -200; min 200; N 40;} + } + } + receiveFile /path/toSCONE; + sendFile /path/fromSCONE; + updateFreq 5; + outputFile /path/tallies_from_SCONE; + outputFormat asciiMATLAB; + fieldPaths (/path1/tempField /path2/densityField); + fields (temperature density); + maxTemp 2000; + maxIt 100; + } + Visualiser ---------- @@ -1689,14 +1746,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 From 7e4fec5d09cda2f45d8f345be7e68407f33ed8a6 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Sat, 14 Mar 2026 23:52:00 +0000 Subject: [PATCH 06/15] Fix to mg majorant --- .../mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 index d874eca34..682e2ae94 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 @@ -486,7 +486,7 @@ subroutine initMajorant(self, loud, maxTemp, scaleDensity) ! TODO: Update should there be a temperature model developed for MG XSs ! Allocate majorant - allocate (self % majorant(self % nG)) + if (.not. allocated(self % majorant)) allocate (self % majorant(self % nG)) ! Loop over energy groups do g = 1,self % nG From c1dc2a479dbd6ec4a1a459d6469d9bc6216f9494 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Mon, 16 Mar 2026 16:44:53 +0000 Subject: [PATCH 07/15] Fix to Brunner lattice --- InputFiles/Benchmarks/Transients/LatticeBrunner | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/InputFiles/Benchmarks/Transients/LatticeBrunner b/InputFiles/Benchmarks/Transients/LatticeBrunner index 0e82a9d08..c357d10db 100644 --- a/InputFiles/Benchmarks/Transients/LatticeBrunner +++ b/InputFiles/Benchmarks/Transients/LatticeBrunner @@ -8,12 +8,12 @@ type kineticPhysicsPackage; -pop 5000000; +pop 500000; timeSteps 1; XSdata mg; dataType mg; dt 7.0; -cycles 750; +cycles 250; buffer 100; collisionOperator { neutronMG {type neutronMGstd;}} From ec2fc4ed196c3cab2bb777b2bac1596cf086bd84 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Sat, 21 Mar 2026 22:25:46 +0000 Subject: [PATCH 08/15] Bug fix --- NuclearData/Reactions/reactionMG/fissionMG_class.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 index b37eb625c..c2aeb91ea 100644 --- a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 +++ b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 @@ -18,8 +18,6 @@ module fissionMG_class !! public :: fissionMG_TptrCast - real(defReal), private, parameter :: KAPPA_DEFAULT = 202.27 ! [MeV] - !! !! A special type of MG reaction that contains data related to fission !! From c4ffa53b58098662cce49cd8a51ef72c6cff82bf Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Sat, 21 Mar 2026 22:50:23 +0000 Subject: [PATCH 09/15] Fix for build tests --- Coupling/couplingAdmin_class.f90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Coupling/couplingAdmin_class.f90 b/Coupling/couplingAdmin_class.f90 index acede20d2..b7da176f4 100644 --- a/Coupling/couplingAdmin_class.f90 +++ b/Coupling/couplingAdmin_class.f90 @@ -116,7 +116,7 @@ subroutine init(self, dict) class(couplingAdmin), intent(inout) :: self class(dictionary), intent(in) :: dict integer(shortInt) :: i, nFields, j - logical(defBool) :: found, duplicate + logical(defBool) :: found class(dictionary), pointer :: tallyDict type(outputFile) :: test_out character(100), parameter :: Here ='init (couplingAdmin_class.f90)' @@ -178,10 +178,11 @@ subroutine init(self, dict) call fatalError(Here,'Field name is invalid. See list above.') end if - if (i < nFields) then - duplicate = any(trim(self % fieldNames(i)) == [ (trim(self % fieldNames(j)), j=i+1,nFields) ]) - if (duplicate) call fatalError(Here,'Field names must be unique: '//trim(self % fieldNames(i))) - end if + do j = i+1, nFields + if (trim(self % fieldNames(i)) == trim(self % fieldNames(j))) then + call fatalError(Here,'Field names must be unique: '//trim(self % fieldNames(i))) + end if + end do end do From 5e66d1c2e66ba83039ea75873440d9924931cbb1 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Mon, 23 Mar 2026 18:21:52 +0000 Subject: [PATCH 10/15] 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 6597103130909b28c2508676a17646ede9f20719 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Mon, 20 Apr 2026 11:41:32 +0100 Subject: [PATCH 11/15] Fixes to data handling --- .../PieceConstantFields/pieceConstantField_inter.f90 | 2 +- Geometry/geometry_inter.f90 | 6 +++--- NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 | 5 ++--- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 b/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 index 8d21708e5..3889b05fd 100644 --- a/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 +++ b/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 @@ -117,7 +117,7 @@ subroutine init(self, dict) class(pieceConstantField), intent(inout) :: self class(dictionary), intent(in) :: dict character(pathLen) :: path - class(dictionary), pointer :: tempDict + type(dictionary) :: tempDict if (dict % isPresent('path')) then call dict % get(path, 'path') diff --git a/Geometry/geometry_inter.f90 b/Geometry/geometry_inter.f90 index c473f5194..8e9161e14 100644 --- a/Geometry/geometry_inter.f90 +++ b/Geometry/geometry_inter.f90 @@ -1,7 +1,7 @@ module geometry_inter use numPrecision - use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST, INFINITY, COLL_EV + use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST, INF, COLL_EV use genericProcedures, only : fatalError use dictionary_class, only : dictionary use charMap_class, only : charMap @@ -552,14 +552,14 @@ subroutine phongTrace(self, ray, matIdx, bright, mats, ambient, light) ! Trace until an opaque material is hit do while (any(mats == matIdx)) - dist = INFINITY + dist = INF call self % moveNoBC(ray, dist, event, normal0) matIdx = ray % matIdx ! If distance is infinite or particle collided, not pointing towards anything opaque. ! Hence, short circuit the trace - if ((dist == INFINITY) .or. (event == COLL_EV)) then + if ((dist == INF) .or. (event == COLL_EV)) then lightTrace = .false. exit end if diff --git a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 index 4175c10c3..b329eb30e 100644 --- a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 +++ b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 @@ -393,9 +393,9 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) else ! Update nuclide cache if needed - !!if (E /= nucCache % E_tot .or. temp /= materialCache(self % matIdx) % T_tot) then + if (E /= nucCache % E_tot .or. temp /= materialCache(self % matIdx) % T_tot) then call self % data % updateTotalNucXS(E, nucIdx, kT, rand) - !!end if + end if totNucXS = nucCache % xss % total end if @@ -440,7 +440,6 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) end if - materialCache(self % matIdx) % T_tot = temp return end if From a0723e9beb84a7b2183121d544ece8c34260ba76 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Mon, 20 Apr 2026 12:41:53 +0100 Subject: [PATCH 12/15] Fix to mg integration test --- .../baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 index 6e1625d89..3b00ebb56 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 @@ -86,6 +86,7 @@ subroutine testBaseMgNeutronDatabaseWithP0() p % pRNG => pRNG p % type = P_NEUTRON p % G = 1 + p % isMG = .true. @assertEqual(2.1_defReal, database % getTrackingXS(p, 1, MATERIAL_XS), TOL) ! Test getting Total XS @@ -217,6 +218,7 @@ subroutine testBaseMgNeutronDatabaseWithP1() p % pRNG => pRNG p % type = P_NEUTRON p % G = 1 + p % isMG = .true. @assertEqual(2.1_defReal, database % getTrackingXS(p, 1, MATERIAL_XS), TOL) ! Test getting Total XS From 5af1841358bb9e309bb7049f1f6d6db61ee64d13 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Mon, 1 Jun 2026 22:07:15 +0100 Subject: [PATCH 13/15] Finalising couplingAdmin --- .../neutronCEstd_class.f90 | 32 ++-- CollisionOperator/scatteringKernels_func.f90 | 13 +- Coupling/couplingAdmin_class.f90 | 32 +++- DataStructures/charMap_class.f90 | 2 +- DataStructures/heapQueue_class.f90 | 6 +- .../cartesianField_class.f90 | 8 +- .../thermalScatterInelastic_class.f90 | 6 +- .../Tests/thermalScatteringData_iTest.f90 | 8 +- .../aceDatabase/aceNeutronDatabase_class.f90 | 92 +++++----- .../aceDatabase/aceNeutronNuclide_class.f90 | 170 ++++++++++++------ .../angleLawENDF/tabularAngle_class.f90 | 4 + NuclearData/materialMenu_mod.f90 | 22 +-- .../Tests/baseMgNeutronDatabase_iTest.f90 | 2 + NuclearData/nuclearDataReg_mod.f90 | 2 +- PhysicsPackages/alphaPhysicsPackage_class.f90 | 2 +- PhysicsPackages/eigenPhysicsPackage_class.f90 | 6 +- .../transportOperatorDT_class.f90 | 6 +- .../transportOperatorHT_class.f90 | 8 +- .../transportOperatorST_class.f90 | 2 + 19 files changed, 266 insertions(+), 157 deletions(-) diff --git a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 index 5748a388d..fbe9199e8 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 @@ -273,7 +273,7 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) r = p % rGlobal() do i = 1, n - call fission % sampleOut(mu, phi, E_out, p % E, p % pRNG, lambda) + call fission % sampleOut(mu, phi, E_out, collDat % E, p % pRNG, lambda) ! Skip if a delayed particle is produced in prompt-only mode if (self % neglectDelayed .and. lambda < huge(lambda)) cycle @@ -363,7 +363,7 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) character(100),parameter :: Here = 'elastic (neutronCEstd_class.f90)' ! Assess if thermal scattering data is needed or not - if (self % nuc % needsSabEl(p % E)) collDat % MT = N_N_ThermEL + if (self % nuc % needsSabEl(colLDat % E)) collDat % MT = N_N_ThermEL ! Get reaction reac => uncorrelatedReactionCE_CptrCast( self % xsData % getReaction(collDat % MT, collDat % nucIdx)) @@ -373,7 +373,7 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) collDat % A = self % nuc % getMass() ! Retrieve kT from either material or nuclide - if (self % mat % useTMS(p % E)) then + if (self % mat % useTMS(collDat % E)) then collDat % kT = self % mat % kT else collDat % kT = self % nuc % getkT() @@ -383,7 +383,7 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) ! Check is DBRC is on hasDBRC = self % nuc % hasDBRC() - isFixed = (.not. hasDBRC) .and. (p % E > collDat % kT * self % threshE) & + isFixed = (.not. hasDBRC) .and. (collDat % E > collDat % kT * self % threshE) & & .and. (collDat % A > self % threshA) ! Apply criterion for Free-Gas vs Fixed Target scattering @@ -414,7 +414,7 @@ subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) collDat % MT = self % nuc % invertInelastic(collDat % E, p % pRNG) reac => uncorrelatedReactionCE_CptrCast(self % xsData % getReaction(collDat % MT, collDat % nucIdx)) if (.not.associated(reac)) call fatalError(Here, "Failed to get scattering reaction") - + ! Scatter particle if (reac % inCMFrame()) then collDat % A = self % nuc % getMass() @@ -423,8 +423,8 @@ subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) call self % scatterInLAB(p, collDat, reac) end if - ! Apply weigth change - p % w = p % w * reac % release(p % E) + ! Apply weight change using ingoing collision particle energy + p % w = p % w * reac % release(collDat % E) end subroutine inelastic @@ -456,7 +456,7 @@ subroutine scatterInLAB(self, p, collDat, reac) real(defReal) :: E_out, mu ! Sample scattering angles and post-collision energy - call reac % sampleOut(mu, phi, E_out, p % E, p % pRNG) + call reac % sampleOut(mu, phi, E_out, collDat % E, p % pRNG) ! Update neutron state p % E = E_out @@ -483,10 +483,10 @@ subroutine scatterFromFixed(self, p, collDat, reac) MT = collDat % MT ! Sample mu, phi and outgoing energy - call reac % sampleOut(mu, phi, E_outCM, p % E, p % pRNG) + call reac % sampleOut(mu, phi, E_outCM, collDat % E, p % pRNG) ! Save incident energy - E_out = p % E + E_out = collDat % E if (MT == N_N_elastic) then call asymptoticScatter(E_out, mu, collDat % A) @@ -530,11 +530,11 @@ subroutine scatterFromMoving(self, p, collDat, reac) ! Get neutron direction and velocity dir_pre = p % dirGlobal() - V_n = dir_pre * sqrt(p % E) + V_n = dir_pre * sqrt(collDat % E) ! Sample target velocity with constant XS or with DBRC ! Check energy range - inEnergyRange = ((p % E <= self % DBRCeMax) .and. (self % DBRCeMin <= p % E)) + inEnergyRange = ((collDat % E <= self % DBRCeMax) .and. (self % DBRCeMin <= collDat % E)) ! Check if DBRC is on for this target nuclide hasDBRC = self % nuc % hasDBRC() @@ -548,14 +548,14 @@ subroutine scatterFromMoving(self, p, collDat, reac) if(.not.associated(ceNuc0K)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide') ! Get elastic scattering 0K majorant - maj = self % xsData % getScattMicroMajXS(p % E, kT, A, nucIdx) + maj = self % xsData % getScattMicroMajXS(collDat % E, kT, A, nucIdx) ! Use DBRC to sample target velocity - V_t = targetVelocity_DBRCXS(ceNuc0K, p % E, dir_pre, A, kT, p % pRNG, maj) + V_t = targetVelocity_DBRCXS(ceNuc0K, collDat % E, dir_pre, A, kT, p % pRNG, maj) else ! Constant cross section approximation - V_t = targetVelocity_constXS(p % E, dir_pre, A, kT, p % pRNG) + V_t = targetVelocity_constXS(collDat % E, dir_pre, A, kT, p % pRNG) end if @@ -568,7 +568,7 @@ subroutine scatterFromMoving(self, p, collDat, reac) V_n = V_n / U_n ! Sample mu and phi in CM frame - call reac % sampleOut(mu, phi, dummy, p % E, p % pRNG) + call reac % sampleOut(mu, phi, dummy, collDat % E, p % pRNG) ! Obtain post collision speed V_n = rotateVector(V_n, mu, phi) * U_n diff --git a/CollisionOperator/scatteringKernels_func.f90 b/CollisionOperator/scatteringKernels_func.f90 index c80df4afb..95952c93e 100644 --- a/CollisionOperator/scatteringKernels_func.f90 +++ b/CollisionOperator/scatteringKernels_func.f90 @@ -1,9 +1,10 @@ module scatteringKernels_func use numPrecision - use genericProcedures, only : rotateVector - use RNG_class, only : RNG - use particle_class, only : particle + use universalVariables, only : MINIMUM_ENERGY + use genericProcedures, only : rotateVector + use RNG_class, only : RNG + use particle_class, only : particle ! Nuclear Data use ceNeutronNuclide_inter, only : ceNeutronNuclide @@ -53,7 +54,8 @@ subroutine asymptoticScatter(E,mu,A) mu = (A*mu + 1)*sqrt(E_in/E)* inv_Ap1 ! Correct possible nuclear data shortcomings - if (mu > ONE) mu = ONE + mu = min(mu, ONE) + E = max(E, MINIMUM_ENERGY) end subroutine asymptoticScatter @@ -87,7 +89,8 @@ subroutine asymptoticInelasticScatter(E,mu,E_out,A) mu = mu * sqrt(E_out/E) + sqrt(E_in/E)* inv_Ap1 ! Correct possible nuclear data shortcomings - if (mu > ONE) mu = ONE + mu = min(mu, ONE) + E = max(E, MINIMUM_ENERGY) end subroutine asymptoticInelasticScatter diff --git a/Coupling/couplingAdmin_class.f90 b/Coupling/couplingAdmin_class.f90 index b7da176f4..cf437fbfb 100644 --- a/Coupling/couplingAdmin_class.f90 +++ b/Coupling/couplingAdmin_class.f90 @@ -38,9 +38,11 @@ module couplingAdmin_class !! - output the results of those tallies to an appropriate location !! - know the location of fields provided by other physics and update them in SCONE !! - possibly provide checks on validity of information fed back to SCONE + !! - decide whether to continue coupling during active cycles !! !! Private members: !! isCoupled -> flag for whether coupling is being performed + !! endWhenActive -> flag to end coupling when moving to active cycles !! commFileIn -> path to the file which triggers a SCONE calculation by the driver !! commFileOut -> path to the file which SCONE outputs to the driver on completing a calculation !! outputFile -> path to the file which SCONE outputs containing results to be used by the driver @@ -56,6 +58,7 @@ module couplingAdmin_class !! kill -> clean up !! couple -> performs the main coupling logic !! endCoupling -> deactivates coupling + !! moveToActive -> informs admin of shift to active cycles, ending coupling if needs be !! !! doCoupling -> returns a logical if coupling is happening !! attachTally -> attach coupling tallies to another tallyAdmin used by a physics package @@ -74,6 +77,7 @@ module couplingAdmin_class type, public :: couplingAdmin private logical(defBool) :: isCoupled = .false. + logical(defBool) :: endWhenActive = .true. character(pathLen) :: commFileIn = '' character(pathLen) :: commFileOut = '' character(pathLen) :: outputFile = '' @@ -89,6 +93,7 @@ module couplingAdmin_class procedure :: kill procedure :: couple procedure :: endCoupling + procedure :: moveToActive ! Status checks procedure :: doCoupling @@ -114,7 +119,7 @@ module couplingAdmin_class !! subroutine init(self, dict) class(couplingAdmin), intent(inout) :: self - class(dictionary), intent(in) :: dict + type(dictionary), intent(in) :: dict integer(shortInt) :: i, nFields, j logical(defBool) :: found class(dictionary), pointer :: tallyDict @@ -127,6 +132,9 @@ subroutine init(self, dict) call dict % get(self % commFileIn, 'receiveFile') call dict % get(self % commFileOut, 'sendFile') + ! Read whether to end when active + call dict % getOrDefault(self % endWhenActive, 'endActive', .true.) + ! Read how often to exchange data call dict % get(self % updateFreq, 'updateFreq') if (self % updateFreq < 1) call fatalError(Here,& @@ -195,6 +203,7 @@ subroutine kill(self) class(couplingAdmin), intent(inout) :: self self % isCoupled = .false. + self % endWhenActive = .true. self % updateFreq = 0 self % commFileIn = '' self % commFileOut = '' @@ -279,6 +288,23 @@ subroutine endCoupling(self) end subroutine endCoupling + !! + !! Signals to admin that active cycles have begun + !! Decides whether coupling should end. + !! If it continues, attach the active tally + !! + subroutine moveToActive(self, tallyPtr) + class(couplingAdmin), intent(inout) :: self + type(tallyAdmin), pointer, intent(inout) :: tallyPtr + + if (self % endWhenActive) then + call self % endCoupling() + else + call self % attachTally(tallyPtr) + end if + + end subroutine moveToActive + !! !! Returns whether to perform a physics update based on iteration number !! and update frequency @@ -412,8 +438,8 @@ end subroutine waitForSignal !! Attach coupling tallies to another tally !! subroutine attachTally(self, tallyPtr) - class(couplingAdmin), intent(in) :: self - type(tallyAdmin), pointer, intent(inout) :: tallyPtr + class(couplingAdmin), intent(in) :: self + type(tallyAdmin), pointer, intent(inout) :: tallyPtr if (self % doCoupling()) then call tallyPtr % push(self % tally) diff --git a/DataStructures/charMap_class.f90 b/DataStructures/charMap_class.f90 index b1e14f798..a9b5103fe 100644 --- a/DataStructures/charMap_class.f90 +++ b/DataStructures/charMap_class.f90 @@ -34,7 +34,7 @@ module charMap_class !! Implementation is based on intMap !! !! NOTE: Following structure can be used to loop over entire map - !! But note that you CANNOT modifay entries with "add" inside this loop. + !! But note that you CANNOT modify entries with "add" inside this loop. !! Use "atSet" instead !! it = map % begin() !! do while (it /= map % end()) diff --git a/DataStructures/heapQueue_class.f90 b/DataStructures/heapQueue_class.f90 index 6daf01854..cbdcdf818 100644 --- a/DataStructures/heapQueue_class.f90 +++ b/DataStructures/heapQueue_class.f90 @@ -138,8 +138,10 @@ subroutine replace(self, val) ! We need to consider the case where there is only one child in a node ! (can happen when the size is even). If child is the last element it is the ! node with no sibling - if (child /= self % size .and. self % heap(child) < self % heap(child + 1)) then - child = child + 1 + if (child /= self % size) then + if (self % heap(child) < self % heap(child + 1)) then + child = child + 1 + end if end if ! If the parent is larger than the larger child we are done diff --git a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 index 74a54fe58..fd4335f1f 100644 --- a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 +++ b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 @@ -163,15 +163,21 @@ subroutine init_dict(self, dict) end if ! Size field value array - self % N = product(self % sizeN * self % nMat) + self % N = product(self % sizeN) * self % nMat allocate(self % val(self % N + 1)) ! Read field values for each material idx0 = 0 do i = 1, size(mats) + if (allocated(temp)) deallocate(temp) call dict % get(temp, mats(i)) + ! Ensure size is as expected + if (size(temp) /= product(self % sizeN)) call fatalError(Here,& + 'Invalid input size to field. Should be '//numToChar(self % sizeN)//'. '//& + 'Was: '//numToChar(size(temp))) + ! Flip array up-down for more natural input ! Reshape into rank 2 array tempMap = reshape(temp, [self % sizeN(1), self % sizeN(2) * self % sizeN(3)]) diff --git a/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 b/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 index e569690bc..7829910a4 100644 --- a/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 +++ b/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 @@ -2,6 +2,7 @@ module thermalScatterInelastic_class use numPrecision use endfConstants + use universalVariables, only : MAXIMUM_ENERGY, MINIMUM_ENERGY use genericProcedures, only : binarySearch, fatalError, numToChar, endfInterpolate, searchError use RNG_class, only : RNG use dataDeck_inter, only : dataDeck @@ -295,10 +296,13 @@ subroutine sampleOut(self, mu, phi, E_out, E_in, rand, lambda) end if - if (E_out > 20.0 .or. E_out <= ZERO) then + if (E_out > MAXIMUM_ENERGY .or. E_out <= ZERO) then call fatalError(Here,'Failed to find energy: '//numToChar(E_out)) end if + ! Allow energies lower than 1E-11 to be sampled, but floor to the minimum + E_out = max(E_out, MINIMUM_ENERGY) + ! Sample phi phi = rand % get() * TWO_PI diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 index a16ec7a1c..5143dc950 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 @@ -80,9 +80,9 @@ subroutine test_thermalScatteringData() ! Get nuclides O16 => aceNeutronNuclide_CptrCast( data % getNuclide(1)) - H1 => aceNeutronNuclide_CptrCast( data % getNuclide(2)) + C12 => aceNeutronNuclide_CptrCast( data % getNuclide(2)) H1_2 => aceNeutronNuclide_CptrCast( data % getNuclide(3)) - C12 => aceNeutronNuclide_CptrCast( data % getNuclide(4)) + H1 => aceNeutronNuclide_CptrCast( data % getNuclide(4)) !<><><><><><><><><><><><><><><><><><><><><><><><> ! Test scattering tables @@ -158,8 +158,8 @@ subroutine test_thermalScatteringData() !<><><><><><><><><><><><><><><><><><><><><><><><> ! Test getting XSs ! H-1 - nuc => ceNeutronNuclide_CptrCast(data % getNuclide(2)) - nuclideCache(2) % E_tot = ONE + nuc => ceNeutronNuclide_CptrCast(data % getNuclide(4)) + nuclideCache(4) % E_tot = ONE call nuc % getMicroXSs(microXSs, 1.8E-6_defReal, ZERO, p % pRNG) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index a5a99e86c..377a6ebf2 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -10,6 +10,7 @@ module aceNeutronDatabase_class use RNG_class, only : RNG use charMap_class, only : charMap use intMap_class, only : intMap + use hashFunctions_func, only : FNV_1 ! Nuclear Data Interfaces use nuclearDatabase_inter, only : nuclearDatabase @@ -927,19 +928,21 @@ subroutine init(self, dict, ptr, silent ) type(materialItem), pointer :: mat class(ceNeutronDatabase), pointer :: ptr_ceDatabase type(charMap) :: nucSet + type(dictionary) :: sabDict type(aceCard) :: ACE - type(aceSabCard) :: ACE_Sab1, ACE_Sab2 + type(aceSabCard), dimension(:), allocatable :: ACE_Sab + character(nameLen), dimension(:), allocatable :: fileNames character(pathLen) :: aceLibPath - character(nameLen) :: name, name_file1, name_file2, nucDBRC_temp - integer(shortInt) :: i, j, envFlag, nucIdx, idx, idx1, idx2 + character(nameLen) :: nucDBRC_temp, name + integer(shortInt) :: i, j, envFlag, nucIdx, idx integer(shortInt) :: maxNuc logical(defBool) :: isFissileMat integer(shortInt),dimension(:),allocatable :: nucIdxs, zaidDBRC character(nameLen),dimension(:),allocatable :: nucDBRC real(defReal) :: A, nuckT, eUpSab, eUpSabNuc, & eLowURR, eLowUrrNuc, alpha, & - deltakT, eUpper, eLower, kT, & - temp + deltakT, eUpper, eLower, & + temp, kT real(defReal), dimension(2) :: sabT integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 character(100), parameter :: Here = 'init (aceNeutronDatabase_class.f90)' @@ -960,6 +963,9 @@ subroutine init(self, dict, ptr, silent ) ptr_ceDatabase => ceNeutronDatabase_CptrCast(ptr) if(.not.associated(ptr_ceDatabase)) call fatalError(Here,"Should not happen. WTF?!") + ! Initialise dictionary to store S(a,b) info + call sabDict % init(1) + ! Create list of all nuclides. Loop over materials ! Find maximum number of nuclides: maxNuc do i = 1, mm_nMat() @@ -970,6 +976,11 @@ subroutine init(self, dict, ptr, silent ) do j = 1, size(mat % nuclides) name = self % makeNuclideName(mat % nuclides(j)) call nucSet % add(name, IN_SET) + + ! Add S(alpha,beta) data + if (mat % nuclides(j) % hasSab) then + call sabDict % store(name, mat % nuclides(j) % file_Sab) + end if end do end do @@ -1042,41 +1053,41 @@ subroutine init(self, dict, ptr, silent ) nucIdx = 1 do while (i /= nucSet % end()) - idx1 = index(nucSet % atKey(i),'+') - idx2 = index(nucSet % atKey(i),'#') - if (idx1 /= 0) then - name = trim(nucSet % atKey(i)) - if (idx2 == 0) then - name_file1 = trim(name(idx1+1:nameLen)) - else - name_file1 = trim(name(idx1+1:idx2-1)) - name_file2 = trim(name(idx2+1:nameLen)) - end if - name = name(1:idx1-1) - else - name = nucSet % atKey(i) + name = trim(nucSet % atKey(i)) + + ! Get Sab filenames + if (sabDict % isPresent(name)) then + call sabDict % get(fileNames, name) + + ! Truncate name to remove hash + idx = index(name,'_') + if (idx == 0) call fatalError(Here,'Failed to parse Sab nuclide: '//name) + name = name(1:idx-1) end if if (loud) then call statusMsg("Building: "// trim(name)// " with index: " //numToChar(nucIdx)) - if (idx1 /= 0 .and. idx2 == 0) & - call statusMsg("including S(alpha,beta) tables with file: " //trim(name_file1)) - if (idx1 /= 0 .and. idx2 /= 0) & - call statusMsg("including S(alpha,beta) tables with files: " //trim(name_file1)//' '//trim(name_file2)) + if (allocated(fileNames)) then + call statusMsg("including S(alpha,beta) tables with files: ") + do idx = 1, size(fileNames) + call statusMsg(fileNames(idx)) + end do + end if end if call new_neutronACE(ACE, name) call self % nuclides(nucIdx) % init(ACE, nucIdx, ptr_ceDatabase) ! Initialise S(alpha,beta) tables - if (idx1 /= 0 ) then - call new_moderACE(ACE_Sab1, name_file1) - if (idx2 /= 0) then - call new_moderACE(ACE_Sab2, name_file2) - call self % nuclides(nucIdx) % initSab(ACE_Sab1, ACE_Sab2) - else - call self % nuclides(nucIdx) % initSab(ACE_Sab1) - end if + if (allocated(fileNames)) then + allocate(ACE_Sab(size(fileNames))) + do idx = 1, size(fileNames) + call new_moderACE(ACE_Sab(idx), fileNames(idx)) + end do + call self % nuclides(nucIdx) % initSab(ACE_Sab) + fileNames = '' + deallocate(fileNames) + deallocate(ACE_Sab) end if ! Initialise probability tables @@ -1125,7 +1136,7 @@ subroutine init(self, dict, ptr, silent ) kT = mat % T * kBoltzmannMev if ((kT < sabT(1)) .or. (kT > sabT(2))) call fatalError(Here,& 'Material temperature must be bounded by the provided S(alpha,beta) data. '//& - 'The material temperature is '//numToChar(kT / kBoltzmannMeV)//& + 'The material temperature is '//numToChar(mat % T)//& 'K while the data bounds are '//numToChar(sabT(1) / kBoltzmannMeV)//& 'K and '//numToChar(sabT(2) / kBoltzmannMeV)//'K.') end if @@ -1219,7 +1230,8 @@ function makeNuclideName(self, nuclide) result(name) class(aceNeutronDatabase), intent(in) :: self type(nuclideInfo), intent(in) :: nuclide character(nameLen) :: name - character(:), allocatable :: file + character(:), allocatable :: extension + integer(shortInt) :: i, hash name = trim(nuclide % toChar()) @@ -1228,17 +1240,13 @@ function makeNuclideName(self, nuclide) result(name) ! scattering if (nuclide % hasSab) then - file = trim(nuclide % file_Sab1) - name = trim(name) // '+' // file - deallocate(file) - - ! Attach second Sab file for stochastic mixing - if (nuclide % sabMix) then - file = trim(nuclide % file_Sab2) - name = trim(name) // '#' // file - deallocate(file) - end if + extension = '' + do i = 1, size(nuclide % file_Sab) + extension = trim(extension)//'_'// trim(nuclide % file_Sab(i)) + end do + call FNV_1(extension, hash) + name = trim(name) //'_'// numToChar(hash) end if end function makeNuclideName diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index 6486c987c..4324550cc 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -321,6 +321,24 @@ elemental subroutine kill(self) if (allocated(self % mainData)) deallocate(self % mainData) if (allocated(self % eGrid)) deallocate(self % eGrid) call self % idxMT % kill() + + ! Clean S(alpha, beta) + if (allocated(self % thData)) then + do i = 1, size(self % thData) + call self % thData(i) % kill() + end do + deallocate(self % thData) + self % SabEl = ZERO + self % SabInel = ZERO + self % hasThData = .false. + self % stochasticMixing = .false. + end if + + ! Clean URR + self % urrE = ZERO + call self % probTab % kill() + self % hasProbTab = .false. + self % IFF = 0 end subroutine kill @@ -885,7 +903,7 @@ subroutine init(self, ACE, nucIdx, database) ! Find bottom and top of the MT cross section energy grid bottom = self % MTdata(i) % firstIdx - top = bottom + size(self % MTdata(i) % xs) + top = bottom + size(self % MTdata(i) % xs) - 1 if (MT >= 50 .and. MT <= 91) firstIdxMT4 = min(firstIdxMT4, bottom) ! Build total inelastic cross sections and MT4 from partials @@ -986,67 +1004,85 @@ subroutine initUrr(self, ACE) end subroutine initUrr !! - !! Initialise thermal scattering tables from ACE card + !! Initialise thermal scattering tables from ACE card. + !! If the array has size greater than 1, activates stochastic mixing treatment. !! !! Args: - !! ACE1 [inout] -> ACE S(a,b) card - !! ACE2 [inout] -> Optional second ACE S(a,b) card + !! ACE [inout] -> array of ACE S(a,b) cards !! !! Errors: !! fatalError if the inelastic scattering S(a,b) energy grid starts at a - !! lower energy than the nuclide energy grid - !! - subroutine initSab(self, ACE1, ACE2) - class(aceNeutronNuclide), intent(inout) :: self - class(aceSabCard), intent(inout) :: ACE1 - class(aceSabCard), intent(inout), optional :: ACE2 - real(defReal), dimension(2) :: EBounds - real(defReal) :: T1, T2 - type(thermalData) :: temp + !! lower energy than the nuclide energy grid. + !! + !! fatalError if any of the input temperatures are equal. + !! + subroutine initSab(self, ACE) + class(aceNeutronNuclide), intent(inout) :: self + class(aceSabCard), dimension(:), intent(inout) :: ACE + real(defReal), dimension(2) :: EBounds + real(defReal) :: T + integer(shortInt) :: nFiles + type(thermalData) :: temp + integer(shortInt) :: i, j character(100), parameter :: Here = "initSab (aceNeutronNuclide_class.f90)" - if (present(ACE2)) then - allocate(self % thData(2)) - else - allocate(self % thData(1)) - end if - - ! Initialise S(a,b) class from ACE file - call self % thData(1) % init(ACE1) self % hasThData = .true. + nFiles = size(ACE) + allocate(self % thData(nFiles)) + if (nFiles > 1) self % stochasticMixing = .true. - ! Initialise energy boundaries - self % SabInel = self % thData(1) % getEBounds('inelastic') - self % SabEl = self % thData(1) % getEBounds('elastic') + ! Initialise S(a,b) class from ACE file + do i = 1, nFiles + call self % thData(i) % init(ACE(i)) + + if (i == 1) then + ! Initialise energy boundaries + self % SabInel = self % thData(1) % getEBounds('inelastic') + self % SabEl = self % thData(1) % getEBounds('elastic') - ! Add second S(a,b) file for stochastic mixing - if (present(ACE2)) then + else + ! Ensure energy bounds are conservative + EBounds = self % thData(i) % getEBounds('inelastic') + if (EBounds(1) > self % SabInel(1)) self % SabInel(1) = EBounds(1) + if (EBounds(2) < self % SabInel(2)) self % SabInel(2) = EBounds(2) + + EBounds = self % thData(i) % getEbounds('elastic') + if (EBounds(1) > self % SabEl(1)) self % SabEl(1) = EBounds(1) + if (EBounds(2) < self % SabEl(2)) self % SabEl(2) = EBounds(2) - self % stochasticMixing = .true. - call self % thData(2) % init(ACE2) + end if - ! Ensure energy bounds are conservative - EBounds = self % thData(2) % getEBounds('inelastic') - if (EBounds(1) > self % SabInel(1)) self % SabInel(1) = EBounds(1) - if (EBounds(2) < self % SabInel(2)) self % SabInel(2) = EBounds(2) + end do - EBounds = self % thData(2) % getEbounds('elastic') - if (EBounds(1) > self % SabEl(1)) self % SabEl(1) = EBounds(1) - if (EBounds(2) < self % SabEl(2)) self % SabEl(2) = EBounds(2) + ! Rearrange files to be in order of ascending temperature + if (nFiles > 1) then + do i = 2, nFiles + temp = self % thData(i) + T = temp % getTemperature() + j = i - 1 - ! Identify which data is higher temperature and which is lower - ! 1 should be lower than 2 - swap if necessary - T1 = self % thData(1) % getTemperature() - T2 = self % thData(2) % getTemperature() + do while (j >= 1) + if (self % thData(j) % getTemperature() <= T) exit + self % thData(j+1) = self % thData(j) + j = j - 1 + end do - if (T1 > T2) then - temp = self % thData(1) - self % thData(1) = self % thData(2) - self % thData(2) = temp - end if + self % thData(j+1) = temp + end do + + ! Check to ensure all temperatures are unique + T = self % thData(1) % getTemperature() / kBoltzmannMeV + do i = 2, nFiles + if (abs(T - self % thData(i) % getTemperature()) / kBoltzmannMeV < 1.0E-6) call fatalError(& + Here,'S(alpha,beta) files have been input with the same temperature. '//& + 'Data '//numToChar(i-1)//': '//numToChar(T)//'. Data '//numToChar(i)//& + ': '//numToChar(self % thData(i) % getTemperature()/kBoltzmannMeV)) + T = self % thData(i) % getTemperature() / kBoltzmannMeV + end do end if + ! Check consistency of energy grid if (self % SabInel(1) < self % eGrid(1)) then call fatalError(Here, 'S(a,b) low energy boundary is lower than nuclide first energy point') @@ -1062,6 +1098,11 @@ end subroutine initSab !! Returns the Sab index for later consistent handling !! of angular distributions by storing in a cache !! + !! Errors: + !! fatalError if requested temperature is outside of the possible temperature range + !! + !! fatalError if the linear search over temperatures for interpolation fails + !! subroutine getSabPointer(self, kT, rand, ptr, idx) class(aceNeutronNuclide), intent(in), target :: self real(defReal), intent(in) :: kT @@ -1069,23 +1110,35 @@ subroutine getSabPointer(self, kT, rand, ptr, idx) type(thermalData), pointer, intent(out) :: ptr integer(shortInt), intent(out) :: idx real(defReal) :: kT1, kT2 + real(defReal), dimension(2) :: kTMax + integer(shortInt) :: i character(100), parameter :: Here = "getSabPointer (aceNeutronNuclide_class.f90)" if (self % stochasticMixing) then - kT1 = self % thData(1) % getTemperature() - kT2 = self % thData(2) % getTemperature() + kTMax = self % getSabTBounds() - if ((kT < kT1) .or. (kT > kT2)) call fatalError(Here,& + if ((kT < kTMax(1)) .or. (kT > kTMax(2))) call fatalError(Here,& 'Requested temperature '//numToChar(kT)//' not in temperature bounds: '//& - numToChar(kT1)//' and '//numToChar(kT2)) + numToChar(kTMax(1))//' and '//numToChar(kTMax(2))) + + ! Search through temperatures to find the interpolation range + do i = 1, size(self % thData) - 1 + if (kT <= self % thData(i+1) % getTemperature()) then + kT1 = self % thData(i) % getTemperature() + kT2 = self % thData(i+1) % getTemperature() + + if ((kT2 - kT)/(kT2 - kT1) > rand % get()) then + idx = i + else + idx = i + 1 + end if + ptr => self % thData(idx) + return + end if + end do - if ((kT2 - kT)/(kT2 - kT1) > rand % get()) then - ptr => self % thData(1) - idx = 1 - else - ptr => self % thData(2) - idx = 2 - end if + call fatalError(Here,'Temperature search failed. Requested temperature '//numToChar(kT)//& + ' exceed maximum: '//numToChar(self % thData(i+1) % getTemperature())) else ptr => self % thData(1) idx = 1 @@ -1097,13 +1150,13 @@ end subroutine getSabPointer !! Return the temperature bounds of S(alpha,beta) data !! If only one library, both bounds are the same temperature !! - function getSabTBounds(self) result(kT) + pure function getSabTBounds(self) result(kT) class(aceNeutronNuclide), intent(in) :: self real(defReal), dimension(2) :: kT kT(1) = self % thData(1) % getTemperature() if (self % stochasticMixing) then - kT(2) = self % thData(2) % getTemperature() + kT(2) = self % thData(size(self % thData)) % getTemperature() else kT(2) = kT(1) end if @@ -1115,7 +1168,7 @@ end function getSabTBounds !! !! Prints: !! nucIdx; Size of Energy grid; Maximum and Minumum energy; Mass; Temperature; ACE ZAID; - !! MT numbers used in tranposrt (active MTs) and not used in transport (inactiveMTs) + !! MT numbers used in transport (active MTs) and not used in transport (inactiveMTs) !! !! NOTE: For Now Formatting is horrible. Can be used only for debug !! MT = 18 may appear as inactive MT but is included from FIS block ! @@ -1149,7 +1202,6 @@ subroutine display(self) allMT = size(self % MTdata) print '(A)', "Active MTs: " // numToChar(self % MTdata(1:sMT) % MT) print '(A)', "Inactive MTs: "// numToChar(self % MTdata(sMT+1:allMT) % MT) - print '(A)', "This implementation ignores MT=5 (N,anything) !" end subroutine display diff --git a/NuclearData/emissionENDF/angleLawENDF/tabularAngle_class.f90 b/NuclearData/emissionENDF/angleLawENDF/tabularAngle_class.f90 index fde43798c..a7a802c8e 100644 --- a/NuclearData/emissionENDF/angleLawENDF/tabularAngle_class.f90 +++ b/NuclearData/emissionENDF/angleLawENDF/tabularAngle_class.f90 @@ -96,6 +96,10 @@ function sample(self,E,rand) result (mu) character(100),parameter :: Here='sample (tabularAngle_class.f90)' idx = binarySearch(self % eGrid,E) + if (idx < 1) then + print *, 'FUCK' + print *, E + end if call searchError(idx,Here) eps = (E - self % eGrid(idx)) / (self % eGrid(idx+1) - self % eGrid(idx)) diff --git a/NuclearData/materialMenu_mod.f90 b/NuclearData/materialMenu_mod.f90 index 3231ab5d1..4cd8b13a5 100644 --- a/NuclearData/materialMenu_mod.f90 +++ b/NuclearData/materialMenu_mod.f90 @@ -46,8 +46,7 @@ module materialMenu_mod !! T -> Evaluation number !! hasSab -> Does the nuclide have S(a,b) data? !! sabMix -> Does the nuclide mix S(a,b) data? - !! file_Sab1 -> First (and maybe only) S(a,b) file - !! file_Sab2 -> Second S(a,b) file + !! file_Sab -> Array of S(a,b) file names (generally >1 for stochastic interpolation) !! !! Interface: !! init -> build from a string @@ -58,8 +57,7 @@ module materialMenu_mod integer(shortInt) :: T = -1 logical(defBool) :: hasSab = .false. logical(defBool) :: sabMix = .false. - character(nameLen) :: file_Sab1 - character(nameLen) :: file_Sab2 + character(nameLen), allocatable, dimension(:) :: file_Sab contains procedure :: init => init_nuclideInfo procedure :: toChar => toChar_nuclideInfo @@ -354,17 +352,11 @@ subroutine init_materialItem(self, name, idx, dict) ! Check for stochastic mixing - this will depend on the ! size of the array of files produce call moderDict % get(filenames, keys(i)) - if (size(filenames) == 2) then - self % nuclides(i) % file_Sab1 = filenames(1) - self % nuclides(i) % file_Sab2 = filenames(2) - self % nuclides(i) % sabMix = .true. - elseif (size(filenames) == 1) then - self % nuclides(i) % file_Sab1 = filenames(1) - else - print *,filenames - call fatalError(Here,'Unexpectedly long moder contents. Should be 1 or 2 '//& - 'entries.') - end if + allocate(self % nuclides(i) % file_Sab(size(filenames))) + self % nuclides(i) % file_Sab = filenames + if (size(filenames) > 1) self % nuclides(i) % sabMix = .true. + filenames = '' + deallocate(filenames) end if ! Initialise the nuclides diff --git a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 index 6e1625d89..3b00ebb56 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 @@ -86,6 +86,7 @@ subroutine testBaseMgNeutronDatabaseWithP0() p % pRNG => pRNG p % type = P_NEUTRON p % G = 1 + p % isMG = .true. @assertEqual(2.1_defReal, database % getTrackingXS(p, 1, MATERIAL_XS), TOL) ! Test getting Total XS @@ -217,6 +218,7 @@ subroutine testBaseMgNeutronDatabaseWithP1() p % pRNG => pRNG p % type = P_NEUTRON p % G = 1 + p % isMG = .true. @assertEqual(2.1_defReal, database % getTrackingXS(p, 1, MATERIAL_XS), TOL) ! Test getting Total XS diff --git a/NuclearData/nuclearDataReg_mod.f90 b/NuclearData/nuclearDataReg_mod.f90 index 8b16e96e7..228204356 100644 --- a/NuclearData/nuclearDataReg_mod.f90 +++ b/NuclearData/nuclearDataReg_mod.f90 @@ -291,7 +291,7 @@ subroutine activate(type, name, activeMat, silent) if(idx == 0 ) then call fatalError(Here, trim(name)//' is was not defined. Cannot activate it!') else if(idx < 0) then - call fatalError(Here, '-ve idx from databaseNameMap. Quite immpossible. WTF?') + call fatalError(Here, '-ve idx from databaseNameMap. Quite impossible. WTF?') end if ! Make if it is not already made diff --git a/PhysicsPackages/alphaPhysicsPackage_class.f90 b/PhysicsPackages/alphaPhysicsPackage_class.f90 index f2aa219cc..d21789b19 100644 --- a/PhysicsPackages/alphaPhysicsPackage_class.f90 +++ b/PhysicsPackages/alphaPhysicsPackage_class.f90 @@ -150,7 +150,7 @@ subroutine run(self) call self % generateInitialState() call self % cycles(self % inactiveTally, self % inactiveAtch, self % N_inactive) call self % buildActiveAtch() - call self % couplingInfo % endCoupling() + call self % couplingInfo % moveToActive(self % activeTally) call self % cycles(self % activeTally, self % activeAtch, self % N_active) ! Collect results from other processes diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index 993424c17..2d0f1fa5d 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -151,11 +151,13 @@ subroutine run(self) call self % cycles(self % inactiveTally, self % inactiveAtch, self % N_inactive) - ! Deactivate coupling for active cycles - call self % couplingInfo % endCoupling() + ! Deactivate coupling for active cycles? + call self % couplingInfo % moveToActive(self % activeTally) call self % cycles(self % activeTally, self % activeAtch, self % N_active) + call self % couplingInfo % endCoupling() + ! Collect results from other processes call self % inactiveTally % collectDistributed() call self % activeTally % collectDistributed() diff --git a/TransportOperator/transportOperatorDT_class.f90 b/TransportOperator/transportOperatorDT_class.f90 index 8aefa39d7..345b50c63 100644 --- a/TransportOperator/transportOperatorDT_class.f90 +++ b/TransportOperator/transportOperatorDT_class.f90 @@ -91,12 +91,14 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle is in a region with overlapping cells case(OVERLAP_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index 1bc9c0015..3434db458 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -131,11 +131,13 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle somehow ended in an overlap material case(OVERLAP_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default @@ -251,12 +253,14 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle ended in an overlap material case(OVERLAP_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default diff --git a/TransportOperator/transportOperatorST_class.f90 b/TransportOperator/transportOperatorST_class.f90 index 80a35b7ef..63055bfbf 100644 --- a/TransportOperator/transportOperatorST_class.f90 +++ b/TransportOperator/transportOperatorST_class.f90 @@ -134,11 +134,13 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle is in a region with overlapping cells case(OVERLAP_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default From e528dd4c1c0fece9b3902b52d66a5187b21cb826 Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Tue, 30 Jun 2026 16:18:40 +0100 Subject: [PATCH 14/15] Finalising coupling --- Coupling/couplingAdmin_class.f90 | 22 ++++++++++++++-------- SharedModules/timer_mod.f90 | 11 ++++++----- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/Coupling/couplingAdmin_class.f90 b/Coupling/couplingAdmin_class.f90 index cf437fbfb..012319051 100644 --- a/Coupling/couplingAdmin_class.f90 +++ b/Coupling/couplingAdmin_class.f90 @@ -8,7 +8,7 @@ module couplingAdmin_class use errors_mod, only : fatalError use display_func, only : statusMsg use genericProcedures, only : numToChar - use timer_mod, only : c_sleep + use timer_mod, only : c_usleep use field_inter, only : field use fieldFactory_func, only : new_field @@ -42,7 +42,7 @@ module couplingAdmin_class !! !! Private members: !! isCoupled -> flag for whether coupling is being performed - !! endWhenActive -> flag to end coupling when moving to active cycles + !! endWhenActive -> flag to end coupling when moving to active cycles (false by default) !! commFileIn -> path to the file which triggers a SCONE calculation by the driver !! commFileOut -> path to the file which SCONE outputs to the driver on completing a calculation !! outputFile -> path to the file which SCONE outputs containing results to be used by the driver @@ -77,7 +77,7 @@ module couplingAdmin_class type, public :: couplingAdmin private logical(defBool) :: isCoupled = .false. - logical(defBool) :: endWhenActive = .true. + logical(defBool) :: endWhenActive = .false. character(pathLen) :: commFileIn = '' character(pathLen) :: commFileOut = '' character(pathLen) :: outputFile = '' @@ -133,7 +133,7 @@ subroutine init(self, dict) call dict % get(self % commFileOut, 'sendFile') ! Read whether to end when active - call dict % getOrDefault(self % endWhenActive, 'endActive', .true.) + call dict % getOrDefault(self % endWhenActive, 'endActive', .false.) ! Read how often to exchange data call dict % get(self % updateFreq, 'updateFreq') @@ -203,7 +203,7 @@ subroutine kill(self) class(couplingAdmin), intent(inout) :: self self % isCoupled = .false. - self % endWhenActive = .true. + self % endWhenActive = .false. self % updateFreq = 0 self % commFileIn = '' self % commFileOut = '' @@ -402,7 +402,7 @@ end subroutine signalCalculationOver subroutine waitForSignal(self) class(couplingAdmin), intent(inout) :: self logical(defBool) :: exists - integer(shortInt) :: unit + integer(shortInt) :: unit, ios character(nameLen) :: line character(100), parameter :: Here ='waitForSignal (couplingAdmin_class.f90)' @@ -415,7 +415,13 @@ subroutine waitForSignal(self) ! Read and delete the communication file open(newunit=unit, file=self % commFileIn, status="old") - read(unit, '(A)') line + read(unit, '(A)', iostat=ios) line + + ! Make sure the file isn't partially written: + if (ios /= 0) then + call c_usleep(int(0.005 * 1E6)) + cycle + end if ! Check whether to continue iteration or end coupling if (trim(adjustl(line)) == END_SIGNAL) then @@ -428,7 +434,7 @@ subroutine waitForSignal(self) return end if - call c_sleep(1) + call c_usleep(int(0.005 * 1E6)) end do diff --git a/SharedModules/timer_mod.f90 b/SharedModules/timer_mod.f90 index eac4f16c7..5735115cc 100644 --- a/SharedModules/timer_mod.f90 +++ b/SharedModules/timer_mod.f90 @@ -22,7 +22,7 @@ !! timerReset -> Reset a defined timer !! timerTime -> Get total elapsed time in Seconds !! -!! c_sleep -> pauses SCONE for several seconds +!! c_usleep -> pauses SCONE for several microseconds !! !! module timer_mod @@ -41,7 +41,7 @@ module timer_mod public :: timerStop public :: timerReset public :: timerTime - public :: c_sleep + public :: c_usleep !! !! This derived type allows to measure time in a program @@ -81,11 +81,12 @@ module timer_mod real(defReal),parameter :: GROWTH_RATIO = 1.6_defReal integer(shortInt),parameter :: MIN_SIZE = 5 - ! Interface to C sleep procedure + ! Interface to C usleep procedure + ! Takes microseconds as input interface - subroutine c_sleep(seconds) bind(C, name="sleep") + subroutine c_usleep(useconds) bind(C, name="usleep") import :: c_int - integer(c_int), value :: seconds + integer(c_int), value :: useconds end subroutine end interface From f973b7b594eb36e12f1000b82da52fa04665e2dd Mon Sep 17 00:00:00 2001 From: ChasingNeutrons Date: Tue, 30 Jun 2026 16:21:19 +0100 Subject: [PATCH 15/15] Modified docs --- docs/Input Manual.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index 0bdb45edb..fb069e907 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -1023,6 +1023,8 @@ in this dictionary are: * maxIt: (*optional*, defaults to infinity) the maximum number of MC cycles that will take place before SCONE terminates the coupling. SCONE will continue to run, but it will signal to the external solver that coupling has finished and fields will no longer be updated. +* endActive: (*optional*, defaults to false) decides whether to end coupling during active cycles. This can save + computational time due to the cost of communication, but it may also produce higher variance/bias results. An example of a coupling dictionary is: :: @@ -1042,6 +1044,7 @@ An example of a coupling dictionary is: :: fields (temperature density); maxTemp 2000; maxIt 100; + endActive 1; } Visualiser