Skip to content

Commit ce347f5

Browse files
The9CatThe9Cat
authored
Merge pull request #143 from EXP-code/fieldFixes
Minor updates to field coefficient bases and OutVel
2 parents 7f9b6e4 + f145ea7 commit ce347f5

8 files changed

Lines changed: 538 additions & 113 deletions

File tree

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features
22

33
project(
44
EXP
5-
VERSION "7.8.4"
5+
VERSION "7.8.5"
66
HOMEPAGE_URL https://github.com/EXP-code/EXP
77
LANGUAGES C CXX Fortran)
88

‎expui/FieldBasis.cc‎

Lines changed: 43 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ namespace BasisClasses
3333
"ascl",
3434
"delta",
3535
"lmax",
36+
"mmax",
3637
"nmax",
3738
"model"
3839
};
@@ -54,7 +55,7 @@ namespace BasisClasses
5455

5556
// Allocate coefficient storage
5657
//
57-
nfld = p.size() + 2;
58+
nfld = p.size() + 1; // Add the density to the phase-space return
5859
allocateStore();
5960

6061
// Okay to register
@@ -67,7 +68,7 @@ namespace BasisClasses
6768
//
6869
void FieldBasis::configure()
6970
{
70-
nfld = 2; // Weight and density fields by default
71+
nfld = 1; // Density field by default
7172
lmax = 4;
7273
nmax = 10;
7374
rmin = 1.0e-4;
@@ -137,7 +138,8 @@ namespace BasisClasses
137138

138139
// Compute interpolation functionoid
139140
//
140-
double rmin = r.front(), rmax = r.back();
141+
rmin = r.front();
142+
rmax = r.back();
141143

142144
interp = std::make_shared<Linear1d>(r, d);
143145
densfunc = [this](double r)
@@ -163,12 +165,15 @@ namespace BasisClasses
163165

164166
// Generate the orthogonal function instance
165167
//
166-
ortho = std::make_shared<OrthoFunction>(nmax, densfunc, rmin, rmax, rmapping, dof);
168+
ortho = std::make_shared<OrthoFunction>
169+
(nmax-1, densfunc, rmin, rmax, rmapping, dof);
170+
// ^
171+
// |
172+
// +--- This is the polynomial order, not the rank
167173

168174
// Initialize fieldlabels
169175
//
170176
fieldLabels.clear();
171-
fieldLabels.push_back("weight");
172177
fieldLabels.push_back("density");
173178

174179
// Debug
@@ -204,9 +209,11 @@ namespace BasisClasses
204209

205210
void FieldBasis::initialize()
206211
{
207-
// Remove matched keys
212+
// Check for unmatched keys
208213
//
209-
// for (auto v : valid_keys) current_keys.erase(v);
214+
auto unmatched = YamlCheck(conf, valid_keys);
215+
if (unmatched.size())
216+
throw YamlConfigError("Basis::FieldBasis", "parameter", unmatched, __FILE__, __LINE__);
210217

211218
// Assign values from YAML
212219
//
@@ -215,6 +222,7 @@ namespace BasisClasses
215222
if (conf["model" ]) model = conf["model" ].as<std::string>();
216223
if (conf["nfld" ]) nfld = conf["nfld" ].as<int>();
217224
if (conf["lmax" ]) lmax = conf["lmax" ].as<int>();
225+
if (conf["mmax" ]) lmax = conf["mmax" ].as<int>();
218226
if (conf["nmax" ]) nmax = conf["nmax" ].as<int>();
219227
if (conf["dof" ]) dof = conf["dof" ].as<int>();
220228
if (conf["rmin" ]) rmin = conf["rmin" ].as<double>();
@@ -283,7 +291,7 @@ namespace BasisClasses
283291
// Sanity test dimensions
284292
if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) {
285293
std::ostringstream serr;
286-
serr << "FieldBasis::set_coefs: dimension error! "
294+
serr << "FieldBasis::set_coefs: dimension error for dof=2! "
287295
<< " nfld [" << nfld << "!= " << p->nfld << "]"
288296
<< " mmax [" << lmax << "!= " << p->mmax << "]"
289297
<< " nmax [" << nmax << "!= " << p->nmax << "]";
@@ -298,7 +306,7 @@ namespace BasisClasses
298306
// Sanity test dimensions
299307
if (nfld!=p->nfld || lmax!=p->lmax || nmax!=p->nmax) {
300308
std::ostringstream serr;
301-
serr << "FieldBasis::set_coefs: dimension error! "
309+
serr << "FieldBasis::set_coefs: dimension error for dof=3! "
302310
<< " nfld [" << nfld << "!= " << p->nfld << "]"
303311
<< " lmax [" << lmax << "!= " << p->lmax << "]"
304312
<< " nmax [" << nmax << "!= " << p->nmax << "]";
@@ -312,8 +320,7 @@ namespace BasisClasses
312320
double x, double y, double z,
313321
double u, double v, double w)
314322
{
315-
constexpr std::complex<double> I(0, 1);
316-
constexpr double fac0 = 0.25*M_2_SQRTPI;
323+
constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2*pi)=0.3989422804
317324

318325
int tid = omp_get_thread_num();
319326
PS3 pos{x, y, z}, vel{u, v, w};
@@ -337,45 +344,42 @@ namespace BasisClasses
337344

338345
auto p = (*ortho)(R);
339346

340-
(*coefs[tid])(0, 0, 0) += mass*p(0)*fac0;
341-
342347
for (int m=0; m<=lmax; m++) {
343348

344-
std::complex<double> P = std::exp(-I*(phi*m));
349+
auto P = std::exp(std::complex<double>(0, -phi*m))*fac2;
345350

346351
for (int n=0; n<nmax; n++) {
347352

348-
(*coefs[tid])(1, m, n) += mass*P*p(n);
353+
(*coefs[tid])(0, m, n) += mass*P*p(n);
349354

350355
for (int k=0; k<vec.size(); k++)
351-
(*coefs[tid])(k+2, m, n) += mass*P*p(n)*vec[k];
356+
(*coefs[tid])(k+1, m, n) += mass*P*p(n)*vec[k];
352357
}
353358
}
354359

355360
} else {
356361

357362
auto p = (*ortho)(r);
358363

359-
(*coefs[tid])(0, 0, 0) += mass*p(0);
360-
361364
for (int l=0, lm=0; l<=lmax; l++) {
362365

363366
double s = 1.0;
364367

365368
for (int m=0; m<=l; m++, lm++) {
366369

367370
// Spherical harmonic value
368-
std::complex<double> P =
369-
std::exp(-I*(phi*m))*Ylm(l, m)*plgndr(l, m, cth) * s;
370-
371-
s *= -1.0; // Flip sign for next evaluation
371+
auto P = std::exp(std::complex<double>(0, -phi*m)) *
372+
Ylm(l, m)*plgndr(l, m, cth) * s;
373+
374+
// Flip sign for next evaluation
375+
s *= -1.0;
372376

373377
for (int n=0; n<nmax; n++) {
374378

375-
(*coefs[tid])(1, lm, n) += mass*P*p(n);
379+
(*coefs[tid])(0, lm, n) += mass*P*p(n);
376380

377381
for (int k=0; k<vec.size(); k++)
378-
(*coefs[tid])(k+2, lm, n) += mass*P*p(n)*vec[k];
382+
(*coefs[tid])(k+1, lm, n) += mass*P*p(n)*vec[k];
379383
}
380384
}
381385
}
@@ -407,8 +411,10 @@ namespace BasisClasses
407411
//
408412
if (not coefret) {
409413
if (dof==2) coefret = std::make_shared<CoefClasses::CylFldStruct>();
410-
else coefret = std::make_shared<CoefClasses::SphFldStruct >();
411-
// Add the configuration data
414+
else coefret = std::make_shared<CoefClasses::SphFldStruct>();
415+
416+
// Add the configuration data
417+
//
412418
std::ostringstream sout; sout << node;
413419
coefret->buf = sout.str();
414420
}
@@ -443,18 +449,21 @@ namespace BasisClasses
443449
std::vector<double>
444450
FieldBasis::sph_eval(double r, double costh, double phi)
445451
{
446-
constexpr std::complex<double> I(0, 1);
452+
constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2 pi)
447453

448454
if (dof==2) {
449455
std::vector<double> ret(nfld, 0);
450-
auto p = (*ortho)(r*sqrt(fabs(1.0 - costh*costh)));
456+
auto p = (*ortho)(r);
451457
for (int i=0; i<nfld; i++) {
452458
for (int m=0; m<=lmax; m++) {
453-
double fac = m>0 ? 2.0 : 1.0;
459+
auto P = std::exp(std::complex<double>(0, -phi*m));
460+
461+
std::complex<double> sum = 0.0;
454462
for (int n=0; n<nmax; n++) {
455-
ret[i] += fac *
456-
std::real((*coefs[0])(i, m, n)*exp(I*(phi*m)))*p[n];
463+
sum += (*coefs[0])(i, m, n)*P * p(n) * fac2;
457464
}
465+
466+
ret[i] += std::real(sum);
458467
}
459468
}
460469
return ret;
@@ -478,9 +487,11 @@ namespace BasisClasses
478487
for (int i=0; i<nfld; i++) {
479488
for (int l=0, L=0; l<=lmax; l++) {
480489
for (int m=0; m<=l; m++, L++) {
490+
auto P = std::exp(std::complex<double>(0, -phi*m));
491+
481492
double fac = m>0 ? 2.0 : 1.0;
482493
for (int n=0; n<nmax; n++) {
483-
ret[i] += std::real((*coefs[0])(i, L, n)*exp(I*(phi*m)))*p[n]*
494+
ret[i] += std::real((*coefs[0])(i, L, n)*P)*p(n)*
484495
Ylm(l, m)*plgndr(l, m, costh) * fac;
485496
}
486497
}
@@ -751,7 +762,6 @@ namespace BasisClasses
751762

752763
// Field labels (force field labels added below)
753764
//
754-
fieldLabels.push_back("weight");
755765
fieldLabels.push_back("density");
756766

757767
if (coordinates == Coord::Cylindrical) {

‎exputil/OrthoFunction.cc‎

Lines changed: 19 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@
22
#include <OrthoFunction.H>
33

44
OrthoFunction::OrthoFunction
5-
(int nmax, DensFunc W, double rmin, double rmax, double scl, int dof, bool seg) :
6-
nmax(nmax), W(W), rmin(rmin), rmax(rmax), scale(scl), dof(dof), segment(seg)
5+
(int nmax, DensFunc W, double rmin, double rmax, double scl, int dof) :
6+
nmax(nmax), W(W), rmin(rmin), rmax(rmax), scale(scl), dof(dof)
77
{
88
// Quadrature knots/weights
99
lq = std::make_shared<LegeQuad>(knots);
@@ -27,14 +27,14 @@ void OrthoFunction::generate()
2727
norm.resize(nmax+1);
2828

2929
// Initial values
30-
beta[0] = norm[0] = scalar_prod(0, 0);
31-
alph[0] = scalar_prod(0, 1)/norm[0];
30+
beta(0) = norm(0) = scalar_prod(0, 0);
31+
alph(0) = scalar_prod(0, 1)/norm(0);
3232

3333
// Remaining values
3434
for (int i=1; i<=nmax; i++) {
35-
norm[i] = scalar_prod(i, 0);
36-
alph[i] = scalar_prod(i, 1)/norm[i];
37-
beta[i] = norm[i]/norm[i-1];
35+
norm(i) = scalar_prod(i, 0);
36+
alph(i) = scalar_prod(i, 1)/norm(i);
37+
beta(i) = norm(i)/norm(i-1);
3838
}
3939
}
4040

@@ -45,10 +45,10 @@ double OrthoFunction::scalar_prod(int n, int m)
4545
for (int i=0; i<knots; i++) {
4646
double x = xmin + dx*lq->knot(i);
4747
double r = x_to_r(x);
48-
double y = 2.0*r/scale; if (segment) y = x;
49-
Eigen::VectorXd p = poly_eval(y, n) * W(r);
48+
Eigen::VectorXd p = poly_eval(r, n) * W(r);
49+
5050
ans += dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1) *
51-
p(n) * p(n) * pow(y, m);
51+
p(n) * p(n) * pow(r, m);
5252
}
5353

5454
return ans;
@@ -58,11 +58,11 @@ double OrthoFunction::scalar_prod(int n, int m)
5858
Eigen::VectorXd OrthoFunction::poly_eval(double t, int n)
5959
{
6060
Eigen::VectorXd ret(n+1);
61-
ret[0] = 1.0;
61+
ret(0) = 1.0;
6262
if (n) {
63-
ret[1] = (t - alph[0])*ret[0];
63+
ret(1) = (t - alph(0))*ret(0);
6464
for (int j=1; j<n; j++) {
65-
ret[j+1] = (t - alph[j])*ret[j] - beta[j]*ret[j-1];
65+
ret(j+1) = (t - alph(j))*ret(j) - beta(j)*ret(j-1);
6666
}
6767
}
6868

@@ -81,15 +81,14 @@ Eigen::MatrixXd OrthoFunction::testOrtho()
8181
double x = xmin + dx*lq->knot(i);
8282
double r = x_to_r(x);
8383
double f = dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1);
84-
double y = 2.0*r/scale; if (segment) y = x;
8584

8685
// Evaluate the unnormalized polynomial
87-
Eigen::VectorXd p = poly_eval(y, nmax) * W(r);
86+
Eigen::VectorXd p = poly_eval(r, nmax) * W(r);
8887

8988
// Compute scalar products with the normalizations
9089
for (int n1=0; n1<=nmax; n1++) {
9190
for (int n2=0; n2<=nmax; n2++) {
92-
ret(n1, n2) += f * p(n1) * p(n2) / sqrt(norm[n1]*norm[n2]);
91+
ret(n1, n2) += f * p(n1) * p(n2) / sqrt(norm(n1)*norm(n2));
9392
}
9493
}
9594
}
@@ -111,13 +110,12 @@ void OrthoFunction::dumpOrtho(const std::string& filename)
111110
//
112111
double x = xmin + dx*i/(number-1);
113112
double r = x_to_r(x);
114-
double y = 2.0*r/scale; if (segment) y = x;
115113

116114
// Evaluate polynomial and apply the normalization
117115
//
118-
Eigen::VectorXd p = poly_eval(y, nmax) * W(r);
116+
Eigen::VectorXd p = poly_eval(r, nmax) * W(r);
119117
fout << std::setw(16) << r;
120-
for (int n=0; n<nmax; n++) fout << std::setw(16) << p(n)/sqrt(norm[n]);
118+
for (int n=0; n<=nmax; n++) fout << std::setw(16) << p(n)/sqrt(norm(n));
121119
fout << std::endl;
122120
}
123121
}
@@ -131,13 +129,9 @@ Eigen::VectorXd OrthoFunction::operator()(double r)
131129
// Weight
132130
double w = W(r);
133131

134-
// Choose system;
135-
double y = 2.0*r/scale;
136-
if (segment) y = r_to_x(r);
137-
138132
// Generate normalized orthogonal functions with weighting
139-
auto p = poly_eval(y, nmax);
140-
for (int i=0; i<=nmax; i++) p(i) *= w/sqrt(norm[i]);
133+
auto p = poly_eval(r, nmax);
134+
for (int i=0; i<=nmax; i++) p(i) *= w/sqrt(norm(i));
141135

142136
// The inner product of p(i), p(j) is Kronecker delta(i, j)
143137

‎exputil/mestel.cc‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -158,7 +158,7 @@ double TaperedMestelDisk::get_mass(double R)
158158
for (int i=1; i<=N; i++) {
159159
double rr = rmin + dr*i;
160160
double cur = get_density(rr);
161-
cum += 0.5*dr*(lst + cur) * 2.0*M_PI*rr;
161+
cum += 0.5*dr*(lst*(rr-dr) + cur*rr) * 2.0*M_PI;
162162
lst = cur;
163163

164164
vecR(i) = rr;

0 commit comments

Comments
 (0)