From 00f82b94cdd1cd2655369a4c246c38f91f297c11 Mon Sep 17 00:00:00 2001 From: Wim Haeck Date: Fri, 15 Nov 2024 09:54:37 -0700 Subject: [PATCH] Reading average energy data from gnds files --- .../format/gnds/createReactionProduct.hpp | 32 +++--- .../gnds/createTabulatedAverageEnergy.hpp | 107 ++++++++++++++++++ .../gnds/test/test_verification_functions.hpp | 94 +++++++-------- 3 files changed, 173 insertions(+), 60 deletions(-) create mode 100644 src/dryad/format/gnds/createTabulatedAverageEnergy.hpp diff --git a/src/dryad/format/gnds/createReactionProduct.hpp b/src/dryad/format/gnds/createReactionProduct.hpp index 5ce9e7c..8545e7e 100644 --- a/src/dryad/format/gnds/createReactionProduct.hpp +++ b/src/dryad/format/gnds/createReactionProduct.hpp @@ -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 { @@ -41,12 +42,13 @@ 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 @@ -54,28 +56,32 @@ namespace gnds { 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 diff --git a/src/dryad/format/gnds/createTabulatedAverageEnergy.hpp b/src/dryad/format/gnds/createTabulatedAverageEnergy.hpp new file mode 100644 index 0000000..b374898 --- /dev/null +++ b/src/dryad/format/gnds/createTabulatedAverageEnergy.hpp @@ -0,0 +1,107 @@ +#ifndef NJOY_DRYAD_FORMAT_GNDS_CREATETABULATEDAVERAGEENERGY +#define NJOY_DRYAD_FORMAT_GNDS_CREATETABULATEDAVERAGEENERGY + +// system includes +#include + +// 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 diff --git a/src/dryad/format/gnds/test/test_verification_functions.hpp b/src/dryad/format/gnds/test/test_verification_functions.hpp index 075ebdc..808ea1d 100644 --- a/src/dryad/format/gnds/test/test_verification_functions.hpp +++ b/src/dryad/format/gnds/test/test_verification_functions.hpp @@ -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() ); @@ -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() );