Skip to content

Commit

Permalink
Reading average energy data from gnds files
Browse files Browse the repository at this point in the history
  • Loading branch information
whaeck committed Nov 15, 2024
1 parent 92b65c3 commit 00f82b9
Show file tree
Hide file tree
Showing 3 changed files with 173 additions and 60 deletions.
32 changes: 19 additions & 13 deletions src/dryad/format/gnds/createReactionProduct.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "dryad/format/gnds/createMultiplicity.hpp"
#include "dryad/format/gnds/createTwoBodyDistributionData.hpp"
#include "dryad/format/gnds/createUncorrelatedDistributionData.hpp"
#include "dryad/format/gnds/createTabulatedAverageEnergy.hpp"
#include "dryad/ReactionProduct.hpp"

namespace njoy {
Expand Down Expand Up @@ -41,41 +42,46 @@ namespace gnds {
// create the multiplicity
auto multiplicity = createMultiplicity( product.child( "multiplicity" ), style );

// get distribution
auto distribution = product.child( "distribution" );
if ( distribution ) {
// get distribution data
std::optional< ReactionProduct::DistributionData > distribution = std::nullopt;
auto node = product.child( "distribution" );
if ( node ) {

// get the first node of the requested style and act accordingly
auto node = distribution.find_child_by_attribute( "label", style.c_str() );
node = node.find_child_by_attribute( "label", style.c_str() );
if ( strcmp( node.name(), "angularTwoBody" ) == 0 ) {

// ignore recoil or regions2d distributions for now
auto recoil = node.child( "recoil" );
auto regions2d = node.child( "regions2d" );
if ( ! recoil && ! regions2d ) {

return ReactionProduct( id, multiplicity, createTwoBodyDistributionData( node ) );
distribution = createTwoBodyDistributionData( node );
}
}
if ( strcmp( node.name(), "uncorrelated" ) == 0 ) {
else if ( strcmp( node.name(), "uncorrelated" ) == 0 ) {

// ignore discrete gamma distributions for now
auto discreteGamma = node.child( "energy" ).child( "discreteGamma" );
auto primaryGamma = node.child( "energy" ).child( "primaryGamma" );
if ( ! discreteGamma && ! primaryGamma ) {

return ReactionProduct( id, multiplicity, createUncorrelatedDistributionData( node ) );
distribution = createUncorrelatedDistributionData( node );
}
}
else {
}

// for now, return a basic reaction product
return ReactionProduct( id, multiplicity );
}
// get distribution data
std::optional< TabulatedAverageEnergy > average = std::nullopt;
node = product.child( "averageProductEnergy" );
if ( node ) {

average = createTabulatedAverageEnergy( node );
}

// there is no distribution, return a basic reaction product
return ReactionProduct( id, multiplicity );
return ReactionProduct( id, multiplicity,
std::move( distribution ),
std::move( average ) );
}

} // gnds namespace
Expand Down
107 changes: 107 additions & 0 deletions src/dryad/format/gnds/createTabulatedAverageEnergy.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
#ifndef NJOY_DRYAD_FORMAT_GNDS_CREATETABULATEDAVERAGEENERGY
#define NJOY_DRYAD_FORMAT_GNDS_CREATETABULATEDAVERAGEENERGY

// system includes
#include <vector>

// other includes
#include "pugixml.hpp"
#include "tools/Log.hpp"
#include "dryad/format/gnds/createInterpolationType.hpp"
#include "dryad/format/gnds/readXYs1d.hpp"
#include "dryad/format/gnds/convertEnergies.hpp"
#include "dryad/TabulatedAverageEnergy.hpp"

namespace njoy {
namespace dryad {
namespace format {
namespace gnds {

/**
* @brief Create a TabulatedAverageEnergy from a GNDS average node
*/
static TabulatedAverageEnergy
createTabulatedAverageEnergy( const pugi::xml_node& average,
const std::string& style = "eval" ) {

std::vector< double > energies;
std::vector< double > values;
std::vector< std::size_t > boundaries;
std::vector< InterpolationType > interpolants;

// check that this is a valid average energy node
throwExceptionOnWrongNode( average, "averageProductEnergy" );

auto node = average.find_child_by_attribute( "label", style.c_str() );
if ( strcmp( node.name(), "XYs1d" ) == 0 ) {

// read the average energy data
auto data = readXYs1D( node );

// get the interpolation type
auto interpolant = createInterpolationType( std::get< 6 >( data ) );

// convert units - if necessary
convertEnergies( std::get< 2 >( data ), std::get< 3 >( data ) );

// assign data
energies = std::move( std::get< 2 >( data ) );
values = std::move( std::get< 4 >( data ) );
boundaries.emplace_back( energies.size() - 1 );
interpolants.emplace_back( interpolant );
}
else if ( strcmp( node.name(), "regions1d" ) == 0 ) {

// get the units
auto units = readAxes( node.child( "axes" ) );

// loop over the children of function1ds
pugi::xml_node function1ds = node.child( "function1ds" );
for ( pugi::xml_node xys1d = function1ds.child( "XYs1d" );
xys1d; xys1d = xys1d.next_sibling( "XYs1d" ) ) {

// read the current interpolation region
auto data = readXYs1D( xys1d, units );

// get the interpolation type
auto interpolant = createInterpolationType( std::get< 6 >( data ) );

// convert units - if necessary
convertEnergies( std::get< 2 >( data ), std::get< 3 >( data ) );

// check for duplicate points at interpolation region boundaries
std::size_t offset = 0;
if ( energies.size() > 0 ) {

if ( energies.back() == std::get< 2 >( data ).front() &&
values.back() == std::get< 4 >( data ).front() ) {

offset = 1;
}
}

// grow the data accordingly
energies.insert( energies.end(), std::get< 2 >( data ).begin() + offset, std::get< 2 >( data ).end() );
values.insert( values.end(), std::get< 4 >( data ).begin() + offset, std::get< 4 >( data ).end() );
boundaries.emplace_back( energies.size() - 1 );
interpolants.emplace_back( interpolant );
}
}
else {

Log::error( "Expected either an XYs1d node or regions1d node with XYs1d nodes"
"for average energy data" );
throw std::exception();
}

return TabulatedAverageEnergy(
std::move( energies ), std::move( values ),
std::move( boundaries ), std::move( interpolants ) );
}

} // gnds namespace
} // format namespace
} // dryad namespace
} // njoy namespace

#endif
94 changes: 47 additions & 47 deletions src/dryad/format/gnds/test/test_verification_functions.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -514,32 +514,32 @@ void verifyElectronBremsstrahlungReaction( const Reaction& bremsstrahlung ) {

auto electron = bremsstrahlung.products()[1];
CHECK( id::ParticleID( "e-" ) == electron.identifier() );
// CHECK( true == electron.isLinearised() );
// CHECK( true == electron.hasAverageEnergy() );
// CHECK( false == electron.hasDistributionData() );
// CHECK( true == std::holds_alternative< int >( electron.multiplicity() ) );
// multiplicity = std::get< int >( electron.multiplicity() );
// CHECK( 1 == multiplicity );
// CHECK( std::nullopt != electron.averageEnergy() );
// auto average = electron.averageEnergy().value();
// CHECK( true == average.isLinearised() );
// CHECK( 82 == average.numberPoints() );
// CHECK( 1 == average.numberRegions() );
// CHECK( 82 == average.energies().size() );
// CHECK( 82 == average.values().size() );
// CHECK( 1 == average.boundaries().size() );
// CHECK( 1 == average.interpolants().size() );
// CHECK( 81 == average.boundaries()[0] );
// CHECK( InterpolationType::LinearLinear == average.interpolants()[0] );
// CHECK_THAT( 10. , WithinRel( average.energies()[0] ) );
// CHECK_THAT( 14.4 , WithinRel( average.energies()[1] ) );
// CHECK_THAT( 7.86876E+10, WithinRel( average.energies()[80] ) );
// CHECK_THAT( 1e+11, WithinRel( average.energies()[81] ) );
// CHECK_THAT( 10. - 2.14426 , WithinRel( average.values()[0] ) );
// CHECK_THAT( 14.4 - 2.87181 , WithinRel( average.values()[1] ) );
// CHECK_THAT( 7.86876E+10 - 2.11850E+9, WithinRel( average.values()[80] ) );
// CHECK_THAT( 1e+11 - 2.66810E+9, WithinRel( average.values()[81] ) );
// CHECK( std::nullopt == electron.distributionData() );
CHECK( true == electron.isLinearised() );
CHECK( true == electron.hasAverageEnergy() );
CHECK( false == electron.hasDistributionData() );
CHECK( true == std::holds_alternative< int >( electron.multiplicity() ) );
multiplicity = std::get< int >( electron.multiplicity() );
CHECK( 1 == multiplicity );
CHECK( std::nullopt != electron.averageEnergy() );
auto average = electron.averageEnergy().value();
CHECK( true == average.isLinearised() );
CHECK( 82 == average.numberPoints() );
CHECK( 1 == average.numberRegions() );
CHECK( 82 == average.energies().size() );
CHECK( 82 == average.values().size() );
CHECK( 1 == average.boundaries().size() );
CHECK( 1 == average.interpolants().size() );
CHECK( 81 == average.boundaries()[0] );
CHECK( InterpolationType::LinearLinear == average.interpolants()[0] );
CHECK_THAT( 10. , WithinRel( average.energies()[0] ) );
CHECK_THAT( 14.4 , WithinRel( average.energies()[1] ) );
CHECK_THAT( 7.86876E+10, WithinRel( average.energies()[80] ) );
CHECK_THAT( 1e+11, WithinRel( average.energies()[81] ) );
CHECK_THAT( 10. - 2.14426 , WithinRel( average.values()[0] ) );
CHECK_THAT( 14.4 - 2.87181 , WithinRel( average.values()[1] ) );
CHECK_THAT( 7.86876E+10 - 2.11850E+9, WithinRel( average.values()[80] ) );
CHECK_THAT( 1e+11 - 2.66810E+9, WithinRel( average.values()[81] ) );
CHECK( std::nullopt == electron.distributionData() );

auto hydrogen = bremsstrahlung.products()[2];
CHECK( id::ParticleID( "H" ) == hydrogen.identifier() );
Expand Down Expand Up @@ -575,31 +575,31 @@ void verifyElectronExcitationReaction( const Reaction& subionisation ) {
auto electron = subionisation.products()[0];
CHECK( id::ParticleID( "e-" ) == electron.identifier() );
CHECK( true == electron.isLinearised() );
// CHECK( true == electron.hasAverageEnergy() );
CHECK( true == electron.hasAverageEnergy() );
CHECK( false == electron.hasDistributionData() );
CHECK( true == std::holds_alternative< int >( electron.multiplicity() ) );
auto multiplicity = std::get< int >( electron.multiplicity() );
CHECK( 1 == multiplicity );
// CHECK( std::nullopt != electron.averageEnergy() );
// auto average = electron.averageEnergy().value();
// CHECK( true == average.isLinearised() );
// CHECK( 170 == average.numberPoints() );
// CHECK( 1 == average.numberRegions() );
// CHECK( 170 == average.energies().size() );
// CHECK( 170 == average.values().size() );
// CHECK( 1 == average.boundaries().size() );
// CHECK( 1 == average.interpolants().size() );
// CHECK( 169 == average.boundaries()[0] );
// CHECK( InterpolationType::LinearLinear == average.interpolants()[0] );
// CHECK_THAT( 13.6 , WithinRel( average.energies()[0] ) );
// CHECK_THAT( 13.7740142 , WithinRel( average.energies()[1] ) );
// CHECK_THAT( 9e+10, WithinRel( average.energies()[168] ) );
// CHECK_THAT( 1e+11, WithinRel( average.energies()[169] ) );
// CHECK_THAT( 13.6 - 13.6 , WithinRel( average.values()[0] ) );
// CHECK_THAT( 13.7740142 - 13.6069279, WithinRel( average.values()[1] ) );
// CHECK_THAT( 9e+10 - 21.0777000, WithinRel( average.values()[168] ) );
// CHECK_THAT( 1e+11 - 21.0777000, WithinRel( average.values()[169] ) );
// CHECK( std::nullopt == electron.distributionData() );
CHECK( std::nullopt != electron.averageEnergy() );
auto average = electron.averageEnergy().value();
CHECK( true == average.isLinearised() );
CHECK( 170 == average.numberPoints() );
CHECK( 1 == average.numberRegions() );
CHECK( 170 == average.energies().size() );
CHECK( 170 == average.values().size() );
CHECK( 1 == average.boundaries().size() );
CHECK( 1 == average.interpolants().size() );
CHECK( 169 == average.boundaries()[0] );
CHECK( InterpolationType::LinearLinear == average.interpolants()[0] );
CHECK_THAT( 13.6 , WithinRel( average.energies()[0] ) );
CHECK_THAT( 13.7740142 , WithinRel( average.energies()[1] ) );
CHECK_THAT( 9e+10, WithinRel( average.energies()[168] ) );
CHECK_THAT( 1e+11, WithinRel( average.energies()[169] ) );
CHECK_THAT( 13.6 - 13.6 , WithinRel( average.values()[0] ) );
CHECK_THAT( 13.7740142 - 13.6069279, WithinRel( average.values()[1] ) );
CHECK_THAT( 9e+10 - 21.0777000, WithinRel( average.values()[168] ) );
CHECK_THAT( 1e+11 - 21.0777000, WithinRel( average.values()[169] ) );
CHECK( std::nullopt == electron.distributionData() );

auto hydrogen = subionisation.products()[1];
CHECK( id::ParticleID( "H" ) == hydrogen.identifier() );
Expand Down

0 comments on commit 00f82b9

Please sign in to comment.