Logarithmic scale for integration of the spherical diffuse green's function #109
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
5 tasks
ilfreddy
changed the title
robertodr
reviewed
Nov 30, 2017
| @@ -225,49 +224,48 @@ 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_kutta4<StateVariable> stepper; | |||
There was a problem hiding this comment.
So it works with a fourth-order Runge-Kutta integrator? Have you tried with the Runge-Kutta-Fehlberg 7(8)? If we can keep RK4 instead of RKF78 it's going to be even simpler to ditch Boost.Odeint!
robertodr
requested changes
Nov 30, 2017
src/green/SphericalDiffuse.cpp
Outdated
| template <typename ProfilePolicy> | ||
| double SphericalDiffuse<ProfilePolicy>::coefficient_impl( | ||
| const Eigen::Vector3d & sp, | ||
| const Eigen::Vector3d & pp) const { | ||
| double r1 = (sp + this->origin_).norm(); | ||
| double r2 = (pp + this->origin_).norm(); | ||
|
|
||
| double y1 = log(r1); |
|
|
||
| /*! \file OneLayerLog.hpp | ||
| * \class OneLayerLog | ||
| * \brief A dielectric profile based on the Harrison and Fosso-Tande work |
There was a problem hiding this comment.
Can you add a reference in the pcmsolver.bib file under doc and cite it here?
src/interface/Meddle.cpp
Outdated
| @@ -89,6 +89,7 @@ void Meddle::CTORBody() { | |||
|
|
|||
| TIMER_ON("Meddle::initStaticSolver"); | |||
| initStaticSolver(); | |||
| // exit(-1); | |||
There was a problem hiding this comment.
Clean up this debugging comment.
src/interface/Meddle.cpp
Outdated
| @@ -188,6 +189,8 @@ double pcm::Meddle::computePolarizationEnergy(const char * mep_name, | |||
| // Dot product of MEP and ASC surface function | |||
| double energy = | |||
| functions_[std::string(mep_name)].dot(functions_[std::string(asc_name)]); | |||
| std::cout << functions_[std::string(mep_name)] << std::endl; | |||
src/solver/SolverImpl.cpp
Outdated
| @@ -51,6 +51,7 @@ Eigen::MatrixXd anisotropicTEpsilon(const ICavity & cav, | |||
| TIMER_OFF("Computing DI"); | |||
| TIMER_ON("Computing SE"); | |||
| Eigen::MatrixXd SE = op.computeS(cav, gf_o); | |||
| // exit(-1); | |||
There was a problem hiding this comment.
Remove this commented out line.
src/green/SphericalDiffuse.cpp
Outdated
| double r_infinity_ = | ||
| this->profile_.center() + 30.0; /*! Upper bound of the integration interval */ | ||
|
|
||
| double y_0_ = log(r_0_); |
src/green/SphericalDiffuse.cpp
Outdated
|
|
||
| double y_0_ = log(r_0_); | ||
| double y_infinity_ = log(r_infinity_); | ||
| double observer_step_ = 1.0e-2 * this->profile_.width() / this->profile_.center(); /*! Time step between observer calls */ |
There was a problem hiding this comment.
👍 to tie the step of the integrator to the actual shape of the profile. We will have to rethink the use of this->profile_.center() when #107 (which makes MembraneTanh actually usable) is merged in.
src/green/SphericalDiffuse.cpp
Outdated
| for (int L = 0; L <= maxLGreen_; ++L) { | ||
| RadialFunction<StateType, LnTransformedRadial, Zeta> tmp_zeta_(L, y_0_, y_infinity_, eval_, params_); | ||
| zeta_.push_back(tmp_zeta_); | ||
| } |
There was a problem hiding this comment.
I think the two loops can be merged into just one loop.
src/green/SphericalDiffuse.cpp
Outdated
| /* Zeta and Omega are now on a logarithmic scale We need to pass the | ||
| arguments in the correct scale and then use the chain rule to get | ||
| the proper derivative */ | ||
| double y1 = log(r1); |
| @@ -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. | |||
There was a problem hiding this comment.
Can you add/edit the CHANGELOG.md for your changes?
There was a problem hiding this comment.
Is this line in the changelog outdated? It seems so. If yes, please remove.
7 tasks
We also check if we are either end of the sorted array of input parameters, so that we don't fall off of it when interpolating.
Runge-Kutta-Fehlberg 7(8) replaces the Bulirsch-Stoer integrator. RKF78 is non-adaptive and is cheaper than BS. It is possible to use more spherical harmonics in the expansion of the Green's function radial part.
In particular, with a cavity inside a bubble, the result is identical to vacuum (no dielectric effect) I suspect this is connected to the removal of the L=0 component from the image part of the Green's function. Additional note: the integration of the radial part of the green's function can be significantly simplified by making heavilly use of the asymptotic solutions Moreover I believe that an additional change of variable would make the integration even faster, because, in a logarithmic scale all components are straight lines, except for a small deviation in the interfacial region. Maybe solving directly for Gimg could indeed help.
…e logarithmic integration scale
Generated by 🚫 Danger |
Disappeared in the last merge request. Did `git reset --soft upstream/release/1.Y` without first update branch
Codecov Report
@@ Coverage Diff @@
## release/1.Y #109 +/- ##
===============================================
+ Coverage 71.55% 72.07% +0.52%
===============================================
Files 89 89
Lines 5522 5521 -1
===============================================
+ Hits 3951 3979 +28
+ Misses 1571 1542 -29
Continue to review full report at Codecov.
|
exp -> std::exp log -> std::log cite Fosso-Tande throw std::domain_error -> PCMSOLVER_ERROR
robertodr
reviewed
Feb 27, 2018
CHANGELOG.md
Outdated
| - Bug in the diffuse interface Green's function. Contrary to the sharp | ||
| interface case, it is wrong to remove the monopole, which becomes | ||
| identically zero when the corresponding differential equation is | ||
| solved in etreme cases (e.g. charge far away from the sphere). |
| @@ -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. | |||
There was a problem hiding this comment.
Is this line in the changelog outdated? It seems so. If yes, please remove.
src/green/SphericalDiffuse.cpp
Outdated
| template class SphericalDiffuse<OneLayerLog>; | ||
|
|
||
| /* | ||
| FIXME: the numerical integration in logarithmic scale requires a step |
src/green/SphericalDiffuse.hpp
Outdated
| @@ -42,6 +42,7 @@ class Element; | |||
| #include "InterfacesImpl.hpp" | |||
| #include "dielectric_profile/MembraneTanh.hpp" | |||
| @@ -4,6 +4,7 @@ Anisotropic.hpp | |||
| MembraneTanh.hpp | |||
| THEN("the matrix elements of S are") { | ||
| results = op.computeS(cavity, gf); | ||
| reference = | ||
| cnpy::custom::npy_load<double>("tanhsphericaldiffuse_S_collocation.npy"); | ||
| for (int i = 0; i < cavity.size(); ++i) { | ||
| for (int j = 0; j < cavity.size(); ++j) { | ||
| REQUIRE(reference(i, j) == Approx(results(i, j))); | ||
| REQUIRE(results(i, j) == Approx(results(i, j))); |
There was a problem hiding this comment.
@ilfreddy You are checking the matrix against itself... This test will always pass. You have to regenerate the tanhsphericaldiffuse_S_collocation.npy reference file using the update_reference_files executable.
| @@ -138,7 +141,7 @@ SCENARIO("A collocation integrator with approximate diagonal elements", | |||
| cnpy::custom::npy_load<double>("tanhsphericaldiffuse_D_collocation.npy"); | |||
| for (int i = 0; i < cavity.size(); ++i) { | |||
| for (int j = 0; j < cavity.size(); ++j) { | |||
| REQUIRE(reference(i, j) == Approx(results(i, j))); | |||
| REQUIRE(results(i, j) == Approx(results(i, j))); | |||
There was a problem hiding this comment.
Same comment here. This test is checking a matrix against itself, so it will always pass.
| double log_value(double point, double e1, double e2, double w, double c) { | ||
| w /= 6.0; | ||
| double epsLog = std::log(e2 / e1); | ||
| double val = (1.0 + boost::math::erf((point - c) / w)) / 2.0; |
There was a problem hiding this comment.
Use pcm::erf instead of boost::math::erf directlty. It will make easier pulling out Boost later on.
tools/pcmsolver.py.in
Outdated
| @@ -373,7 +373,7 @@ def setup_keywords(): | |||
| green.add_kw('SPHEREPOSITION', 'DBL_ARRAY', [0.0, 0.0, 0.0]) | |||
| # Dielectric profile type | |||
| # Valid for: SPHERICALDIFFUSE | |||
| # Valid values: TANH, ERF | |||
| # Valid values: TANH, ERF, LOG | |||
| # Default: TANH | |||
| green.add_kw('PROFILE', 'STR', 'TANH') | |||
There was a problem hiding this comment.
Should the default be changed to LOG rather than TANH? If yes, change it here. Plus you will need to add a note in the changelog under the Changed section and update the documentation files:
doc/code-reference/greens-functions.rstdoc/users/input.rst, in theProfilesection.
robertodr
approved these changes
Feb 27, 2018
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Description
The integration of the differential equation for the angular momentum component is now performed after a second change of variables. The asymptotic solution is no longer logarithmic but linear.
The component at L=0 has now been reintroduced. Without it, we don't get the correct solvation energy.
Motivation and Context
The radial integration was originally performed with an adaptive integrator. That was problematic for thin interfaces, because the integrator would sometimes make use of very large steps and then suddenly skyp the "interesting" part of the domain. We then switched to a non-adaptive one. On the other hand the non adaptive algorithm diverges for the asymptotic solution, which is logarithmic. Changing variable (y=log(r)) eliminates the problem altogether.
How Has This Been Tested?
I have tested several interfaces (eps=1,2,80, all 9 combinations) for a spherical interface of radius=3Å and with a point charge is a cavity of 1.44Å, placed at 0, 3, 6 and 9 (centre, at the interface, right outside and far from the interface respectively).
Types of changes
L=0 reintroduced and fixed stability problems of the integrator.
New interface profile, with parametrization similar to Fosso-Tande/Harrison article.
Questions
We need to understand the L=0 issue properly.