Skip to content

Updates to IC generation routines in utils/ICs #107

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
The9Cat wants to merge 13 commits into from
Closed

Updates to IC generation routines in utils/ICs #107

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
cf49265
Added ring generation routine, position quick start to gensph, and mo…
Jan 17, 2025
2ee5f3b
Use additive machine constants to set epsilon in Legendre functions o…
Jan 17, 2025
51298ca
Update particle numbers to exactly fill each radial shell
Jan 17, 2025
37cbf88
Stagger radial bins from mass end points; add diagnostic mass ring pr…
Jan 17, 2025
7c20e98
Merge remote-tracking branch 'origin/FixEPS' into ICupdate
Jan 17, 2025
7285737
Added Fibonacci spiral mapped to the sphere for uniform coverage of s…
Jan 18, 2025
846e4f4
Fixed error in velocity assignment at fixed radius
Jan 18, 2025
87666df
Updates to gensph for quiet start variations
Jan 21, 2025
8539dc5
Make goldenRatio a const rather than a constexpr
Jan 21, 2025
4fe5c44
Print mismatched cache variables
Jan 24, 2025
95d9a99
Add quiet-start method
Jan 24, 2025
ee3fa9a
Report cache name tagged for recomputation
Feb 4, 2025
72fc23b
Removed the angle-only quick start for the sphere, which leaves too m…
Feb 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 16 additions & 26 deletions exputil/EXPmath.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,12 +61,12 @@ namespace AltMath
return ans;
}

#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10

double bessj(int n,double x)
{
const double ACC = 40.0;
const double BIGNO = 1.0e10;
const double BIGNI = 1.0e-10;

int j,jsum,m;
double ax,bj,bjm,bjp,sum,tox,ans;

Expand Down Expand Up @@ -114,10 +114,6 @@ namespace AltMath
return x < 0.0 && n%2 == 1 ? -ans : ans;
}

#undef ACC
#undef BIGNO
#undef BIGNI

double bessk0(double x)
{
double y,ans;
Expand Down Expand Up @@ -174,12 +170,12 @@ namespace AltMath
return ans;
}

#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10

double bessi(int n,double x)
{
const double ACC = 40.0;
const double BIGNO = 1.0e10;
const double BIGNI = 1.0e-10;

int j;
double bi,bim,bip,tox,ans;
double bessi0(double );
Expand Down Expand Up @@ -207,10 +203,6 @@ namespace AltMath
}
}

#undef ACC
#undef BIGNO
#undef BIGNI

double bessi0(double x)
{
double ax,ans;
Expand Down Expand Up @@ -253,12 +245,12 @@ namespace AltMath
return x < 0.0 ? -ans : ans;
}

#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10

double jn_sph(int n, double x)
{
const double ACC = 40.0;
const double BIGNO = 1.0e10;
const double BIGNI = 1.0e-10;

int j,m;
double ax,bj,bjm,bjp,tox,ans;
double jn_sph0(double x),jn_sph1(double x),jn_sphm1(double x);
Expand Down Expand Up @@ -306,13 +298,10 @@ namespace AltMath
return x < 0.0 && n%2 == 1 ? -ans : ans;
}

#undef ACC
#undef BIGNO
#undef BIGNI

#define TOL 1.0e-7
double jn_sph0(double x)
{
const double TOL = 1.0e-7;

if (fabs(x)<TOL)
return 1.0-x*x/6.0;
else
Expand All @@ -321,6 +310,8 @@ namespace AltMath

double jn_sph1(double x)
{
const double TOL = 1.0e-7;

if (fabs(x)<TOL)
return x/3.0;
else
Expand All @@ -331,7 +322,6 @@ namespace AltMath
{
return cos(x)/x;
}
#undef TOL

double check_nu(double nu)
{
Expand Down
22 changes: 17 additions & 5 deletions exputil/EmpCyl2d.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -780,19 +780,31 @@ bool EmpCyl2d::ReadH5Cache()
auto checkInt = [&file](int value, std::string name)
{
int v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
if (value == v) return true; return false;
if (value == v) return true;
if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
<< name << " expected " << value << " but found "
<< v << std::endl;
return false;
};

auto checkDbl = [&file](double value, std::string name)
{
double v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
if (fabs(value - v) < 1.0e-16) return true; return false;
if (fabs(value - v) < 1.0e-16) return true;
if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
<< name << " expected " << value << " but found "
<< v << std::endl;
return false;
};

auto checkStr = [&file](std::string value, std::string name)
{
std::string v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
if (value.compare(v)==0) return true; return false;
if (value.compare(v)==0) return true;
if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
<< name << " expected " << value << " but found "
<< v << std::endl;
return false;
};

//
Expand All @@ -803,8 +815,8 @@ bool EmpCyl2d::ReadH5Cache()
if (not checkStr(Version, "Version")) return false;
} else {
if (myid==0)
std::cout << "---- EmpCyl2d::ReadH5Cache: "
<< "recomputing cache for HighFive API change"
std::cout << "---- EmpCyl2d::ReadH5Cache: " << "<" << cache_name_2d
<< "> recomputing cache for HighFive API change"
<< std::endl;
return false;
}
Expand Down
40 changes: 29 additions & 11 deletions exputil/massmodel_dist.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,14 +49,14 @@
#include <massmodel.H>
#include <interp.H>

#define OFFSET 1.0e-3
#define OFFTOL 1.2
const double OFFSET=1.0e-3;
const double OFFTOL=1.2;

extern double gint_0(double a, double b, std::function<double(double)> f, int NGauss);
extern double gint_2(double a, double b, std::function<double(double)> f, int NGauss);

#define TSTEP 1.0e-8
#define NGauss 96
const double TSTEP=1.0e-8;
const int NGauss=96;

static int DIVERGE=0;

Expand Down Expand Up @@ -109,6 +109,10 @@ void SphericalModelTable::setup_df(int NUM, double RA)
rhoQy.resize(num);
rhoQy2.resize(num);

rhoQx. setZero();
rhoQy. setZero();
rhoQy2.setZero();

for (int i=0; i<num; i++) {
x = density.x[i];
rhoQx[i] = pot.y[i];
Expand Down Expand Up @@ -138,14 +142,20 @@ void SphericalModelTable::setup_df(int NUM, double RA)
dfc.Q .resize(NUM);
dfc.fQ .resize(NUM);
dfc.ffQ.resize(NUM);

dfc.Q .setZero();
dfc.fQ .setZero();
dfc.ffQ.setZero();

dfc.num = NUM;
dfc.ra2 = ra2;

Qmax = get_pot(pot.x[pot.num-1]);
Qmin = get_pot(pot.x[0]);
dQ = (Qmax - Qmin)/(double)(dfc.num-1);

foffset = -std::numeric_limits<double>::max();
// foffset = -std::numeric_limits<double>::max();
foffset = -1.0e42;
dfc.Q[dfc.num-1] = Qmax;
dfc.ffQ[dfc.num-1] = 0.0;
fac = 1.0/(sqrt(8.0)*M_PI*M_PI);
Expand Down Expand Up @@ -184,6 +194,13 @@ void SphericalModelTable::setup_df(int NUM, double RA)
df.ffQ .resize(NUM);
df.fQ2 .resize(NUM);
df.ffQ2.resize(NUM);

df.Q .setZero();
df.fQ .setZero();
df.ffQ .setZero();
df.fQ2 .setZero();
df.ffQ2.setZero();

df.num = NUM;
df.ra2 = ra2;

Expand All @@ -194,7 +211,8 @@ void SphericalModelTable::setup_df(int NUM, double RA)
df.Q[df.num-1] = Qmax;
df.ffQ[df.num-1] = 0.0;
fac = 1.0/(sqrt(8.0)*M_PI*M_PI);
foffset = -std::numeric_limits<double>::max();
// foffset = -std::numeric_limits<double>::max();
foffset = -1.0e42;
for (int i=df.num-2; i>=0; i--) {
df.Q[i] = df.Q[i+1] - dQ;
Q = df.Q[i];
Expand Down Expand Up @@ -233,7 +251,7 @@ void SphericalModelTable::setup_df(int NUM, double RA)

dist_defined = true;

debug_fdist();
// debug_fdist();
}


Expand Down Expand Up @@ -294,7 +312,7 @@ double SphericalModelTable::distf(double E, double L)

if (!dist_defined) bomb("distribution function not defined");

double d, g;
double d=0.0, g=0.0;

if (chebyN) {

Expand Down Expand Up @@ -341,7 +359,7 @@ double SphericalModelTable::dfde(double E, double L)

if (!dist_defined) bomb("distribution function not defined");

double d, g, h, d1;
double d=0, g=0, h=0, d1=0;

if (chebyN) {

Expand Down Expand Up @@ -399,7 +417,7 @@ double SphericalModelTable::dfdl(double E, double L)

if (!dist_defined) bomb("distribution function not defined");

double d, g, h, d1;
double d=0, g=0, h=0, d1=0;

if (chebyN) {

Expand Down Expand Up @@ -451,7 +469,7 @@ double SphericalModelTable::d2fde2(double E, double L)
{
if (!dist_defined) bomb("distribution function not defined");

double d, g, h, k, d2;
double d=0, g=0, h=0, k=0, d2=0;

if (chebyN) {

Expand Down
Loading