Skip to content

Commit e1bb501

Browse files
amlewis1578amlewis1578
committed
Merge branch 'feature/lrf3' into 'develop'
feature lrf3 See merge request njoy/dryad!105
2 parents 7381572 + 455f523 commit e1bb501

File tree

19 files changed

+95406
-55
lines changed

19 files changed

+95406
-55
lines changed

‎cmake/unit_testing.cmake‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -164,6 +164,7 @@ add_cpp_test( dryad.format.endf.createProductIdentifier
164164
add_cpp_test( dryad.format.endf.createTargetIdentifier dryad/format/endf/createTargetIdentifier.test.cpp )
165165
add_cpp_test( dryad.format.endf.createInteractionType dryad/format/endf/createInteractionType.test.cpp )
166166
add_cpp_test( dryad.format.endf.resonances.createTabulatedRadius dryad/format/endf/resonances/createTabulatedRadius.test.cpp )
167+
add_cpp_test( dryad.format.endf.resonances.lrf3.createSpinGroups dryad/format/endf/resonances/lrf3/createSpinGroups.test.cpp )
167168
add_cpp_test( dryad.format.endf.resonances.lrf7.createBoundaryCondition dryad/format/endf/resonances/lrf7/createBoundaryCondition.test.cpp )
168169
add_cpp_test( dryad.format.endf.resonances.lrf7.createFormalism dryad/format/endf/resonances/lrf7/createFormalism.test.cpp )
169170
add_cpp_test( dryad.format.endf.resonances.lrf7.createKinematics dryad/format/endf/resonances/lrf7/createKinematics.test.cpp )

‎python/test/dryad/resonances/Test_ResonanceTable.py‎

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -253,14 +253,6 @@ def test_failures( self ) :
253253
energies = [ 1., 2., 3., 4. ],
254254
amplitudes = [ [ 11., 12., 13., 14. ], [ 21., 22., 23., 24. ] ] )
255255

256-
# the energies are not unique
257-
with self.assertRaises( Exception ) :
258-
259-
table = ResonanceTable( channels = [ ChannelID( 'n,U235->n,U235{0,1/2,1/2+}' ),
260-
ChannelID( 'n,U235->n,U235_e1{0,1/2,1/2+}' ) ],
261-
energies = [ 1., 2., 2., 4. ],
262-
amplitudes = [ [ 11., 12., 13., 14. ], [ 21., 22., 23., 24. ] ] )
263-
264256
# the number of energies and number of amplitudes is inconsistent
265257
with self.assertRaises( Exception ) :
266258

‎src/njoy/dryad/format/endf/resonances/createChannelRadii.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ namespace resonances {
2020
* @brief Create the channel radii based on ENDF information
2121
*
2222
* @param[in] naps the channel radius option as given in the ENDF file
23-
* @param[in] nro the energy dependent scattering radius (if defined)
23+
* @param[in] nro the energy dependent scattering radius (if defined, given in fm)
2424
* @param[in] ap the l-dependent scattering radius (if defined, given in fm)
2525
* @param[in] awr the atomic weight ratio as given in the ENDF file
2626
*/
@@ -47,7 +47,7 @@ namespace resonances {
4747
case 2 : return dryad::resonances::ChannelRadii( ap, nro.value() );
4848
default : {
4949

50-
Log::error( "Encountered unknown value for NAPS = {}", naps );
50+
Log::error( "Encountered unknown value for NAPS = {} with NRO = {}", naps, nro.has_value() );
5151
throw std::exception();
5252
}
5353
}
@@ -62,7 +62,7 @@ namespace resonances {
6262
case 1 : return dryad::resonances::ChannelRadii( ap );
6363
default : {
6464

65-
Log::error( "Encountered unknown value for NAPS = {}", naps );
65+
Log::error( "Encountered unknown value for NAPS = {} with NRO = {}", naps, nro.has_value() );
6666
throw std::exception();
6767
}
6868
}

‎src/njoy/dryad/format/endf/resonances/createResonanceParameters.hpp‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,8 @@
77
// other includes
88
#include "tools/Log.hpp"
99
#include "njoy/dryad/resonances/ResonanceParameters.hpp"
10+
#include "njoy/dryad/format/endf/resonances/createTabulatedRadius.hpp"
11+
#include "njoy/dryad/format/endf/resonances/lrf3/createCompoundSystem.hpp"
1012
#include "njoy/dryad/format/endf/resonances/lrf7/createCompoundSystem.hpp"
1113
#include "ENDFtk/section/2/151.hpp"
1214

@@ -35,11 +37,24 @@ namespace resonances {
3537
double lower = range.lowerEnergy();
3638
double upper = range.upperEnergy();
3739

40+
auto naps = range.scatteringRadiusCalculationOption();
41+
std::optional< dryad::resonances::TabulatedRadius > nro = std::nullopt;
42+
if ( range.scatteringRadius().has_value() ) {
43+
44+
nro = createTabulatedRadius( range.scatteringRadius().value() );
45+
}
46+
3847
if ( range.type() == 1 ) {
3948

4049
Log::info( "Reading resolved resonance region between {} and {} eV", lower, upper );
4150
switch ( range.representation() ) {
4251

52+
case 3 : {
53+
54+
decltype(auto) parameters = std::get< njoy::ENDFtk::section::Type<2,151>::ReichMoore >( range.parameters() );
55+
resolved.emplace_back( lrf3::createCompoundSystem( projectile, target, lower, upper, naps, nro, parameters ) );
56+
break;
57+
}
4358
case 7 : {
4459

4560
decltype(auto) parameters = std::get< njoy::ENDFtk::section::Type<2,151>::RMatrixLimited >( range.parameters() );

‎src/njoy/dryad/format/endf/resonances/lrf3/createChannelData.hpp‎

Lines changed: 180 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
#ifndef NJOY_DRYAD_FORMAT_ENDF_RESONANCES_LRF3_CREATECHANNELDATA
2+
#define NJOY_DRYAD_FORMAT_ENDF_RESONANCES_LRF3_CREATECHANNELDATA
3+
4+
// system includes
5+
#include <algorithm>
6+
#include <vector>
7+
8+
// other includes
9+
#include "tools/Log.hpp"
10+
#include "njoy/dryad/id/ReactionID.hpp"
11+
#include "njoy/dryad/id/ChannelID.hpp"
12+
#include "njoy/dryad/resonances/SpinGroup.hpp"
13+
#include "njoy/dryad/format/createVector.hpp"
14+
#include "njoy/dryad/format/endf/resonances/createChannelRadii.hpp"
15+
#include "ENDFtk/section/2/151.hpp"
16+
17+
namespace njoy {
18+
namespace dryad {
19+
namespace format {
20+
namespace endf {
21+
namespace resonances {
22+
namespace lrf3 {
23+
24+
inline dryad::resonances::ChannelQuantumNumbers
25+
retrieveQuantumNumber( unsigned int l, double j,
26+
std::vector< dryad::resonances::ChannelQuantumNumbers >& available ) {
27+
28+
auto find = [l,j] ( auto&& numbers ) {
29+
30+
return numbers.orbitalAngularMomentum() == l &&
31+
numbers.totalAngularMomentum() == std::abs( j );
32+
};
33+
34+
auto first = std::find_if( available.begin(), available.end(), find );
35+
if ( first != available.end() ) {
36+
37+
auto second = std::find_if( first, available.end(), find );
38+
if ( second != available.end() ) {
39+
40+
if ( j > 0 ) {
41+
42+
first = second;
43+
}
44+
}
45+
46+
auto numbers = *first;
47+
available.erase( first );
48+
return numbers;
49+
}
50+
else {
51+
52+
throw std::runtime_error( "None of the expected spin groups has l = "
53+
+ std::to_string( l ) + " and J = "
54+
+ std::to_string( std::abs( j ) ) );
55+
}
56+
}
57+
58+
/**
59+
* @brief Create the channel data for Reich-Moore data for a given l value
60+
*
61+
* @param[in] projectile the projectile identifier
62+
* @param[in] target the target identifier
63+
* @param[in] incident the incident particle pair
64+
* @param[in] naps the channel radius option as given in the ENDF file
65+
* @param[in] nro the energy dependent scattering radius (if defined, given in fm)
66+
* @param[in] ap the l-dependent scattering radius (if defined, given in fm)
67+
* @param[in,out] available the quantum number combinations that are still available
68+
* @param[in] endfReichMooreLValue the parsed ENDF Reich-Moore l-value data
69+
*/
70+
inline auto createChannelData(
71+
const id::ParticleID& projectile,
72+
const id::ParticleID& target,
73+
const dryad::resonances::ParticlePair& incident,
74+
int naps,
75+
const std::optional< dryad::resonances::TabulatedRadius >& nro,
76+
double ap,
77+
std::vector< dryad::resonances::ChannelQuantumNumbers >& available,
78+
const ENDFtk::section::Type< 2, 151 >::ReichMooreLValue& endfReichMooreLValue ) {
79+
80+
std::vector< dryad::resonances::SpinGroup::ChannelData > channel_data;
81+
82+
unsigned int l = endfReichMooreLValue.orbitalMomentum();
83+
double apl = endfReichMooreLValue.lDependentScatteringRadius() * constants::deca;
84+
double awri = endfReichMooreLValue.atomicWeightRatio();
85+
86+
// create the channel radii (check for l-dependent radius)
87+
dryad::resonances::ChannelRadii radii = createChannelRadii( naps, nro, apl != 0. ? apl : ap, awri );
88+
89+
// get all possible total angular momentum values
90+
std::vector< double > jvalues = createVector( endfReichMooreLValue.spinValues() );
91+
std::sort( jvalues.begin(), jvalues.end() );
92+
jvalues.erase( std::unique( jvalues.begin(), jvalues.end() ), jvalues.end() );
93+
94+
// lambda to calculate reduced width
95+
auto reduced_width = [] ( auto&& width, auto&& penetrability ) {
96+
97+
int sign = width >= 0. ? +1 : -1;
98+
return sign * std::sqrt( std::abs( width ) / 2. / penetrability );
99+
};
100+
101+
// go over each one to create channel data
102+
for ( double j : jvalues ) {
103+
104+
// elastic channel
105+
id::ChannelID elastic_id( id::ReactionID( projectile, target, 2 ),
106+
retrieveQuantumNumber( l, j, available ) );
107+
dryad::resonances::Channel elastic( elastic_id, incident, incident, 0., std::nullopt, radii );
108+
109+
// collect level energies and widths
110+
std::vector< double > energies;
111+
std::vector< double > elastic_widths;
112+
std::vector< double > capture_widths;
113+
std::vector< double > fission1_widths;
114+
std::vector< double > fission2_widths;
115+
116+
// go over all resonances
117+
for ( unsigned int i = 0; i < endfReichMooreLValue.numberResonances(); ++i ) {
118+
119+
if ( endfReichMooreLValue.spinValues()[i] == j ) {
120+
121+
energies.emplace_back( endfReichMooreLValue.resonanceEnergies()[i] );
122+
elastic_widths.emplace_back( reduced_width( endfReichMooreLValue.neutronWidths()[i],
123+
elastic.penetrability( energies.back() ) ) );
124+
capture_widths.emplace_back( reduced_width( endfReichMooreLValue.gammaWidths()[i], 1. ) );
125+
fission1_widths.emplace_back( reduced_width( endfReichMooreLValue.firstFissionWidths()[i], 1. ) );
126+
fission2_widths.emplace_back( reduced_width( endfReichMooreLValue.secondFissionWidths()[i], 1. ) );
127+
}
128+
}
129+
130+
// treat elastic
131+
channel_data.emplace_back( std::move( elastic ), dryad::resonances::ResonanceTable{ { elastic_id }, energies, std::move( elastic_widths ) } );
132+
133+
// treat capture
134+
dryad::resonances::ChannelQuantumNumbers other( l, 0, std::abs( j ), l%2 == 0 ? +1 : -1 );
135+
id::ChannelID capture_id( id::ReactionID( projectile, target, 102 ), other );
136+
dryad::resonances::ParticlePair capture_pair( { id::ParticleID::photon(), 0., 0., +1 },
137+
{ capture_id.reaction().residual().value(), 0., 0., +1 } );
138+
dryad::resonances::Channel capture( capture_id, incident, capture_pair, 0., std::nullopt, radii );
139+
channel_data.emplace_back( std::move( capture ), dryad::resonances::ResonanceTable{ { capture_id }, energies, std::move( capture_widths ) } );
140+
141+
// check for fission
142+
auto is_non_zero = [] ( auto&& value ) { return value != 0.; };
143+
bool has_fission1 = std::any_of( fission1_widths.begin(), fission1_widths.end(), is_non_zero );
144+
bool has_fission2 = std::any_of( fission2_widths.begin(), fission2_widths.end(), is_non_zero );
145+
if ( has_fission1 && has_fission2 ) {
146+
147+
id::ChannelID fission1_id( id::ReactionID( projectile, target, 18 ), other, 0 );
148+
id::ChannelID fission2_id( id::ReactionID( projectile, target, 18 ), other, 1 );
149+
dryad::resonances::Channel fission1( fission1_id, incident, std::nullopt, 0., std::nullopt, radii );
150+
dryad::resonances::Channel fission2( fission2_id, incident, std::nullopt, 0., std::nullopt, radii );
151+
channel_data.emplace_back( std::move( fission1 ), dryad::resonances::ResonanceTable{ { fission1_id }, energies, std::move( fission1_widths ) } );
152+
channel_data.emplace_back( std::move( fission2 ), dryad::resonances::ResonanceTable{ { fission2_id }, energies, std::move( fission2_widths ) } );
153+
154+
}
155+
else if ( has_fission1 || has_fission2 ) {
156+
157+
id::ChannelID fission_id( id::ReactionID( projectile, target, 18 ), other );
158+
dryad::resonances::Channel fission( fission_id, incident, std::nullopt, 0., std::nullopt, radii );
159+
if ( has_fission1 ) {
160+
161+
channel_data.emplace_back( std::move( fission ), dryad::resonances::ResonanceTable{ { fission_id }, energies, std::move( fission1_widths ) } );
162+
}
163+
else {
164+
165+
channel_data.emplace_back( std::move( fission ), dryad::resonances::ResonanceTable{ { fission_id }, energies, std::move( fission2_widths ) } );
166+
}
167+
}
168+
}
169+
170+
return channel_data;
171+
}
172+
173+
} // lrf3 namespace
174+
} // resonances namespace
175+
} // endf namespace
176+
} // format namespace
177+
} // dryad namespace
178+
} // njoy namespace
179+
180+
#endif

‎src/njoy/dryad/format/endf/resonances/lrf3/createCompoundSystem.hpp‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
#ifndef NJOY_DRYAD_FORMAT_ENDF_RESONANCES_LRF3_CREATECOMPOUNDSYSTEM
2+
#define NJOY_DRYAD_FORMAT_ENDF_RESONANCES_LRF3_CREATECOMPOUNDSYSTEM
3+
4+
// system includes
5+
6+
// other includes
7+
#include "tools/Log.hpp"
8+
#include "njoy/dryad/resonances/CompoundSystem.hpp"
9+
#include "njoy/dryad/format/endf/resonances/lrf3/createSpinGroups.hpp"
10+
#include "ENDFtk/section/2/151.hpp"
11+
12+
namespace njoy {
13+
namespace dryad {
14+
namespace format {
15+
namespace endf {
16+
namespace resonances {
17+
namespace lrf3 {
18+
19+
/**
20+
* @brief Create the compound system for LRF3 resonance parameters
21+
*
22+
* @param[in] projectile the projectile identifier
23+
* @param[in] target the target identifier
24+
* @param[in] lower the lower energy limit
25+
* @param[in] upper the upper energy limit
26+
* @param[in] naps the channel radius option as given in the ENDF file
27+
* @param[in] nro the energy dependent scattering radius (if defined, given in fm)
28+
* @param[in] endf the parsed ENDF LRF3 data
29+
*/
30+
inline auto createCompoundSystem( const id::ParticleID& projectile,
31+
const id::ParticleID& target,
32+
double lower,
33+
double upper,
34+
int naps,
35+
const std::optional< dryad::resonances::TabulatedRadius >& nro,
36+
const ENDFtk::section::Type< 2, 151 >::ReichMoore& endf ) {
37+
38+
auto groups = lrf3::createSpinGroups( projectile, target, naps, nro, endf );
39+
40+
return dryad::resonances::CompoundSystem( lower, upper, std::move( groups ) );
41+
}
42+
43+
} // lrf3 namespace
44+
} // resonances namespace
45+
} // endf namespace
46+
} // format namespace
47+
} // dryad namespace
48+
} // njoy namespace
49+
50+
#endif

0 commit comments

Comments
 (0)