Adding zero and first moment integration to the AnalyticalIntegrator

njoy · Oct 1, 2024 · 5808c54 · 5808c54
1 parent 62f0c1a
commit 5808c54
Show file tree

Hide file tree

Showing 6 changed files with 108 additions and 5 deletions.
diff --git a/src/scion/integration.hpp b/src/scion/integration.hpp
@@ -4,5 +4,6 @@
 #include "scion/integration/LinearLogarithmic.hpp"
 #include "scion/integration/LogarithmicLinear.hpp"
 #include "scion/integration/LogarithmicLogarithmic.hpp"
+#include "scion/integration/FirstMomentLinearLinear.hpp"
 #include "scion/integration/GaussLegendre.hpp"
 #include "scion/integration/GaussLobatto.hpp"
diff --git a/src/scion/integration/FirstMomentLinearLinear.hpp b/src/scion/integration/FirstMomentLinearLinear.hpp
@@ -1,5 +1,5 @@
-#ifndef NJOY_SCION_INTEGRATION_LINEARLINEAR
-#define NJOY_SCION_INTEGRATION_LINEARLINEAR
+#ifndef NJOY_SCION_INTEGRATION_FIRSTMOMENTLINEARLINEAR
+#define NJOY_SCION_INTEGRATION_FIRSTMOMENTLINEARLINEAR
 
 // system includes
 
@@ -40,7 +40,7 @@ namespace integration {
       auto slope = ( yRight - yLeft ) / delta / 3.;
       auto constant = 0.5 * ( xRight * yLeft - xLeft * yRight ) / delta;
       return xRight * xRight * ( slope * xRight + constant )
-             - xLeft * xLeft * ( slope * yLeft + constant );
+             - xLeft * xLeft * ( slope * xLeft + constant );
     }
 
   public:

diff --git a/src/scion/math/AnalyticalIntegrator.hpp b/src/scion/math/AnalyticalIntegrator.hpp
@@ -14,7 +14,7 @@ namespace math {
 
   /**
    *  @class
-   *  @brief Analytical integration of data over a number of successive
+   *  @brief Analytical integration of tabulated data over a number of successive
    *         integration intervals
    */
   template < typename X >
@@ -53,6 +53,8 @@ namespace math {
     }
 
     #include "scion/math/AnalyticalIntegrator/src/call.hpp"
+    #include "scion/math/AnalyticalIntegrator/src/zerothMoment.hpp"
+    #include "scion/math/AnalyticalIntegrator/src/firstMoment.hpp"
   };
 
 } // math namespace

diff --git a/src/scion/math/AnalyticalIntegrator/src/firstMoment.hpp b/src/scion/math/AnalyticalIntegrator/src/firstMoment.hpp
@@ -0,0 +1,9 @@
+template < typename Y,
+           typename I = decltype( std::declval< X >() * std::declval< X >() * std::declval< Y >() ) >
+std::vector< I > firstMoment( const InterpolationTable< X, Y >& table ) const {
+
+  auto interpolator = interpolation::LinearLinear();
+  auto integrator = integration::FirstMomentLinearLinear();
+
+  return this->integrate( interpolator, integrator, table );
+}
diff --git a/src/scion/math/AnalyticalIntegrator/src/zerothMoment.hpp b/src/scion/math/AnalyticalIntegrator/src/zerothMoment.hpp
@@ -0,0 +1,6 @@
+template < typename Y,
+           typename I = decltype( std::declval< X >() * std::declval< Y >() ) >
+std::vector< I > zerothMoment( const InterpolationTable< X, Y >& table ) const {
+
+  return this->operator()( table );
+}
diff --git a/src/scion/math/AnalyticalIntegrator/test/AnalyticalIntegrator.test.cpp b/src/scion/math/AnalyticalIntegrator/test/AnalyticalIntegrator.test.cpp
@@ -53,7 +53,7 @@ SCENARIO( "AnalyticalIntegrator" ) {
     }
   }
 
-  GIVEN( "tabulated data" ) {
+  GIVEN( "tabulated data to be integrated using the call operator" ) {
 
     std::vector< double > x = { 1., 2., 3., 4. };
     std::vector< double > y = { 4., 3., 2., 1. };
@@ -73,6 +73,20 @@ SCENARIO( "AnalyticalIntegrator" ) {
         CHECK_THAT( 3.5, WithinRel( result[0] ) );
         CHECK_THAT( 2.5, WithinRel( result[1] ) );
         CHECK_THAT( 1.5, WithinRel( result[2] ) );
+
+        result = integrator.zerothMoment( table );
+
+        CHECK( 3 == result.size() );
+        CHECK_THAT( 3.5, WithinRel( result[0] ) );
+        CHECK_THAT( 2.5, WithinRel( result[1] ) );
+        CHECK_THAT( 1.5, WithinRel( result[2] ) );
+
+        result = integrator.firstMoment( table );
+
+        CHECK( 3 == result.size() );
+        CHECK_THAT( 5. + 1./6., WithinRel( result[0] ) );
+        CHECK_THAT( 6. + 1./6., WithinRel( result[1] ) );
+        CHECK_THAT( 5. + 1./6., WithinRel( result[2] ) );
       } // THEN
     } // WHEN
 
@@ -89,6 +103,18 @@ SCENARIO( "AnalyticalIntegrator" ) {
         CHECK( 2 == result.size() );
         CHECK_THAT( 3., WithinRel( result[0] ) );
         CHECK_THAT( 2., WithinRel( result[1] ) );
+
+        result = integrator.zerothMoment( table );
+
+        CHECK( 2 == result.size() );
+        CHECK_THAT( 3., WithinRel( result[0] ) );
+        CHECK_THAT( 2., WithinRel( result[1] ) );
+
+        result = integrator.firstMoment( table );
+
+        CHECK( 2 == result.size() );
+        CHECK_THAT( 5. + 11./12., WithinRel( result[0] ) );
+        CHECK_THAT( 5. + 11./12., WithinRel( result[1] ) );
       } // THEN
     } // WHEN
 
@@ -108,6 +134,24 @@ SCENARIO( "AnalyticalIntegrator" ) {
         CHECK_THAT( 3.5, WithinRel( result[2] ) );
         CHECK_THAT( 2.5, WithinRel( result[3] ) );
         CHECK_THAT( 1.5, WithinRel( result[4] ) );
+
+        result = integrator.zerothMoment( table );
+
+        CHECK( 5 == result.size() );
+        CHECK_THAT( 0.0, WithinRel( result[0] ) );
+        CHECK_THAT( 0.0, WithinRel( result[1] ) );
+        CHECK_THAT( 3.5, WithinRel( result[2] ) );
+        CHECK_THAT( 2.5, WithinRel( result[3] ) );
+        CHECK_THAT( 1.5, WithinRel( result[4] ) );
+
+        result = integrator.firstMoment( table );
+
+        CHECK( 5 == result.size() );
+        CHECK_THAT( 0.0       , WithinRel( result[0] ) );
+        CHECK_THAT( 0.0       , WithinRel( result[1] ) );
+        CHECK_THAT( 5. + 1./6., WithinRel( result[2] ) );
+        CHECK_THAT( 6. + 1./6., WithinRel( result[3] ) );
+        CHECK_THAT( 5. + 1./6., WithinRel( result[4] ) );
       } // THEN
     } // WHEN
 
@@ -127,6 +171,24 @@ SCENARIO( "AnalyticalIntegrator" ) {
         CHECK_THAT( 1.5, WithinRel( result[2] ) );
         CHECK_THAT( 0.0, WithinRel( result[3] ) );
         CHECK_THAT( 0.0, WithinRel( result[4] ) );
+
+        result = integrator.zerothMoment( table );
+
+        CHECK( 5 == result.size() );
+        CHECK_THAT( 3.5, WithinRel( result[0] ) );
+        CHECK_THAT( 2.5, WithinRel( result[1] ) );
+        CHECK_THAT( 1.5, WithinRel( result[2] ) );
+        CHECK_THAT( 0.0, WithinRel( result[3] ) );
+        CHECK_THAT( 0.0, WithinRel( result[4] ) );
+
+        result = integrator.firstMoment( table );
+
+        CHECK( 5 == result.size() );
+        CHECK_THAT( 5. + 1./6., WithinRel( result[0] ) );
+        CHECK_THAT( 6. + 1./6., WithinRel( result[1] ) );
+        CHECK_THAT( 5. + 1./6., WithinRel( result[2] ) );
+        CHECK_THAT( 0.0       , WithinRel( result[3] ) );
+        CHECK_THAT( 0.0       , WithinRel( result[4] ) );
       } // THEN
     } // WHEN
 
@@ -148,6 +210,29 @@ SCENARIO( "AnalyticalIntegrator" ) {
         CHECK_THAT( 0.62, WithinRel( result[4] ) );
         CHECK_THAT( 2.5, WithinRel( result[5] ) );
         CHECK_THAT( 1.5, WithinRel( result[6] ) );
+
+        result = integrator.zerothMoment( table );
+
+        CHECK( 7 == result.size() );
+        CHECK_THAT( 0.78, WithinRel( result[0] ) );
+        CHECK_THAT( 0.74, WithinRel( result[1] ) );
+        CHECK_THAT( 0.7, WithinRel( result[2] ) );
+        CHECK_THAT( 0.66, WithinRel( result[3] ) );
+        CHECK_THAT( 0.62, WithinRel( result[4] ) );
+        CHECK_THAT( 2.5, WithinRel( result[5] ) );
+        CHECK_THAT( 1.5, WithinRel( result[6] ) );
+
+        result = integrator.firstMoment( table );
+
+        CHECK( 7 == result.size() );
+        CHECK_THAT( 5.144 / 6., WithinRel( result[0] ) );
+        CHECK_THAT( 5.768 / 6., WithinRel( result[1] ) );
+        CHECK_THAT( 6.296 / 6., WithinRel( result[2] ) );
+        CHECK_THAT( 6.728 / 6., WithinRel( result[3] ) );
+        CHECK_THAT( 7.064 / 6., WithinRel( result[4] ) );
+        CHECK_THAT( 6. + 1./6., WithinRel( result[5] ) );
+        CHECK_THAT( 5. + 1./6., WithinRel( result[6] ) );
+
       } // THEN
     } // WHEN
   } // GIVEN