Skip to content

Commit

Permalink
Merge pull request #61 from EXP-code/TaperedMestel
Browse files Browse the repository at this point in the history 
Implement a biorthongal basis target for the inner/outer tapered Mestel disk
  • Loading branch information
@The9Cat
The9Cat authored Mar 1, 2024
2 parents ebab458 + 9daee16 commit 071fd19
Show file tree
Hide file tree
Showing 40 changed files with 1,980 additions and 660 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.21) # Needed for CUDA, MPI, and CTest features

project(
EXP
VERSION "7.7.27"
VERSION "7.7.28"
HOMEPAGE_URL https://github.com/EXP-code/EXP
LANGUAGES C CXX Fortran)

Expand Down
21 changes: 15 additions & 6 deletions coefs/BiorthBasis.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -594,7 +594,7 @@ namespace BasisClasses
}

// Radial grid spacing
double dx = (logxmax - logxmin)/numgrid;
double dx = (logxmax - logxmin)/(numgrid-1);

// Basis storage
Eigen::MatrixXd tabpot, tabden, tabfrc;
Expand Down Expand Up @@ -1394,6 +1394,7 @@ namespace BasisClasses
"numx",
"numy",
"numr",
"NQDHT",
"knots",
"logr",
"model",
Expand All @@ -1405,14 +1406,18 @@ namespace BasisClasses
"NO_M0",
"NO_M1",
"EVEN_M",
"M0_BACK",
"M0_ONLY",
"NO_MONO",
"diskconf",
"ssfrac",
"playback",
"coefMaster",
"Lmax",
"Mmax",
"nmax",
"mmax",
"mlim",
"dof",
"subsamp",
"samplesz",
Expand Down Expand Up @@ -1456,6 +1461,8 @@ namespace BasisClasses
//
try {
if (conf["cmap"]) cmap = conf["cmap"].as<int>();
if (conf["Lmax"]) mmax = conf["Lmax"].as<int>(); // Proxy
if (conf["Mmax"]) mmax = conf["Mmax"].as<int>(); // Proxy
if (conf["mmax"]) mmax = conf["mmax"].as<int>();
if (conf["nmax"]) nmax = conf["nmax"].as<int>();

Expand Down Expand Up @@ -1722,18 +1729,20 @@ namespace BasisClasses
}

// Get the basis fields
//
ortho->get_dens (dend[tid], R, z);
ortho->get_pot (potd[tid], R, z);
ortho->get_rforce (potR[tid], R, z);
ortho->get_zforce (potZ[tid], R, z);

// m loop
//
for (int m=0, moffset=0; m<=mmax; m++) {

if (m==0 and NO_M0) continue;
if (m==1 and NO_M1) continue;
if (EVEN_M and m/2*2 != m) continue;
if (m>0 and M0_only) break;
if (m==0 and NO_M0) { moffset++; continue; }
if (m==1 and NO_M1) { moffset += 2; continue; }
if (EVEN_M and m/2*2 != m) { moffset += 2; continue; }
if (m>0 and M0_only) break;

if (m==0) {
for (int n=std::max<int>(0, N1); n<=std::min<int>(nmax-1, N2); n++) {
Expand Down Expand Up @@ -1850,7 +1859,7 @@ namespace BasisClasses
}

// Radial grid spacing
double dx = (logxmax - logxmin)/numgrid;
double dx = (logxmax - logxmin)/(numgrid-1);

// Basis storage
Eigen::MatrixXd tabpot, tabden, tabrfc;
Expand Down
2 changes: 1 addition & 1 deletion doc/exp.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ PROJECT_NAME = "EXP"
# could be handy for archiving the generated documentation or if some version
# control system is used.

PROJECT_NUMBER = 7.7.27
PROJECT_NUMBER = 7.7.28

# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewer a
Expand Down
46 changes: 25 additions & 21 deletions exputil/BiorthCyl.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf)
mmax = conf["Lmax"].as<int>();

if (conf["nmaxfid"]) nmaxfid = conf["nmaxfid"].as<int>();
else nmaxfid = 40;
else nmaxfid = 256;

if (conf["nmax"]) nmax = conf["nmax"].as<int>();
else nmax = nmaxfid;
Expand Down Expand Up @@ -75,11 +75,8 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf)
if (conf["scale"]) scale = conf["scale"].as<double>();
else scale = 0.01;

if (conf["acyltbl"]) acyltbl = conf["acyltbl"].as<double>();
else acyltbl = 0.6;

if (conf["cachename"]) cachename = conf["cachename"].as<std::string>();
else cachename = default_cache;
else throw std::runtime_error("BiorthCyl: you must specify a cachename");

// Add output directory and runtag
cachename = outdir + cachename;
Expand All @@ -101,6 +98,9 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf)

if (conf["EVEN_M"]) EVEN_M = conf["EVEN_M"].as<bool>();
else EVEN_M = false;

if (conf["diskconf"]) diskconf = conf["diskconf"];
else throw std::runtime_error("BiorthCyl: you must specify the diskconf stanza");
}
catch (YAML::Exception & error) {
if (myid==0) std::cout << "Error parsing parameters in BiorthCyl: "
Expand All @@ -116,6 +116,11 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf)

geometry = "cylinder";
forceID = "BiorthCyl";

// Hardwired for testing. Could be parameters.
//
logr = false;
biorth = "bess";

initialize();

Expand Down Expand Up @@ -158,15 +163,9 @@ void BiorthCyl::initialize()

void BiorthCyl::create_tables()
{
// Hardwired for testing. Could be variables.
//
bool logr = false;
std::string target("expon");
std::string biorth("bess");

emp = EmpCyl2d(mmax, nmaxfid, knots, numr,
rcylmin*scale, rcylmax*scale, acyltbl*scale, scale,
cmapR, logr, target, biorth);
emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr,
rcylmin*scale, rcylmax*scale,
scale, cmapR, logr, diskconf, biorth);

if (conf["basis"]) emp.basisTest(true);

Expand Down Expand Up @@ -476,7 +475,6 @@ void BiorthCyl::WriteH5Params(HighFive::File& file)
file.createAttribute<double> ("rcylmin", HighFive::DataSpace::From(rcylmin)).write(rcylmin);
file.createAttribute<double> ("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax);
file.createAttribute<double> ("scale", HighFive::DataSpace::From(scale)).write(scale);
file.createAttribute<double> ("acyltbl", HighFive::DataSpace::From(acyltbl)).write(acyltbl);
file.createAttribute<int> ("cmapR", HighFive::DataSpace::From(cmapR)).write(cmapR);
file.createAttribute<int> ("cmapZ", HighFive::DataSpace::From(cmapZ)).write(cmapZ);
}
Expand Down Expand Up @@ -597,7 +595,7 @@ bool BiorthCyl::ReadH5Cache()
double v; HighFive::Attribute vv = h5file.getAttribute(name); vv.read(v);
if (fabs(value - v) < 1.0e-16) return true;
else {
if (myid==0) std::cout << "ReadH5Cache(): wanted "
if (myid==0) std::cout << "---- BiortyCyl::ReadH5Cache: wanted "
<< name << "=" << value
<< " but found "
<< name << "=" << v << std::endl;
Expand Down Expand Up @@ -635,7 +633,6 @@ bool BiorthCyl::ReadH5Cache()
if (not checkDbl(rcylmin, "rcylmin")) return false;
if (not checkDbl(rcylmax, "rcylmax")) return false;
if (not checkDbl(scale, "scale")) return false;
if (not checkDbl(acyltbl, "acyltbl")) return false;
if (not checkInt(cmapR, "cmapR")) return false;
if (not checkInt(cmapZ, "cmapZ")) return false;

Expand All @@ -649,6 +646,14 @@ bool BiorthCyl::ReadH5Cache()
std::cerr << "---- BiorthCyl::ReadH5Cache: "
<< "read <" << cachename << ">" << std::endl;

// Create the basis instance
//
emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr,
rcylmin*scale, rcylmax*scale,
scale, cmapR, logr, diskconf, biorth);

if (conf["basis"]) emp.basisTest(true);

return true;

} catch (HighFive::Exception& err) {
Expand Down Expand Up @@ -755,7 +760,6 @@ YAML::Node BiorthCyl::getHeader(const std::string& cachefile)
node["rcylmin"] = getDbl("rcylmin");
node["rcylmax"] = getDbl("rcylmax");
node["scale"] = getDbl("scale");
node["acyltbl"] = getDbl("acyltbl");
node["cmapR"] = getInt("cmapR");
node["cmapZ"] = getInt("cmapZ");
}
Expand Down Expand Up @@ -797,9 +801,9 @@ BiorthCyl::cacheInfo(const std::string& cachefile, bool verbose)
std::vector<Eigen::MatrixXd> BiorthCyl::orthoCheck()
{
if (not emp()) {
emp = EmpCyl2d(mmax, nmaxfid, knots, numr,
rcylmin*scale, rcylmax*scale, acyltbl*scale, scale,
cmapR, false, "expon", "bess");
emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr,
rcylmin*scale, rcylmax*scale, scale,
cmapR, logr, diskconf, biorth);
}

return emp.orthoCheck();
Expand Down
2 changes: 1 addition & 1 deletion exputil/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ set(PHASE_SRC phase.cc ensemble.cc io_ensemble.cc move_ensemble.cc
rotcurv.cc needle.cc hubble.cc quad.cc)
set(SPECFUNC_SRC gammln.cc bessel.cc OrthoPoly.cc CauchyPV.cc) # modbessel.cc
set(INTERP_SRC Spline.cc SplintE.cc Vodd2.cc Vlocate.cc levsurf.cc Interp1d.cc Cheby1d.cc MonotCubicInterpolator.cc)
set(MASSMODEL_SRC massmodel.cc massmodel_dist.cc embedded.cc isothermal.cc realize_model.cc GenPoly.cc
set(MASSMODEL_SRC massmodel.cc massmodel_dist.cc embedded.cc isothermal.cc realize_model.cc GenPoly.cc mestel.cc
toomre.cc exponential.cc)
set(ORBIT_SRC orbit.cc orbit_trans.cc FindOrb.cc)

Expand Down
Loading

0 comments on commit 071fd19

Please sign in to comment.