Skip to content

Commit f6f2f3f

Browse files
author
Martin D. Weinberg
committed
Add inner/outer boundary restrictions for SL solution to prevent underflow while computing basis elements with large harmonic orders
1 parent f7ec660 commit f6f2f3f

2 files changed

Lines changed: 100 additions & 14 deletions

File tree

‎exputil/SLGridMP2.cc‎

Lines changed: 88 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -120,6 +120,13 @@ double sphdens(double r)
120120
return 4.0*M_PI * model->get_density(r);
121121
}
122122

123+
// Use entire grid for for l<Lswitch
124+
int SLGridSph::Lswitch = 32;
125+
126+
// For l>=Lswitch: rmin=rmap/rfac, rmax=rmap*rfac with rfac=pow(10,
127+
// Lalpha/l)
128+
double SLGridSph::Lalpha = 100.0;
129+
123130

124131
void SLGridSph::bomb(string oops)
125132
{
@@ -1109,10 +1116,25 @@ void SLGridSph::compute_table(struct TableSph* table, int l)
11091116
if (myid==0) VERBOSE = SLEDGE_VERBOSE;
11101117
#endif
11111118

1112-
cons[6] = rmin;
1113-
cons[7] = rmax;
1119+
// Inner radial boundary heuristic
1120+
//
1121+
int Nlo = 0, Nhi = numr;
1122+
if (l>Lswitch) {
1123+
double Rfac = std::pow(10.0, Lalpha/l);
1124+
double Rmin = std::max(rmin, rmap/Rfac);
1125+
double Rmax = std::min(rmax, rmap*Rfac);
1126+
1127+
// Find index boundaries in r grid
1128+
auto lower = std::lower_bound(r.data(), r.data() + r.size(), Rmin);
1129+
Nlo = std::distance(r.data(), lower);
1130+
auto upper = std::upper_bound(r.data(), r.data() + r.size(), Rmax);
1131+
Nhi = std::distance(r.data(), upper);
1132+
}
1133+
1134+
cons[6] = r[Nlo];
1135+
cons[7] = r[Nhi-1];
11141136
L2 = l*(l+1);
1115-
NUM = numr;
1137+
NUM = Nhi - Nlo;
11161138
N = nmax;
11171139

11181140
// integer iflag[nmax], invec[nmax+3];
@@ -1164,7 +1186,7 @@ void SLGridSph::compute_table(struct TableSph* table, int l)
11641186
//
11651187
// Output mesh
11661188
//
1167-
for (int i=0; i<NUM; i++) xef[i] = r[i];
1189+
for (int i=0; i<NUM; i++) xef[i] = r[i+Nlo];
11681190

11691191
//
11701192
// Open file for output.
@@ -1197,6 +1219,7 @@ void SLGridSph::compute_table(struct TableSph* table, int l)
11971219

11981220
std::cout << std::endl
11991221
<< "Tolerance errors in Sturm-Liouville solver for l=" << l
1222+
<< ", rmin=" << cons[6] << ", rmax=" << cons[7]
12001223
<< std::endl << std::endl;
12011224

12021225
std::cout << std::setw(15) << "order"
@@ -1216,6 +1239,32 @@ void SLGridSph::compute_table(struct TableSph* table, int l)
12161239
}
12171240
std::cout << std::endl;
12181241

1242+
// Additional diagnostic output for a failed l-value computation
1243+
//
1244+
if (VERBOSE) {
1245+
1246+
for (int i=0; i<N; i++) {
1247+
1248+
if (iflag[i] > -10) {
1249+
std::cout << "# order=" << i
1250+
<< " error: " << sledge_error(iflag[i])
1251+
<< std::endl << "#" << std::endl;
1252+
1253+
std::cout << std::setw(14) << "x"
1254+
<< std::setw(25) << "u(x)"
1255+
<< std::setw(25) << "(pu`)(x)"
1256+
<< std::endl;
1257+
int k = NUM*i;
1258+
for (int j=0; j<NUM; j++) {
1259+
std::cout << std::setw(25) << xef[j]
1260+
<< std::setw(25) << ef[j+k]
1261+
<< std::setw(25) << pdef[j+k]
1262+
<< std::endl;
1263+
}
1264+
}
1265+
}
1266+
}
1267+
12191268
sout << std::endl
12201269
<< "SLGridSph found " << bad
12211270
<< " tolerance errors in computing SL solutions." << std::endl
@@ -1273,6 +1322,7 @@ void SLGridSph::compute_table(struct TableSph* table, int l)
12731322
for (int i=0; i<N; i++) table->ev[i] = ev[i];
12741323

12751324
table->ef.resize(N, numr);
1325+
table->ef.setZero();
12761326

12771327
// Choose sign conventions for the ef table
12781328
//
@@ -1282,9 +1332,11 @@ void SLGridSph::compute_table(struct TableSph* table, int l)
12821332
if (ef[j*NUM+nfid]<0.0) sgn(j) = -1;
12831333
}
12841334

1285-
for (int i=0; i<numr; i++) {
1335+
// Pack offset eigen function values
1336+
//
1337+
for (int i=0; i<NUM; i++) {
12861338
for(int j=0; j<N; j++)
1287-
table->ef(j, i) = ef[j*NUM+i] * sgn(j);
1339+
table->ef(j, i+Nlo) = ef[j*NUM+i] * sgn(j);
12881340
}
12891341

12901342
table->l = l;
@@ -1374,11 +1426,26 @@ void SLGridSph::compute_table_worker(void)
13741426
if (tbdbg)
13751427
std::cout << "Worker " << mpi_myid << ": ordered to compute l = " << L << "" << std::endl;
13761428

1429+
// Inner radial boundary heuristic
1430+
//
1431+
int Nlo = 0, Nhi = numr;
1432+
if (L>Lswitch) {
1433+
double Rfac = std::pow(10.0, Lalpha/L);
1434+
double Rmin = std::max(rmin, rmap/Rfac);
1435+
double Rmax = std::min(rmax, rmap*Rfac);
1436+
1437+
// Find index boundaries in r grid
1438+
auto lower = std::lower_bound(r.data(), r.data() + r.size(), Rmin);
1439+
Nlo = std::distance(r.data(), lower);
1440+
auto upper = std::upper_bound(r.data(), r.data() + r.size(), Rmax);
1441+
Nhi = std::distance(r.data(), upper);
1442+
}
1443+
13771444
cons[0] = cons[1] = cons[2] = cons[3] = cons[4] = cons[5] = 0.0;
1378-
cons[6] = rmin;
1379-
cons[7] = rmax;
1445+
cons[6] = r[Nlo];
1446+
cons[7] = r[Nhi-1];
13801447
L2 = L*(L+1);
1381-
NUM = numr;
1448+
NUM = Nhi - Nlo;
13821449
N = nmax;
13831450

13841451
// integer iflag[nmax], invec[nmax+3];
@@ -1422,7 +1489,7 @@ void SLGridSph::compute_table_worker(void)
14221489
//
14231490
// Output mesh
14241491
//
1425-
for (int i=0; i<NUM; i++) xef[i] = r[i];
1492+
for (int i=0; i<NUM; i++) xef[i] = r[i+Nlo];
14261493

14271494
//
14281495
// Open file for output.
@@ -1454,6 +1521,7 @@ void SLGridSph::compute_table_worker(void)
14541521

14551522
std::cout << std::endl
14561523
<< "Tolerance errors in Sturm-Liouville solver for l=" << L
1524+
<< ", rmin=" << cons[6] << ", rmax=" << cons[7]
14571525
<< std::endl << std::endl;
14581526

14591527
std::cout << std::setw(15) << "order"
@@ -1523,10 +1591,16 @@ void SLGridSph::compute_table_worker(void)
15231591
if (ef[j*NUM+nfid]<0.0) sgn(j) = -1;
15241592
}
15251593

1594+
// Allocate table
1595+
//
15261596
table.ef.resize(N, numr);
1527-
for (int i=0; i<numr; i++) {
1528-
for (int j=0; j<N; j++)
1529-
table.ef(j, i) = ef[j*NUM+i] * sgn(j);
1597+
table.ef.setZero();
1598+
1599+
// Pack offset eigen function values
1600+
//
1601+
for (int i=0; i<NUM; i++) {
1602+
for(int j=0; j<N; j++)
1603+
table.ef(j, i+Nlo) = ef[j*NUM+i] * sgn(j);
15301604
}
15311605

15321606
table.l = L;
@@ -2833,7 +2907,7 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY)
28332907
<< std::setw( 5) << iflag[i]
28342908
<< std::endl;
28352909

2836-
if (VERBOSE) {
2910+
if (VERBOSE and iflag[i] != 0) {
28372911

28382912
if (iflag[i] > -10) {
28392913
cout << std::setw(14) << "x"

‎include/SLGridMP2.H‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,18 @@ private:
9494
//! Cache versioning
9595
inline static const std::string Version = "1.0";
9696

97+
//@{
98+
/** Rmin heuristic for grid construction to prevent underflow issues
99+
in Sturm-Liouville solver: sledge */
100+
101+
//! Use entire grid for l<Lswitch
102+
static int Lswitch;
103+
104+
//! For l>=Lswitch:
105+
//! rmin=rmap/rfac, rmax=rmap*rfac with rfac=pow(10, Lalpha/l)
106+
static double Lalpha;
107+
//@}
108+
97109
public:
98110

99111
//! Flag for MPI enabled (default: 0=off)

0 commit comments

Comments
 (0)