Adding reading tabulated form factors (coherent real and imaginary pa…

…rts from gnds)

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 <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/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

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 )

src/dryad/format/gnds/test/createTabulatedFormFactorFromNodes.test.cpp

@@ -0,0 +1,61 @@
+// include Catch2
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+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] ) );
+}

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() );