diff --git a/src/dryad/format/gnds/createTabulatedFormFactorFromNodes.hpp b/src/dryad/format/gnds/createTabulatedFormFactorFromNodes.hpp new file mode 100644 index 0000000..c9cfba6 --- /dev/null +++ b/src/dryad/format/gnds/createTabulatedFormFactorFromNodes.hpp @@ -0,0 +1,102 @@ +#ifndef NJOY_DRYAD_FORMAT_GNDS_CREATETABULATEDFORMFACTORFROMNODES +#define NJOY_DRYAD_FORMAT_GNDS_CREATETABULATEDFORMFACTORFROMNODES + +// 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/TabulatedFormFactor.hpp" + +namespace njoy { +namespace dryad { +namespace format { +namespace gnds { + + /** + * @brief Create a TabulatedFormFactor from a GNDS form factor node + */ + static TabulatedFormFactor + createTabulatedFormFactorFromNodes( const pugi::xml_node& node ) { + + std::vector< double > energies; + std::vector< double > values; + std::vector< std::size_t > boundaries; + std::vector< InterpolationType > interpolants; + + if ( strcmp( node.name(), "XYs1d" ) == 0 ) { + + // read the form factor 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 tabulated form factor data data" ); + throw std::exception(); + } + + return TabulatedFormFactor( + 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/CMakeLists.txt b/src/dryad/format/gnds/test/CMakeLists.txt index 8e5c154..a9de5b9 100644 --- a/src/dryad/format/gnds/test/CMakeLists.txt +++ b/src/dryad/format/gnds/test/CMakeLists.txt @@ -14,6 +14,8 @@ add_cpp_test( format.gnds.readConstant1d readLegendre.test.cpp ) add_cpp_test( format.gnds.readXYs1d readXYs1d.test.cpp ) add_cpp_test( format.gnds.readLegendre readLegendre.test.cpp ) +add_cpp_test( format.gnds.createTabulatedFormFactorFromNodes createTabulatedFormFactorFromNodes.test.cpp ) + add_cpp_test( format.gnds.createTabulatedCrossSection createTabulatedCrossSection.test.cpp ) add_cpp_test( format.gnds.createMultiplicity createMultiplicity.test.cpp ) add_cpp_test( format.gnds.createQValue createQValue.test.cpp ) diff --git a/src/dryad/format/gnds/test/createTabulatedFormFactorFromNodes.test.cpp b/src/dryad/format/gnds/test/createTabulatedFormFactorFromNodes.test.cpp new file mode 100644 index 0000000..6703568 --- /dev/null +++ b/src/dryad/format/gnds/test/createTabulatedFormFactorFromNodes.test.cpp @@ -0,0 +1,61 @@ +// include Catch2 +#include +#include +using Catch::Matchers::WithinRel; + +// what we are testing +#include "dryad/format/gnds/createTabulatedFormFactorFromNodes.hpp" + +// other includes +#include "pugixml.hpp" + +// convenience typedefs +using namespace njoy::dryad; + +void verifyChunk( const TabulatedFormFactor& ); + +SCENARIO( "createTabulatedFormFactorFromNodes" ) { + + GIVEN( "GNDS average energy node from photoatomic data" ) { + + pugi::xml_document document; + pugi::xml_parse_result result = document.load_file( "photoat-001_H_000.endf.gnds.xml" ); + pugi::xml_node node = document.child( "reactionSuite" ).child( "reactions" ). + find_child_by_attribute( "reaction", "ENDF_MT", "502" ). + child( "doubleDifferentialCrossSection" ). + child( "coherentPhotonScattering" ). + child( "realAnomalousFactor" ). + first_child(); + + WHEN( "a single average energy node is given" ) { + + THEN( "it can be converted" ) { + + auto chunk = format::gnds::createTabulatedFormFactorFromNodes( node ); + + verifyChunk( chunk ); + } // THEN + } // WHEN + } // GIVEN +} // SCENARIO + +void verifyChunk( const TabulatedFormFactor& chunk ) { + + CHECK( true == chunk.isLinearised() ); + CHECK( 297 == chunk.numberPoints() ); + CHECK( 1 == chunk.numberRegions() ); + CHECK( 297 == chunk.energies().size() ); + CHECK( 297 == chunk.values().size() ); + CHECK( 1 == chunk.boundaries().size() ); + CHECK( 1 == chunk.interpolants().size() ); + CHECK( 296 == chunk.boundaries()[0] ); + CHECK( InterpolationType::LinearLinear == chunk.interpolants()[0] ); + CHECK_THAT( 1. , WithinRel( chunk.energies()[0] ) ); + CHECK_THAT( 2. , WithinRel( chunk.energies()[1] ) ); + CHECK_THAT( 9549925.86, WithinRel( chunk.energies()[295] ) ); + CHECK_THAT( 1e+7 , WithinRel( chunk.energies()[296] ) ); + CHECK_THAT( -1.00260813, WithinRel( chunk.values()[0] ) ); + CHECK_THAT( -1.01054501, WithinRel( chunk.values()[1] ) ); + CHECK_THAT( 2.8024E-11, WithinRel( chunk.values()[295] ) ); + CHECK_THAT( 0. , WithinRel( chunk.values()[296] ) ); +} diff --git a/src/dryad/format/gnds/test/test_verification_functions.hpp b/src/dryad/format/gnds/test/test_verification_functions.hpp index 808ea1d..2c02f19 100644 --- a/src/dryad/format/gnds/test/test_verification_functions.hpp +++ b/src/dryad/format/gnds/test/test_verification_functions.hpp @@ -783,8 +783,8 @@ void verifyPhotonCoherentReaction( const Reaction& coherent ) { auto multiplicity = std::get< int >( gamma.multiplicity() ); CHECK( 1 == multiplicity ); -// CHECK( std::nullopt == gamma.averageEnergy() ); -// + CHECK( std::nullopt == gamma.averageEnergy() ); + // CHECK( std::nullopt != gamma.distributionData() ); // CHECK( true == std::holds_alternative< CoherentDistributionData >( gamma.distributionData().value() ) ); // auto data = std::get< CoherentDistributionData >( gamma.distributionData().value() ); @@ -891,8 +891,8 @@ void verifyPhotonIncoherentReaction( const Reaction& incoherent ) { auto multiplicity = std::get< int >( gamma.multiplicity() ); CHECK( 1 == multiplicity ); -// CHECK( std::nullopt == gamma.averageEnergy() ); -// + CHECK( std::nullopt == gamma.averageEnergy() ); + // CHECK( std::nullopt != gamma.distributionData() ); // CHECK( true == std::holds_alternative< IncoherentDistributionData >( gamma.distributionData().value() ) ); // auto data = std::get< IncoherentDistributionData >( gamma.distributionData().value() );