Skip to content

Commit 74bbe79

Browse files
author
Roberto Di Remigio
committed
Change integrator: from Bulirsch-Stoer to Runge-Kutta-Fehlberg(7,8)
1 parent 0f0462b commit 74bbe79

File tree

8 files changed

+421
-59
lines changed

8 files changed

+421
-59
lines changed

‎CMakeLists.txt‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -223,6 +223,3 @@ execute_process(
223223
include(BuildInfo)
224224
# Printout useful info
225225
include(ConfigInfo)
226-
227-
add_executable(plot_green_spherical.x ${PROJECT_SOURCE_DIR}/src/plot_green_spherical.cpp)
228-
target_link_libraries(plot_green_spherical.x green bi_operators utils ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY} ${Boost_TIMER_LIBRARY} ${Boost_CHRONO_LIBRARY})

‎src/bin/CMakeLists.txt‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,3 +3,9 @@
33

44
#add_executable(debug_wavcav.x debug_wavcav.cpp)
55
#target_link_libraries(debug_wavcav.x solver green cavity wavcav pwl wem utils)
6+
7+
add_executable(plot_green_spherical.x plot_green_spherical.cpp)
8+
target_link_libraries(plot_green_spherical.x green bi_operators utils ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY} ${Boost_TIMER_LIBRARY} ${Boost_CHRONO_LIBRARY})
9+
10+
add_executable(check_Coulomb_coefficient.x check_Coulomb_coefficient.cpp)
11+
target_link_libraries(check_Coulomb_coefficient.x green bi_operators utils ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY} ${Boost_TIMER_LIBRARY} ${Boost_CHRONO_LIBRARY})

‎src/bin/check_Coulomb_coefficient.cpp‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
#include <iostream>
2+
#include <ostream>
3+
#include <fstream>
4+
5+
#include <Eigen/Dense>
6+
7+
#include "DerivativeTypes.hpp"
8+
#include "UniformDielectric.hpp"
9+
#include "TanhSphericalDiffuse.hpp"
10+
11+
int main()
12+
{
13+
int nPoints = 10000;
14+
double epsInside = 80.0;
15+
double epsOutside = 1.0;
16+
Eigen::Vector3d sphereCenter;
17+
sphereCenter << 0.0, 0.0, 0.0;
18+
double sphereRadius = 100.0;
19+
double width = 10.0;
20+
21+
Eigen::Vector3d source;
22+
source << 0.0, 0.0, 0.0;
23+
Eigen::Vector3d probe;
24+
double delta = 0.05;
25+
26+
double zMin = 90.0;
27+
double zMax = 200.0;
28+
double step = (zMax - zMin) / nPoints;
29+
30+
TanhSphericalDiffuse gf(epsInside, epsOutside, width, sphereRadius);
31+
// Plot Green's function, image potential and separation coefficient
32+
std::ofstream out;
33+
out.open("checkCoulomb_spherical.log");
34+
out << "#" << '\t' << "Distance" << '\t' << "coefficient" << std::endl;
35+
out.precision(16);
36+
for (int i = 0; i < nPoints; ++i) {
37+
source(2) = zMin + i*step;
38+
probe = source;
39+
probe(2) +=delta;
40+
out << '\t' << source(2) << '\t' << gf.Coulomb(source, probe) << std::endl;
41+
}
42+
out.close();
43+
44+
UniformDielectric<AD_directional> gf_inside(epsInside);
45+
out.open("checkCoulomb_uniform_inside.log");
46+
out << "#" << '\t' << "Distance" << '\t' << "gf_value" << std::endl;
47+
out.precision(16);
48+
for (int i = 0; i < nPoints; ++i) {
49+
source(2) = zMin + i*step;
50+
probe = source;
51+
probe(2) +=delta;
52+
out << '\t' << source(2) << '\t' << gf_inside.function(source, probe) << std::endl;
53+
}
54+
out.close();
55+
56+
UniformDielectric<AD_directional> gf_outside(epsOutside);
57+
out.open("checkCoulomb_uniform_outside.log");
58+
out << "#" << '\t' << "Distance" << '\t' << "gf_value" << std::endl;
59+
out.precision(16);
60+
for (int i = 0; i < nPoints; ++i) {
61+
source(2) = zMin + i*step;
62+
probe = source;
63+
probe(2) +=delta;
64+
out << '\t' << source(2) << '\t' << gf_outside.function(source, probe) << std::endl;
65+
}
66+
out.close();
67+
}

‎src/bin/plot_green_spherical.cpp‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
#include <iostream>
2+
#include <ostream>
3+
#include <fstream>
4+
5+
#include <Eigen/Dense>
6+
7+
#include "DerivativeTypes.hpp"
8+
#include "UniformDielectric.hpp"
9+
#include "TanhSphericalDiffuse.hpp"
10+
11+
int main()
12+
{
13+
int nPoints = 10000;
14+
double epsInside = 80.0;
15+
double epsOutside = 1.0;
16+
Eigen::Vector3d sphereCenter;
17+
sphereCenter << 0.0, 0.0, 0.0;
18+
double sphereRadius = 100.0;
19+
double width = 10.0;
20+
21+
Eigen::Vector3d source;
22+
source << 0.5, 0.0, 0.0;
23+
Eigen::Vector3d probe;
24+
probe << 50.0, 0.0, 0.0;
25+
26+
double zMin = 0.0;
27+
double zMax = 200.0;
28+
double step = (zMax - zMin) / nPoints;
29+
30+
TanhSphericalDiffuse gf(epsInside, epsOutside, width, sphereRadius);
31+
// Plot Green's function, image potential and separation coefficient
32+
std::ofstream out;
33+
out.open("gf_spherical.log");
34+
out << "#" << '\t' << "Distance" << '\t' << "gf_value" << '\t' << "image" << '\t' << "coefficient" << std::endl;
35+
out.precision(16);
36+
for (int i = 0; i < nPoints; ++i) {
37+
probe(2) = zMin + i*step;
38+
out << '\t' << probe(2) << '\t' << gf.function(source, probe) << '\t' << gf.imagePotential(source, probe) << '\t' << gf.Coulomb(source, probe) << std::endl;
39+
}
40+
out.close();
41+
42+
UniformDielectric<AD_directional> gf_inside(epsInside);
43+
out.open("gf_uniform_inside.log");
44+
out << "#" << '\t' << "Distance" << '\t' << "gf_value" << std::endl;
45+
out.precision(16);
46+
for (int i = 0; i < nPoints; ++i) {
47+
probe(2) = zMin + i*step;
48+
out << '\t' << probe(2) << '\t' << gf_inside.function(source, probe) << std::endl;
49+
}
50+
out.close();
51+
52+
UniformDielectric<AD_directional> gf_outside(epsOutside);
53+
out.open("gf_uniform_outside.log");
54+
out << "#" << '\t' << "Distance" << '\t' << "gf_value" << std::endl;
55+
out.precision(16);
56+
for (int i = 0; i < nPoints; ++i) {
57+
probe(2) = zMin + i*step;
58+
out << '\t' << probe(2) << '\t' << gf_outside.function(source, probe) << std::endl;
59+
}
60+
out.close();
61+
}

‎src/green/GreensFunction.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -110,7 +110,7 @@ class GreensFunction: public IGreensFunction
110110
* Notice that this method returns the directional derivative with respect
111111
* to the probe point.
112112
*
113-
* \param[in] normal_p2 the normal vector to p1
113+
* \param[in] normal_p2 the normal vector to p2
114114
* \param[in] p1 first point
115115
* \param[in] p2 second point
116116
*/

‎src/green/SphericalDiffuse.hpp‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -181,6 +181,8 @@ class SphericalDiffuse : public GreensFunction<Numerical, ProfilePolicy>
181181
friend std::ostream & operator<<(std::ostream & os, SphericalDiffuse & gf) {
182182
return gf.printObject(os);
183183
}
184+
double Coulomb(const Eigen::Vector3d & source, const Eigen::Vector3d & probe) const;
185+
double imagePotential(const Eigen::Vector3d & source, const Eigen::Vector3d & probe) const;
184186
private:
185187
/*! Evaluates the Green's function given a pair of points
186188
*
@@ -222,7 +224,7 @@ class SphericalDiffuse : public GreensFunction<Numerical, ProfilePolicy>
222224
/*! This calculates all the components needed to evaluate the Green's function */
223225
void initSphericalDiffuse();
224226

225-
/**@{ Parameters and functions for the calculation of the Green's function, including Coulomb singularity */
227+
/**@{ Parameters for the calculation of the Green's function, including Coulomb singularity */
226228
/*! Maximum angular momentum in the final summation over Legendre polynomials to obtain G */
227229
int maxLGreen_ = 30;
228230
/*! \brief First independent radial solution, used to build Green's function.
@@ -235,7 +237,7 @@ class SphericalDiffuse : public GreensFunction<Numerical, ProfilePolicy>
235237
std::vector<RadialFunction> omega_;
236238
/**@}*/
237239

238-
/**@{ Parameters and functions for the calculation of the Coulomb singularity separation coefficient */
240+
/**@{ Parameters for the calculation of the Coulomb singularity separation coefficient */
239241
/*! Maximum angular momentum to obtain C(r, r'), needed to separate the Coulomb singularity */
240242
int maxLC_ = maxLGreen_ + 30;
241243
/*! \brief First independent radial solution, used to build coefficient.
@@ -250,6 +252,9 @@ class SphericalDiffuse : public GreensFunction<Numerical, ProfilePolicy>
250252

251253
double coefficient(double r1, double r2) const;
252254
double functionSummation(int L, double r1, double r2, double cos_gamma, double Cr12) const;
255+
256+
Eigen::Vector3d coefficientGradient(const Eigen::Vector3d & p1, const Eigen::Vector3d & p2) const;
257+
Eigen::Vector3d functionSummationGradient(int L, const Eigen::Vector3d & p1, const Eigen::Vector3d & p2, double Cr12) const;
253258
};
254259

255260
#endif // SPHERICALDIFFUSE_HPP

0 commit comments

Comments
 (0)