[WIP] Use Runge-Kutta-Fehlberg 7(8) integrator in SphericalDiffuse

CHANGELOG.md

@@ -20,6 +20,8 @@
 - The uppercased contents of the `.pcm` input file are written to a temporary
   file, instead of overwriting the user provided file. The temporary file is
   removed after it has been parsed. Fixes #91 as noted by @ilfreddy.
+- Use Runge-Kutta-Fehlberg 7(8) ODE solver to integrate the radial equation
+  in the spherical diffuse Green's function class.
 
 ## [Version 1.1.10] - 2017-03-27
 

Dangerfile

@@ -20,7 +20,9 @@ end
 commit_lint.check disable: [:subject_cap, :subject_period]
 
 # Check code style with clang-format
-code_style_validation.check file_extensions: ['.hpp', '.cpp', '.h'], ignore_file_patterns: [/^external\//]
+code_style_validation.check validator: 'clang-format-3.9',
+                            file_extensions: ['.hpp', '.cpp', '.h'],
+                            ignore_file_patterns: [/^external\//]
 
 # ------------------------------------------------------------------------------
 # What changed?

Gemfile

@@ -4,4 +4,4 @@ source "https://rubygems.org"
 gem "danger"
 gem "danger-commit_lint"
 gem "danger-lgtm"
-gem 'danger-code_style_validation', '0.1.0' , :git => 'https://github.com/flix-tech/danger-code_style_validation.git'
+gem 'danger-code_style_validation', :git => 'https://github.com/robertodr/danger-code_style_validation.git', :branch => 'pass-exe'

Gemfile.lock

@@ -1,6 +1,7 @@
 GIT
-  remote: https://github.com/flix-tech/danger-code_style_validation.git
-  revision: 05925775fa1ba9cd3eda09437c8c7efba31ef59e
+  remote: https://github.com/robertodr/danger-code_style_validation.git
+  revision: 49cb68348ff8e5440cdeb3512f85a2c602a9db83
+  branch: pass-exe
   specs:
     danger-code_style_validation (0.1.0)
       danger-plugin-api (~> 1.0)
@@ -59,7 +60,7 @@ PLATFORMS
 
 DEPENDENCIES
   danger
-  danger-code_style_validation (= 0.1.0)!
+  danger-code_style_validation!
   danger-commit_lint
   danger-lgtm
 

src/green/InterfacesImpl.hpp

@@ -225,15 +225,15 @@ template <typename StateVariable, typename ODESystem> class Zeta __final {
    */
   void compute(const ProfileEvaluator & eval, const IntegratorParameters & parms) {
     namespace odeint = boost::numeric::odeint;
-    odeint::bulirsch_stoer_dense_out<StateVariable> stepper(
-        parms.eps_abs_, parms.eps_rel_, parms.factor_x_, parms.factor_dxdt_);
+    odeint::runge_kutta_fehlberg78<StateVariable> stepper;
+
     ODESystem system(eval, L_);
     // Holds the initial conditions
     StateVariable init_zeta(2);
     // Set initial conditions
     init_zeta[0] = L_ * std::log(r_0_);
     init_zeta[1] = L_ / r_0_;
-    odeint::integrate_adaptive(
+    odeint::integrate_const(
         stepper,
         system,
         init_zeta,
@@ -329,16 +329,16 @@ template <typename StateVariable, typename ODESystem> class Omega __final {
    */
   void compute(const ProfileEvaluator & eval, const IntegratorParameters & parms) {
     namespace odeint = boost::numeric::odeint;
-    odeint::bulirsch_stoer_dense_out<StateVariable> stepper(
-        parms.eps_abs_, parms.eps_rel_, parms.factor_x_, parms.factor_dxdt_);
+    odeint::runge_kutta_fehlberg78<StateVariable> stepper;
+
     ODESystem system(eval, L_);
     // Holds the initial conditions
     StateVariable init_omega(2);
     // Set initial conditions
     init_omega[0] = -(L_ + 1) * std::log(r_infinity_);
     init_omega[1] = -(L_ + 1) / r_infinity_;
     // Notice that we integrate BACKWARDS, so we pass -step to integrate_adaptive
-    boost::numeric::odeint::integrate_adaptive(
+    odeint::integrate_const(
         stepper,
         system,
         init_omega,

src/green/SphericalDiffuse.cpp

@@ -39,9 +39,9 @@
 #include "GreenData.hpp"
 #include "GreensFunction.hpp"
 #include "cavity/Element.hpp"
-#include "utils/MathUtils.hpp"
 #include "dielectric_profile/ProfileTypes.hpp"
 #include "utils/ForId.hpp"
+#include "utils/MathUtils.hpp"
 
 namespace pcm {
 namespace green {
@@ -272,10 +272,10 @@ void SphericalDiffuse<ProfilePolicy>::initSphericalDiffuse() {
   double eps_rel_ = 1.0e-06; /*! Relative tolerance level */
   double factor_x_ = 0.0;    /*! Weight of the state      */
   double factor_dxdt_ = 0.0; /*! Weight of the state derivative */
-  double r_0_ = 0.5;         /*! Lower bound of the integration interval */
+  double r_0_ = 0.1;         /*! Lower bound of the integration interval */
   double r_infinity_ =
       this->profile_.center() + 200.0; /*! Upper bound of the integration interval */
-  double observer_step_ = 1.0e-03;     /*! Time step between observer calls */
+  double observer_step_ = 1.0e-02;     /*! Time step between observer calls */
   IntegratorParameters params_(eps_abs_,
                                eps_rel_,
                                factor_x_,

src/utils/MathUtils.hpp

@@ -292,14 +292,26 @@ inline double splineInterpolation(const double point,
                                   const std::vector<double> & grid,
                                   const std::vector<double> & function) {
   // Find nearest points on grid to the arbitrary point given
-  size_t index = std::distance(grid.begin(),
-                               std::lower_bound(grid.begin(), grid.end(), point)) -
-                 1;
+  int index =
+      std::distance(grid.begin(), std::lower_bound(grid.begin(), grid.end(), point));
+
+  int imax = grid.size() - 1;
+  if (index <= 0)
+    index = 1;
+  if (index >= imax - 1)
+    index = imax - 2;
 
   // Parameters for the interpolating spline
-  Eigen::Vector3d x, y;
-  x << grid[index - 1], grid[index], grid[index + 1];
-  y << function[index - 1], function[index], function[index + 1];
+  Eigen::Vector4d x = (Eigen::Vector4d() << grid[index - 1],
+                       grid[index],
+                       grid[index + 1],
+                       grid[index + 2])
+                          .finished();
+  Eigen::Vector4d y = (Eigen::Vector4d() << function[index - 1],
+                       function[index],
+                       function[index + 1],
+                       function[index + 2])
+                          .finished();
   SplineFunction s(x, y);
 
   return s(point);

tests/bi_operators/bi_operators_collocation.cpp

@@ -24,17 +24,17 @@
 #include "catch.hpp"
 
 #include <cmath>
-#include <iostream>
 #include <iomanip>
+#include <iostream>
 #include <limits>
 
 #include <Eigen/Core>
 
+#include "TestingMolecules.hpp"
+#include "bi_operators/Collocation.hpp"
 #include "cavity/Element.hpp"
 #include "cavity/GePolCavity.hpp"
-#include "bi_operators/Collocation.hpp"
 #include "green/SphericalDiffuse.hpp"
-#include "TestingMolecules.hpp"
 #include "green/UniformDielectric.hpp"
 #include "green/Vacuum.hpp"
 #include "utils/MathUtils.hpp"

tests/green/green_spherical_diffuse.cpp

@@ -43,7 +43,7 @@ using dielectric_profile::OneLayerTanh;
 SCENARIO("Evaluation of the spherical diffuse Green's function and its derivatives",
          "[green][green_spherical_diffuse]") {
   GIVEN("A permittivity profile modelled by the hyperbolic tangent function") {
-    int maxL = 3;
+    int maxL = 10;
     // High dielectric constant inside
     double eps1 = 80.0;
     // Low dielectric constant outside
@@ -79,7 +79,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
         REQUIRE(value == Approx(gf_value));
       }
       AND_THEN("the value of the Green's function outside the droplet is") {
-        double value = 0.50004567416494572;
+        double value = 0.49987476890456245;
         double gf_value = gf.kernelS(source2, probe2);
         INFO("ref_value = " << std::setprecision(
                                    std::numeric_limits<long double>::digits10)
@@ -103,7 +103,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
       }
       AND_THEN("the value of the Green's function directional derivative wrt the "
                "probe point outside the droplet is") {
-        double derProbe = -0.29005549287308696;
+        double derProbe = -0.289952813654903441;
         double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2);
         INFO("ref_derProbe = "
              << std::setprecision(std::numeric_limits<long double>::digits10)
@@ -127,7 +127,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
       }
       AND_THEN("the value of the Green's function directional derivative wrt the "
                "source point outside the droplet is") {
-        double derSource = 0.28879776627577236;
+        double derSource = 0.288676184370117994;
         double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2);
         INFO("ref_derSource = "
              << std::setprecision(std::numeric_limits<long double>::digits10)
@@ -155,7 +155,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
         REQUIRE(value == Approx(gf_value));
       }
       AND_THEN("the value of the Green's function outside the droplet is") {
-        double value = 0.5000329900631173;
+        double value = 0.499902162204405809;
         double gf_value = gf.kernelS(source2, probe2);
         INFO("ref_value = " << std::setprecision(
                                    std::numeric_limits<long double>::digits10)
@@ -179,7 +179,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
       }
       AND_THEN("the value of the Green's function directional derivative wrt the "
                "probe point outside the droplet is") {
-        double derProbe = -0.289999709125188243;
+        double derProbe = -0.289953098676354326;
         double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2);
         INFO("ref_derProbe = "
              << std::setprecision(std::numeric_limits<long double>::digits10)
@@ -203,7 +203,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
       }
       AND_THEN("the value of the Green's function directional derivative wrt the "
                "source point outside the droplet is") {
-        double derSource = 0.288657584958107449;
+        double derSource = 0.288675899050294671;
         double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2);
         INFO("ref_derSource = "
              << std::setprecision(std::numeric_limits<long double>::digits10)
@@ -254,7 +254,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
         REQUIRE(value == Approx(gf_value));
       }
       AND_THEN("the value of the Green's function outside the droplet is") {
-        double value = 0.49991229576650942;
+        double value = 0.499874677668854961;
         double gf_value = gf.kernelS(source2, probe2);
         INFO("ref_value = " << std::setprecision(
                                    std::numeric_limits<long double>::digits10)
@@ -278,7 +278,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
       }
       AND_THEN("the value of the Green's function directional derivative wrt the "
                "probe point outside the droplet is") {
-        double derProbe = -0.28997553534915177;
+        double derProbe = -0.289951551978862021;
         double gf_derProbe = gf.derivativeProbe(probeNormal2, source2, probe2);
         INFO("ref_derProbe = "
              << std::setprecision(std::numeric_limits<long double>::digits10)
@@ -302,7 +302,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
       }
       AND_THEN("the value of the Green's function directional derivative wrt the "
                "source point outside the droplet is") {
-        double derSource = 0.28871292628573908;
+        double derSource = 0.288674990983894819;
         double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2);
         INFO("ref_derSource = "
              << std::setprecision(std::numeric_limits<long double>::digits10)
@@ -379,7 +379,7 @@ SCENARIO("Evaluation of the spherical diffuse Green's function and its derivativ
       }
       AND_THEN("the value of the Green's function directional derivative wrt the "
                "source point outside the droplet is") {
-        double derSource = 0.28869402146192158;
+        double derSource = 0.288675321960529807;
         double gf_derSource = gf.derivativeSource(sourceNormal2, source2, probe2);
         INFO("ref_derSource = "
              << std::setprecision(std::numeric_limits<long double>::digits10)

tests/iefpcm/iefpcm_diffuse-gepol-point.cpp

@@ -28,12 +28,13 @@
 #include <Eigen/Core>
 
 #include "bi_operators/Collocation.hpp"
-#include "green/DerivativeTypes.hpp"
 #include "cavity/GePolCavity.hpp"
-#include "utils/Molecule.hpp"
+#include "green/DerivativeTypes.hpp"
+#include "green/SphericalDiffuse.hpp"
 #include "green/Vacuum.hpp"
 #include "solver/IEFSolver.hpp"
-#include "green/SphericalDiffuse.hpp"
+#include "utils/Molecule.hpp"
+
 #include "TestingMolecules.hpp"
 
 using namespace pcm;
@@ -74,7 +75,9 @@ SCENARIO("Test solver for the IEFPCM for a point charge in a spherical diffuse "
       double minRadius = 100.0;
       GePolCavity cavity(point, area, probeRadius, minRadius);
 
-      SphericalDiffuse<> gf_o(eps1, eps2, width, center, Eigen::Vector3d::Zero(), 3);
+      int maxL = 5;
+      SphericalDiffuse<> gf_o(
+          eps1, eps2, width, center, Eigen::Vector3d::Zero(), maxL);
       IEFSolver solver(symm);
       solver.buildSystemMatrix(cavity, gf_i, gf_o, op);
       int size = cavity.size();
@@ -105,7 +108,8 @@ SCENARIO("Test solver for the IEFPCM for a point charge in a spherical diffuse "
       double minRadius = 100.0;
       GePolCavity cavity(point, area, probeRadius, minRadius);
 
-      SphericalDiffuse<> gf_o(eps1, eps2, width, center, origin, 3);
+      int maxL = 5;
+      SphericalDiffuse<> gf_o(eps1, eps2, width, center, origin, maxL);
 
       IEFSolver solver(symm);
       solver.buildSystemMatrix(cavity, gf_i, gf_o, op);