Skip to content

Commit f30746b

Browse files
robertodrrobertodr
committed
Calculate powers using a loop and not pow. Should help avoid overflows
1 parent 3e4073c commit f30746b

File tree

2 files changed

+182
-43
lines changed

2 files changed

+182
-43
lines changed

‎src/bin/plot_green_spherical.cpp‎

Lines changed: 160 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -7,24 +7,56 @@
77
#include <Eigen/Core>
88

99
#include "CollocationIntegrator.hpp"
10+
#include "DerivativeTypes.hpp"
1011
#include "GeneralPurpose.hpp"
1112
#include "LoggerInterface.hpp"
1213
#include "OneLayerTanh.hpp"
14+
#include "PhysicalConstants.hpp"
1315
#include "SphericalDiffuse.hpp"
1416
#include "UniformDielectric.hpp"
1517

18+
/*! \brief Generate uniform grid of points in interval
19+
* \param[in] min_val lower bound of interval
20+
* \param[in] max_val upper bound of interval
21+
* \param[in] nPoints number of points in interval
22+
*/
1623
std::vector<double> generate_grid(double, double, size_t);
24+
25+
/*! \brief Plot Green's function
26+
* The dielectric sphere is centered in (25.0, 25.0, 0.0) with radius 10.0
27+
* and width 5.0 The dielectric profile is parametrized by a tanh sigmoid
28+
* with epsInside = 78.39 (water) and epsOutside = 1.0 (vacuum)
29+
* The maximum value of angular momentum is set to 150
30+
* All distances in Angstrom.
31+
* The source charge is at (26.0, 25.0, 0.0) i.e. inside the dielectric sphere
32+
*/
1733
void case1();
34+
35+
/*! \brief Plot Green's function
36+
* The dielectric sphere is centered in (25.0, 25.0, 0.0) with radius 10.0
37+
* and width 5.0 The dielectric profile is parametrized by a tanh sigmoid
38+
* with epsInside = 78.39 (water) and epsOutside = 1.0 (vacuum)
39+
* The maximum value of angular momentum is set to 150
40+
* All distances in Angstrom.
41+
* The source charge is at (1.0, 1.0, 0.0) i.e. outside the dielectric sphere
42+
*/
1843
void case2();
44+
45+
/*! \brief Plot Green's function
46+
* The dielectric sphere is centered in (25.0, 25.0, 0.0) with radius 10.0
47+
* and width 5.0 The dielectric profile is parametrized by a tanh sigmoid
48+
* with epsInside = 78.39 (water) and epsOutside = 1.0 (vacuum)
49+
* The maximum value of angular momentum is set to 150
50+
* All distances in Angstrom.
51+
* The source charge is at (15.0, 25.0, 0.0) i.e. in the transition layer
52+
*/
1953
void case3();
20-
void case4();
2154

2255
int main()
2356
{
2457
case1();
25-
// case2();
26-
// case3();
27-
// case4();
58+
case2();
59+
case3();
2860

2961
return EXIT_SUCCESS;
3062
}
@@ -42,76 +74,167 @@ std::vector<double> generate_grid(double min_val, double max_val, size_t nPoints
4274
void case1()
4375
{
4476
size_t nPoints = 100;
45-
double epsInside = 80.0;
46-
double epsOutside = 2.0;
77+
double epsInside = 78.39;
78+
double epsOutside = 1.0;
79+
double angToBohr = angstromToBohr(2010);
4780
Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
48-
double sphereRadius = 10.0;
49-
double width = 5.0;
81+
sphereCenter << 25.0, 25.0, 0.0;
82+
sphereCenter *= angToBohr;
83+
double sphereRadius = 10.0 * angToBohr;
84+
double width = 5.0 * angToBohr;
5085
int maxL = 150;
5186

5287
Eigen::Vector3d source = Eigen::Vector3d::Zero();
88+
source << 26.0, 25.0, 0.0;
89+
source *= angToBohr;
5390
Eigen::Vector3d probe = Eigen::Vector3d::Zero();
5491

55-
double min = -50.0;
56-
double max = 50.0;
92+
// Conversion is done later, so that x, y values are printed in Angstrom
93+
double min = 0.0;
94+
double max = 40.0;
5795
std::vector<double> grid = generate_grid(min, max, nPoints);
5896
double delta = 0.1;
5997

6098
std::ofstream out;
6199

62100
SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
63101
LOG(gf);
102+
gf.toFile("CASE1");
103+
UniformDielectric<AD_directional, CollocationIntegrator> gf_i(epsInside);
104+
UniformDielectric<AD_directional, CollocationIntegrator> gf_o(epsOutside);
64105

65106
out.open("gf_spherical_CASE1.dat");
66107
out.precision(16);
67108
for (size_t i = 0; i < nPoints; ++i) {
68-
source << 1.0, 1.0, 0.0;
69109
for (size_t j = 0; j < nPoints; ++j) {
70-
probe << (grid[i] + delta), (grid[j] + delta), 0.0;
110+
probe << (grid[i] + delta), (grid[j] + delta), delta;
111+
probe *= angToBohr;
71112
out << '\t' << i
72-
<< '\t' << j
73-
<< '\t' << probe(0)
74-
<< '\t' << probe(1)
113+
<< '\t' << j
114+
<< '\t' << grid[i] + delta
115+
<< '\t' << grid[j] + delta
75116
<< '\t' << gf.Coulomb(source, probe)
76117
<< '\t' << gf.imagePotential(source, probe)
77118
<< '\t' << gf.kernelS(source, probe)
78-
<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe)) << std::endl;
119+
<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe))
120+
<< '\t' << gf_i.kernelS(source, probe)
121+
<< '\t' << gf_o.kernelS(source, probe)
122+
<< std::endl;
79123
}
124+
out << std::endl;
80125
}
81126
out.close();
82127
LOG_TIME;
128+
}
83129

84-
/*
85-
UniformDielectric<AD_directional, CollocationIntegrator> gf_inside(epsInside);
86-
out.open("gf_uniform_inside_CASE1.dat");
87-
out << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
88-
out << "gf_value" << '\t' << "derivativeProbe" << std::endl;
130+
void case2()
131+
{
132+
size_t nPoints = 100;
133+
double epsInside = 78.39;
134+
double epsOutside = 1.0;
135+
double angToBohr = angstromToBohr(2010);
136+
Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
137+
sphereCenter << 25.0, 25.0, 0.0;
138+
sphereCenter *= angToBohr;
139+
double sphereRadius = 10.0 * angToBohr;
140+
double width = 5.0 * angToBohr;
141+
int maxL = 150;
142+
143+
Eigen::Vector3d source = Eigen::Vector3d::Zero();
144+
source << 1.0, 1.0, 0.0;
145+
source *= angToBohr;
146+
Eigen::Vector3d probe = Eigen::Vector3d::Zero();
147+
148+
// Conversion is done later, so that x, y values are printed in Angstrom
149+
double min = 0.0;
150+
double max = 40.0;
151+
std::vector<double> grid = generate_grid(min, max, nPoints);
152+
double delta = 0.1;
153+
154+
std::ofstream out;
155+
156+
SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
157+
LOG(gf);
158+
gf.toFile("CASE2");
159+
UniformDielectric<AD_directional, CollocationIntegrator> gf_i(epsInside);
160+
UniformDielectric<AD_directional, CollocationIntegrator> gf_o(epsOutside);
161+
162+
out.open("gf_spherical_CASE2.dat");
163+
out.precision(16);
89164
for (size_t i = 0; i < nPoints; ++i) {
90165
for (size_t j = 0; j < nPoints; ++j) {
91-
probe << x[i], 0.0, z[j];
92-
out << '\t' << x[i]
93-
<< '\t' << z[j]
94-
<< '\t' << gf_inside.kernelS(source, probe)
95-
<< '\t' << gf_inside.derivativeProbe(Eigen::Vector3d::UnitZ(), source, probe) << std::endl;
166+
probe << (grid[i] + delta), (grid[j] + delta), delta;
167+
probe *= angToBohr;
168+
out << '\t' << i
169+
<< '\t' << j
170+
<< '\t' << grid[i] + delta
171+
<< '\t' << grid[j] + delta
172+
<< '\t' << gf.Coulomb(source, probe)
173+
<< '\t' << gf.imagePotential(source, probe)
174+
<< '\t' << gf.kernelS(source, probe)
175+
<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe))
176+
<< '\t' << gf_i.kernelS(source, probe)
177+
<< '\t' << gf_o.kernelS(source, probe)
178+
<< std::endl;
96179
}
180+
out << std::endl;
97181
}
98-
out.precision(16);
99182
out.close();
183+
LOG_TIME;
184+
}
185+
186+
void case3()
187+
{
188+
size_t nPoints = 100;
189+
double epsInside = 78.39;
190+
double epsOutside = 1.0;
191+
double angToBohr = angstromToBohr(2010);
192+
Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
193+
sphereCenter << 25.0, 25.0, 0.0;
194+
sphereCenter *= angToBohr;
195+
double sphereRadius = 10.0 * angToBohr;
196+
double width = 5.0 * angToBohr;
197+
int maxL = 150;
198+
199+
Eigen::Vector3d source = Eigen::Vector3d::Zero();
200+
source << 15.0, 25.0, 0.0;
201+
source *= angToBohr;
202+
Eigen::Vector3d probe = Eigen::Vector3d::Zero();
203+
204+
// Conversion is done later, so that x, y values are printed in Angstrom
205+
double min = 0.0;
206+
double max = 40.0;
207+
std::vector<double> grid = generate_grid(min, max, nPoints);
208+
double delta = 0.1;
209+
210+
std::ofstream out;
211+
212+
SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
213+
LOG(gf);
214+
gf.toFile("CASE3");
215+
UniformDielectric<AD_directional, CollocationIntegrator> gf_i(epsInside);
216+
UniformDielectric<AD_directional, CollocationIntegrator> gf_o(epsOutside);
100217

101-
UniformDielectric<AD_directional, CollocationIntegrator> gf_outside(epsOutside);
102-
out.open("gf_uniform_outside_CASE1.dat");
103-
out << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
104-
out << "gf_value" << '\t' << "derivativeProbe" << std::endl;
218+
out.open("gf_spherical_CASE3.dat");
105219
out.precision(16);
106220
for (size_t i = 0; i < nPoints; ++i) {
107221
for (size_t j = 0; j < nPoints; ++j) {
108-
probe << x[i], 0.0, z[j];
109-
out << '\t' << x[i]
110-
<< '\t' << z[j]
111-
<< '\t' << gf_outside.kernelS(source, probe)
112-
<< '\t' << gf_outside.derivativeProbe(Eigen::Vector3d::UnitZ(), source, probe) << std::endl;
222+
probe << (grid[i] + delta), (grid[j] + delta), delta;
223+
probe *= angToBohr;
224+
out << '\t' << i
225+
<< '\t' << j
226+
<< '\t' << grid[i] + delta
227+
<< '\t' << grid[j] + delta
228+
<< '\t' << gf.Coulomb(source, probe)
229+
<< '\t' << gf.imagePotential(source, probe)
230+
<< '\t' << gf.kernelS(source, probe)
231+
<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe))
232+
<< '\t' << gf_i.kernelS(source, probe)
233+
<< '\t' << gf_o.kernelS(source, probe)
234+
<< std::endl;
113235
}
236+
out << std::endl;
114237
}
115238
out.close();
116-
*/
239+
LOG_TIME;
117240
}

‎src/green/SphericalDiffuse.hpp‎

Lines changed: 22 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -207,6 +207,16 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
207207
std::tuple<double, double> epsilon(const Eigen::Vector3d & point) const {
208208
return this->profile_((point + this->origin_).norm());
209209
}
210+
void toFile(const std::string & prefix = "") {
211+
std::string tmp;
212+
prefix.empty() ? tmp = prefix : tmp = prefix + "-";
213+
writeToFile(zetaC_, tmp + "zetaC.dat");
214+
writeToFile(omegaC_, tmp + "omegaC.dat");
215+
for (int L = 1; L <= maxLGreen_; ++L) {
216+
writeToFile(zeta_[L], tmp + "zeta_" + std::to_string(L) + ".dat");
217+
writeToFile(omega_[L], tmp + "omega_" + std::to_string(L) + ".dat");
218+
}
219+
}
210220
EIGEN_MAKE_ALIGNED_OPERATOR_NEW /* See http://eigen.tuxfamily.org/dox/group__TopicStructHavingEigenMembers.html */
211221
private:
212222
using StateType = interfaces::StateType;
@@ -273,14 +283,12 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
273283
LOG("Computing first radial solution L = " + std::to_string(maxLC_));
274284
timerON("computeZeta for coefficient");
275285
zetaC_ = RadialFunction<StateType, LnTransformedRadial, Zeta>(maxLC_, r_0_, r_infinity_, eval_, params_);
276-
writeToFile(zetaC_, "zetaC.dat");
277286
timerOFF("computeZeta for coefficient");
278287
LOG("DONE: Computing first radial solution L = " + std::to_string(maxLC_));
279288

280289
LOG("Computing second radial solution L = " + std::to_string(maxLC_));
281290
timerON("computeOmega for coefficient");
282291
omegaC_ = RadialFunction<StateType, LnTransformedRadial, Omega>(maxLC_, r_0_, r_infinity_, eval_, params_);
283-
writeToFile(omegaC_, "omegaC.dat");
284292
timerOFF("computeOmega for coefficient");
285293
LOG("Computing second radial solution L = " + std::to_string(maxLC_));
286294
LOG("DONE: Computing coefficient for the separation of the Coulomb singularity");
@@ -368,10 +376,14 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
368376
double gr12 = 0.0;
369377
if (r1 < r2) {
370378
gr12 = std::exp(zeta1 - zeta2) * (2*L +1) / denominator;
371-
gr12 = (gr12 - std::pow(r1/r2, L) / (r2 * Cr12) ) * pl_x ;
379+
double f_L = r1 / r2;
380+
for (int i = 1; i < L; ++i) { f_L *= r1 / r2; }
381+
gr12 = (gr12 - f_L / (r2 * Cr12) ) * pl_x ;
372382
} else {
373383
gr12 = std::exp(omega1 - omega2) * (2*L +1) / denominator;
374-
gr12 = (gr12 - std::pow(r2/r1, L) / (r1 * Cr12) ) * pl_x ;
384+
double f_L = r2 / r1;
385+
for (int i = 1; i < L; ++i) { f_L *= r2 / r1; }
386+
gr12 = (gr12 - f_L / (r1 * Cr12) ) * pl_x ;
375387
}
376388

377389
return gr12;
@@ -421,11 +433,15 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
421433
double denominator = (d_zeta2 - d_omega2) * std::pow(r2, 2) * eps_r2;
422434

423435
if (r1 < r2) {
436+
double f_L = r1 / r2;
437+
for (int i = 1; i < maxLC_; ++i) { f_L *= r1 / r2; }
424438
tmp = std::exp(zeta1 - zeta2) * (2*maxLC_ +1) / denominator;
425-
coeff = std::pow(r1/r2, maxLC_) / (tmp * r2);
439+
coeff = f_L / (tmp * r2);
426440
} else {
441+
double f_L = r2 / r1;
442+
for (int i = 1; i < maxLC_; ++i) { f_L *= r2 / r1; }
427443
tmp = std::exp(omega1 - omega2) * (2*maxLC_ +1) / denominator;
428-
coeff = std::pow(r2/r1, maxLC_) / (tmp * r1);
444+
coeff = f_L / (tmp * r1);
429445
}
430446

431447
return coeff;

0 commit comments

Comments
 (0)