From 9c62018beb8fe75f0cd91a0f292b5d93ce170f47 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 20 Dec 2023 17:07:13 -0500 Subject: [PATCH 01/65] Implement inner/outer tapered Mestel disk - Added "zang" type disk - Added parameter vector rather than single parameter for disk density targets - Updated "testeof" routine to use "zang" type disk --- exputil/BiorthCyl.cc | 30 +++++-- exputil/EmpCyl2d.cc | 183 ++++++++++++++++++++++++++++++++++--------- include/BiorthCyl.H | 5 +- include/EmpCyl2d.H | 28 ++++--- src/FlatDisk.cc | 4 + utils/SL/EOF2d.cc | 72 +++++++++++------ 6 files changed, 240 insertions(+), 82 deletions(-) diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index 24e90f634..e4305e9bf 100644 --- a/exputil/BiorthCyl.cc +++ b/exputil/BiorthCyl.cc @@ -78,6 +78,18 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) if (conf["acyltbl"]) acyltbl = conf["acyltbl"].as(); else acyltbl = 0.6; + if (conf["acylcut"]) acylcut = conf["acylcut"].as(); + else acylcut = 10.0; + + if (conf["Ninner"]) Ninner = conf["Ninner"].as(); + else Ninner = 2.0; + + if (conf["Mouter"]) Mouter = conf["Mouter"].as(); + else Mouter = 4.0; + + if (conf["disktype"]) disktype = conf["disktype"].as(); + else disktype = "expon"; + if (conf["cachename"]) cachename = conf["cachename"].as(); else cachename = default_cache; @@ -161,12 +173,13 @@ 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); + EmpCyl2d::ParamVec par {acyltbl*scale, acylcut, Ninner, Mouter}; + + emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, + rcylmin*scale, rcylmax*scale, par, scale, + cmapR, logr, disktype, biorth); if (conf["basis"]) emp.basisTest(true); @@ -477,6 +490,9 @@ void BiorthCyl::WriteH5Params(HighFive::File& file) file.createAttribute ("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax); file.createAttribute ("scale", HighFive::DataSpace::From(scale)).write(scale); file.createAttribute ("acyltbl", HighFive::DataSpace::From(acyltbl)).write(acyltbl); + file.createAttribute ("acylcut", HighFive::DataSpace::From(acylcut)).write(acylcut); + file.createAttribute ("Ninner", HighFive::DataSpace::From(Ninner)).write(Ninner); + file.createAttribute ("Mouter", HighFive::DataSpace::From(Mouter)).write(Mouter); file.createAttribute ("cmapR", HighFive::DataSpace::From(cmapR)).write(cmapR); file.createAttribute ("cmapZ", HighFive::DataSpace::From(cmapZ)).write(cmapZ); } @@ -797,8 +813,10 @@ BiorthCyl::cacheInfo(const std::string& cachefile, bool verbose) std::vector BiorthCyl::orthoCheck() { if (not emp()) { - emp = EmpCyl2d(mmax, nmaxfid, knots, numr, - rcylmin*scale, rcylmax*scale, acyltbl*scale, scale, + EmpCyl2d::ParamVec par {acyltbl*scale, acyltbl*scale*10.0, 2.0, 2.0}; + + emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, + rcylmin*scale, rcylmax*scale, par, scale, cmapR, false, "expon", "bess"); } diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 3c2152076..31c182d4c 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -336,12 +336,12 @@ class EmpCyl2d::ExponCyl : public EmpCyl2d::ModelCyl private: - double sigma0; + double acyl, sigma0; public: - ExponCyl(double scl) - { acyl = scl; sigma0 = 0.5/(M_PI*acyl*acyl); id = "expon"; } + ExponCyl(const EmpCyl2d::ParamVec& par) + { acyl = par[0]; sigma0 = 0.5/(M_PI*acyl*acyl); id = "expon"; } double pot(double r) { double y = 0.5 * r / acyl; @@ -365,9 +365,14 @@ class EmpCyl2d::ExponCyl : public EmpCyl2d::ModelCyl class EmpCyl2d::KuzminCyl : public EmpCyl2d::ModelCyl { +private: + + double acyl; + public: - KuzminCyl(double scl) { acyl = scl; id = "kuzmin"; } + KuzminCyl(const EmpCyl2d::ParamVec& par) + { acyl = par[0]; id = "kuzmin"; } double pot(double R) { double a2 = acyl * acyl; @@ -392,9 +397,14 @@ class EmpCyl2d::KuzminCyl : public EmpCyl2d::ModelCyl class EmpCyl2d::MestelCyl : public EmpCyl2d::ModelCyl { +private: + + double acyl; + public: - MestelCyl(double scl) { acyl = scl; id = "mestel"; } + MestelCyl(const EmpCyl2d::ParamVec& par) + { acyl = par[0]; id = "mestel"; } double pot(double R) { return M_PI/(2.0*acyl)*log(0.5*R/acyl); @@ -418,8 +428,78 @@ class EmpCyl2d::MestelCyl : public EmpCyl2d::ModelCyl }; +class EmpCyl2d::ZangCyl : public EmpCyl2d::ModelCyl +{ + +private: + //! Parameters + double a, c, N, M; + + //! Normalization + double norm = 1.0; + + //! Jacobian for inf to unit-segment mapping + double drdz(double x) + { + double ra = x/a, faca = pow(x/a, N); + return a * pow(ra, 1.0 - N) * (faca + 1.0)*(faca + 1.0)/N; + } + + //! Normalize the surface density + void compute_norm() + { + const int grid = 4000; + double cur = 0.0; + + norm = 1.0; + + double dz = 1.0/grid, z = 0.0, last = 0.0; + + for (int i=1; i -EmpCyl2d::createModel(const std::string type, double acyl) +EmpCyl2d::createModel(const std::string type, const EmpCyl2d::ParamVec& par) { // Convert ID string to lower case // @@ -429,26 +509,37 @@ EmpCyl2d::createModel(const std::string type, double acyl) if (data.find("kuzmin") != std::string::npos) { if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making a Kuzmin disk" << std::endl; - return std::make_shared(acyl); + std::cout << "---- EmpCyl2d::ModelCyl: Making a Kuzmin disk with A=" + << par[0] << std::endl; + return std::make_shared(par); } if (data.find("mestel") != std::string::npos) { if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making a finite Mestel disk" << std::endl; - return std::make_shared(acyl); + std::cout << "---- EmpCyl2d::ModelCyl: Making a finite Mestel disk " + << "with A=" << par[0] << std::endl; + return std::make_shared(par); + } + + if (data.find("zang") != std::string::npos) { + if (myid==0) + std::cout << "---- EmpCyl2d::ModelCyl: Making a double-tapered Zang " + "disk with A=" << par[0] << std::endl; + return std::make_shared(par); } if (data.find("expon") != std::string::npos) { if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk" << std::endl; - return std::make_shared(acyl); + std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk " + << "with A=" << par[0] << std::endl; + return std::make_shared(par); } // Default if nothing else matches if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk [Default]" << std::endl; - return std::make_shared(acyl); + std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk [Default]" + << " with A=" << par[0] << std::endl; + return std::make_shared(par); } @@ -516,47 +607,51 @@ const std::string EmpCyl2d::default_cache_name = ".eof_cache_2d"; // The main constructor // -EmpCyl2d::EmpCyl2d(int mmax, int nmax, int knots, int numr, - double rmin, double rmax, double A, double scale, - bool cmap, bool logr, +EmpCyl2d::EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, + double rmin, double rmax, const EmpCyl2d::ParamVec& par, + double scale, bool cmap, bool logr, const std::string model, const std::string biorth, const std::string cache) : - mmax(mmax), nmax(nmax), knots(knots), numr(numr), - rmin(rmin), rmax(rmax), acyl(A), scale(scale), cmap(cmap), logr(logr), + mmax(mmax), nmaxfid(nmaxfid), nmax(nmax), knots(knots), numr(numr), + rmin(rmin), rmax(rmax), par(par), scale(scale), cmap(cmap), logr(logr), model(model), biorth(biorth), cache_name_2d(cache) { if (cache_name_2d.size()==0) cache_name_2d = default_cache_name; - disk = createModel(model, acyl); - basis = Basis2d::createBasis(mmax, nmax, rmax, biorth); + disk = createModel(model, par); + basis = Basis2d::createBasis(mmax, nmaxfid, rmax, biorth); basis_test = false; if (not ReadH5Cache()) create_tables(); configured = true; + + nmax = std::min(nmax, nmaxfid); } -EmpCyl2d::EmpCyl2d(int mmax, int nmax, int knots, int numr, - double rmin, double rmax, double A, double scale, - bool cmap, bool logr, +EmpCyl2d::EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, + double rmin, double rmax, const EmpCyl2d::ParamVec& par, + double scale, bool cmap, bool logr, std::shared_ptr disk, const std::string biorth, const std::string cache) : - mmax(mmax), nmax(nmax), knots(knots), numr(numr), - rmin(rmin), rmax(rmax), acyl(A), scale(scale), cmap(cmap), logr(logr), + mmax(mmax), nmaxfid(nmaxfid), nmax(nmax), knots(knots), numr(numr), + rmin(rmin), rmax(rmax), par(par), scale(scale), cmap(cmap), logr(logr), disk(disk), biorth(biorth), cache_name_2d(cache) { if (cache_name_2d.size()==0) cache_name_2d = default_cache_name; model = disk->ID(); - basis = Basis2d::createBasis(mmax, nmax, rmax, biorth); + basis = Basis2d::createBasis(mmax, nmaxfid, rmax, biorth); basis_test = false; if (not ReadH5Cache()) create_tables(); configured = true; + + nmax = std::min(nmax, nmaxfid); } @@ -571,7 +666,7 @@ void EmpCyl2d::create_tables() dpot_array.resize(mmax+1); rot_matrix.resize(mmax+1); - Eigen::MatrixXd D(nmax, nmax); + Eigen::MatrixXd D(nmaxfid, nmaxfid); for (int m=0; m<=mmax; m++) { @@ -587,8 +682,8 @@ void EmpCyl2d::create_tables() double rr = map.xi_to_r(xx); double fac = lw.weight(k) * rr / map.d_xi_to_r(xx) * (ximax - ximin); - for (int j=0; jdens(rr) * basis->potl(m, j, rr) * basis->potl(m, l, rr) / sqrt(basis->norm(j, m)*basis->norm(l, m)); @@ -613,7 +708,7 @@ void EmpCyl2d::create_tables() } double dr = (lrmax - lrmin)/(numr - 1), r; - Eigen::VectorXd pot(nmax), den(nmax), dph(nmax); + Eigen::VectorXd pot(nmaxfid), den(nmaxfid), dph(nmaxfid); potl_array[m].resize(numr, nmax); dens_array[m].resize(numr, nmax); @@ -627,7 +722,7 @@ void EmpCyl2d::create_tables() if (logr) r = exp(r); if (m==0) xgrid[i] = r; - for (int n=0; npotl(m, n, r) / sqrt(basis->norm(n, m)); den(n) = basis->dens(m, n, r) / sqrt(basis->norm(n, m)); dph(n) = basis->dpot(m, n, r) / sqrt(basis->norm(n, m)); @@ -655,7 +750,7 @@ void EmpCyl2d::writeBasis(int M, const std::string& filename) if (out) { out << std::endl << "# EOF basis grid for M=" << M << ":" << std::endl; - for (int i=0; i ("mmax", HighFive::DataSpace::From(mmax)). write(mmax); - file.createAttribute ("nmax", HighFive::DataSpace::From(nmax)). write(nmax); + file.createAttribute ("nmaxfid", HighFive::DataSpace::From(nmaxfid)). write(nmaxfid); + file.createAttribute ("nmax", HighFive::DataSpace::From(nmax)).write(nmax); file.createAttribute ("numr", HighFive::DataSpace::From(numr)). write(numr); file.createAttribute ("knots", HighFive::DataSpace::From(knots)). write(knots); file.createAttribute ("ilogr", HighFive::DataSpace::From(ilogr)). write(ilogr); @@ -770,7 +866,7 @@ void EmpCyl2d::WriteH5Cache() file.createAttribute ("rmin", HighFive::DataSpace::From(rmin)). write(rmin); file.createAttribute ("rmax", HighFive::DataSpace::From(rmax)). write(rmax); file.createAttribute ("scale", HighFive::DataSpace::From(scale)). write(scale); - file.createAttribute ("acyl", HighFive::DataSpace::From(acyl)). write(acyl); + file.createAttribute ("params", HighFive::DataSpace::From(par)).write(par); file.createAttribute("model", HighFive::DataSpace::From(model)). write(model); file.createAttribute("biorth", HighFive::DataSpace::From(biorth)).write(biorth); @@ -824,6 +920,15 @@ bool EmpCyl2d::ReadH5Cache() if (fabs(value - v) < 1.0e-16) return true; return false; }; + auto checkArr = [&file](ParamVec& value, std::string name) + { + ParamVec v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); + for (int i=0; i 1.0e-16) return false; + } + return true; + }; + auto checkStr = [&file](std::string value, std::string name) { std::string v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); @@ -836,7 +941,7 @@ bool EmpCyl2d::ReadH5Cache() if (cmap) icmap = 1; if (not checkInt(mmax, "mmax")) return false; - if (not checkInt(nmax, "nmax")) return false; + if (not checkInt(nmaxfid, "nmaxfid")) return false; if (not checkInt(nmax, "nmax")) return false; if (not checkInt(numr, "numr")) return false; if (not checkInt(knots, "knots")) return false; @@ -845,9 +950,9 @@ bool EmpCyl2d::ReadH5Cache() if (not checkDbl(rmin, "rmin")) return false; if (not checkDbl(rmax, "rmax")) return false; if (not checkDbl(scale, "scale")) return false; - if (not checkDbl(acyl, "acyl")) return false; + if (not checkArr(par, "params")) return false; if (not checkStr(model, "model")) return false; - if (not checkStr(biorth, "biorth")) return false; + if (not checkStr(biorth, "biorth")) return false; // Arrays // @@ -1019,7 +1124,7 @@ void EmpCyl2d::checkCoefs() if (myid) return; Mapping map(scale, cmap); - auto disk = createModel(model, acyl); + auto disk = createModel(model, par); LegeQuad lw(knots); Eigen::VectorXd coefs(nmax), coef0(nmax); diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index 790201aa4..eda0a08ad 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -46,7 +46,7 @@ protected: std::string geometry, forceID; int mmax, nmax, numr, nmaxfid, mmin, mlim, nmin, nlim, knots, NQDHT; - double rcylmin, rcylmax, scale, acyltbl; + double rcylmin, rcylmax, scale, acyltbl, acylcut, Ninner, Mouter; bool EVEN_M, verbose; @@ -82,6 +82,9 @@ protected: //! Default cache file name const std::string default_cache = ".biorth_cache"; + //! Density target name + std::string disktype; + //! Active cache file name std::string cachename; diff --git a/include/EmpCyl2d.H b/include/EmpCyl2d.H index 8a12c6892..b220d9538 100644 --- a/include/EmpCyl2d.H +++ b/include/EmpCyl2d.H @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -26,12 +27,15 @@ class EmpCyl2d { public: + + //! Datatype for passing model parameters + using ParamVec = std::array; + //! A two-dimensional disk model for computing the EOF class ModelCyl { protected: - double acyl; std::string id; public: @@ -45,9 +49,10 @@ public: protected: - double acyl, rmin, rmax, scale; + ParamVec par; + double rmin, rmax, scale; double xmin, xmax, dxi; - int mmax, nmax, numr, knots; + int mmax, nmaxfid, numr, knots, nmax; std::string model, biorth; bool cmap, logr, basis_test, configured; @@ -65,7 +70,7 @@ protected: //! Basis creation factory static std::shared_ptr - createBasis(int mmax, int nmax, double rmax, const std::string& type); + createBasis(int mmax, int nmaxfid, double rmax, const std::string& type); //@{ //! Required basis-function members @@ -83,11 +88,12 @@ protected: std::shared_ptr basis; static std::shared_ptr - createModel(const std::string type, double A); + createModel(const std::string type, const ParamVec& params); class ExponCyl; class MestelCyl; class KuzminCyl; + class ZangCyl; //! Map the radius class Mapping @@ -142,15 +148,15 @@ public: EmpCyl2d() : configured(false) {} //! Constructor - EmpCyl2d(int mmax, int nmax, int knots, int numr, - double rmin, double rmax, double A, double scale, - bool cmap, bool logr, const std::string model, + EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, + double rmin, double rmax, const ParamVec& P, + double scale, bool cmap, bool logr, const std::string model, const std::string biorth, const std::string cache=""); //! Constructor with user-supplied target model - EmpCyl2d(int mmax, int nmax, int knots, int numr, - double rmin, double rmax, double A, double scale, - bool cmap, bool logr, std::shared_ptr disk, + EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, + double rmin, double rmax, const ParamVec& P, + double scale, bool cmap, bool logr, std::shared_ptr disk, const std::string biorth, const std::string cache=""); //! Use underlying basis for testing? Default: false. diff --git a/src/FlatDisk.cc b/src/FlatDisk.cc index 6875ef853..4a727ebe8 100644 --- a/src/FlatDisk.cc +++ b/src/FlatDisk.cc @@ -13,6 +13,9 @@ FlatDisk::valid_keys = { "rcylmin", "rcylmax", "acyltbl", + "acylcut", + "Ninner", + "Nouter", "mmax", "numx", "numy", @@ -20,6 +23,7 @@ FlatDisk::valid_keys = { "logr", "model", "biorth", + "disktype", "cachename" }; diff --git a/utils/SL/EOF2d.cc b/utils/SL/EOF2d.cc index 32b207c01..915c279f6 100644 --- a/utils/SL/EOF2d.cc +++ b/utils/SL/EOF2d.cc @@ -14,16 +14,16 @@ int main(int argc, char** argv) { bool logr = false, cmap = false, ortho = false, plane = false; - int numr, mmax, nmax, knots, M, N, nradial, nout; - double A, rmin, rmax; + int numr, mmax, nmaxfid, nmax, knots, M, N, nradial, nout; + double A, rmin, rmax, O, ZN, ZM, rout; std::string filename, type, biorth; // Parse command line // cxxopts::Options options(argv[0], "Computes an EOF two-dimensional disk basis from the Clutton-Brock basis for\n" - "one of the Kuzmin, finite Mestel or Exponential disk targets. [The Mestel\n" - "disk will work very poorly because the Clutton-Brock basis has infinite\n" + "one of the Kuzmin, finite Mestel, Zang or Exponential disk targets. [The\n" + "Mestel disk will work very poorly because the Clutton-Brock basis has infinite\n" "support and looks nothing like the Mestel disk profile.] The new basis, the\n" "orthgonogality matrix and the rotation matrices may be written to files.\n"); options.add_options() @@ -44,31 +44,41 @@ int main(int argc, char** argv) ("totforce", "Compute the total radial force") ("M,harmonic", "Aximuthal harmonic m=0,1,2,3,...", cxxopts::value(M)->default_value("0")) - ("N,norder", "Default number of knots", + ("N,nsize", "Default radial grid size", cxxopts::value(N)->default_value("256")) - ("n,nradial", "Radial order for vertical potential output", + ("n,nradial", "Radial order for vertical potential output", cxxopts::value(nradial)->default_value("0")) - ("A,length", "characteristic disk scale length", + ("A,length", "characteristic disk scale length", cxxopts::value(A)->default_value("1.0")) - ("mmax", "maximum number of angular harmonics in the expansion", + ("O,outer", "outer disk truncation", + cxxopts::value(O)->default_value("10.0")) + ("1,TaperN", "Inner taper exponent", + cxxopts::value(ZN)->default_value("1.0")) + ("2,TaperM", "Outer taper exponent", + cxxopts::value(ZM)->default_value("4.0")) + ("mmax", "maximum number of angular harmonics in the expansion", cxxopts::value(mmax)->default_value("4")) - ("nmax", "maximum number of radial harmonics in the expansion", + ("nmaxfid", "maximum number of radial basis harmonics for EOF construction", + cxxopts::value(nmaxfid)->default_value("64")) + ("nmax", "maximum number of radial harmonics in the expansion", cxxopts::value(nmax)->default_value("10")) - ("numr", "radial knots for the SL grid", + ("numr", "radial knots for the SL grid", cxxopts::value(numr)->default_value("4000")) - ("r,rmin", "minimum radius for the SL grid", + ("r,rmin", "minimum radius for the SL grid", cxxopts::value(rmin)->default_value("0.0001")) - ("R,rmax", "maximum radius for the SL grid", + ("R,rmax", "maximum radius for the SL grid", cxxopts::value(rmax)->default_value("20.0")) - ("knots", "Number of Legendre integration knots", + ("knots", "Number of Legendre integration knots", cxxopts::value(knots)->default_value("200")) - ("nout", "number of points in the output grid per side", + ("rout", "Outer radius for evaluation", + cxxopts::value(rout)->default_value("10.0")) + ("nout", "number of points in the output grid per side", cxxopts::value(nout)->default_value("40")) - ("type", "Target model type (kuzmin, mestel, expon)", + ("type", "Target model type (kuzmin, mestel, zang, expon)", cxxopts::value(type)->default_value("expon")) - ("biorth", "Biorthogonal type (cb, bess)", + ("biorth", "Biorthogonal type (cb, bess)", cxxopts::value(biorth)->default_value("bess")) - ("o,filename", "Output filename", + ("o,filename", "Output filename", cxxopts::value(filename)->default_value("testeof")) ; @@ -101,9 +111,13 @@ int main(int argc, char** argv) // if (vm.count("cmap")) cmap = true; + // Parameter vector for the EmpCyl2d models + // + EmpCyl2d::ParamVec par {A, O, ZN, ZM}; + // Make the class instance // - EmpCyl2d emp(mmax, nmax, knots, numr, rmin, rmax, A, 1.0, cmap, logr, + EmpCyl2d emp(mmax, nmaxfid, nmax, knots, numr, rmin, rmax, par, A, cmap, logr, type, biorth); if (vm.count("basis")) emp.basisTest(true); @@ -134,8 +148,8 @@ int main(int argc, char** argv) // Define some representative limits // - double Rmax = 4.0*A; - double Zmax = 4.0*A; + double Rmax = rout; + double Zmax = rout; // Grid spacing // @@ -207,9 +221,16 @@ int main(int argc, char** argv) for (int i=0; i Date: Sun, 24 Dec 2023 18:17:09 -0500 Subject: [PATCH 02/65] Added background to tapered Mestel model --- exputil/EmpCyl2d.cc | 33 +++++++++++++++++++++++++++------ utils/SL/EOF2d.cc | 18 +++++++++--------- 2 files changed, 36 insertions(+), 15 deletions(-) diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 31c182d4c..c32e5b95b 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -435,9 +435,15 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::ModelCyl //! Parameters double a, c, N, M; + //! Softening factor + double asoft = 1.0e-8; + //! Normalization double norm = 1.0; + //! Ignore inner cut-off for N<0.05 + bool Inner = true; + //! Jacobian for inf to unit-segment mapping double drdz(double x) { @@ -477,23 +483,38 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::ModelCyl c = par[1]; N = par[2]; M = par[3]; + + if (N<0.05) { + // Exponent is now for mapping only + N = 1.0; + Inner = false; + } + + // Normalization compute_norm(); } - //! Potential (not implemented) + //! Potential (infinite Mestel) double pot(double R) { - return 0.0; + return -2.0*M_PI*norm/a/R; } - //! Potential gradient (not implemented) + //! Potential gradient (infinite Mestel) double dpot(double R) { - return 0.0; + return 2.0*M_PI*norm/a/(R*R); } //! Surface density double dens(double R) { - double faca = pow(R/a, N), facc = pow(R/c, M); - return norm * a/(R+1.0e-8) * faca/(faca + 1.0)/(facc + 1.0); + // Outer cut-out + double facc = pow(R/c, M); + double ret = norm * a/(R+asoft) /(facc + 1.0); + // Inner cut-out + if (Inner) { + double faca = pow(R/a, N); + ret *= faca/(faca + 1.0); + } + return ret; } }; diff --git a/utils/SL/EOF2d.cc b/utils/SL/EOF2d.cc index 915c279f6..8db73dc11 100644 --- a/utils/SL/EOF2d.cc +++ b/utils/SL/EOF2d.cc @@ -131,6 +131,12 @@ int main(int argc, char** argv) if (vm.count("trans")) emp.writeTrans(M, filename + ".trans"); if (vm.count("ortho")) emp.orthoCheck(M, filename + ".ortho"); + // Get field type + // + PotRZ::Field F = PotRZ::Field::potential; + if (vm.count("rforce")) F = PotRZ::Field::rforce; + if (vm.count("zforce")) F = PotRZ::Field::zforce; + emp.checkCoefs(); if (vm.count("vertical")) { @@ -156,12 +162,6 @@ int main(int argc, char** argv) double dR = Rmax/(nout - 1); double dz = Zmax/(nout - 1); - // Get field type - // - PotRZ::Field F = PotRZ::Field::potential; - if (vm.count("rforce")) F = PotRZ::Field::rforce; - if (vm.count("zforce")) F = PotRZ::Field::zforce; - // Potential instance with radially sensitive convergence parameters // PotRZ pot(rmax, N, M); @@ -278,7 +278,7 @@ int main(int argc, char** argv) double Rmax = rout; double dR = Rmax/(nout-1); - Eigen::MatrixXd outP(nmax, nout); + Eigen::MatrixXd outF(nmax, nout); for (int n=0; n Date: Tue, 6 Feb 2024 17:34:47 -0500 Subject: [PATCH 03/65] Updates for fixed m=0 background; e.g. the Zang disk [no ci] --- include/BiorthCyl.H | 8 ++++++++ include/EmpCyl2d.H | 7 +++++++ src/FlatDisk.H | 5 +++++ src/PolarBasis.H | 11 ++++++++++- src/PolarBasis.cc | 17 +++++++++++++---- 5 files changed, 43 insertions(+), 5 deletions(-) diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index eda0a08ad..dd5144acb 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -201,6 +201,14 @@ public: //! Evaluate all orders in matrices; for n-body void get_pot(Eigen::MatrixXd& Vc, Eigen::MatrixXd& Vs, double r, double z); + //! Background evaluation + virtual void background(double r, double z, + double & p, double & dr, double & dz) + { + emp.background(r, p, dr); + dz = 0.0; + } + //! Get the table bounds double getRtable() { return rcylmax*scale; } diff --git a/include/EmpCyl2d.H b/include/EmpCyl2d.H index b220d9538..e3d39761e 100644 --- a/include/EmpCyl2d.H +++ b/include/EmpCyl2d.H @@ -204,6 +204,13 @@ public: //! Check for configuration bool operator()() { return configured; } + + //! Background evaluation + void background(double r, double & p, double & dr) + { + p = disk->pot(r); + dr = disk->dpot(r); + } }; diff --git a/src/FlatDisk.H b/src/FlatDisk.H index ad8ae58df..27c1d9c5d 100644 --- a/src/FlatDisk.H +++ b/src/FlatDisk.H @@ -74,6 +74,11 @@ private: } #endif + //! Background evaluation + virtual void get_pot_background(double r, double z, + double & p, double & dr, double & dz) + { ortho->background(r, z, p, dr, dz); } + void get_dens(double r, double z, Eigen::MatrixXd& p, int tid); void get_potl_dens(double r, double z, diff --git a/src/PolarBasis.H b/src/PolarBasis.H index aecffe886..55fe546ed 100644 --- a/src/PolarBasis.H +++ b/src/PolarBasis.H @@ -35,6 +35,7 @@ class MixtureBasis; @param NO_M1 true to omint the m=1 harmonic (default: false) @param EVEN_M true to include only even m harmonics (default: false) @param M0_ONLY true to include only the m=0 harmonic (default: false) + @param M0_BACK true includes fixed m=0 harmonic (default: false) @param ssfrac set > 0.0 to compute a fraction of particles only @param playback true to replay from a coefficient file @param coefMaster true to have only the root node read and distribute the coefficients @@ -376,6 +377,11 @@ protected: #endif + //! Background evaluation + virtual void get_pot_background(double r, double z, + double & p, double & dr, double & dz) + { p = dr = dz = 0.0; } + //! Thread method for accerlation compuation virtual void * determine_acceleration_and_potential_thread(void * arg); @@ -394,7 +400,7 @@ protected: //! Flag self_consitency bool self_consistent; - //! Flag fixed monopole + //! Omit monopole bool NO_M0; //! Flag drop m=1 term @@ -406,6 +412,9 @@ protected: //! Use axisymmetric terms bool M0_only; + //! Flag fixed monopole + bool M0_back; + //! Expected coefficient values for RMS noise computation Eigen::VectorXd meanC; diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index cee056eab..2c0f09827 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -40,6 +40,7 @@ PolarBasis::valid_keys = { "NO_M1", "EVEN_M", "M0_ONLY", + "M0_BACK", "mlim", "ssfrac", "playback", @@ -69,6 +70,7 @@ PolarBasis::PolarBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : NO_M1 = false; EVEN_M = false; M0_only = false; + M0_back = false; ssfrac = 0.0; subset = false; coefMaster = true; @@ -104,6 +106,7 @@ PolarBasis::PolarBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (conf["NO_M1"]) NO_M1 = conf["NO_M1"]. as(); if (conf["EVEN_M"]) EVEN_M = conf["EVEN_M"]. as(); if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); + if (conf["M0_BACK"]) M0_back = conf["M0_BACK"].as(); if (conf["ssfrac"]) { ssfrac = conf["ssfrac"].as(); @@ -1419,11 +1422,17 @@ void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) // Add m=0 force contribution // if (not NO_M0) { - get_pot_coefs_safe(0, *expcoef[0], p, drc, dzc, potd[id], dpotR[id], dpotZ[id]); - potl = mfac*norm0 * p; - potr = mfac*norm0 * drc; - potz = mfac*norm0 * dzc; + if (M0_back) { + get_pot_background(r, zz, p, drc, dzc); + } else { + get_pot_coefs_safe(0, *expcoef[0], p, drc, dzc, + potd[id], dpotR[id], dpotZ[id]); + + potl = mfac*norm0 * p; + potr = mfac*norm0 * drc; + potz = mfac*norm0 * dzc; + } } // Asymmetric terms? From 458783f88317be47f9dfabeb545905bef6c2293e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 8 Feb 2024 19:31:48 -0500 Subject: [PATCH 04/65] Added a draft IC generator for the doubly tapered Mestel disk --- exputil/EmpCyl2d.cc | 119 +++++++++++-------------- include/mestel.H | 159 ++++++++++++++++++++++++++++++--- utils/ICs/CMakeLists.txt | 5 +- utils/ICs/ZangICs.cc | 184 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 383 insertions(+), 84 deletions(-) create mode 100644 utils/ICs/ZangICs.cc diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index c32e5b95b..0ef3379c9 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -399,121 +399,100 @@ class EmpCyl2d::MestelCyl : public EmpCyl2d::ModelCyl { private: - double acyl; + double vrot, rot; public: MestelCyl(const EmpCyl2d::ParamVec& par) - { acyl = par[0]; id = "mestel"; } + { vrot = par[0]; rot = vrot*vrot; id = "mestel"; } - double pot(double R) { - return M_PI/(2.0*acyl)*log(0.5*R/acyl); + virtual double pot(double R) { + if (R>0.0) return rot*log(R); + else throw std::runtime_error("MestelCyl::pot: R<=0"); } - double dpot(double R) { - double a2 = acyl * acyl; - double fac = sqrt(1.0 + R*R/a2); - return M_PI/(2.0*acyl*R); + virtual double dpot(double R) { + if (R>0.0) return rot/R; + else throw std::runtime_error("MestelCyl::dpot: R<=0"); } - double dens(double R) { - if (R>acyl) - return 0.0; - else - return 4.0*M_PI/(2.0*M_PI*acyl*R)*acos(R/acyl); - // ^ - // | - // This 4pi from Poisson's eqn + virtual double dens(double R) { + if (R>0.0) return rot/(2.0*M_PI*R); + else throw std::runtime_error("MestelCyl::dens: R<=0"); } }; -class EmpCyl2d::ZangCyl : public EmpCyl2d::ModelCyl +class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl { private: //! Parameters - double a, c, N, M; + double vr, mu, nu, ri; //! Softening factor double asoft = 1.0e-8; - //! Normalization - double norm = 1.0; - //! Ignore inner cut-off for N<0.05 bool Inner = true; - //! Jacobian for inf to unit-segment mapping - double drdz(double x) + //! Taper factor + double Tifac; + + //! Inner taper function + double Tinner(double Jp) { - double ra = x/a, faca = pow(x/a, N); - return a * pow(ra, 1.0 - N) * (faca + 1.0)*(faca + 1.0)/N; + double fac = pow(Jp, nu); + return fac/(Tifac + fac); } - //! Normalize the surface density - void compute_norm() + //! Outer taper function + double Touter(double Jp) { - const int grid = 4000; - double cur = 0.0; - - norm = 1.0; - - double dz = 1.0/grid, z = 0.0, last = 0.0; + return 1.0/(1.0 + pow(Jp/vr, mu)); + } - for (int i=1; i + +//! Infiinte Mestel Disk class MestelDisk : public AxiSymModel { -private: +protected: + double vrot; double rot; double rmin; double rmax; @@ -14,21 +18,30 @@ private: public: + //! Constructor MestelDisk(double VROT = 1.0, - double RMIN = 1.0e-6, double RMAX = 1.0e6) - { rot = VROT*VROT; rmin = RMIN; rmax = RMAX; dist_defined = true; - dim=2; ModelID = "MestelDisk"; } + double RMIN = 1.0e-6, double RMAX = 1.0e6) + { + vrot = VROT; + rot = VROT*VROT; + rmin = RMIN; + rmax = RMAX; + dim = 2; + + dist_defined = true; + ModelID = "MestelDisk"; + } - // Required member functions - + //@{ + //! Required member functions double get_mass(const double r) { if (r>0.0) return rot*r; else bomb("radius cannot be zero!"); return 0; } - double get_density(const double r) { + virtual double get_density(const double r) { if (r>0.0) return rot/(2.0*M_PI*r); else bomb("radius cannot be zero!"); return 0; @@ -56,43 +69,163 @@ public: if (r>0.0) {ur = rot*log(r); dur = rot/r;} else bomb("radius cannot be zero!"); } + //@} - // Addiional member functions + //@{ + //! Addiional member functions double get_min_radius(void) { return rmin; } double get_max_radius(void) { return rmax; } - void setup_df(double sigma) { + virtual void setup_df(double sigma) { sig2 = sigma*sigma; q = rot/sig2 - 1.0; F = rot/(4.0*M_PI) / ( sqrt(M_PI) * exp(lgamma(0.5*(q+1.0)) + (2.0 + q)*log(sigma) + 0.5*q*log(2.0)) ); } - double distf(double E, double L) { + virtual double distf(double E, double L) { L = fabs(L); if (L==0.0) return 0.0; return F*pow(L, q) * exp(-E/sig2); } - double dfde(double E, double L) { + virtual double dfde(double E, double L) { L = fabs(L); if (L==0.0) return 0.0; return -F*pow(L, q) * exp(-E/sig2)/sig2; } - double dfdl(double E, double L) { + virtual double dfdl(double E, double L) { int sgn=1; if (L<0) {sgn=-1; L *= sgn;} if (L==0.0) return 0.0; return q*F*pow(L, q-1.0) * exp(-E/sig2) * sgn; } - double d2fde2(double E, double L) { + virtual double d2fde2(double E, double L) { L = fabs(L); if (L<=0.0) return 0.0; return F*pow(L, q) * exp(-E/sig2)/sig2/sig2; } + //@} + +}; + + +//! Doubly tapered Mestel Disk +class TaperedMestelDisk : public MestelDisk +{ +protected: + + //! Taper parameters + double nu, mu; + + //! Inner radius + double Ri; + + //! Taper factor + double Tifac; + + //! Inner taper function + double Tinner(double Jp) + { + double fac = pow(Jp, nu); + return fac/(Tifac + fac); + } + + //! Outer taper function + double Touter(double Jp) + { + return 1.0/(1.0 + pow(Jp/vrot, mu)); + } + + //! Deriv of inner taper function + double dTinner(double Jp) + { + double fac = pow(Jp, nu); + double fac2 = Tifac + fac; + return Tifac*nu/Jp/fac2; + } + + //! Deriv of outer taper function + double dTouter(double Jp) + { + double fac = pow(Jp/vrot, mu); + double fac2 = 1.0 + fac; + return -nu*fac/Jp/fac2; + } + + +public: + + //! Constructor + TaperedMestelDisk(double nu, double mu, double Ri, + double VROT = 1.0, + double RMIN = 1.0e-6, double RMAX = 1.0e6) : + nu(nu), mu(mu), Ri(Ri), MestelDisk(VROT, RMIN, RMAX) + { + Tifac = pow(Ri*vrot, nu); + ModelID = "TaperedMestelDisk"; + } + + //! Overloaded density function + double get_density(const double r) { + if (r>0.0) return rot/(2.0*M_PI*r) * Tinner(r) * Touter(r); + else bomb("radius cannot be zero!"); + return 0; + } + + double get_mass_ratio(double Rmin, double Rmax) + { + double dr = 0.01*Ri; + int N = std::floor((Rmax - Rmin)/dr) + 1; + dr = (Rmax - Rmin)/N; + double lst = get_density(Rmin); + double ret = 0.0; + for (int i=1; i<=N; i++) { + double cur = get_density(Rmin+dr*i); + ret += 0.5*dr*(lst + cur); + lst = cur; + } + + return ret / (get_mass(Rmax) - get_mass(Rmin)); + } + + //@{ + //! Overloaded DF member functions + + double distf(double E, double L) + { + L = fabs(L); + if (L==0.0) return 0.0; + return F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2); + } + + double dfde(double E, double L) + { + L = fabs(L); + if (L==0.0) return 0.0; + return -F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2)/sig2; + } + + double dfdl(double E, double L) + { + int sgn=1; + if (L<0) {sgn=-1; L *= sgn;} + if (L==0.0) return 0.0; + double Tfac = pow(L, q)*Tinner(L)*Touter(L); + double dL = Tfac*(q/L + dTinner(L) + dTouter(L)); + return F* dL * exp(-E/sig2) * sgn; + } + + double d2fde2(double E, double L) + { + L = fabs(L); + if (L<=0.0) return 0.0; + return F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2)/sig2/sig2; + } + //@} }; diff --git a/utils/ICs/CMakeLists.txt b/utils/ICs/CMakeLists.txt index 950ed7179..f10de0787 100644 --- a/utils/ICs/CMakeLists.txt +++ b/utils/ICs/CMakeLists.txt @@ -1,7 +1,8 @@ set(bin_PROGRAMS shrinkics gensph gendisk gendisk2d gsphere pstmod empinfo empdump eofbasis eofcomp testcoefs testcoefs2 testdeval - forcetest hdf52accel forcetest2 cylcache modelfit test2d addsphmod cubeics) + forcetest hdf52accel forcetest2 cylcache modelfit test2d addsphmod + cubeics zangics) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp exputil ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) @@ -87,6 +88,8 @@ add_executable(addsphmod addsphmod.cc AddSpheres.cc) add_executable(cubeics cubeICs.cc) +add_executable(zangics ZangICs.cc) + foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) target_include_directories(${program} PUBLIC ${common_INCLUDE}) diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc new file mode 100644 index 000000000..a226fba68 --- /dev/null +++ b/utils/ICs/ZangICs.cc @@ -0,0 +1,184 @@ +/* + A tapered Mestel disk IC generator +*/ + // C++/STL headers +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +int +main(int ac, char **av) +{ + //===================== + // Begin option parsing + //===================== + + int N; // Number of particles + double mu, nu, Ri; // Taper paramters + double Rmin, Rmax; // Radial range + double sigma; // Velocity dispersion + std::string bodyfile; // Output file + unsigned seed; // Will be inialized by /dev/random if + // not set on the command line + + cxxopts::Options options(av[0], "Ideal tapered Mestel IC generator"); + + options.add_options() + ("h,help", "Print this help message") + ("z,zerovel", "Zero the mean velocity") + ("N,number", "Number of particles to generate", + cxxopts::value(N)->default_value("100000")) + ("n,nu", "Inner taper exponent (0 for no taper)", + cxxopts::value(nu)->default_value("2.0")) + ("m,mu", "Outer taper exponent (0 for no taper)", + cxxopts::value(mu)->default_value("2.0")) + ("i,Ri", "Inner radius for taper", + cxxopts::value(Ri)->default_value("0.1")) + ("r,Rmin", "Inner radius for model", + cxxopts::value(Rmin)->default_value("0.01")) + ("R,Rmax", "Outer radius for model", + cxxopts::value(Rmax)->default_value("10.0")) + ("S,sigma", "Radial velocity dispersion", + cxxopts::value(sigma)) + ("s,seed", "Random number seed. Default: use /dev/random", + cxxopts::value(seed)) + ("o,file", "Output body file", + cxxopts::value(bodyfile)->default_value("cube.bods")) + ; + + cxxopts::ParseResult vm; + + try { + vm = options.parse(ac, av); + } catch (cxxopts::OptionException& e) { + std::cout << "Option error: " << e.what() << std::endl; + exit(-1); + } + + // Print help message and exit + // + if (vm.count("help")) { + std::cout << options.help() << std::endl << std::endl; + return 1; + } + + // Set from /dev/random if not specified + if (vm.count("seed")==0) { + seed = std::random_device{}(); + } + + // Make sure N>0 + if (N<=0) { + std::cerr << av[0] << ": you must requiest at least one body" + << std::endl; + } + + // Open the output file + // + std::ofstream out(bodyfile); + if (not out) { + std::string msg(av[0]); + msg += ": output file <" + bodyfile + "> can not be opened"; + throw std::runtime_error(msg); + } + + // Create the model + // + auto model = std::make_shared(nu, mu, Ri); + + // Create an orbit grid + // + SphericalOrbit orb(model); + + double Ktol = 0.01; + double Kmin = Ktol, Kmax = 1.0 - Ktol; + + double Emin = 0.5*Rmin*model->get_dpot(Rmin) + model->get_pot(Emin); + double Emax = 0.5*Rmax*model->get_dpot(Rmax) + model->get_pot(Rmax); + + // Scan to find the peak df + // + const int num = 100; + double peak = 0.0; + double dE = (Emax - Emin)/num, dK = (1.0 - 2.0*Ktol)/num; + for (int i=0; i<=num; i++) { + double E = Emin + dE*i; + for (int j=0; j<=num; j++) { + double K = Kmin + dK*j; + orb.new_orbit(E, K); + double F = model->distf(E, orb.get_action(2)) / orb.get_freq(1); + peak = std::max(peak, F); + } + } + + // Header + // + out << std::setw(10) << N << std::setw(10) << 0 << std::setw(10) << 0 + << std::endl << std::setprecision(10); + + std::mt19937 gen(seed); + std::uniform_real_distribution<> uniform(0.0, 1.0); + + // Save the position and velocity vectors + std::vector> pos(N), vel(N); + + int itmax = 10000; + int over = 0; + + // Generation loop + for (int n=0; ndistf(E, orb.get_action(2)) / orb.get_freq(1); + if (F/peak > R) break; + } + + if (j==itmax) over++; + + double J = orb.get_action(2); + double T = 2.0*M_PI/orb.get_freq(1)*uniform(gen); + double r = orb.get_angle(6, T); + double w1 = orb.get_angle(1, T); + double phi = 2.0*M_PI*uniform(gen) + orb.get_angle(7, T); + + double vt = J/r; + double vr = sqrt(fabs(2.0*E - model->get_pot(r)) - J*J/(r*r)); + if (w1 > M_PI) vr *= -1.0; + + pos[n][0] = r*cos(phi); + pos[n][1] = r*sin(phi); + pos[n][2] = 0.0; + vel[n][0] = vr*cos(phi) - vt*sin(phi); + vel[n][1] = vr*sin(phi) + vt*cos(phi); + vel[n][2] = 0.0; + } + + std::cout << "** " << over << " particles failed iteration" << std::endl; + + for (int n=0; n Date: Thu, 8 Feb 2024 19:57:00 -0500 Subject: [PATCH 05/65] Compute cumulative mass for Zang disk --- include/mestel.H | 50 ++++++++++++++++++++++++++++++++------------ utils/ICs/ZangICs.cc | 29 ++++++++++++++++++------- 2 files changed, 58 insertions(+), 21 deletions(-) diff --git a/include/mestel.H b/include/mestel.H index df13da9e3..b99866b15 100644 --- a/include/mestel.H +++ b/include/mestel.H @@ -35,7 +35,7 @@ public: //@{ //! Required member functions - double get_mass(const double r) { + virtual double get_mass(const double r) { if (r>0.0) return rot*r; else bomb("radius cannot be zero!"); return 0; @@ -116,8 +116,11 @@ public: //! Doubly tapered Mestel Disk class TaperedMestelDisk : public MestelDisk { +private: + + std::shared_ptr interp; + protected: - //! Taper parameters double nu, mu; @@ -176,20 +179,41 @@ public: return 0; } - double get_mass_ratio(double Rmin, double Rmax) + double get_mass(double R) { - double dr = 0.01*Ri; - int N = std::floor((Rmax - Rmin)/dr) + 1; - dr = (Rmax - Rmin)/N; - double lst = get_density(Rmin); - double ret = 0.0; - for (int i=1; i<=N; i++) { - double cur = get_density(Rmin+dr*i); - ret += 0.5*dr*(lst + cur); - lst = cur; + if (interp) { + // Compute a grid spacing + double dr = 0.01*Ri; + int N = std::floor((rmax - rmin)/dr) + 1; + dr = (rmax - rmin)/N; + + // Storage for table + Eigen::VectorXd vecR(N+1), vecM(N+1); + + double lst = get_density(rmin); + double cum = 0.0; + + // Make the table + // + vecR(0) = rmin; + vecM(0) = 0.0; + + for (int i=1; i<=N; i++) { + double rr = rmin + dr*i; + double cur = get_density(rr); + cum += 0.5*dr*(lst + cur); + lst = cur; + + vecR(i) = rr; + vecM(i) = cum; + } + + // Make the interpolation function + // + interp = std::make_shared(vecR, vecM); } - return ret / (get_mass(Rmax) - get_mass(Rmin)); + return interp->eval(R); } //@{ diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index a226fba68..71c3dc8ab 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -142,6 +142,7 @@ main(int ac, char **av) int j; for (j=0; jget_pot(r)) - J*J/(r*r)); - if (w1 > M_PI) vr *= -1.0; + if (w1 > M_PI) vr *= -1.0; // Branch of radial motion + + // Convert from polar to Cartesian + // pos[n][0] = r*cos(phi); pos[n][1] = r*sin(phi); pos[n][2] = 0.0; + vel[n][0] = vr*cos(phi) - vt*sin(phi); vel[n][1] = vr*sin(phi) + vt*cos(phi); vel[n][2] = 0.0; } - std::cout << "** " << over << " particles failed iteration" << std::endl; + // Compute the particle mass + // + double mass = (model->get_mass(Rmax) - model->get_mass(Rmin))/N; + + std::cout << "** " << over << " particles failed iteration" << std::endl + << "** Particle mass=" << mass << std::endl; + + out << std::setw(8) << N << std::setw(8) << 0 << std::setw(8) << 0 + << std::endl; for (int n=0; n Date: Fri, 9 Feb 2024 16:38:58 -0500 Subject: [PATCH 06/65] Use std::function rather than C-style --- exputil/rombe2d.cc | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/exputil/rombe2d.cc b/exputil/rombe2d.cc index d1bc80c46..fb4e56420 100644 --- a/exputil/rombe2d.cc +++ b/exputil/rombe2d.cc @@ -1,8 +1,7 @@ +#include #include #include -using namespace std; - /* * Computes integral using Simpson's rule with Romberg's * correction with a supplied integrand routine @@ -23,7 +22,7 @@ using namespace std; #define NROMB 12 #define DMACH 1.0e-8 -double rombe2(double a, double b, double (*f) (double), int n) +double rombe2(double a, double b, std::function f, int n) { double t[NROMB],xmesh,d,ends,x2,x4; int nmax,num,i,j,k,index; @@ -32,9 +31,9 @@ double rombe2(double a, double b, double (*f) (double), int n) xmesh = (b-a)/nmax; d = (b-a)/2.0; - ends = (*f)(a) + (*f)(b); + ends = f(a) + f(b); num = nmax/2; - x2 = (*f)(a+xmesh*num); + x2 = f(a+xmesh*num); t[0] = d*(ends+4.0*x2)/3.0; for(i=2; i<=n; i++) { @@ -43,7 +42,7 @@ double rombe2(double a, double b, double (*f) (double), int n) index = num/2; for(j=1; j<=nmax/num; j++) { - x4 = x4 + (*f)(a+xmesh*index); + x4 = x4 + f(a+xmesh*index); index = index + num; } t[i-1] = d*(ends+2.0*x2+4.0*x4)/3.0; From 4bd86b5b57a96a36c998f47c080421a256b3c2f3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 16:39:28 -0500 Subject: [PATCH 07/65] Flip the sign of the boolean test --- include/mestel.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/mestel.H b/include/mestel.H index b99866b15..4be059ef6 100644 --- a/include/mestel.H +++ b/include/mestel.H @@ -181,7 +181,7 @@ public: double get_mass(double R) { - if (interp) { + if (not interp) { // Compute a grid spacing double dr = 0.01*Ri; int N = std::floor((rmax - rmin)/dr) + 1; From ed20fab697f1a4b220f71cefc043baf32341aebf Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 16:40:07 -0500 Subject: [PATCH 08/65] Some trivial clean-up only --- exputil/zbrent.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/exputil/zbrent.cc b/exputil/zbrent.cc index 475e1266e..ee085c988 100644 --- a/exputil/zbrent.cc +++ b/exputil/zbrent.cc @@ -1,3 +1,4 @@ +#include #include #include #include @@ -7,10 +8,10 @@ #include #define ITMAX 100 -// #define EPS 3.0e-8 #define EPS 1.0e-16 -double zbrent(std::function func, double x1, double x2, double tol) +double zbrent(std::function func, + double x1, double x2, double tol) { int iter; double a=x1,b=x2,c,d,e,min1,min2; @@ -18,7 +19,7 @@ double zbrent(std::function func, double x1, double x2, double t if (fb*fa > 0.0) { std::ostringstream str; - str << "Root must be bracketed in ZBRENT:" + str << "root must be bracketed in ZBRENT:" << " f(" << a << ")=" << fa << ", f(" << b << ")=" << fb; throw std::runtime_error(str.str()); From bae633b368a8ff95cb0195642e482cbddefdb692 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 16:40:27 -0500 Subject: [PATCH 09/65] Add 2T/W test --- utils/ICs/ZangICs.cc | 68 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 53 insertions(+), 15 deletions(-) diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index 71c3dc8ab..9be26f849 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -12,9 +12,12 @@ #include #include +#include + #include -#include +#include // Progress bar +#include // Option parsing int main(int ac, char **av) @@ -92,18 +95,32 @@ main(int ac, char **av) throw std::runtime_error(msg); } + SphericalOrbit::ZFRAC=0.3; // TEST + // Create the model // - auto model = std::make_shared(nu, mu, Ri); + auto model = std::make_shared(nu, mu, Ri, 1.0, 0.1*Rmin, 10.0*Rmax); + + // Progress bar + // + std::shared_ptr progress; + int nomp = 1; +#pragma omp parallel + { + nomp = omp_get_num_threads(); + if (omp_get_thread_num()==0) + progress = std::make_shared(N/nomp); + } // Create an orbit grid // - SphericalOrbit orb(model); - + std::vector> orb(nomp); + for (auto & v : orb) v = std::make_shared(model); + double Ktol = 0.01; double Kmin = Ktol, Kmax = 1.0 - Ktol; - double Emin = 0.5*Rmin*model->get_dpot(Rmin) + model->get_pot(Emin); + double Emin = 0.5*Rmin*model->get_dpot(Rmin) + model->get_pot(Rmin); double Emax = 0.5*Rmax*model->get_dpot(Rmax) + model->get_pot(Rmax); // Scan to find the peak df @@ -115,8 +132,8 @@ main(int ac, char **av) double E = Emin + dE*i; for (int j=0; j<=num; j++) { double K = Kmin + dK*j; - orb.new_orbit(E, K); - double F = model->distf(E, orb.get_action(2)) / orb.get_freq(1); + orb[0]->new_orbit(E, K); + double F = model->distf(E, orb[0]->get_action(2)) / orb[0]->get_freq(1); peak = std::max(peak, F); } } @@ -133,11 +150,15 @@ main(int ac, char **av) std::vector> pos(N), vel(N); int itmax = 10000; + int tid = 0; int over = 0; // Generation loop + // +#pragma omp parallel for reduction(+:over) for (int n=0; ndistf(E, orb.get_action(2)) / orb.get_freq(1); + orb[tid]->new_orbit(E, K); + double F = model->distf(E, orb[tid]->get_action(2)) / orb[tid]->get_freq(1); if (F/peak > R) break; } if (j==itmax) over++; - double J = orb.get_action(2); - double T = 2.0*M_PI/orb.get_freq(1)*uniform(gen); - double r = orb.get_angle(6, T); - double w1 = orb.get_angle(1, T); - double phi = 2.0*M_PI*uniform(gen) + orb.get_angle(7, T); + double J = orb[tid]->get_action(2); + double T = 2.0*M_PI/orb[tid]->get_freq(1)*uniform(gen); + double r = orb[tid]->get_angle(6, T); + double w1 = orb[tid]->get_angle(1, T); + double phi = 2.0*M_PI*uniform(gen) + orb[tid]->get_angle(7, T); double vt = J/r; double vr = sqrt(fabs(2.0*E - model->get_pot(r)) - J*J/(r*r)); @@ -174,7 +195,11 @@ main(int ac, char **av) vel[n][0] = vr*cos(phi) - vt*sin(phi); vel[n][1] = vr*sin(phi) + vt*cos(phi); vel[n][2] = 0.0; + + // Print progress bar + if (progress) ++(*progress); } + std::cout << std::endl << "** Main loop complete" << std::endl; // Compute the particle mass // @@ -186,12 +211,25 @@ main(int ac, char **av) out << std::setw(8) << N << std::setw(8) << 0 << std::setw(8) << 0 << std::endl; + double ektot = 0.0, clausius = 0.0; for (int n=0; nget_dpot(r)*r; } + std::cout << "2T/VC=" << ektot/clausius << std::endl; + return 0; } From 87b9775ba22686aae0fdc87e692b9f208cf1c9e5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 16:41:07 -0500 Subject: [PATCH 10/65] Move to functors with lambdas rather than global function pointers which are not multithread friendly --- exputil/orbit.cc | 22 ++- exputil/orbit_trans.cc | 371 +++++++++++++++++++++++------------------ include/orbit.H | 7 + 3 files changed, 225 insertions(+), 175 deletions(-) diff --git a/exputil/orbit.cc b/exputil/orbit.cc index 1d08a6618..b80091ce4 100644 --- a/exputil/orbit.cc +++ b/exputil/orbit.cc @@ -12,23 +12,21 @@ #include #include -void RegularOrbit::bomb(const char *s) { - std::cerr << "ERROR from " << OrbitID << ": " << s << std::endl; -#ifdef DEBUG - abort(); -#endif - exit(-1); +void RegularOrbit::bomb(const char *s) +{ + std::ostringstream msg; + msg << "ERROR from " << OrbitID << ": " << s; + throw std::runtime_error(msg.str()); } void RegularOrbit::bomb(const string& s) { - cerr << "ERROR from " << OrbitID << ": " << s << endl; -#ifdef DEBUG - abort(); -#endif - exit(-1); + std::ostringstream msg; + msg << "ERROR from " << OrbitID << ": " << s; + throw std::runtime_error(msg.str()); } -void RegularOrbit::warn(const string& fct, const string& msg) { +void RegularOrbit::warn(const string& fct, const string& msg) +{ cerr << fct << ": from " << OrbitID << ": " << msg << endl; } diff --git a/exputil/orbit_trans.cc b/exputil/orbit_trans.cc index db777cfff..71e52bb79 100644 --- a/exputil/orbit_trans.cc +++ b/exputil/orbit_trans.cc @@ -12,16 +12,52 @@ #include #include - // Global variables -static AxiSymModPtr mm; -static double EE, JJ, KK; -static double ap, am, sp, sm; - +int SphericalOrbit::Nseg = 40; // Number of trial segments for root finding double SphericalOrbit::tol = 1.0e-8; // root finder tolerance double SphericalOrbit::ZFRAC = 0.001; // fraction below minimum grid for tp search double SphericalOrbit::RMAXF = 3.0; double SphericalOrbit::tolnr = 1.0e-10; // r_apo location using Newton-Raphson refinement +std::tuple SphericalOrbit::search +(std::function func, double rmin, double rmax) +{ + bool use_log = false; + if (rmin > 0.0) { + use_log = true; + rmin = log(rmin); + rmax = log(rmax); + } + + double dx = (rmax - rmin)/Nseg; + bool cond = false; + + // Loop variables + double xend = rmin, xbeg, ylst; + double ycur = func(use_log ? exp(xend) : xend); + + // Now, increment the end point and test for sign change in the + // interval + for (int n=1; n<=Nseg; n++) { + xbeg = xend; + ylst = ycur; + xend += dx; + ycur = func(use_log ? exp(xend) : xend); + + // Found an interval, return it + if (ycur*ylst <= 0.0) { + cond = true; + break; + } + } + + if (use_log) { // Scale back to linear + xbeg = exp(xbeg); + xend = exp(xend); + } + // Done! + return std::tuple(xbeg, xend, cond); +} + void SphericalOrbit::compute_freq(void) { const string rname("compute_freq"); @@ -30,86 +66,123 @@ void SphericalOrbit::compute_freq(void) double tmp2; #endif int i; - double Ecirc(double), denom(double); - - mm = model; - EE = energy; - KK = kappa; // Find turning points - + // double xmin, xmax; xmin = ZFRAC*model->get_min_radius(); xmax = model->get_max_radius(); - if ( Ecirc(xmin)*Ecirc(xmax) < 0.0) { - try { - r_circ = zbrent(Ecirc, xmin, xmax, tol); - } - catch (const std::exception& error) { - std::ostringstream mesg; - mesg << "SphericalOrbit::compute_freq: " << error.what() << std::endl - << "SphericalOrbit::compute_freq @ r_circ: model=" - << model->ModelID << endl - << " energy=" << energy << endl - << " kappa=" << kappa << endl - << " xmin=" << xmin << endl - << " xmax=" << xmax << endl; - throw std::runtime_error(mesg.str()); + // Functor whose zero locates radius of circular orbit with energy EE + // + auto Ecirc = [&](double r) + { + double ur, dudr; + model->get_pot_dpot(r, ur, dudr); + return energy - 0.5*r*dudr - ur; + }; + + // New strategy: break into Nseg segments to find good starting + // points + // + { + auto [xbeg, xend, cond] = search(Ecirc, xmin, xmax); + + if ( Ecirc(xbeg)*Ecirc(xend) < 0.0) { + try { + r_circ = zbrent(Ecirc, xbeg, xend, tol); + } + catch (const std::exception& error) { + std::ostringstream mesg; + mesg << "SphericalOrbit::compute_freq: " << error.what() << std::endl + << "SphericalOrbit::compute_freq @ r_circ: model=" + << model->ModelID << std::endl + << "** energy=" << energy << std::endl + << "** kappa=" << kappa << std::endl + << "** xmin=" << xbeg << std::endl + << "** xmax=" << xend << std::endl + << "** f(xmin)=" << Ecirc(xbeg) << std::endl + << "** f(xmax)=" << Ecirc(xend) << std::endl + << "** cond=" << std::boolalpha << cond << std::endl; + throw std::runtime_error(mesg.str()); + } } - } - else { // Circular radius outside mass distribution - r_circ = -0.5*model->get_mass(model->get_max_radius())/EE; - if (r_circ < model->get_min_radius()) { - r_circ = model->get_min_radius(); + else { // Circular radius outside mass distribution + r_circ = -0.5*model->get_mass(model->get_max_radius())/energy; + if (r_circ < model->get_min_radius()) { + r_circ = model->get_min_radius(); + } + std::string message("SphericalOrbit::compute_freq warning: Circular radius outside mass distribution"); + throw std::runtime_error(message); } - std::string message("SphericalOrbit::compute_freq warning: Circular radius outside mass distribution"); - throw std::runtime_error(message); } dudr = model->get_dpot(r_circ); if (dudr>0.0) jmax = sqrt(r_circ*r_circ*r_circ*dudr); else jmax = 0.0; - JJ = jmax*kappa; + // Functor whose zero locates turning points for orbit (EE,JJ) + // + auto denom = [&](double r) + { + double ur = model->get_pot(r); + double JJ = jmax*kappa; + return 2.0*(energy-ur)*r*r - JJ*JJ; + }; xmin = r_circ; xmax = model->get_max_radius(); #ifdef DEBUG bool used_asymp = false; #endif - if ( denom(xmin)*denom(xmax) < 0.0) { - try { - r_apo = zbrent(denom, xmin, xmax, tol); - } - catch (const std::exception& error) { - std::ostringstream mesg; - mesg << "SphericalOrbit::compute_freq: " << error.what() << std::endl - << "SphericalOrbit::compute_freq @ r_apo[1]: model=" - << model->ModelID << endl - << " energy=" << energy << endl - << " kappa=" << kappa << endl - << " xmin=" << xmin << endl - << " xmax=" << xmax << endl; - throw std::runtime_error(mesg.str()); + + { + auto [xbeg, xend, cond] = search(denom, xmin, xmax); + + if ( denom(xbeg)*denom(xend) < 0.0) { + try { + r_apo = zbrent(denom, xbeg, xend, tol); + } + catch (const std::exception& error) { + std::ostringstream mesg; + mesg << "SphericalOrbit::compute_freq: " << error.what() << std::endl + << "SphericalOrbit::compute_freq @ r_apo[1]: model=" + << model->ModelID << std::endl + << "** energy=" << energy << std::endl + << "** kappa=" << kappa << std::endl + << "** energy=" << energy << std::endl + << "** kappa=" << kappa << std::endl + << "** J=" << kappa*jmax << std::endl + << "** rmin=" << xbeg << std::endl + << "** rmax=" << xend << std::endl + << "** cond=" << std::boolalpha << cond << std::endl; + + throw std::runtime_error(mesg.str()); + } } } + if (r_circ < model->get_max_radius()) { - if (denom(xmin)*denom(RMAXF*xmax) < 0.0) { + + auto [xbeg, xend, cond] = search(denom, r_circ, RMAXF*xmax); + + if (denom(xbeg)*denom(xend) < 0.0) { try { - r_apo = zbrent(denom, xmin, RMAXF*xmax, tol); + r_apo = zbrent(denom, xbeg, xend, tol); } catch (const std::exception& error) { std::ostringstream mesg; mesg << "SphericalOrbit::compute_freq: " << error.what() << std::endl << "SphericalOrbit::compute_freq @ r_apo[2]: model=" - << model->ModelID << endl - << " energy=" << energy << endl - << " pot(max)=" << model->get_pot(xmax) << endl - << " kappa=" << kappa << endl - << " xmin=" << xmin << endl - << " xmax=" << xmax*RMAXF << endl - << " rc check=" << energy - model->get_pot(r_circ) - jmax*jmax/(2.0*r_circ*r_circ) << endl; + << model->ModelID << std::endl + << "** energy=" << energy << std::endl + << "** pot(max)=" << model->get_pot(xend) << std::endl + << "** kappa=" << kappa << std::endl + << "** J=" << kappa*jmax << std::endl + << "** xmin=" << xbeg << std::endl + << "** xmax=" << xend << std::endl + << "** cond=" << std::boolalpha << cond << std::endl + << "** rc check=" << energy - model->get_pot(r_circ) - jmax*jmax/(2.0*r_circ*r_circ) << std::endl; throw std::runtime_error(mesg.str()); } } else { @@ -117,12 +190,13 @@ void SphericalOrbit::compute_freq(void) } } else { // Circular radius outside mass distribution - double r0 = -0.5*model->get_mass(model->get_max_radius())/EE; - r_apo = r0 + sqrt(r0*r0 - 0.5*JJ*JJ/EE); + double JJ = kappa*jmax; + double r0 = -0.5*model->get_mass(model->get_max_radius())/energy; + r_apo = r0 + sqrt(r0*r0 - 0.5*JJ*JJ/energy); // Newton-Raphson refinement (slow . . .) double f, df; for (int i=0; i<200; i++) { - f = 2.0*(EE - model->get_pot(r_apo)) - JJ*JJ/(r_apo*r_apo); + f = 2.0*(energy - model->get_pot(r_apo)) - JJ*JJ/(r_apo*r_apo); if (fabs(f)get_dpot(r_apo) + JJ*JJ/(r_apo*r_apo*r_apo); r_apo += -f/df; @@ -138,33 +212,44 @@ void SphericalOrbit::compute_freq(void) #endif } - if (denom(ZFRAC*model->get_min_radius())*denom(r_circ) >= 0.0) { - r_peri = 0.99*ZFRAC*model->get_min_radius(); - } - else { - try { - r_peri = zbrent(denom, ZFRAC*model->get_min_radius(), r_circ, tol); - } - catch (const std::exception& error) { - std::ostringstream mesg; - mesg << "SphericalOrbit::compute_freq: " << error.what() << std::endl - << "SphericalOrbit::compute_freq @ r_peri: model=" << model->ModelID - << " energy=" << energy - << " kappa=" << kappa - << " xmin=" << ZFRAC*model->get_min_radius() - << " xmax=" << r_circ << endl; - throw std::runtime_error(mesg.str()); + { + auto [xbeg, xend, cond] = search(denom, ZFRAC*model->get_min_radius(), r_circ); + + if (not cond) { + // Assume peri is inside the range + r_peri = 0.99*ZFRAC*model->get_min_radius(); + + } else { + // Find the peri from the bounded interval + try { + r_peri = zbrent(denom, 0.999*xbeg, 1.001*xend, tol); + } + catch (const std::exception& error) { + std::ostringstream mesg; + mesg << "SphericalOrbit::compute_freq: " << error.what() << std::endl + << "SphericalOrbit::compute_freq @ r_peri: model=" << model->ModelID + << std::endl + << "** energy=" << energy << std::endl + << "** kappa=" << kappa << std::endl + << "** J=" << kappa*jmax << std::endl + << "** rmin=" << xbeg << std::endl + << "** rmax=" << xend << std::endl + << "** cond=" << std::boolalpha << cond << std::endl; + + throw std::runtime_error(mesg.str()); + } } } // Ok, prepare to do freq. integrals - + // ap = 0.5*(r_apo + r_peri); am = 0.5*(r_apo - r_peri); sp = ap/(r_apo*r_peri); sm = am/(r_apo*r_peri); // Do using centered rectangle technique + // accum0 = 0.0; accum1 = 0.0; accum2 = 0.0; @@ -251,26 +336,8 @@ void SphericalOrbit::compute_freq(void) #undef FRECS #undef tol -// Function to iteratively locate radius of circular orbit with energy EE -double Ecirc(double r) -{ - double ur, dudr, dif; - - mm->get_pot_dpot(r, ur, dudr); - dif = EE - 0.5*r*dudr - ur; - return(dif); -} - -// Function to iteratively locate turning points for orbit (EE,JJ) - -double denom(double r) -{ - double ur = mm->get_pot(r); - return 2.0*(EE-ur)*r*r - JJ*JJ; -} - // for testing---frequencies in the epicyclic approx - +// void SphericalOrbit::compute_freq_epi(void) { double d2udr2; @@ -285,14 +352,13 @@ void SphericalOrbit::compute_freq_epi(void) } extern -double rombe2(double a, double b, double (*f) (double), int n); +double rombe2(double a, double b, std::function, int n); double dtp, dtm; void SphericalOrbit::compute_angles(void) { double accum1,accum2,r, s, t, sl, tl; - double fw1(double t), ff(double t); int i; l1s = l2s = 0; @@ -313,11 +379,9 @@ void SphericalOrbit::compute_angles(void) } + double JJ = kappa * jmax; + if (freq_defined) { - mm = model; - EE = energy; - KK = kappa; - JJ = jmax*kappa; ap = 0.5*(r_apo + r_peri); am = 0.5*(r_apo - r_peri); sp = ap/(r_apo*r_peri); @@ -326,6 +390,43 @@ void SphericalOrbit::compute_angles(void) else compute_freq(); + + const double tol = 1.0e-8; // tolerance for radicand + + auto fw1 = [&](double t) + { + double r = ap + am*sin(t); + double ur = model->get_pot(r); + double radcand = 2.0*(energy-ur)-JJ*JJ/(r*r); + + if (radcand < tol) { + // values near turning points + double sgn = 1.0; + if (t<0) sgn = -1.0; + r = ap + sgn*am; + double dudr = model->get_dpot(r); + return sqrt( am*sgn/(dudr-JJ*JJ/(r*r*r)) ); + } + return am*cos(t)/sqrt(radcand); + }; + + auto ff = [&](double t) + { + double s = sp + sm*sin(t); + double ur = model->get_pot(1.0/s); + double radcand = 2.0*(energy-ur) - JJ*JJ*s*s; + + if (radcand < tol) { + // values near turning points + double sgn = 1.0; + if (t<0) sgn = -1.0; + s = sp + sgn*sm; + double dudr = model->get_dpot(1.0/s); + return sqrt( -sm*s*s*sgn/(dudr-JJ*JJ*s*s*s) ); + } + return sm*cos(t)/sqrt(radcand); + }; + accum1 = 0.0; accum2 = 0.0; tl = -0.5*M_PI; @@ -333,9 +434,9 @@ void SphericalOrbit::compute_angles(void) for (int i=0; iget_pot(r); - radcand = 2.0*(EE-ur)-JJ*JJ/(r*r); - - if (radcand < tol) { - // values near turning points - if (t<0) - sgn = -1.0; - else - sgn = 1.0; - r = ap+sgn*am; - dudr = mm->get_dpot(r); - return sqrt( am*sgn/(dudr-JJ*JJ/(r*r*r)) ); - } - return am*cos(t)/sqrt(radcand); -} - -double ff(double t) -{ - double s, ur, dudr, radcand, sgn; - - s = sp + sm*sin(t); - ur = mm->get_pot(1.0/s); - radcand = 2.0*(EE-ur) - JJ*JJ*s*s; - - if (radcand < tol) { - // values near turning points - if (t<0) - sgn = -1.0; - else - sgn = 1.0; - s = sp+sgn*sm; - dudr = mm->get_dpot(1.0/s); - return sqrt( -sm*s*s*sgn/(dudr-JJ*JJ*s*s*s) ); - } - return sm*cos(t)/sqrt(radcand); -} - // // epicyclic test // void SphericalOrbit::compute_angles_epi(void) { double r,t,a,a2,fac,ur; - double Jcirc(double r); int i; angle_grid.t. resize(2, recs); @@ -435,17 +491,11 @@ void SphericalOrbit::compute_angles_epi(void) Gkn = LegeQuad(recs); } - if (freq_defined) { - mm = model; - EE = energy; - KK = kappa; - } - else - compute_freq_epi(); + if (not freq_defined) compute_freq_epi(); - JJ = jmax*kappa; - ur = mm->get_pot(r_circ); - a2 = 2.0*(EE-0.5*JJ*JJ/(r_circ*r_circ)-ur)/(freq[0]*freq[0]); + double JJ = jmax*kappa; + ur = model->get_pot(r_circ); + a2 = 2.0*(energy-0.5*JJ*JJ/(r_circ*r_circ)-ur)/(freq[0]*freq[0]); a = sqrt(a2); for (int i=0; iget_dpot(r); - return (JJ*JJ - r*r*r*dudr); -} #undef tol diff --git a/include/orbit.H b/include/orbit.H index 13af2fbbd..c99bf0fd3 100644 --- a/include/orbit.H +++ b/include/orbit.H @@ -1,6 +1,7 @@ #ifndef _orbit_H #define _orbit_H +#include #include #include #include @@ -59,6 +60,7 @@ private: double r_peri; double r_apo; double r_circ; + double ap, am, sp, sm; std::shared_ptr model; std::shared_ptr biorth; @@ -82,8 +84,13 @@ private: GaussQuad Gkn; + std::tuple search + (std::function func, double rmin, double rmax); + public: + static int Nseg; // Number of trial segments for root finding + static double RMAXF; // Factors of model rmax for upper // turning point search [default=3.0] From 0c0650a70f1bec1dcaa577514bf84a2a1ccb2ad1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 17:22:18 -0500 Subject: [PATCH 11/65] Split Mestel implentation from header definition --- exputil/CMakeLists.txt | 2 +- exputil/orbit.cc | 31 +++--- include/massmodel.H | 3 +- include/mestel.H | 213 +++++++---------------------------------- utils/ICs/ZangICs.cc | 6 +- 5 files changed, 58 insertions(+), 197 deletions(-) diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index 295836cca..87a25087a 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -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) diff --git a/exputil/orbit.cc b/exputil/orbit.cc index b80091ce4..7e5637308 100644 --- a/exputil/orbit.cc +++ b/exputil/orbit.cc @@ -47,10 +47,10 @@ SphericalOrbit::SphericalOrbit(void) { model = 0; - model_defined = false; + model_defined = false; action_defined = false; - freq_defined = false; - angle_defined = false; + freq_defined = false; + angle_defined = false; biorth_defined = false; set_numerical_params(); @@ -71,11 +71,11 @@ SphericalOrbit::SphericalOrbit(AxiSymModPtr Model) r_circ = 0.0; model = Model; - model_defined = true; + model_defined = true; action_defined = false; - freq_defined = false; - angle_defined = false; + freq_defined = false; + angle_defined = false; biorth_defined = false; set_numerical_params(); @@ -118,10 +118,10 @@ SphericalOrbit::SphericalOrbit(const SphericalOrbit &t) model = t.model; - model_defined = t.model_defined; + model_defined = t.model_defined; action_defined = t.action_defined; - freq_defined = t.freq_defined; - angle_defined = t.angle_defined; + freq_defined = t.freq_defined; + angle_defined = t.angle_defined; } SphericalOrbit &SphericalOrbit::operator=(const SphericalOrbit &t) @@ -154,8 +154,9 @@ void SphericalOrbit::new_orbit(double Energy, double Kappa, double Beta) bomb(ost.str().c_str()); } -// if (Energy < model->get_pot(model->get_min_radius()) || -// Energy > model->get_pot(model->get_max_radius())) { + // if (Energy < model->get_pot(model->get_min_radius()) || + // Energy > model->get_pot(model->get_max_radius())) { + if (guard && Energy < model->get_pot(model->get_min_radius())) { ostringstream ost; ost << "[new orbit] illegal value of Energy: " << Energy; @@ -166,12 +167,12 @@ void SphericalOrbit::new_orbit(double Energy, double Kappa, double Beta) energy = Energy; - kappa = Kappa; - beta = Beta; + kappa = Kappa; + beta = Beta; action_defined = false; - freq_defined = false; - angle_defined = false; + freq_defined = false; + angle_defined = false; biorth_defined = false; } diff --git a/include/massmodel.H b/include/massmodel.H index 4b717c30d..66c789b1c 100644 --- a/include/massmodel.H +++ b/include/massmodel.H @@ -97,7 +97,8 @@ public: int dof() { return dim; } //! Exception handling using std::exception - void bomb(const char *s) { + void bomb(const char *s) + { std::ostringstream sout; sout << "MassModel error from " << ModelID << ": " << s; throw std::runtime_error(sout.str()); diff --git a/include/mestel.H b/include/mestel.H index 4be059ef6..58bb00af4 100644 --- a/include/mestel.H +++ b/include/mestel.H @@ -19,56 +19,21 @@ protected: public: //! Constructor - MestelDisk(double VROT = 1.0, - double RMIN = 1.0e-6, double RMAX = 1.0e6) - { - vrot = VROT; - rot = VROT*VROT; - rmin = RMIN; - rmax = RMAX; - dim = 2; - - dist_defined = true; - ModelID = "MestelDisk"; - } - + MestelDisk(double VROT = 1.0, double RMIN = 1.0e-6, double RMAX = 1.0e6); //@{ //! Required member functions - virtual double get_mass(const double r) { - if (r>0.0) return rot*r; - else bomb("radius cannot be zero!"); - return 0; - } - - virtual double get_density(const double r) { - if (r>0.0) return rot/(2.0*M_PI*r); - else bomb("radius cannot be zero!"); - return 0; - } - - double get_pot(const double r) { - if (r>0.0) return rot*log(r); - else bomb("radius cannot be zero!"); - return 0; - } - - double get_dpot(const double r) { - if (r>0.0) return rot/r; - else bomb("radius cannot be zero!"); - return 0; - } - - double get_dpot2(const double r) { - if (r>0.0) return -rot/(r*r); - else bomb("radius cannot be zero!"); - return 0; - } + virtual double get_mass(const double r); + + virtual double get_density(const double r); + + double get_pot(const double r); + + double get_dpot(const double r); + + double get_dpot2(const double r); - void get_pot_dpot(const double r, double &ur, double &dur) { - if (r>0.0) {ur = rot*log(r); dur = rot/r;} - else bomb("radius cannot be zero!"); - } + void get_pot_dpot(const double r, double &ur, double &dur); //@} //@{ @@ -77,37 +42,15 @@ public: double get_min_radius(void) { return rmin; } double get_max_radius(void) { return rmax; } - virtual void setup_df(double sigma) { - sig2 = sigma*sigma; - q = rot/sig2 - 1.0; - F = rot/(4.0*M_PI) / ( sqrt(M_PI) * - exp(lgamma(0.5*(q+1.0)) + (2.0 + q)*log(sigma) + 0.5*q*log(2.0)) ); - } - - virtual double distf(double E, double L) { - L = fabs(L); - if (L==0.0) return 0.0; - return F*pow(L, q) * exp(-E/sig2); - } - - virtual double dfde(double E, double L) { - L = fabs(L); - if (L==0.0) return 0.0; - return -F*pow(L, q) * exp(-E/sig2)/sig2; - } - - virtual double dfdl(double E, double L) { - int sgn=1; - if (L<0) {sgn=-1; L *= sgn;} - if (L==0.0) return 0.0; - return q*F*pow(L, q-1.0) * exp(-E/sig2) * sgn; - } - - virtual double d2fde2(double E, double L) { - L = fabs(L); - if (L<=0.0) return 0.0; - return F*pow(L, q) * exp(-E/sig2)/sig2/sig2; - } + virtual void setup_df(double sigma); + + virtual double distf(double E, double L); + + virtual double dfde(double E, double L); + + virtual double dfdl(double E, double L); + + virtual double d2fde2(double E, double L); //@} }; @@ -131,128 +74,42 @@ protected: double Tifac; //! Inner taper function - double Tinner(double Jp) - { - double fac = pow(Jp, nu); - return fac/(Tifac + fac); - } + double Tinner(double Jp); //! Outer taper function - double Touter(double Jp) - { - return 1.0/(1.0 + pow(Jp/vrot, mu)); - } + double Touter(double Jp); //! Deriv of inner taper function - double dTinner(double Jp) - { - double fac = pow(Jp, nu); - double fac2 = Tifac + fac; - return Tifac*nu/Jp/fac2; - } + double dTinner(double Jp); //! Deriv of outer taper function - double dTouter(double Jp) - { - double fac = pow(Jp/vrot, mu); - double fac2 = 1.0 + fac; - return -nu*fac/Jp/fac2; - } - + double dTouter(double Jp); public: //! Constructor TaperedMestelDisk(double nu, double mu, double Ri, double VROT = 1.0, - double RMIN = 1.0e-6, double RMAX = 1.0e6) : - nu(nu), mu(mu), Ri(Ri), MestelDisk(VROT, RMIN, RMAX) - { - Tifac = pow(Ri*vrot, nu); - ModelID = "TaperedMestelDisk"; - } + double RMIN = 1.0e-6, double RMAX = 1.0e6); //! Overloaded density function - double get_density(const double r) { - if (r>0.0) return rot/(2.0*M_PI*r) * Tinner(r) * Touter(r); - else bomb("radius cannot be zero!"); - return 0; - } - - double get_mass(double R) - { - if (not interp) { - // Compute a grid spacing - double dr = 0.01*Ri; - int N = std::floor((rmax - rmin)/dr) + 1; - dr = (rmax - rmin)/N; - - // Storage for table - Eigen::VectorXd vecR(N+1), vecM(N+1); - - double lst = get_density(rmin); - double cum = 0.0; - - // Make the table - // - vecR(0) = rmin; - vecM(0) = 0.0; - - for (int i=1; i<=N; i++) { - double rr = rmin + dr*i; - double cur = get_density(rr); - cum += 0.5*dr*(lst + cur); - lst = cur; - - vecR(i) = rr; - vecM(i) = cum; - } - - // Make the interpolation function - // - interp = std::make_shared(vecR, vecM); - } - - return interp->eval(R); - } + double get_density(const double r); + + //! Overloaded cumulative mass function + double get_mass(double R); //@{ //! Overloaded DF member functions - double distf(double E, double L) - { - L = fabs(L); - if (L==0.0) return 0.0; - return F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2); - } - - double dfde(double E, double L) - { - L = fabs(L); - if (L==0.0) return 0.0; - return -F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2)/sig2; - } - - double dfdl(double E, double L) - { - int sgn=1; - if (L<0) {sgn=-1; L *= sgn;} - if (L==0.0) return 0.0; - double Tfac = pow(L, q)*Tinner(L)*Touter(L); - double dL = Tfac*(q/L + dTinner(L) + dTouter(L)); - return F* dL * exp(-E/sig2) * sgn; - } + double distf(double E, double L); + + double dfde(double E, double L); + + double dfdl(double E, double L); - double d2fde2(double E, double L) - { - L = fabs(L); - if (L<=0.0) return 0.0; - return F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2)/sig2/sig2; - } + double d2fde2(double E, double L); //@} }; - - #endif diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index 9be26f849..bf35d5576 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -99,7 +99,9 @@ main(int ac, char **av) // Create the model // - auto model = std::make_shared(nu, mu, Ri, 1.0, 0.1*Rmin, 10.0*Rmax); + auto model = std::make_shared(nu, mu, Ri, 1.0, + 0.1*Rmin, 10.0*Rmax); + model->setup_df(sigma); // Progress bar // @@ -149,7 +151,7 @@ main(int ac, char **av) // Save the position and velocity vectors std::vector> pos(N), vel(N); - int itmax = 10000; + int itmax = 100000; int tid = 0; int over = 0; From acbea9cbb352f636120bb7f7792828100c4ecdd2 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 17:34:10 -0500 Subject: [PATCH 12/65] A bit of clean up --- exputil/mestel.cc | 204 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 204 insertions(+) create mode 100644 exputil/mestel.cc diff --git a/exputil/mestel.cc b/exputil/mestel.cc new file mode 100644 index 000000000..5256abad3 --- /dev/null +++ b/exputil/mestel.cc @@ -0,0 +1,204 @@ +#include + +// Implementation + +MestelDisk::MestelDisk(double VROT, double RMIN, double RMAX) +{ + vrot = VROT; + rot = VROT*VROT; + rmin = RMIN; + rmax = RMAX; + dim = 2; + + setup_df(1.0); // Default to sigma=1.0 + ModelID = "MestelDisk"; +} + + +double MestelDisk::get_mass(const double r) +{ + if (r>0.0) return rot*r; + else bomb("radius cannot be zero!"); + return 0; +} + +double MestelDisk::get_density(const double r) +{ + if (r>0.0) return rot/(2.0*M_PI*r); + else bomb("radius cannot be zero!"); + return 0; +} + +double MestelDisk::get_pot(const double r) +{ + if (r>0.0) return rot*log(r); + else bomb("radius cannot be zero!"); + return 0; +} + +double MestelDisk::get_dpot(const double r) +{ + if (r>0.0) return rot/r; + else bomb("radius cannot be zero!"); + return 0; +} + +double MestelDisk::get_dpot2(const double r) +{ + if (r>0.0) return -rot/(r*r); + else bomb("radius cannot be zero!"); + return 0; +} + +void MestelDisk::get_pot_dpot(const double r, double &ur, double &dur) +{ + if (r>0.0) {ur = rot*log(r); dur = rot/r;} + else bomb("radius cannot be zero!"); +} + +void MestelDisk::setup_df(double sigma) +{ + sig2 = sigma*sigma; + q = rot/sig2 - 1.0; + F = rot/(4.0*M_PI) / + ( sqrt(M_PI) * + exp(lgamma(0.5*(q+1.0)) + (2.0 + q)*log(sigma) + 0.5*q*log(2.0)) ); + // + dist_defined = true; +} + +double MestelDisk::distf(double E, double L) +{ + L = fabs(L); + if (L==0.0) return 0.0; + return F*pow(L, q) * exp(-E/sig2); +} + +double MestelDisk::dfde(double E, double L) +{ + L = fabs(L); + if (L==0.0) return 0.0; + return -F*pow(L, q) * exp(-E/sig2)/sig2; +} + +double MestelDisk::dfdl(double E, double L) +{ + int sgn = 1; + if (L<0) {sgn=-1; L *= sgn;} + if (L==0.0) return 0.0; + return q*F*pow(L, q-1.0) * exp(-E/sig2) * sgn; +} + +double MestelDisk::d2fde2(double E, double L) +{ + L = fabs(L); + if (L<=0.0) return 0.0; + return F*pow(L, q) * exp(-E/sig2)/sig2/sig2; +} + +double TaperedMestelDisk::Tinner(double Jp) +{ + double fac = pow(Jp, nu); + return fac/(Tifac + fac); +} + +double TaperedMestelDisk::Touter(double Jp) +{ + return 1.0/(1.0 + pow(Jp/vrot, mu)); +} + +double TaperedMestelDisk::dTinner(double Jp) +{ + double fac = pow(Jp, nu); + double fac2 = Tifac + fac; + return Tifac*nu/Jp/fac2; +} + +double TaperedMestelDisk::dTouter(double Jp) +{ + double fac = pow(Jp/vrot, mu); + double fac2 = 1.0 + fac; + return -nu*fac/Jp/fac2; +} + +TaperedMestelDisk::TaperedMestelDisk(double nu, double mu, double Ri, + double VROT, double RMIN, double RMAX) + : nu(nu), mu(mu), Ri(Ri), MestelDisk(VROT, RMIN, RMAX) +{ + Tifac = pow(Ri*vrot, nu); + ModelID = "TaperedMestelDisk"; +} + +double TaperedMestelDisk::get_density(const double r) { + if (r>0.0) return rot/(2.0*M_PI*r) * Tinner(r) * Touter(r); + else bomb("radius cannot be zero!"); + return 0; +} + +double TaperedMestelDisk::get_mass(double R) +{ + if (not interp) { + // Compute a grid spacing + double dr = 0.01*Ri; + int N = std::floor((rmax - rmin)/dr) + 1; + dr = (rmax - rmin)/N; + + // Storage for table + Eigen::VectorXd vecR(N+1), vecM(N+1); + + double lst = get_density(rmin); + double cum = 0.0; + + // Make the table + // + vecR(0) = rmin; + vecM(0) = 0.0; + + for (int i=1; i<=N; i++) { + double rr = rmin + dr*i; + double cur = get_density(rr); + cum += 0.5*dr*(lst + cur); + lst = cur; + + vecR(i) = rr; + vecM(i) = cum; + } + + // Make the interpolation function + // + interp = std::make_shared(vecR, vecM); + } + + return interp->eval(R); +} + +double TaperedMestelDisk::distf(double E, double L) +{ + L = fabs(L); + if (L==0.0) return 0.0; + return F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2); +} + +double TaperedMestelDisk::dfde(double E, double L) +{ + L = fabs(L); + if (L==0.0) return 0.0; + return -F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2)/sig2; +} + +double TaperedMestelDisk::dfdl(double E, double L) +{ + int sgn=1; + if (L<0) {sgn=-1; L *= sgn;} + if (L==0.0) return 0.0; + double Tfac = pow(L, q)*Tinner(L)*Touter(L); + double dL = Tfac*(q/L + dTinner(L) + dTouter(L)); + return F* dL * exp(-E/sig2) * sgn; +} + +double TaperedMestelDisk::d2fde2(double E, double L) +{ + L = fabs(L); + if (L<=0.0) return 0.0; + return F*pow(L, q)* Tinner(L)*Touter(L) * exp(-E/sig2)/sig2/sig2; +} From dc7e9c5a69264ddc2eb9401f11a2c1d893e6ee46 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 17:58:11 -0500 Subject: [PATCH 13/65] Fixed a typo in tangential velocity; seems to make an equilibrium model now --- utils/ICs/ZangICs.cc | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index bf35d5576..947de5386 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -52,11 +52,11 @@ main(int ac, char **av) ("R,Rmax", "Outer radius for model", cxxopts::value(Rmax)->default_value("10.0")) ("S,sigma", "Radial velocity dispersion", - cxxopts::value(sigma)) + cxxopts::value(sigma)->default_value("1.0")) ("s,seed", "Random number seed. Default: use /dev/random", cxxopts::value(seed)) ("o,file", "Output body file", - cxxopts::value(bodyfile)->default_value("cube.bods")) + cxxopts::value(bodyfile)->default_value("zang.bods")) ; cxxopts::ParseResult vm; @@ -135,7 +135,7 @@ main(int ac, char **av) for (int j=0; j<=num; j++) { double K = Kmin + dK*j; orb[0]->new_orbit(E, K); - double F = model->distf(E, orb[0]->get_action(2)) / orb[0]->get_freq(1); + double F = model->distf(E, orb[0]->get_action(1)) / orb[0]->get_freq(0); peak = std::max(peak, F); } } @@ -171,20 +171,23 @@ main(int ac, char **av) R = uniform(gen); orb[tid]->new_orbit(E, K); - double F = model->distf(E, orb[tid]->get_action(2)) / orb[tid]->get_freq(1); + double F = model->distf(E, orb[tid]->get_action(1)) / orb[tid]->get_freq(0); if (F/peak > R) break; } if (j==itmax) over++; - double J = orb[tid]->get_action(2); - double T = 2.0*M_PI/orb[tid]->get_freq(1)*uniform(gen); + double J = orb[tid]->get_action(1); + double T = 2.0*M_PI/orb[tid]->get_freq(0)*uniform(gen); double r = orb[tid]->get_angle(6, T); double w1 = orb[tid]->get_angle(1, T); double phi = 2.0*M_PI*uniform(gen) + orb[tid]->get_angle(7, T); double vt = J/r; - double vr = sqrt(fabs(2.0*E - model->get_pot(r)) - J*J/(r*r)); + double vr = sqrt(fabs(2.0*(E - model->get_pot(r)) - J*J/(r*r))); + // ^ + // | + // For sanity-----+ if (w1 > M_PI) vr *= -1.0; // Branch of radial motion @@ -231,7 +234,7 @@ main(int ac, char **av) clausius += mass*model->get_dpot(r)*r; } - std::cout << "2T/VC=" << ektot/clausius << std::endl; + std::cout << "** 2T/VC=" << ektot/clausius << std::endl; return 0; } From 53ca45ba831f60a9b657f6155f1bf0b1083a5c0a Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 Feb 2024 18:16:04 -0500 Subject: [PATCH 14/65] Fixed up progress bar and a few output issues; more comments --- utils/ICs/ZangICs.cc | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index 947de5386..a17505dea 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -110,8 +110,9 @@ main(int ac, char **av) #pragma omp parallel { nomp = omp_get_num_threads(); - if (omp_get_thread_num()==0) + if (omp_get_thread_num()==0) { progress = std::make_shared(N/nomp); + } } // Create an orbit grid @@ -149,21 +150,25 @@ main(int ac, char **av) std::uniform_real_distribution<> uniform(0.0, 1.0); // Save the position and velocity vectors + // std::vector> pos(N), vel(N); + // Maximum number rejection-method iterations int itmax = 100000; - int tid = 0; + + // Track number of iteration overflows int over = 0; - // Generation loop + // Generation loop with OpenMP // #pragma omp parallel for reduction(+:over) for (int n=0; nget_mass(Rmax) - model->get_mass(Rmin))/N; - std::cout << "** " << over << " particles failed iteration" << std::endl + std::cout << "** " << over << " particles failed acceptance" << std::endl << "** Particle mass=" << mass << std::endl; out << std::setw(8) << N << std::setw(8) << 0 << std::setw(8) << 0 From 2ad932e72bc782aa1402e9f24d919a484ae1a2d6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 10 Feb 2024 13:06:47 -0500 Subject: [PATCH 15/65] Changed the FlatDisk API so that disk-model configuration is a YAML map of names and parameters rather than a string name and fixed array of doubles --- exputil/BiorthCyl.cc | 18 +-- exputil/EmpCyl2d.cc | 284 +++++++++++++++++++++++++++++++---------- exputil/orbit_trans.cc | 6 +- include/BiorthCyl.H | 2 +- include/EmpCyl2d.H | 20 +-- src/FlatDisk.cc | 4 +- utils/ICs/ZangICs.cc | 54 ++++++-- utils/SL/EOF2d.cc | 12 +- 8 files changed, 294 insertions(+), 106 deletions(-) diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index e4305e9bf..d8ec02b37 100644 --- a/exputil/BiorthCyl.cc +++ b/exputil/BiorthCyl.cc @@ -87,9 +87,6 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) if (conf["Mouter"]) Mouter = conf["Mouter"].as(); else Mouter = 4.0; - if (conf["disktype"]) disktype = conf["disktype"].as(); - else disktype = "expon"; - if (conf["cachename"]) cachename = conf["cachename"].as(); else cachename = default_cache; @@ -113,6 +110,9 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) if (conf["EVEN_M"]) EVEN_M = conf["EVEN_M"].as(); 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: " @@ -175,11 +175,9 @@ void BiorthCyl::create_tables() bool logr = false; std::string biorth("bess"); - EmpCyl2d::ParamVec par {acyltbl*scale, acylcut, Ninner, Mouter}; - emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, - rcylmin*scale, rcylmax*scale, par, scale, - cmapR, logr, disktype, biorth); + rcylmin*scale, rcylmax*scale, + scale, cmapR, logr, diskconf, biorth); if (conf["basis"]) emp.basisTest(true); @@ -813,11 +811,9 @@ BiorthCyl::cacheInfo(const std::string& cachefile, bool verbose) std::vector BiorthCyl::orthoCheck() { if (not emp()) { - EmpCyl2d::ParamVec par {acyltbl*scale, acyltbl*scale*10.0, 2.0, 2.0}; - emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, - rcylmin*scale, rcylmax*scale, par, scale, - cmapR, false, "expon", "bess"); + rcylmin*scale, rcylmax*scale, scale, + cmapR, false, diskconf, "bess"); } return emp.orthoCheck(); diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 0ef3379c9..4d8bff39e 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -338,10 +338,38 @@ class EmpCyl2d::ExponCyl : public EmpCyl2d::ModelCyl double acyl, sigma0; + // Assign values from YAML + // + void parse(const YAML::Node& conf) + { + try { + if (conf["acyl"]) + acyl = conf["acyl"].as(); + else + acyl = 1.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::ExponCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + MPI_Finalize(); + exit(-1); + } + } + public: - ExponCyl(const EmpCyl2d::ParamVec& par) - { acyl = par[0]; sigma0 = 0.5/(M_PI*acyl*acyl); id = "expon"; } + ExponCyl(const YAML::Node& par) + { + parse(par); + sigma0 = 0.5/(M_PI*acyl*acyl); + id = "expon"; + } double pot(double r) { double y = 0.5 * r / acyl; @@ -369,10 +397,37 @@ class EmpCyl2d::KuzminCyl : public EmpCyl2d::ModelCyl double acyl; + // Assign values from YAML + // + void parse(const YAML::Node& conf) + { + try { + if (conf["acyl"]) + acyl = conf["acyl"].as(); + else + acyl = 1.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::KuzminCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + MPI_Finalize(); + exit(-1); + } + } + public: - KuzminCyl(const EmpCyl2d::ParamVec& par) - { acyl = par[0]; id = "kuzmin"; } + KuzminCyl(const YAML::Node& par) + { + parse(par); + id = "kuzmin"; + } double pot(double R) { double a2 = acyl * acyl; @@ -397,14 +452,42 @@ class EmpCyl2d::KuzminCyl : public EmpCyl2d::ModelCyl class EmpCyl2d::MestelCyl : public EmpCyl2d::ModelCyl { -private: +protected: double vrot, rot; + // Assign values from YAML + // + void parse(const YAML::Node& conf) + { + try { + if (conf["vrot"]) + vrot = conf["vrot"].as(); + else + vrot = 1.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::MestelCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + MPI_Finalize(); + exit(-1); + } + } + public: - MestelCyl(const EmpCyl2d::ParamVec& par) - { vrot = par[0]; rot = vrot*vrot; id = "mestel"; } + MestelCyl(const YAML::Node& par) + { + parse(par); + rot = vrot*vrot; + id = "mestel"; + } virtual double pot(double R) { if (R>0.0) return rot*log(R); @@ -468,18 +551,57 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl return -nu*fac/Jp/fac2; } +protected: + + //! Assign values from YAML + void parse(const YAML::Node& conf) + { + try { + if (conf["vrot"]) + vrot = conf["vrot"].as(); + else + vrot = 1.0; + + if (conf["Ninner"]) + nu = conf["Ninner"].as(); + else + nu = 2.0; + + if (conf["Mouter"]) + mu = conf["Mouter"].as(); + else + mu = 2.0; + + if (conf["Ri"]) + ri = conf["Ri"].as(); + else + ri = 1.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::ZangCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + MPI_Finalize(); + exit(-1); + } + } public: //! Constructor - ZangCyl(const EmpCyl2d::ParamVec& par) : MestelCyl(par) + ZangCyl(const YAML::Node& par) : MestelCyl(par) { + // Parse the YAML + parse(par); + // Assign the id id = "zang"; - vr = par[0]; - nu = par[1]; - mu = par[2]; - ri = par[3]; + // Cache taper factor Tifac = pow(ri*vr, nu); if (nu<0.05) { @@ -499,47 +621,76 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl }; std::shared_ptr -EmpCyl2d::createModel(const std::string type, const EmpCyl2d::ParamVec& par) +EmpCyl2d::createModel() { - // Convert ID string to lower case - // - std::string data(type); - std::transform(data.begin(), data.end(), data.begin(), - [](unsigned char c){ return std::tolower(c); }); + try { + // Get model name + // + std::string name; + if (Params["name"]) + name = Params["name"].as(); + else + throw std::runtime_error("EmpCyl2d::createModel: the 'diskconf' stanza must specify the model 'name'"); - if (data.find("kuzmin") != std::string::npos) { - if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making a Kuzmin disk with A=" - << par[0] << std::endl; - return std::make_shared(par); + // Convert ID string to lower case + // + std::transform(name.begin(), name.end(), name.begin(), + [](unsigned char c){ return std::tolower(c); }); + + if (name.find("kuzmin") != std::string::npos) { + if (myid==0) { + std::cout << "---- EmpCyl2d::ModelCyl: Making a Kuzmin disk"; + if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + std::cout << std::endl; + } + return std::make_shared(Params["parameters"]); } - if (data.find("mestel") != std::string::npos) { - if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making a finite Mestel disk " - << "with A=" << par[0] << std::endl; - return std::make_shared(par); + if (name.find("mestel") != std::string::npos) { + if (myid==0) { + std::cout << "---- EmpCyl2d::ModelCyl: Making a finite Mestel disk"; + if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + std::cout << std::endl; + } + return std::make_shared(Params["parameters"]); } - if (data.find("zang") != std::string::npos) { - if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making a double-tapered Zang " - "disk with A=" << par[0] << std::endl; - return std::make_shared(par); + if (name.find("zang") != std::string::npos) { + if (myid==0) { + std::cout << "---- EmpCyl2d::ModelCyl: Making a double-tapered Zang"; + if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + std::cout << std::endl; + } + return std::make_shared(Params["parameters"]); } - if (data.find("expon") != std::string::npos) { - if (myid==0) - std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk " - << "with A=" << par[0] << std::endl; - return std::make_shared(par); + if (name.find("expon") != std::string::npos) { + if (myid==0) { + std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk"; + if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + std::cout << std::endl; + } + return std::make_shared(Params["parameters"]); } // Default if nothing else matches if (myid==0) std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk [Default]" - << " with A=" << par[0] << std::endl; - return std::make_shared(par); + << std::endl; + return std::make_shared(Params["parameters"]); + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::modelCreate: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << Params << std::endl + << std::string(60, '-') << std::endl; + MPI_Finalize(); + exit(-1); + } } @@ -607,18 +758,18 @@ const std::string EmpCyl2d::default_cache_name = ".eof_cache_2d"; // The main constructor // -EmpCyl2d::EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, - double rmin, double rmax, const EmpCyl2d::ParamVec& par, - double scale, bool cmap, bool logr, - const std::string model, const std::string biorth, - const std::string cache) : +EmpCyl2d::EmpCyl2d +(int mmax, int nmaxfid, int nmax, int knots, int numr, + double rmin, double rmax, double scale, bool cmap, bool logr, + const YAML::Node& par, + const std::string biorth, const std::string cache) : mmax(mmax), nmaxfid(nmaxfid), nmax(nmax), knots(knots), numr(numr), - rmin(rmin), rmax(rmax), par(par), scale(scale), cmap(cmap), logr(logr), - model(model), biorth(biorth), cache_name_2d(cache) + rmin(rmin), rmax(rmax), scale(scale), cmap(cmap), logr(logr), + Params(par), biorth(biorth), cache_name_2d(cache) { if (cache_name_2d.size()==0) cache_name_2d = default_cache_name; - disk = createModel(model, par); + disk = createModel(); basis = Basis2d::createBasis(mmax, nmaxfid, rmax, biorth); basis_test = false; @@ -631,13 +782,13 @@ EmpCyl2d::EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, } -EmpCyl2d::EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, - double rmin, double rmax, const EmpCyl2d::ParamVec& par, - double scale, bool cmap, bool logr, - std::shared_ptr disk, - const std::string biorth, const std::string cache) : +EmpCyl2d::EmpCyl2d +(int mmax, int nmaxfid, int nmax, int knots, int numr, + double rmin, double rmax, double scale, bool cmap, bool logr, + std::shared_ptr disk, + const std::string biorth, const std::string cache) : mmax(mmax), nmaxfid(nmaxfid), nmax(nmax), knots(knots), numr(numr), - rmin(rmin), rmax(rmax), par(par), scale(scale), cmap(cmap), logr(logr), + rmin(rmin), rmax(rmax), scale(scale), cmap(cmap), logr(logr), disk(disk), biorth(biorth), cache_name_2d(cache) { if (cache_name_2d.size()==0) cache_name_2d = default_cache_name; @@ -854,6 +1005,10 @@ void EmpCyl2d::WriteH5Cache() if (logr) ilogr = 1; if (cmap) icmap = 1; + // Serialize the config and make a string + YAML::Emitter y; y << Params; + std::string params(y.c_str()); + // Parameters // file.createAttribute ("mmax", HighFive::DataSpace::From(mmax)). write(mmax); @@ -866,7 +1021,7 @@ void EmpCyl2d::WriteH5Cache() file.createAttribute ("rmin", HighFive::DataSpace::From(rmin)). write(rmin); file.createAttribute ("rmax", HighFive::DataSpace::From(rmax)). write(rmax); file.createAttribute ("scale", HighFive::DataSpace::From(scale)). write(scale); - file.createAttribute ("params", HighFive::DataSpace::From(par)).write(par); + file.createAttribute ("params", HighFive::DataSpace::From(params)).write(params); file.createAttribute("model", HighFive::DataSpace::From(model)). write(model); file.createAttribute("biorth", HighFive::DataSpace::From(biorth)).write(biorth); @@ -920,21 +1075,18 @@ bool EmpCyl2d::ReadH5Cache() if (fabs(value - v) < 1.0e-16) return true; return false; }; - auto checkArr = [&file](ParamVec& value, std::string name) - { - ParamVec v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); - for (int i=0; i 1.0e-16) return false; - } - return true; - }; - 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; }; + // + + // Serialize the config and make a string for checking + YAML::Emitter y; y << Params; + std::string params(y.c_str()); + // Workaround for lack of HighFive boolean support int ilogr = 0, icmap = 0; if (logr) ilogr = 1; @@ -950,7 +1102,7 @@ bool EmpCyl2d::ReadH5Cache() if (not checkDbl(rmin, "rmin")) return false; if (not checkDbl(rmax, "rmax")) return false; if (not checkDbl(scale, "scale")) return false; - if (not checkArr(par, "params")) return false; + if (not checkStr(params, "params")) return false; if (not checkStr(model, "model")) return false; if (not checkStr(biorth, "biorth")) return false; @@ -1124,7 +1276,7 @@ void EmpCyl2d::checkCoefs() if (myid) return; Mapping map(scale, cmap); - auto disk = createModel(model, par); + auto disk = createModel(); LegeQuad lw(knots); Eigen::VectorXd coefs(nmax), coef0(nmax); diff --git a/exputil/orbit_trans.cc b/exputil/orbit_trans.cc index 71e52bb79..35edbd32f 100644 --- a/exputil/orbit_trans.cc +++ b/exputil/orbit_trans.cc @@ -112,8 +112,10 @@ void SphericalOrbit::compute_freq(void) if (r_circ < model->get_min_radius()) { r_circ = model->get_min_radius(); } - std::string message("SphericalOrbit::compute_freq warning: Circular radius outside mass distribution"); - throw std::runtime_error(message); + /* + std::string message("SphericalOrbit::compute_freq warning: Circular radius outside mass distribution"); + throw std::runtime_error(message); + */ } } diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index dd5144acb..110cc053f 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -41,7 +41,7 @@ class BiorthCyl protected: - YAML::Node conf; + YAML::Node conf, diskconf; std::string geometry, forceID; diff --git a/include/EmpCyl2d.H b/include/EmpCyl2d.H index e3d39761e..2549e29a9 100644 --- a/include/EmpCyl2d.H +++ b/include/EmpCyl2d.H @@ -12,6 +12,7 @@ #include #include +#include // Parameters parsing #include // For config macros // For reading and writing cache file @@ -28,15 +29,13 @@ class EmpCyl2d { public: - //! Datatype for passing model parameters - using ParamVec = std::array; - //! A two-dimensional disk model for computing the EOF class ModelCyl { protected: std::string id; + virtual void parse(const YAML::Node&) = 0; public: @@ -49,7 +48,9 @@ public: protected: - ParamVec par; + //! Contains parameter database + YAML::Node Params; + double rmin, rmax, scale; double xmin, xmax, dxi; int mmax, nmaxfid, numr, knots, nmax; @@ -87,8 +88,7 @@ protected: std::shared_ptr disk; std::shared_ptr basis; - static std::shared_ptr - createModel(const std::string type, const ParamVec& params); + std::shared_ptr createModel(); class ExponCyl; class MestelCyl; @@ -149,14 +149,14 @@ public: //! Constructor EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, - double rmin, double rmax, const ParamVec& P, - double scale, bool cmap, bool logr, const std::string model, + double rmin, double rmax, double scale, bool cmap, bool logr, + const YAML::Node& P, const std::string biorth, const std::string cache=""); //! Constructor with user-supplied target model EmpCyl2d(int mmax, int nmaxfid, int nmax, int knots, int numr, - double rmin, double rmax, const ParamVec& P, - double scale, bool cmap, bool logr, std::shared_ptr disk, + double rmin, double rmax, double scale, bool cmap, bool logr, + std::shared_ptr disk, const std::string biorth, const std::string cache=""); //! Use underlying basis for testing? Default: false. diff --git a/src/FlatDisk.cc b/src/FlatDisk.cc index 4a727ebe8..3a238c10f 100644 --- a/src/FlatDisk.cc +++ b/src/FlatDisk.cc @@ -15,7 +15,7 @@ FlatDisk::valid_keys = { "acyltbl", "acylcut", "Ninner", - "Nouter", + "Mouter", "mmax", "numx", "numy", @@ -23,7 +23,7 @@ FlatDisk::valid_keys = { "logr", "model", "biorth", - "disktype", + "diskconf", "cachename" }; diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index a17505dea..f83ec0540 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -39,6 +39,7 @@ main(int ac, char **av) options.add_options() ("h,help", "Print this help message") ("z,zerovel", "Zero the mean velocity") + ("d,debug", "Print debug grid") ("N,number", "Number of particles to generate", cxxopts::value(N)->default_value("100000")) ("n,nu", "Inner taper exponent (0 for no taper)", @@ -128,19 +129,43 @@ main(int ac, char **av) // Scan to find the peak df // - const int num = 100; + const int numE = 400; + const int numK = 20; + Eigen::VectorXd cumE(numE+1), cumF(numE+1), topF(numE+1); double peak = 0.0; - double dE = (Emax - Emin)/num, dK = (1.0 - 2.0*Ktol)/num; - for (int i=0; i<=num; i++) { + double dE = (Emax - Emin)/numE, dK = (1.0 - 2.0*Ktol)/numK; + for (int i=0; i<=numE; i++) { double E = Emin + dE*i; - for (int j=0; j<=num; j++) { + cumE(i) = E; + cumF(i) = 0.0; + topF(i) = 0.0; + for (int j=0; j<=numK; j++) { double K = Kmin + dK*j; orb[0]->new_orbit(E, K); double F = model->distf(E, orb[0]->get_action(1)) / orb[0]->get_freq(0); peak = std::max(peak, F); + topF(i) = std::max(topF(i), F); + cumF(i) += F * orb[0]->Jmax()/orb[0]->get_freq(0); } } + if (peak <= 0.0) { + throw std::runtime_error(std::string(av[0]) + ": peak DF is zero!"); + } + + // Improve the acceptance rejection by computing the cumulative + // distribution + // + for (int i=1; i<=numE; i++) cumF[i] += cumF[i-1]; + + if (cumF[numE] <= 0.0) { + throw std::runtime_error(std::string(av[0]) + ": no mass on cum DF grid!"); + } + + for (int i=0; i<=numE; i++) cumF[i] /= cumF[numE]; + + Linear1d Ecum(cumF, cumE), Ftop(cumE, topF); + // Header // out << std::setw(10) << N << std::setw(10) << 0 << std::setw(10) << 0 @@ -167,17 +192,17 @@ main(int ac, char **av) int tid = omp_get_thread_num(); // Loop variables - double E, K, R; + double E, F, K; int j; for (j=0; jnew_orbit(E, K); double F = model->distf(E, orb[tid]->get_action(1)) / orb[tid]->get_freq(0); - if (F/peak > R) break; + if (F/peak > uniform(gen)) break; } if (j==itmax) over++; @@ -241,5 +266,18 @@ main(int ac, char **av) std::cout << "** 2T/VC=" << ektot/clausius << std::endl; + if (vm.count("debug")) { + std::cout << std::endl << "Peak per energy" << std::endl; + for (int i=0; i<=numE; i++) { + std::cout << std::setw(6) << i + << std::setw(18) << cumE(i) + << std::setw(18) << cumF(i) + << std::setw(18) << topF(i) + << std::endl; + } + std::cout << std::endl; + } + + return 0; } diff --git a/utils/SL/EOF2d.cc b/utils/SL/EOF2d.cc index 8db73dc11..42b5f8377 100644 --- a/utils/SL/EOF2d.cc +++ b/utils/SL/EOF2d.cc @@ -16,7 +16,7 @@ int main(int argc, char** argv) bool logr = false, cmap = false, ortho = false, plane = false; int numr, mmax, nmaxfid, nmax, knots, M, N, nradial, nout; double A, rmin, rmax, O, ZN, ZM, rout; - std::string filename, type, biorth; + std::string filename, config, biorth; // Parse command line // @@ -74,8 +74,8 @@ int main(int argc, char** argv) cxxopts::value(rout)->default_value("10.0")) ("nout", "number of points in the output grid per side", cxxopts::value(nout)->default_value("40")) - ("type", "Target model type (kuzmin, mestel, zang, expon)", - cxxopts::value(type)->default_value("expon")) + ("config", "Target model config (kuzmin, mestel, zang, expon)", + cxxopts::value(config)->default_value("{name: expon, parameters: {acyl: 0.01}}")) ("biorth", "Biorthogonal type (cb, bess)", cxxopts::value(biorth)->default_value("bess")) ("o,filename", "Output filename", @@ -113,12 +113,12 @@ int main(int argc, char** argv) // Parameter vector for the EmpCyl2d models // - EmpCyl2d::ParamVec par {A, O, ZN, ZM}; + YAML::Node par(config); // Make the class instance // - EmpCyl2d emp(mmax, nmaxfid, nmax, knots, numr, rmin, rmax, par, A, cmap, logr, - type, biorth); + EmpCyl2d emp(mmax, nmaxfid, nmax, knots, numr, rmin, rmax, A, cmap, logr, + par, biorth); if (vm.count("basis")) emp.basisTest(true); if (vm.count("debug")) QDHT::debug = true; From 7d9a293d6bfa36dd6933250e4a7f3022c1692602 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 10 Feb 2024 22:50:51 -0500 Subject: [PATCH 16/65] Correction to HighFive write of params string [no ci] --- exputil/EmpCyl2d.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 4d8bff39e..3e4f41b8e 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -1021,7 +1021,7 @@ void EmpCyl2d::WriteH5Cache() file.createAttribute ("rmin", HighFive::DataSpace::From(rmin)). write(rmin); file.createAttribute ("rmax", HighFive::DataSpace::From(rmax)). write(rmax); file.createAttribute ("scale", HighFive::DataSpace::From(scale)). write(scale); - file.createAttribute ("params", HighFive::DataSpace::From(params)).write(params); + file.createAttribute("params", HighFive::DataSpace::From(params)).write(params); file.createAttribute("model", HighFive::DataSpace::From(model)). write(model); file.createAttribute("biorth", HighFive::DataSpace::From(biorth)).write(biorth); From 9b68ac92bdd74c6cc759e45151838cf8b72e960d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 11 Feb 2024 15:01:36 -0500 Subject: [PATCH 17/65] Return force not potential gradient [no ci] --- src/FlatDisk.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FlatDisk.H b/src/FlatDisk.H index 27c1d9c5d..a3e1fe182 100644 --- a/src/FlatDisk.H +++ b/src/FlatDisk.H @@ -77,7 +77,7 @@ private: //! Background evaluation virtual void get_pot_background(double r, double z, double & p, double & dr, double & dz) - { ortho->background(r, z, p, dr, dz); } + { ortho->background(r, z, p, dr, dz); dr *= -1.0; dz *= -1.0; } void get_dens(double r, double z, Eigen::MatrixXd& p, int tid); From ea2fb4f688e285c2e0bdc76d2320c2a65a7dedc4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 11 Feb 2024 15:02:14 -0500 Subject: [PATCH 18/65] Fix prefactor for force from background [no ci] --- src/PolarBasis.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index 2c0f09827..384697626 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -1425,6 +1425,10 @@ void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) if (M0_back) { get_pot_background(r, zz, p, drc, dzc); + + potl = mfac * p; + potr = mfac * drc; + potz = mfac * dzc; } else { get_pot_coefs_safe(0, *expcoef[0], p, drc, dzc, potd[id], dpotR[id], dpotZ[id]); From 9cbc467e72e56ca734c9e59e139ff42c1bc17d74 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 11 Feb 2024 15:02:50 -0500 Subject: [PATCH 19/65] Edit out double header in body output [no ci] --- utils/ICs/ZangICs.cc | 5 ----- 1 file changed, 5 deletions(-) diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index f83ec0540..fd29cc9c1 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -166,11 +166,6 @@ main(int ac, char **av) Linear1d Ecum(cumF, cumE), Ftop(cumE, topF); - // Header - // - out << std::setw(10) << N << std::setw(10) << 0 << std::setw(10) << 0 - << std::endl << std::setprecision(10); - std::mt19937 gen(seed); std::uniform_real_distribution<> uniform(0.0, 1.0); From f3b5134b6244fcdbddb06a5ebd120bd4d4d50258 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 11 Feb 2024 15:03:33 -0500 Subject: [PATCH 20/65] Remove unused parameters and serialize YAML in HDF5 [no ci] --- exputil/BiorthCyl.cc | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index d8ec02b37..5b5a8db37 100644 --- a/exputil/BiorthCyl.cc +++ b/exputil/BiorthCyl.cc @@ -81,12 +81,6 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) if (conf["acylcut"]) acylcut = conf["acylcut"].as(); else acylcut = 10.0; - if (conf["Ninner"]) Ninner = conf["Ninner"].as(); - else Ninner = 2.0; - - if (conf["Mouter"]) Mouter = conf["Mouter"].as(); - else Mouter = 4.0; - if (conf["cachename"]) cachename = conf["cachename"].as(); else cachename = default_cache; @@ -489,8 +483,6 @@ void BiorthCyl::WriteH5Params(HighFive::File& file) file.createAttribute ("scale", HighFive::DataSpace::From(scale)).write(scale); file.createAttribute ("acyltbl", HighFive::DataSpace::From(acyltbl)).write(acyltbl); file.createAttribute ("acylcut", HighFive::DataSpace::From(acylcut)).write(acylcut); - file.createAttribute ("Ninner", HighFive::DataSpace::From(Ninner)).write(Ninner); - file.createAttribute ("Mouter", HighFive::DataSpace::From(Mouter)).write(Mouter); file.createAttribute ("cmapR", HighFive::DataSpace::From(cmapR)).write(cmapR); file.createAttribute ("cmapZ", HighFive::DataSpace::From(cmapZ)).write(cmapZ); } @@ -663,6 +655,19 @@ bool BiorthCyl::ReadH5Cache() std::cerr << "---- BiorthCyl::ReadH5Cache: " << "read <" << cachename << ">" << std::endl; + // Hardwired for testing. Could be variables. + // + bool logr = false; + std::string biorth("bess"); + + // 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) { From fa663099986a6ec5784b14146c8d716c5656f5e9 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 11 Feb 2024 15:04:10 -0500 Subject: [PATCH 21/65] White-space change only [no ci] --- exputil/EmpCyl2d.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 3e4f41b8e..621339aee 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -620,8 +620,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl }; -std::shared_ptr -EmpCyl2d::createModel() +std::shared_ptr EmpCyl2d::createModel() { try { // Get model name From 0f3b64e77cb4b0f306930837bf4777e320c26e65 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 11 Feb 2024 15:13:55 -0500 Subject: [PATCH 22/65] Added spacer to log output [no ci] --- exputil/BiorthCyl.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index 5b5a8db37..2a06104a4 100644 --- a/exputil/BiorthCyl.cc +++ b/exputil/BiorthCyl.cc @@ -603,7 +603,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; From 9d9948ebf02b1a67a7e68580bb42382b8c20f9d5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 11 Feb 2024 17:34:56 -0500 Subject: [PATCH 23/65] Compute grid counters independent of VERBOSE level [no ci] --- src/multistep.cc | 59 ++++++++++++++++++++++-------------------------- 1 file changed, 27 insertions(+), 32 deletions(-) diff --git a/src/multistep.cc b/src/multistep.cc index d24af5119..c68f840fd 100644 --- a/src/multistep.cc +++ b/src/multistep.cc @@ -35,16 +35,16 @@ struct thrd_pass_sync // // Count offgrid particles in the threads // -static map< Component*, vector > offlo1, offhi1; -static vector< double > mindt1; -static vector< double > maxdt1; -static vector< double > adjtm1; -static vector< double > adjtm2; -static vector< unsigned> numsw; -static vector< unsigned> numtt; +static map< Component*, std::vector > offlo1, offhi1; +static std::vector< double > mindt1; +static std::vector< double > maxdt1; +static std::vector< double > adjtm1; +static std::vector< double > adjtm2; +static std::vector< unsigned> numsw; +static std::vector< unsigned> numtt; // Type counter -static vector< vector< vector > > tmdt; +static std::vector< std::vector< std::vector > > tmdt; // // The threaded routine @@ -225,10 +225,8 @@ void * adjust_multistep_level_thread(void *ptr) } } - if (VERBOSE>0) { - offlo1[c][id] += offlo; - offhi1[c][id] += offhi; - } + offlo1[c][id] += offlo; + offhi1[c][id] += offhi; finish0 = std::chrono::high_resolution_clock::now(); std::chrono::duration duration = finish0 - start0; @@ -388,22 +386,19 @@ void adjust_multistep_level() // // Preliminary data structure and thread creation // - mindt1 = vector< double > (nthrds, 1.0e20); - maxdt1 = vector< double > (nthrds, -1.0e20); - adjtm1 = vector< double > (nthrds, 0.0); - adjtm2 = vector< double > (nthrds, 0.0); - numsw = vector< unsigned > (nthrds, 0); - numtt = vector< unsigned > (nthrds, 0); - - if (VERBOSE>0) { + mindt1 = std::vector< double > (nthrds, 1.0e20); + maxdt1 = std::vector< double > (nthrds, -1.0e20); + adjtm1 = std::vector< double > (nthrds, 0.0); + adjtm2 = std::vector< double > (nthrds, 0.0); + numsw = std::vector< unsigned > (nthrds, 0); + numtt = std::vector< unsigned > (nthrds, 0); - if (offhi1.size()==0 || mdrft==Mstep) { + if (offhi1.size()==0 || mdrft==Mstep) { - for (auto c : comp->components) { - for (int n=0; n(nthrds, 0); - offlo1[c] = vector(nthrds, 0); - } + for (auto c : comp->components) { + for (int n=0; n(nthrds, 0); + offlo1[c] = std::vector(nthrds, 0); } } } @@ -428,10 +423,10 @@ void adjust_multistep_level() if (tmdt.size() == 0) { - tmdt = vector< vector< vector > >(nthrds); + tmdt = std::vector< std::vector< std::vector > >(nthrds); for (int n=0; n >(multistep+1); - for (int k=0; k<=multistep; k++) tmdt[n][k] = vector(mdtDim); + tmdt[n] = std::vector< std::vector >(multistep+1); + for (int k=0; k<=multistep; k++) tmdt[n][k] = std::vector(mdtDim); } } @@ -641,9 +636,9 @@ void initialize_multistep() for (int n=1; n<=multistep; n++) mintvl.push_back(mintvl.back()/2); // Set up the step-level bool array - mactive.push_back(vector(multistep+1, true)); + mactive.push_back(std::vector(multistep+1, true)); for (int ms=1; ms<=Mstep; ms++) - mactive.push_back(vector(multistep+1, false)); + mactive.push_back(std::vector(multistep+1, false)); // Find and save the active levels at each step for (int ms=1; ms<=Mstep; ms++) { @@ -653,7 +648,7 @@ void initialize_multistep() } } - mfirst = vector(Mstep+1); + mfirst = std::vector(Mstep+1); // Lowest active level at each step for (int ms=0; ms<=Mstep; ms++) { for (int mlevel=0; mlevel<=multistep; mlevel++) { From 62b7a75059f0e079a03b2b8f8aebbc23e46e6a7a Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 12 Feb 2024 18:54:32 -0500 Subject: [PATCH 24/65] Update tapered Mestel disk parameters for consistency with the literature [no ci] --- exputil/EmpCyl2d.cc | 16 +++++++++++----- exputil/SLGridMP2.cc | 14 +++++++++++--- exputil/mestel.cc | 7 ++++--- include/mestel.H | 13 ++++++++----- utils/ICs/ZangICs.cc | 36 +++++++++++++++++++++++++++++------- 5 files changed, 63 insertions(+), 23 deletions(-) diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 621339aee..3557160dc 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -511,7 +511,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl private: //! Parameters - double vr, mu, nu, ri; + double vr, mu, nu, ri, ro; //! Softening factor double asoft = 1.0e-8; @@ -519,8 +519,8 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl //! Ignore inner cut-off for N<0.05 bool Inner = true; - //! Taper factor - double Tifac; + //! Taper factors + double Tifac, Tofac; //! Inner taper function double Tinner(double Jp) @@ -546,7 +546,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl //! Deriv of outer taper function double dTouter(double Jp) { - double fac = pow(Jp/vr, mu); + double fac = pow(Jp/Tofac, mu); double fac2 = 1.0 + fac; return -nu*fac/Jp/fac2; } @@ -576,6 +576,11 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl ri = conf["Ri"].as(); else ri = 1.0; + + if (conf["Ro"]) + ro = conf["Ro"].as(); + else + ro = 10.0; } catch (YAML::Exception & error) { if (myid==0) @@ -601,8 +606,9 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl // Assign the id id = "zang"; - // Cache taper factor + // Cache taper factors Tifac = pow(ri*vr, nu); + Tofac = ro*vr; if (nu<0.05) { // Exponent is now for mapping only diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 37ed197f9..299716660 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -2258,10 +2258,17 @@ bool SLGridSph::ReadH5Cache(void) if (not checkInt(cmap, "cmap")) return false; if (not checkDbl(rmin, "rmin")) return false; if (not checkDbl(rmax, "rmax")) return false; - if (not checkDbl(rmap, "rmapping")) return false; if (not checkInt(diverge, "diverge")) return false; if (not checkDbl(dfac, "dfac")) return false; + // Backward compatibility for old 'scale' key word + // + if (h5file.hasAttribute("scale")) { + if (not checkDbl(rmap, "scale")) return false; + } else { + if (not checkDbl(rmap, "rmapping")) return false; + } + // Harmonic order // auto harmonic = h5file.getGroup("Harmonic"); @@ -2291,8 +2298,9 @@ bool SLGridSph::ReadH5Cache(void) } catch (HighFive::Exception& err) { if (myid==0) std::cerr << "---- SLGridSph::ReadH5Cache: " - << "error opening <" << sph_cache_name - << "> as HDF5 basis cache" << std::endl; + << "error reading <" << sph_cache_name << ">" << std::endl + << "---- SLGridSph::ReadH5Cache: HDF5 error is <" << err.what() + << ">" << std::endl; } return false; diff --git a/exputil/mestel.cc b/exputil/mestel.cc index 5256abad3..6fc6d80eb 100644 --- a/exputil/mestel.cc +++ b/exputil/mestel.cc @@ -104,7 +104,7 @@ double TaperedMestelDisk::Tinner(double Jp) double TaperedMestelDisk::Touter(double Jp) { - return 1.0/(1.0 + pow(Jp/vrot, mu)); + return 1.0/(1.0 + pow(Jp/Tofac, mu)); } double TaperedMestelDisk::dTinner(double Jp) @@ -121,11 +121,12 @@ double TaperedMestelDisk::dTouter(double Jp) return -nu*fac/Jp/fac2; } -TaperedMestelDisk::TaperedMestelDisk(double nu, double mu, double Ri, +TaperedMestelDisk::TaperedMestelDisk(double nu, double mu, double Ri, double Ro, double VROT, double RMIN, double RMAX) - : nu(nu), mu(mu), Ri(Ri), MestelDisk(VROT, RMIN, RMAX) + : nu(nu), mu(mu), Ri(Ri), Ro(Ro), MestelDisk(VROT, RMIN, RMAX) { Tifac = pow(Ri*vrot, nu); + Tofac = Ro*vrot; ModelID = "TaperedMestelDisk"; } diff --git a/include/mestel.H b/include/mestel.H index 58bb00af4..41f896ae7 100644 --- a/include/mestel.H +++ b/include/mestel.H @@ -64,14 +64,17 @@ private: std::shared_ptr interp; protected: - //! Taper parameters + //! Taper exponents double nu, mu; - //! Inner radius + //! Inner taper radius double Ri; - //! Taper factor - double Tifac; + //! Outer taper radius + double Ro; + + //! Taper factors + double Tifac, Tofac; //! Inner taper function double Tinner(double Jp); @@ -88,7 +91,7 @@ protected: public: //! Constructor - TaperedMestelDisk(double nu, double mu, double Ri, + TaperedMestelDisk(double nu, double mu, double Ri, double Ro, double VROT = 1.0, double RMIN = 1.0e-6, double RMAX = 1.0e6); diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index fd29cc9c1..dc17d5ca4 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -27,7 +27,7 @@ main(int ac, char **av) //===================== int N; // Number of particles - double mu, nu, Ri; // Taper paramters + double mu, nu, Ri, Ro; // Taper paramters double Rmin, Rmax; // Radial range double sigma; // Velocity dispersion std::string bodyfile; // Output file @@ -47,16 +47,18 @@ main(int ac, char **av) ("m,mu", "Outer taper exponent (0 for no taper)", cxxopts::value(mu)->default_value("2.0")) ("i,Ri", "Inner radius for taper", - cxxopts::value(Ri)->default_value("0.1")) + cxxopts::value(Ri)->default_value("1.0")) + ("o,Ro", "Outer radius for taper", + cxxopts::value(Ro)->default_value("10.0")) ("r,Rmin", "Inner radius for model", cxxopts::value(Rmin)->default_value("0.01")) ("R,Rmax", "Outer radius for model", - cxxopts::value(Rmax)->default_value("10.0")) + cxxopts::value(Rmax)->default_value("100.0")) ("S,sigma", "Radial velocity dispersion", cxxopts::value(sigma)->default_value("1.0")) ("s,seed", "Random number seed. Default: use /dev/random", cxxopts::value(seed)) - ("o,file", "Output body file", + ("f,file", "Output body file", cxxopts::value(bodyfile)->default_value("zang.bods")) ; @@ -100,8 +102,8 @@ main(int ac, char **av) // Create the model // - auto model = std::make_shared(nu, mu, Ri, 1.0, - 0.1*Rmin, 10.0*Rmax); + auto model = std::make_shared(nu, mu, Ri, Ro, + 1.0, 0.1*Rmin, 10.0*Rmax); model->setup_df(sigma); // Progress bar @@ -179,6 +181,10 @@ main(int ac, char **av) // Track number of iteration overflows int over = 0; + // Velocity zeroing + // + std::vector> zerovel(nomp); + // Generation loop with OpenMP // #pragma omp parallel for reduction(+:over) @@ -226,6 +232,10 @@ main(int ac, char **av) vel[n][1] = vr*sin(phi) + vt*cos(phi); vel[n][2] = 0.0; + // Accumulate mean velocity + // + for (int k=0; k<3; k++) zerovel[tid][k] += vel[n][k]; + // Print progress bar if (tid==0) ++(*progress); } @@ -235,6 +245,18 @@ main(int ac, char **av) // double mass = (model->get_mass(Rmax) - model->get_mass(Rmin))/N; + // Reduce the mean velocity + // + for (int n=1; n Date: Tue, 13 Feb 2024 11:44:18 -0500 Subject: [PATCH 25/65] Fixed mistake in IC generator; removed unused parameters from 2d basis parsing and cache [no ic] --- coefs/BiorthBasis.cc | 2 ++ exputil/BiorthCyl.cc | 27 ++++++--------------------- exputil/mestel.cc | 2 +- include/BiorthCyl.H | 4 ++-- include/mestel.H | 14 +++++++------- 5 files changed, 18 insertions(+), 31 deletions(-) diff --git a/coefs/BiorthBasis.cc b/coefs/BiorthBasis.cc index 18427eb08..d72cbb85c 100644 --- a/coefs/BiorthBasis.cc +++ b/coefs/BiorthBasis.cc @@ -1405,7 +1405,9 @@ namespace BasisClasses "NO_M0", "NO_M1", "EVEN_M", + "M0_BACK", "M0_ONLY", + "diskconf", "ssfrac", "playback", "coefMaster", diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index 2a06104a4..3bba20384 100644 --- a/exputil/BiorthCyl.cc +++ b/exputil/BiorthCyl.cc @@ -75,12 +75,6 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) if (conf["scale"]) scale = conf["scale"].as(); else scale = 0.01; - if (conf["acyltbl"]) acyltbl = conf["acyltbl"].as(); - else acyltbl = 0.6; - - if (conf["acylcut"]) acylcut = conf["acylcut"].as(); - else acylcut = 10.0; - if (conf["cachename"]) cachename = conf["cachename"].as(); else cachename = default_cache; @@ -122,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(); @@ -164,11 +163,6 @@ void BiorthCyl::initialize() void BiorthCyl::create_tables() { - // Hardwired for testing. Could be variables. - // - bool logr = false; - std::string biorth("bess"); - emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, rcylmin*scale, rcylmax*scale, scale, cmapR, logr, diskconf, biorth); @@ -481,8 +475,6 @@ void BiorthCyl::WriteH5Params(HighFive::File& file) file.createAttribute ("rcylmin", HighFive::DataSpace::From(rcylmin)).write(rcylmin); file.createAttribute ("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax); file.createAttribute ("scale", HighFive::DataSpace::From(scale)).write(scale); - file.createAttribute ("acyltbl", HighFive::DataSpace::From(acyltbl)).write(acyltbl); - file.createAttribute ("acylcut", HighFive::DataSpace::From(acylcut)).write(acylcut); file.createAttribute ("cmapR", HighFive::DataSpace::From(cmapR)).write(cmapR); file.createAttribute ("cmapZ", HighFive::DataSpace::From(cmapZ)).write(cmapZ); } @@ -641,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; @@ -655,11 +646,6 @@ bool BiorthCyl::ReadH5Cache() std::cerr << "---- BiorthCyl::ReadH5Cache: " << "read <" << cachename << ">" << std::endl; - // Hardwired for testing. Could be variables. - // - bool logr = false; - std::string biorth("bess"); - // Create the basis instance // emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, @@ -774,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"); } @@ -818,7 +803,7 @@ std::vector BiorthCyl::orthoCheck() if (not emp()) { emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, rcylmin*scale, rcylmax*scale, scale, - cmapR, false, diskconf, "bess"); + cmapR, logr, diskconf, biorth); } return emp.orthoCheck(); diff --git a/exputil/mestel.cc b/exputil/mestel.cc index 6fc6d80eb..d13207ef9 100644 --- a/exputil/mestel.cc +++ b/exputil/mestel.cc @@ -158,7 +158,7 @@ double TaperedMestelDisk::get_mass(double R) for (int i=1; i<=N; i++) { double rr = rmin + dr*i; double cur = get_density(rr); - cum += 0.5*dr*(lst + cur); + cum += 0.5*dr*(lst + cur) * 2.0*M_PI*rr; lst = cur; vecR(i) = rr; diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index 110cc053f..d767ab42e 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -43,12 +43,12 @@ protected: YAML::Node conf, diskconf; - std::string geometry, forceID; + std::string geometry, biorth, forceID; int mmax, nmax, numr, nmaxfid, mmin, mlim, nmin, nlim, knots, NQDHT; double rcylmin, rcylmax, scale, acyltbl, acylcut, Ninner, Mouter; - bool EVEN_M, verbose; + bool EVEN_M, verbose, logr; //@{ //! Grid parameters diff --git a/include/mestel.H b/include/mestel.H index 41f896ae7..fa5541a22 100644 --- a/include/mestel.H +++ b/include/mestel.H @@ -23,17 +23,17 @@ public: //@{ //! Required member functions - virtual double get_mass(const double r); + virtual double get_mass(const double R); - virtual double get_density(const double r); + virtual double get_density(const double R); - double get_pot(const double r); + double get_pot(const double R); - double get_dpot(const double r); + double get_dpot(const double R); - double get_dpot2(const double r); + double get_dpot2(const double R); - void get_pot_dpot(const double r, double &ur, double &dur); + void get_pot_dpot(const double R, double &ur, double &dur); //@} //@{ @@ -96,7 +96,7 @@ public: double RMIN = 1.0e-6, double RMAX = 1.0e6); //! Overloaded density function - double get_density(const double r); + double get_density(const double R); //! Overloaded cumulative mass function double get_mass(double R); From c82ef1e585870c75cc45f6e47003ef755b6edef1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 13 Feb 2024 18:13:27 -0500 Subject: [PATCH 26/65] Enforce a specified cachename parameter [no ci] --- exputil/BiorthCyl.cc | 2 +- exputil/EmpCylSL.cc | 6 ++++-- exputil/SLGridMP2.cc | 4 ++-- include/BiorthCyl.H | 3 --- include/EmpCylSL.H | 3 --- src/Cylinder.cc | 6 ++---- src/FlatDisk.cc | 7 ------- src/PolarBasis.cc | 2 +- src/Sphere.cc | 8 ++++---- 9 files changed, 14 insertions(+), 27 deletions(-) diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index 3bba20384..f4f2abd1b 100644 --- a/exputil/BiorthCyl.cc +++ b/exputil/BiorthCyl.cc @@ -76,7 +76,7 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) else scale = 0.01; if (conf["cachename"]) cachename = conf["cachename"].as(); - else cachename = default_cache; + else throw std::runtime_error("BiorthCyl: you must specify a cachename"); // Add output directory and runtag cachename = outdir + cachename; diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 55fa4a0c0..90840b972 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -99,7 +99,6 @@ EmpCylSL::EmpCylSL() Nodd = 0; minSNR = std::numeric_limits::max(); maxSNR = 0.0; - cachefile = default_cache; // Initial allocation; set to false // @@ -142,7 +141,7 @@ EmpCylSL::EmpCylSL(int nmax, int lmax, int mmax, int nord, { // Use default name? if (cachename.size()) cachefile = cachename; - else cachefile = default_cache; + else throw std::runtime_error("EmpCylSL: you must specify a cachename"); // Sanity check if (lmax <= mmax) { @@ -227,6 +226,7 @@ void EmpCylSL::reset(int numr, int lmax, int mmax, int nord, { // Reset cache file name if (cachename.size()) cachefile = cachename; + else throw std::runtime_error("EmpCylSL:reset: you must specify a cachename"); // Option sanity check if (lmax <= mmax) { @@ -825,6 +825,8 @@ int EmpCylSL::cache_grid(int readwrite, string cachename) // Option to reset cache file name // if (cachename.size()) cachefile = cachename; + else throw std::runtime_error("EmpCylSL::cache_grid: you must specify a cachename"); + if (readwrite) { diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 299716660..c8b5e5d88 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1884,7 +1884,7 @@ SLGridSph::SLGridSph(std::string modelname, else model_file_name = default_model; if (cachename.size()) sph_cache_name = cachename; - else sph_cache_name = default_cache; + else throw std::runtime_error("SLGridSph: you must specify a cachename"); mpi_buf = 0; model = SphModTblPtr(new SphericalModelTable(model_file_name, DIVERGE, DFAC)); @@ -1907,7 +1907,7 @@ SLGridSph::SLGridSph(std::shared_ptr mod, dfac = 1; if (cachename.size()) sph_cache_name = cachename; - else sph_cache_name = default_cache; + else throw std::runtime_error("SLGridSph: you must specify a cachename"); initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, RMAP); } diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index d767ab42e..e6912e840 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -79,9 +79,6 @@ protected: const std::vector>& mat, Eigen::MatrixXd& ret, bool anti_symmetric=false); - //! Default cache file name - const std::string default_cache = ".biorth_cache"; - //! Density target name std::string disktype; diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 9c33136ff..55fa787df 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -225,9 +225,6 @@ protected: void setup_eof_grid(void); void parityCheck(const std::string& prefix); - //! The default cache file name - const std::string default_cache = ".eof.cache.file"; - //! Write HDF5 cache void WriteH5Cache(); diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 6fd7b6126..d68bb447f 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -201,10 +201,8 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; - // EOF default file name override. Default uses runtag suffix as - // above. Override file must exist if explicitly specified. - // - if (cachename.size()==0) cachename = outdir + ".eof.cache." + runtag; + if (cachename.size()==0) + throw std::runtime_error("EmpCylSL: you must specify a cachename"); // Make the empirical orthogonal basis instance // diff --git a/src/FlatDisk.cc b/src/FlatDisk.cc index 3a238c10f..2d2dc7345 100644 --- a/src/FlatDisk.cc +++ b/src/FlatDisk.cc @@ -12,10 +12,6 @@ FlatDisk::valid_keys = { "numr", "rcylmin", "rcylmax", - "acyltbl", - "acylcut", - "Ninner", - "Mouter", "mmax", "numx", "numy", @@ -34,7 +30,6 @@ FlatDisk::FlatDisk(Component* c0, const YAML::Node& conf, MixtureBasis* m) : // Defaults knots = 40; numr = 2000; - acyltbl = 1.0; rcylmin = 0.0; rcylmax = 10.0; logr = false; @@ -50,7 +45,6 @@ FlatDisk::FlatDisk(Component* c0, const YAML::Node& conf, MixtureBasis* m) : std::cout << "---- FlatDisk parameters: " << std::endl << sep << "lmax=" << Lmax << std::endl << sep << "nmax=" << nmax - << std::endl << sep << "acyltbl=" << acyltbl << std::endl << sep << "rcylmin=" << rcylmin << std::endl << sep << "rcylmax=" << rcylmax << std::endl << sep << "logr=" << std::boolalpha << logr @@ -76,7 +70,6 @@ void FlatDisk::initialize() // Assign values from YAML // try { - if (conf["acyltbl"]) acyltbl = conf["acyltbl"].as(); if (conf["rcylmin"]) rcylmin = conf["rcylmin"].as(); if (conf["rcylmax"]) rcylmax = conf["rcylmax"].as(); if (conf["mmax"]) mmax = conf["mmax"].as(); diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index 384697626..f060f0a91 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -1445,7 +1445,7 @@ void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) // m loop // ------ - for (int m=1, moffset=1; m<=Mmax; m++, moffset+=2) { + for (int m=1, moffset=1; m<=std::min(mlim, Mmax); m++, moffset+=2) { // Skip m=1 terms? // diff --git a/src/Sphere.cc b/src/Sphere.cc index 51ffe5fa1..bd81633c3 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -36,7 +36,6 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : diverge = 0; dfac = 1.0; model_file = "SLGridSph.model"; - cache_file = "SLGridSph.cache"; tnext = dtime = 0.0; recompute = false; plummer = true; @@ -54,7 +53,7 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : if (numprocs>1) SLGridSph::mpi = 1; std::string modelname = model_file; - std::string cachename = outdir + cache_file + "." + runtag; + std::string cachename = outdir + cache_file; id += ", model=" + modelname; @@ -114,6 +113,7 @@ void Sphere::initialize() if (conf["dfac"]) dfac = conf["dfac"].as(); if (conf["modelname"]) model_file = conf["modelname"].as(); if (conf["cachename"]) cache_file = conf["cachename"].as(); + else throw std::runtime_error("Sphere: you must specify a cachname"); if (conf["dtime"]) dtime = conf["dtime"].as(); if (conf["logr"]) logr = conf["logr"].as(); if (conf["plummer"]) plummer = conf["plummer"].as(); @@ -324,7 +324,7 @@ void Sphere::make_model_bin() // Regenerate Sturm-Liouville grid // - std::string cachename = outdir + cache_file + "." + runtag; + std::string cachename = outdir + cache_file; ortho = std::make_shared(mod, Lmax, nmax, numR, Rmin, Rmax, false, 1, 1.0, cachename); // Test for basis consistency (will generate an exception if maximum @@ -470,7 +470,7 @@ void Sphere::make_model_plummer() // Regenerate Sturm-Liouville grid // - std::string cachename = outdir + cache_file + "." + runtag; + std::string cachename = outdir + cache_file; ortho = std::make_shared(mod, Lmax, nmax, numr, Rmin, Rmax, false, 1, 1.0, cachename); // Test for basis consistency (will generate an exception if maximum From 2a5ae647a34a025261048221534920ffc3f772bc Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 15 Feb 2024 13:22:44 -0500 Subject: [PATCH 27/65] Added orthgonality check [no ci] --- utils/SL/EOF2d.cc | 116 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 107 insertions(+), 9 deletions(-) diff --git a/utils/SL/EOF2d.cc b/utils/SL/EOF2d.cc index 42b5f8377..9db5a1e35 100644 --- a/utils/SL/EOF2d.cc +++ b/utils/SL/EOF2d.cc @@ -6,16 +6,18 @@ #include #include +#include + #include // Hankel computation for potential #include // 2d empirical basis +#include // Gauss-Legendre quadrature #include - int main(int argc, char** argv) { bool logr = false, cmap = false, ortho = false, plane = false; int numr, mmax, nmaxfid, nmax, knots, M, N, nradial, nout; - double A, rmin, rmax, O, ZN, ZM, rout; + double A, rmin, rmax, rout; std::string filename, config, biorth; // Parse command line @@ -42,6 +44,7 @@ int main(int argc, char** argv) ("debug", "Check unitarity in QDHT") ("full", "Use full transform rather than grid evaluation") ("totforce", "Compute the total radial force") + ("bartest", "Check expansion against a known quadrupole") ("M,harmonic", "Aximuthal harmonic m=0,1,2,3,...", cxxopts::value(M)->default_value("0")) ("N,nsize", "Default radial grid size", @@ -50,12 +53,6 @@ int main(int argc, char** argv) cxxopts::value(nradial)->default_value("0")) ("A,length", "characteristic disk scale length", cxxopts::value(A)->default_value("1.0")) - ("O,outer", "outer disk truncation", - cxxopts::value(O)->default_value("10.0")) - ("1,TaperN", "Inner taper exponent", - cxxopts::value(ZN)->default_value("1.0")) - ("2,TaperM", "Outer taper exponent", - cxxopts::value(ZM)->default_value("4.0")) ("mmax", "maximum number of angular harmonics in the expansion", cxxopts::value(mmax)->default_value("4")) ("nmaxfid", "maximum number of radial basis harmonics for EOF construction", @@ -113,7 +110,7 @@ int main(int argc, char** argv) // Parameter vector for the EmpCyl2d models // - YAML::Node par(config); + YAML::Node par = YAML::Load(config); // Make the class instance // @@ -235,6 +232,8 @@ int main(int argc, char** argv) } } } + // END: vertical + if (vm.count("plane")) { @@ -305,6 +304,7 @@ int main(int argc, char** argv) } } } + // END: plane if (vm.count("Sk")) { @@ -339,6 +339,104 @@ int main(int argc, char** argv) out << std::endl; } } + // END: Sk + + if (vm.count("bartest")) { + double a = 1.0, b = 2.0*A; + + auto pot = [a, b](double r) + { + return a*r*r*pow(1.0 + r/b, -5.0); + }; + + auto rho = [a, b](double r) + { + double x = r/b; + return -6.0*a*(1.0 - 3.0*x + x*x)*pow(1.0 + x, -7.0); + }; + + auto r_to_x = [b](double r) + { + double y = r/b; + return y/(1.0 + y); + }; + + auto x_to_r = [b](double x) + { + return b*x/(1.0 - x); + }; + + auto drdx = [b](double x) + { + double d = 1.0 - x; + return b/(d*d); + }; + + // Let's set up for quadrature + // + double xmin = r_to_x(rmin), xmax = r_to_x(rmax); + + const int num = 400; + LegeQuad lq(num); + + // TEST 1: compute inner product against potential + // TEST 2: compute inner product against density + // + Eigen::VectorXd coef1(nmax), coef2(nmax); + Eigen::MatrixXd ortho(nmax, nmax); + + coef1.setZero(); + coef2.setZero(); + ortho.setZero(); + + for (int i=0; i Date: Thu, 15 Feb 2024 13:25:35 -0500 Subject: [PATCH 28/65] Fixed derivatives of taper functions (not currently used but there was a mistake nonetheless) [no ci] --- exputil/EmpCyl2d.cc | 23 ++++++++++------------- exputil/mestel.cc | 6 +++--- include/EmpCyl2d.H | 7 +++++++ 3 files changed, 20 insertions(+), 16 deletions(-) diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 3557160dc..2e54e3f66 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -511,7 +511,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl private: //! Parameters - double vr, mu, nu, ri, ro; + double mu, nu, ri, ro; //! Softening factor double asoft = 1.0e-8; @@ -532,7 +532,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl //! Outer taper function double Touter(double Jp) { - return 1.0/(1.0 + pow(Jp/vr, mu)); + return 1.0/(1.0 + pow(Jp/Tofac, mu)); } //! Deriv of inner taper function @@ -540,7 +540,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl { double fac = pow(Jp, nu); double fac2 = Tifac + fac; - return Tifac*nu/Jp/fac2; + return nu*fac/Jp/(fac2*fac2); } //! Deriv of outer taper function @@ -548,7 +548,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl { double fac = pow(Jp/Tofac, mu); double fac2 = 1.0 + fac; - return -nu*fac/Jp/fac2; + return -mu*fac/Jp/(fac2*fac2); } protected: @@ -557,11 +557,6 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl void parse(const YAML::Node& conf) { try { - if (conf["vrot"]) - vrot = conf["vrot"].as(); - else - vrot = 1.0; - if (conf["Ninner"]) nu = conf["Ninner"].as(); else @@ -603,12 +598,13 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl { // Parse the YAML parse(par); + // Assign the id id = "zang"; // Cache taper factors - Tifac = pow(ri*vr, nu); - Tofac = ro*vr; + Tifac = pow(ri*vrot, nu); + Tofac = ro*vrot; if (nu<0.05) { // Exponent is now for mapping only @@ -635,7 +631,8 @@ std::shared_ptr EmpCyl2d::createModel() if (Params["name"]) name = Params["name"].as(); else - throw std::runtime_error("EmpCyl2d::createModel: the 'diskconf' stanza must specify the model 'name'"); + throw std::runtime_error("EmpCyl2d::createModel: the 'diskconf' " + "stanza must specify the model 'name'"); // Convert ID string to lower case // @@ -906,7 +903,7 @@ void EmpCyl2d::writeBasis(int M, const std::string& filename) if (out) { out << std::endl << "# EOF basis grid for M=" << M << ":" << std::endl; - for (int i=0; ipot(r); dr = disk->dpot(r); } + + std::tuple background(double r) + { + return {disk->pot(r), disk->dpot(r), disk->dens(r)}; + } + //@} }; From 6bc779dfaefd649416d1a9249706823b04944bac Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 15 Feb 2024 17:22:13 -0500 Subject: [PATCH 29/65] Replace Lmax by Mmax in header signatures for clarity [no ci] --- src/PolarBasis.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/PolarBasis.H b/src/PolarBasis.H index 55fe546ed..73ac4f432 100644 --- a/src/PolarBasis.H +++ b/src/PolarBasis.H @@ -266,9 +266,9 @@ protected: std::vector> T_covr; - void resize_coefs(int nmax, int Lmax, int N, int gridSize, int stride, + void resize_coefs(int nmax, int Mmax, int N, int gridSize, int stride, int sampT, bool pca, bool pcaeof, bool subsamp); - void resize_acc (int Lmax, int Nthread); + void resize_acc (int Mmax, int Nthread); }; double get_coef(int m, int n, char c) From b8c26f8a0b49976f818ac0df6b3d3b64b6fb7b55 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 15 Feb 2024 17:22:36 -0500 Subject: [PATCH 30/65] White space and comments only [no ci] --- src/PolarBasis.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index f060f0a91..9677d5cdb 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -35,6 +35,7 @@ PolarBasis::valid_keys = { "scale", "rmin", "rmax", + "Mmax", "self_consistent", "NO_M0", "NO_M1", @@ -108,12 +109,16 @@ PolarBasis::PolarBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); if (conf["M0_BACK"]) M0_back = conf["M0_BACK"].as(); + if (conf["ssfrac"]) { ssfrac = conf["ssfrac"].as(); // Check for sane value if (ssfrac>0.0 && ssfrac<1.0) subset = true; } + // AxisymmetricBasis only knows about Lmax and Mmax. These + // make sure that they are equivalent. Mmax is preferred. + // if (conf["Lmax"] and not conf["Mmax"]) Mmax = Lmax; if (conf["Mmax"] and not conf["Lmax"]) Lmax = Mmax; if (Lmax != Mmax) Lmax = Mmax; From 55b212c07d7abaa85073098dd3eecb0c27232d1a Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 15 Feb 2024 17:24:02 -0500 Subject: [PATCH 31/65] Simplify logic for substep skip [no ci] --- src/OutCoef.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutCoef.cc b/src/OutCoef.cc index 07c9a1a23..4522fabab 100644 --- a/src/OutCoef.cc +++ b/src/OutCoef.cc @@ -101,7 +101,7 @@ void OutCoef::Run(int n, int mstep, bool last) // Skip this sub step // - if (mstep < std::numeric_limits::max() and mstep % nintsub != 0) return; + if (mstep % nintsub != 0) return; // Check for repeat time // From 8ffef0985c97048062a110f855938112e6722fe5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 15 Feb 2024 17:24:26 -0500 Subject: [PATCH 32/65] Comment change only [no ci] --- src/AxisymmetricBasis.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AxisymmetricBasis.H b/src/AxisymmetricBasis.H index a1b8138d1..ce456d554 100644 --- a/src/AxisymmetricBasis.H +++ b/src/AxisymmetricBasis.H @@ -9,7 +9,7 @@ @param Lmax is the maximum spherical harmonic order - @param Mmax is the maximum spherical harmonic order + @param Mmax is the maximum cylindrical harmonic order @param nmax is the maximum radial order From 10347258d7fc841a264c4761ca199f61b61f6f10 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 15 Feb 2024 17:25:12 -0500 Subject: [PATCH 33/65] Allow Lmax and Mmax as proxies for mmax for consistency with PolarBasis in the N-body code [no ci] --- coefs/BiorthBasis.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/coefs/BiorthBasis.cc b/coefs/BiorthBasis.cc index d72cbb85c..8e357d77d 100644 --- a/coefs/BiorthBasis.cc +++ b/coefs/BiorthBasis.cc @@ -1458,6 +1458,8 @@ namespace BasisClasses // try { if (conf["cmap"]) cmap = conf["cmap"].as(); + if (conf["Lmax"]) mmax = conf["Lmax"].as(); // Proxy + if (conf["Mmax"]) mmax = conf["Mmax"].as(); // Proxy if (conf["mmax"]) mmax = conf["mmax"].as(); if (conf["nmax"]) nmax = conf["nmax"].as(); From d1d435188740b70fd3715ec65882ececfcba5c9e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 16 Feb 2024 11:43:54 -0500 Subject: [PATCH 34/65] Comment change only [no ci] --- exputil/EmpCyl2d.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 2e54e3f66..9247bd0d3 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -513,7 +513,7 @@ class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl //! Parameters double mu, nu, ri, ro; - //! Softening factor + //! Softening factor (not currently used) double asoft = 1.0e-8; //! Ignore inner cut-off for N<0.05 From af43ca0584eca26e802b2366c3dfbb736ba8ac3f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 16 Feb 2024 11:45:16 -0500 Subject: [PATCH 35/65] Update energy grid defaults for DF [no ci] --- utils/ICs/ZangICs.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index dc17d5ca4..75b92219b 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -49,11 +49,11 @@ main(int ac, char **av) ("i,Ri", "Inner radius for taper", cxxopts::value(Ri)->default_value("1.0")) ("o,Ro", "Outer radius for taper", - cxxopts::value(Ro)->default_value("10.0")) + cxxopts::value(Ro)->default_value("20.0")) ("r,Rmin", "Inner radius for model", - cxxopts::value(Rmin)->default_value("0.01")) + cxxopts::value(Rmin)->default_value("0.001")) ("R,Rmax", "Outer radius for model", - cxxopts::value(Rmax)->default_value("100.0")) + cxxopts::value(Rmax)->default_value("50.0")) ("S,sigma", "Radial velocity dispersion", cxxopts::value(sigma)->default_value("1.0")) ("s,seed", "Random number seed. Default: use /dev/random", @@ -103,7 +103,7 @@ main(int ac, char **av) // Create the model // auto model = std::make_shared(nu, mu, Ri, Ro, - 1.0, 0.1*Rmin, 10.0*Rmax); + 1.0, Rmin, Rmax); model->setup_df(sigma); // Progress bar @@ -131,8 +131,8 @@ main(int ac, char **av) // Scan to find the peak df // - const int numE = 400; - const int numK = 20; + const int numE = 800; + const int numK = 40; Eigen::VectorXd cumE(numE+1), cumF(numE+1), topF(numE+1); double peak = 0.0; double dE = (Emax - Emin)/numE, dK = (1.0 - 2.0*Ktol)/numK; From f28404695a582c316b099fa08c0e0d286f3f6235 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 16 Feb 2024 11:45:50 -0500 Subject: [PATCH 36/65] Use the C++17 structured binding rather than reference return for background() function [no ci] --- include/BiorthCyl.H | 7 +++---- src/FlatDisk.H | 9 ++++++--- src/PolarBasis.H | 5 ++--- src/PolarBasis.cc | 4 ++-- 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index e6912e840..da3cd51f5 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -199,11 +199,10 @@ public: void get_pot(Eigen::MatrixXd& Vc, Eigen::MatrixXd& Vs, double r, double z); //! Background evaluation - virtual void background(double r, double z, - double & p, double & dr, double & dz) + virtual std::tuple background(double r, double z) { - emp.background(r, p, dr); - dz = 0.0; + auto [p, dr, d] = emp.background(r); + return {p, dr, 0.0}; } //! Get the table bounds diff --git a/src/FlatDisk.H b/src/FlatDisk.H index a3e1fe182..64be763a4 100644 --- a/src/FlatDisk.H +++ b/src/FlatDisk.H @@ -75,9 +75,12 @@ private: #endif //! Background evaluation - virtual void get_pot_background(double r, double z, - double & p, double & dr, double & dz) - { ortho->background(r, z, p, dr, dz); dr *= -1.0; dz *= -1.0; } + virtual std::tuple + get_pot_background(double r, double z) + { + auto [p, dr, dz] = ortho->background(r, z); + return {p, -dr, -dz}; + } void get_dens(double r, double z, Eigen::MatrixXd& p, int tid); diff --git a/src/PolarBasis.H b/src/PolarBasis.H index 73ac4f432..ae43d4af8 100644 --- a/src/PolarBasis.H +++ b/src/PolarBasis.H @@ -378,9 +378,8 @@ protected: #endif //! Background evaluation - virtual void get_pot_background(double r, double z, - double & p, double & dr, double & dz) - { p = dr = dz = 0.0; } + virtual std::tuple + get_pot_background(double r, double z) { return {0.0, 0.0, 0.0}; } //! Thread method for accerlation compuation virtual void * determine_acceleration_and_potential_thread(void * arg); diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index 9677d5cdb..f7fa3ff60 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -726,7 +726,7 @@ void PolarBasis::determine_coefficients_particles(void) for (auto & v : expcoef0) { for (auto & u : v) u->setZero(); } use1 = 0; - if (multistep==0) used = 0; + if (multistep==0 or (mstep==0 and mlevel==multistep)) used = 0; #ifdef DEBUG cout << "Process " << myid @@ -1429,7 +1429,7 @@ void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) if (not NO_M0) { if (M0_back) { - get_pot_background(r, zz, p, drc, dzc); + auto [p, drc, dzc] = get_pot_background(r, zz); potl = mfac * p; potr = mfac * drc; From d0c1f8a0039eb05737e06f5ce896aeee2712c6bb Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 16 Feb 2024 13:55:09 -0500 Subject: [PATCH 37/65] Allow a monopole override option for PolarBasis [no ci] --- src/PolarBasis.H | 6 +++++- src/PolarBasis.cc | 11 +++++++++-- utils/ICs/ZangICs.cc | 3 ++- 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/PolarBasis.H b/src/PolarBasis.H index ae43d4af8..caab18f0d 100644 --- a/src/PolarBasis.H +++ b/src/PolarBasis.H @@ -32,7 +32,8 @@ class MixtureBasis; @param rmax is the maximum radius of the basis @param self_consistent true for updating the the coefficients at time steps @param NO_M0 true to omit the m=0 harmonic (default: false) - @param NO_M1 true to omint the m=1 harmonic (default: false) + @param NO_M1 true to omit the m=1 harmonic (default: false) + @param NO_MONO true to omit the monopole off-grid force algorithm (default: false) @param EVEN_M true to include only even m harmonics (default: false) @param M0_ONLY true to include only the m=0 harmonic (default: false) @param M0_BACK true includes fixed m=0 harmonic (default: false) @@ -402,6 +403,9 @@ protected: //! Omit monopole bool NO_M0; + //! Omit off-grid force algorithm + bool NO_MONO; + //! Flag drop m=1 term bool NO_M1; diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index f7fa3ff60..4e7926c4f 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -42,6 +42,7 @@ PolarBasis::valid_keys = { "EVEN_M", "M0_ONLY", "M0_BACK", + "NO_MONO", "mlim", "ssfrac", "playback", @@ -72,6 +73,7 @@ PolarBasis::PolarBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : EVEN_M = false; M0_only = false; M0_back = false; + NO_MONO = false; ssfrac = 0.0; subset = false; coefMaster = true; @@ -108,7 +110,7 @@ PolarBasis::PolarBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (conf["EVEN_M"]) EVEN_M = conf["EVEN_M"]. as(); if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); if (conf["M0_BACK"]) M0_back = conf["M0_BACK"].as(); - + if (conf["NO_MONO"]) NO_MONO = conf["NO_MONO"].as(); if (conf["ssfrac"]) { ssfrac = conf["ssfrac"].as(); @@ -1404,7 +1406,9 @@ void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) constexpr double midpt = ratmin + 0.5*(1.0 - ratmin); constexpr double rsmth = 0.5*(1.0 - ratmin)/maxerf; - if (ratio >= 1.0) { + if (NO_MONO) { + ratio = 0.0; + } else if (ratio >= 1.0) { frac = 0.0; cfrac = 1.0; } else if (ratio > ratmin) { @@ -1414,6 +1418,9 @@ void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) frac = 1.0; } + // TEST ratio override + ratio = 0.0; + // Ongrid contribution // if (ratio < 1.0) { diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index 75b92219b..bece0fce8 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -127,7 +127,8 @@ main(int ac, char **av) double Kmin = Ktol, Kmax = 1.0 - Ktol; double Emin = 0.5*Rmin*model->get_dpot(Rmin) + model->get_pot(Rmin); - double Emax = 0.5*Rmax*model->get_dpot(Rmax) + model->get_pot(Rmax); + // double Emax = 0.5*Rmax*model->get_dpot(Rmax) + model->get_pot(Rmax); + double Emax = model->get_pot(Rmax); // Scan to find the peak df // From d18209d6e97fa73fd29221eb42443165d7ef8d28 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 17 Feb 2024 11:36:22 -0500 Subject: [PATCH 38/65] Preliminary implementation of background arrays for flat disk [no ci] --- src/cudaBiorthCyl.cu | 39 +++++++++++++++++++++++++- src/cudaPolarBasis.cu | 65 +++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 100 insertions(+), 4 deletions(-) diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index 9be032270..82b668d5a 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -60,7 +60,9 @@ void BiorthCyl::initialize_cuda { // Number of texture arrays // - size_t ndim = (mmax+1)*nmax; + size_t ndim1 = (mmax+1)*nmax; + size_t ndim2 = 2; + size_t ndim = ndim1 + ndim2; // Interpolation data array // @@ -146,6 +148,41 @@ void BiorthCyl::initialize_cuda } // harmonic subspace loop + // Add background arrays + // + std::vector> vv(ndim2); + for (auto & v : vv) v.resize(numr); + + double dx0 = (xmax - xmin)/(numr - 1); + + for (int i=0; i P, dArray I, cuFP_t ratio = sqrt( (R2 + zz*zz)/rmax2 ); cuFP_t mfactor = 1.0, frac = 1.0, cfrac = 0.0; - if (ratio >= 1.0) { + if (plrNoMono) { + ratio = 0.0; + } else if (ratio >= 1.0) { frac = 0.0; cfrac = 1.0; } else if (ratio > ratmin) { @@ -664,6 +686,43 @@ forceKernelPlr6(dArray P, dArray I, if (plrNO_M1 and mm==1 ) continue; if (plrEVEN_M and (mm/2)*2 != mm) continue; + if (plrM0back and mm==0 ) { + int ndim1 = (mmax+1)*nmax; + + cuFP_t xx = X*plrDxi/plrDx0; + int ind = floor(xx); + + if (ind<0) ind = 0; + if (ind>plrNumr-2) ind = plrNumr - 2; + + cuFP_t a = (cuFP_t)(ind+1) - xx; + cuFP_t b = 1.0 - a; + + // Do the interpolation for the prefactor potential + // +#if cuREAL == 4 + cuFP_t pp0 = tex1D(tex._v[ndim1+0], ind ); + cuFP_t pp1 = tex1D(tex._v[ndim1+0], ind+1); + cuFP_t dp0 = -tex1D(tex._v[ndim1+1], ind ); + cuFP_t dp1 = -tex1D(tex._v[ndim1+1], ind+1); + +#else + cuFP_t pp0 = int2_as_double(tex1D(tex._v[ndim1+0], ind )); + cuFP_t pp1 = int2_as_double(tex1D(tex._v[ndim1+0], ind+1)); + cuFP_t dp0 = -int2_as_double(tex1D(tex._v[ndim1+1], ind )); + cuFP_t dp1 = -int2_as_double(tex1D(tex._v[ndim1+1], ind+1)); +#endif + if (xx<=0.0) { + pp += pp0; + fr += dp0; + } else { + pp += a*pp0 + b*pp1; + pp += a*dp0 + b*dp1; + } + + continue; + } + for (int n=0; n Date: Sat, 17 Feb 2024 11:55:45 -0500 Subject: [PATCH 39/65] Add NO_MONO as allowed flag for compatbility with EXP n-body [no ci] --- coefs/BiorthBasis.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/coefs/BiorthBasis.cc b/coefs/BiorthBasis.cc index 8e357d77d..a9e49f98a 100644 --- a/coefs/BiorthBasis.cc +++ b/coefs/BiorthBasis.cc @@ -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; @@ -1407,6 +1407,7 @@ namespace BasisClasses "EVEN_M", "M0_BACK", "M0_ONLY", + "NO_MONO", "diskconf", "ssfrac", "playback", @@ -1854,7 +1855,7 @@ namespace BasisClasses } // Radial grid spacing - double dx = (logxmax - logxmin)/numgrid; + double dx = (logxmax - logxmin)/(numgrid-1); // Basis storage Eigen::MatrixXd tabpot, tabden, tabrfc; From 7f1d85d4e48bc17448082b4e4f0f8df9714025ca Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 17 Feb 2024 11:56:07 -0500 Subject: [PATCH 40/65] A few minor cuda implementation gaffs [no ci] --- src/cudaBiorthCyl.cu | 53 +++++++++++++++++++++++++++++++------------ src/cudaPolarBasis.cu | 9 +++++++- 2 files changed, 46 insertions(+), 16 deletions(-) diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index 82b668d5a..1a2864ac5 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -150,37 +150,60 @@ void BiorthCyl::initialize_cuda // Add background arrays // - std::vector> vv(ndim2); - for (auto & v : vv) v.resize(numr); + std::vector> tt(ndim2); + for (auto & v : tt) v.resize(numr); double dx0 = (xmax - xmin)/(numr - 1); for (int i=0; i(); +#else + cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(); +#endif + + // Define the texture parameters + // + cudaTextureDesc texDesc1; - memset(&texDesc[n], 0, sizeof(cudaTextureDesc)); - texDesc[n].addressMode[n] = cudaAddressModeClamp; - texDesc[n].filterMode = cudaFilterModePoint; - texDesc[n].readMode = cudaReadModeElementType; - texDesc[n].normalizedCoords = 0; + memset(&texDesc1, 0, sizeof(cudaTextureDesc)); + texDesc1.addressMode[0] = cudaAddressModeClamp; + texDesc1.filterMode = cudaFilterModePoint; + texDesc1.readMode = cudaReadModeElementType; + texDesc1.normalizedCoords = 0; + // Size of interpolation array + // + size_t tsize = numr*sizeof(cuFP_t); + + // Make the textures + // + for (int n=ndim1; n P, dArray I, fr += dp0; } else { pp += a*pp0 + b*pp1; - pp += a*dp0 + b*dp1; + fr += a*dp0 + b*dp1; } +#ifdef TEMP_DEBUG + if (npart < 10) { + printf("background: %12.4e [%12.4e %12.4e] [%12.4e %12.4e]\n", + R, pp, log(R), fr, -1.0/R); + } +#endif continue; } From 24d844f1c1306b134a08aff487945956df43b70b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 17 Feb 2024 12:12:38 -0500 Subject: [PATCH 41/65] Improve texture counting logic for BiorthCyl and fix check [no ci] --- src/cudaBiorthCyl.cu | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index 1a2864ac5..cbe29d718 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -186,11 +186,11 @@ void BiorthCyl::initialize_cuda // Make the textures // - for (int n=ndim1; n Date: Sat, 17 Feb 2024 12:23:16 -0500 Subject: [PATCH 42/65] Background calculation needs to be the flat disk force routine [no ci] --- src/cudaBiorthCyl.cu | 10 ++++- src/cudaPolarBasis.cu | 90 ++++++++++++++++++++++--------------------- 2 files changed, 54 insertions(+), 46 deletions(-) diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index cbe29d718..4d6f1ebe4 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -131,16 +131,21 @@ void BiorthCyl::initialize_cuda cuda_safe_call(cudaMemcpy3D(©Params), __FILE__, __LINE__, "Error in copying 3d pitched array"); + // Specify the texture + // cudaResourceDesc resDesc; memset(&resDesc, 0, sizeof(cudaResourceDesc)); resDesc.resType = cudaResourceTypeArray; resDesc.res.array.array = cuArray[k]; + // Create texture object + // cuda_safe_call (cudaCreateTextureObject(&tex[k], &resDesc, &texDesc, NULL), __FILE__, __LINE__, "Failure in 2d texture creation"); + // Advance to next array k++; } // radial order loop @@ -192,7 +197,7 @@ void BiorthCyl::initialize_cuda cuda_safe_call(cudaMemcpyToArray(cuArray[k], 0, 0, &tt[n], tsize, cudaMemcpyHostToDevice), __FILE__, __LINE__, "copy texture to array"); - // Specify texture + // Specify the texture // cudaResourceDesc resDesc; @@ -204,10 +209,11 @@ void BiorthCyl::initialize_cuda // cuda_safe_call(cudaCreateTextureObject(&tex[k], &resDesc, &texDesc1, NULL), __FILE__, __LINE__, "create texture object"); + // Advance to next array k++; } - assert(k == ndim); + assert(k == ndim); // Sanity check cuda_safe_call(cudaFree(d_Interp), __FILE__, __LINE__, "Failure freeing device memory"); diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index 404bb5065..ab18a21ad 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -687,49 +687,6 @@ forceKernelPlr6(dArray P, dArray I, if (plrNO_M1 and mm==1 ) continue; if (plrEVEN_M and (mm/2)*2 != mm) continue; - if (plrM0back and mm==0 ) { - int ndim1 = (mmax+1)*nmax; - - cuFP_t xx = X*plrDxi/plrDx0; - int ind = floor(xx); - - if (ind<0) ind = 0; - if (ind>plrNumr-2) ind = plrNumr - 2; - - cuFP_t a = (cuFP_t)(ind+1) - xx; - cuFP_t b = 1.0 - a; - - // Do the interpolation for the prefactor potential - // -#if cuREAL == 4 - cuFP_t pp0 = tex1D(tex._v[ndim1+0], ind ); - cuFP_t pp1 = tex1D(tex._v[ndim1+0], ind+1); - cuFP_t dp0 = -tex1D(tex._v[ndim1+1], ind ); - cuFP_t dp1 = -tex1D(tex._v[ndim1+1], ind+1); - -#else - cuFP_t pp0 = int2_as_double(tex1D(tex._v[ndim1+0], ind )); - cuFP_t pp1 = int2_as_double(tex1D(tex._v[ndim1+0], ind+1)); - cuFP_t dp0 = -int2_as_double(tex1D(tex._v[ndim1+1], ind )); - cuFP_t dp1 = -int2_as_double(tex1D(tex._v[ndim1+1], ind+1)); -#endif - if (xx<=0.0) { - pp += pp0; - fr += dp0; - } else { - pp += a*pp0 + b*pp1; - fr += a*dp0 + b*dp1; - } - -#ifdef TEMP_DEBUG - if (npart < 10) { - printf("background: %12.4e [%12.4e %12.4e] [%12.4e %12.4e]\n", - R, pp, log(R), fr, -1.0/R); - } -#endif - continue; - } - for (int n=0; n P, dArray I, cuFP_t ratio = sqrt( (R2 + zz*zz)/rmax2 ); cuFP_t mfactor = 1.0, frac = 1.0, cfrac = 0.0; - if (ratio >= 1.0) { + if (plrNoMono) { + ratio = 0.0; + } if (ratio >= 1.0) { frac = 0.0; cfrac = 1.0; } else if (ratio > ratmin) { @@ -1207,6 +1166,49 @@ forceKernelPlr3(dArray P, dArray I, if (plrNO_M1 and mm==1 ) continue; if (plrEVEN_M and (mm/2)*2 != mm) continue; + if (plrM0back and mm==0 ) { + int ndim1 = (mmax+1)*nmax; + + cuFP_t xx = X*plrDxi/plrDx0; + int ind = floor(xx); + + if (ind<0) ind = 0; + if (ind>plrNumr-2) ind = plrNumr - 2; + + cuFP_t a = (cuFP_t)(ind+1) - xx; + cuFP_t b = 1.0 - a; + + // Do the interpolation for the prefactor potential + // +#if cuREAL == 4 + cuFP_t pp0 = tex1D(tex._v[ndim1+0], ind ); + cuFP_t pp1 = tex1D(tex._v[ndim1+0], ind+1); + cuFP_t dp0 = -tex1D(tex._v[ndim1+1], ind ); + cuFP_t dp1 = -tex1D(tex._v[ndim1+1], ind+1); + +#else + cuFP_t pp0 = int2_as_double(tex1D(tex._v[ndim1+0], ind )); + cuFP_t pp1 = int2_as_double(tex1D(tex._v[ndim1+0], ind+1)); + cuFP_t dp0 = -int2_as_double(tex1D(tex._v[ndim1+1], ind )); + cuFP_t dp1 = -int2_as_double(tex1D(tex._v[ndim1+1], ind+1)); +#endif + if (xx<=0.0) { + pp += pp0; + fr += dp0; + } else { + pp += a*pp0 + b*pp1; + fr += a*dp0 + b*dp1; + } + +#ifdef TEMP_DEBUG + if (npart < 10) { + printf("background: %8d %12.4e [%12.4e %12.4e] [%12.4e %12.4e] [%5.4f %5.4f]\n", + i, R, pp, log(R), fr, -1.0/R, a, b); + } +#endif + continue; + } + if (mm) norm = norm1; else norm = norm0; From fc333cf8a8684a6543ccfe8bd53079b18796c2aa Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 17 Feb 2024 12:36:32 -0500 Subject: [PATCH 43/65] Typo on cuda boolean constant assignement [no ci] --- src/cudaPolarBasis.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index ab18a21ad..0f832bf74 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -248,7 +248,7 @@ void PolarBasis::initialize_mapping_constants() cuda_safe_call(cudaMemcpyToSymbol(plrM0only, &M0_only, sizeof(bool), size_t(0), cudaMemcpyHostToDevice), __FILE__, __LINE__, "Error copying plrM0only"); - cuda_safe_call(cudaMemcpyToSymbol(plrM0only, &M0_back, sizeof(bool), size_t(0), cudaMemcpyHostToDevice), + cuda_safe_call(cudaMemcpyToSymbol(plrM0back, &M0_back, sizeof(bool), size_t(0), cudaMemcpyHostToDevice), __FILE__, __LINE__, "Error copying plrM0back"); cuda_safe_call(cudaMemcpyToSymbol(plrNoMono, &NO_MONO, sizeof(bool), size_t(0), cudaMemcpyHostToDevice), From 8390f5f03cb304000bf0423ab6cd5f71cd05bf9c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 17 Feb 2024 12:38:48 -0500 Subject: [PATCH 44/65] Add numr to cuda constant verbose output [no ci] --- src/cudaPolarBasis.cu | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index 0f832bf74..fd79d5b52 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -82,6 +82,7 @@ void testConstantsPlr() printf(" Dx0 = %f\n", plrDx0 ); printf(" Numx = %d\n", plrNumx ); printf(" Numy = %d\n", plrNumy ); + printf(" Numr = %d\n", plrNumr ); printf(" CmapR = %d\n", plrCmapR ); printf(" CmapZ = %d\n", plrCmapZ ); printf(" Orient = %d\n", plrOrient); From a79437ad6ee5f5ae2046d17020ca915db5608f7d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 18 Feb 2024 11:09:39 -0500 Subject: [PATCH 45/65] Move texture arrays to class data; fix typo in array reference [no ci] --- exputil/cudaSLGridMP2.cu | 3 - include/BiorthCyl.H | 8 ++ include/SLGridMP2.H | 8 ++ src/cudaBiorthCyl.cu | 164 +++++++++++++++++++++++++++++++-------- 4 files changed, 147 insertions(+), 36 deletions(-) diff --git a/exputil/cudaSLGridMP2.cu b/exputil/cudaSLGridMP2.cu index 58202f19f..780fc34eb 100644 --- a/exputil/cudaSLGridMP2.cu +++ b/exputil/cudaSLGridMP2.cu @@ -49,9 +49,6 @@ thrust::host_vector returnTestSph return f_d; } -static std::vector resDesc; -static std::vector texDesc; - struct Element { int l; double a; diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index da3cd51f5..e61c1a9ce 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -100,6 +100,14 @@ protected: //! Read the HDF5 cache virtual bool ReadH5Cache(); +#if HAVE_LIBCUDA==1 + //@{ + //! Texture sturctures + std::vector resDesc; + std::vector texDesc; + //@} +#endif + public: //! Flag for MPI enabled (default: 0=off) diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 2a04a4c44..b2184e822 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -238,6 +238,14 @@ private: //! Read HDF5 cache bool ReadH5Cache(); +#if HAVE_LIBCUDA==1 + //@{ + //! Texture sturctures + std::vector resDesc; + std::vector texDesc; + //@} +#endif + public: //! Flag for MPI enabled (default: 0=off) diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index 4d6f1ebe4..fe91c4e0c 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -53,6 +53,39 @@ thrust::host_vector returnTestBioCyl return f_d; } +__global__ +void testFetchBioCyl1(dArray T, dArray f, int l) +{ + const int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i(T._v[l], i); +#else + f._v[i] = int2_as_double(tex1D(T._v[l], i)); +#endif + } +} + +thrust::host_vector returnTestBioCyl1 + (thrust::host_vector& tex, int l, unsigned int N) +{ + thrust::device_vector t_d = tex; + + unsigned int gridSize = N/BLOCK_SIZE; + if (N > gridSize*BLOCK_SIZE) gridSize++; + + thrust::device_vector f_d(N); + + testFetchBioCyl1<<>>(toKernel(t_d), toKernel(f_d), l); + + cudaDeviceSynchronize(); + + return f_d; +} + +static std::vector resDesc; +static std::vector texDesc; + void BiorthCyl::initialize_cuda (std::vector& cuArray, thrust::host_vector& tex @@ -64,6 +97,14 @@ void BiorthCyl::initialize_cuda size_t ndim2 = 2; size_t ndim = ndim1 + ndim2; + // Cuda resouce descriptions + // + resDesc.resize(ndim); + + // Specify texture object parameters + // + texDesc.resize(ndim); + // Interpolation data array // cuArray.resize(ndim); @@ -73,15 +114,6 @@ void BiorthCyl::initialize_cuda tex.resize(ndim); thrust::fill(tex.begin(), tex.end(), 0); - cudaTextureDesc texDesc; - - memset(&texDesc, 0, sizeof(texDesc)); - texDesc.readMode = cudaReadModeElementType; - texDesc.filterMode = cudaFilterModePoint; - texDesc.addressMode[0] = cudaAddressModeClamp; - texDesc.addressMode[1] = cudaAddressModeClamp; - texDesc.addressMode[2] = cudaAddressModeClamp; - // Temporary storage // cuFP_t *d_Interp; @@ -96,6 +128,13 @@ void BiorthCyl::initialize_cuda for (size_t n=0; n(); #endif - // Define the texture parameters - // - cudaTextureDesc texDesc1; - - memset(&texDesc1, 0, sizeof(cudaTextureDesc)); - texDesc1.addressMode[0] = cudaAddressModeClamp; - texDesc1.filterMode = cudaFilterModePoint; - texDesc1.readMode = cudaReadModeElementType; - texDesc1.normalizedCoords = 0; - // Size of interpolation array // size_t tsize = numr*sizeof(cuFP_t); @@ -193,21 +220,31 @@ void BiorthCyl::initialize_cuda // for (int n=0; n::epsilon(); struct Element { @@ -288,4 +325,65 @@ void BiorthCyl::initialize_cuda << "**Hi [1] : " << hi1->first << " (" << hi1->second.m << ", " << hi1->second.k << ", " << hi1->second.a << ", " << hi1->second.b << ", " << hi1->second.a - hi1->second.b << ")" << std::endl << "**" << std::endl; } + + if (false and myid==0) { + + constexpr cuFP_t tol = 10.0*std::numeric_limits::epsilon(); + + struct Element { + int l; + int k; + double a; + double b; + }; + + std::multimap compare; + unsigned tot = 0, bad = 0; + + for (int L=0; L<2; L++) { + + std::cout << "**HOST** Texture 1D compare L=" << L << std::endl; + + thrust::host_vector xyg = returnTestBioCyl1(tex, ndim1 + L, numr); + + for (int i=0; i1.0e-18) { + + Element comp = {L, i, a, b}; + compare.insert(std::make_pair(fabs((a - b)/a), comp)); + + if ( fabs((a - b)/a ) > tol) { + std::cout << std::setw( 5) << L << std::setw( 8) << i + << std::setw(20) << a << std::setw(20) << b + << std::setw(20) << (a-b)/a << std::endl; + bad++; + } + } + tot++; + } + } + + std::multimap::iterator beg = compare.begin(); + std::multimap::iterator end = compare.end(); + std::multimap::iterator lo1=beg, lo9=beg, mid=beg, hi9=end, hi1=end; + + std::advance(lo9, 9); + std::advance(mid, compare.size()/2); + std::advance(hi1, -1); + std::advance(hi9, -10); + + std::cout << "**Found " << bad << "/" << tot << " values > eps" << std::endl + << "**Low[1] : " << lo1->first << " (" << lo1->second.l << ", " << lo1->second.k << ", " << lo1->second.a << ", " << lo1->second.b << ", " << lo1->second.a - lo1->second.b << ")" << std::endl + << "**Low[9] : " << lo9->first << " (" << lo9->second.l << ", " << lo9->second.k << ", " << lo9->second.a << ", " << lo9->second.b << ", " << lo9->second.a - lo9->second.b << ")" << std::endl + << "**Middle : " << mid->first << " (" << mid->second.l << ", " << mid->second.k << ", " << mid->second.a << ", " << mid->second.b << ", " << mid->second.a - mid->second.b << ")" << std::endl + << "**Hi [9] : " << hi9->first << " (" << hi9->second.l << ", " << hi9->second.k << ", " << hi9->second.a << ", " << lo1->second.b << ", " << hi9->second.a - hi9->second.b << ")" << std::endl + << "**Hi [1] : " << hi1->first << " (" << hi1->second.l << ", " << hi1->second.k << ", " << hi1->second.a << ", " << hi1->second.b << ", " << hi1->second.a - hi1->second.b << ")" << std::endl + << "**" << std::endl; + } + + } From 0e86a0ef379fe618751a43b7eb2c36cbfc356fdb Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 18 Feb 2024 11:17:03 -0500 Subject: [PATCH 46/65] Make texture description and resource sturctures automatic variables [no ci] --- exputil/cudaSLGridMP2.cu | 51 +++++++++++++++++++---------------- include/BiorthCyl.H | 8 ------ include/SLGridMP2.H | 8 ------ src/cudaBiorthCyl.cu | 57 +++++++++++++++++++--------------------- 4 files changed, 55 insertions(+), 69 deletions(-) diff --git a/exputil/cudaSLGridMP2.cu b/exputil/cudaSLGridMP2.cu index 780fc34eb..8d51ac338 100644 --- a/exputil/cudaSLGridMP2.cu +++ b/exputil/cudaSLGridMP2.cu @@ -83,18 +83,17 @@ void SLGridSph::initialize_cuda(std::vector& cuArray, tex.resize(ndim); thrust::fill(tex.begin(), tex.end(), 0); - // cudaResourceDesc resDesc; - resDesc.resize(ndim); + // std::vector resDesc; + // std::vector texDesc; - // Specify texture object parameters - // - texDesc.resize(ndim); - memset(&texDesc[0], 0, sizeof(cudaTextureDesc)); - texDesc[0].addressMode[0] = cudaAddressModeClamp; - texDesc[0].filterMode = cudaFilterModePoint; - texDesc[0].readMode = cudaReadModeElementType; - texDesc[0].normalizedCoords = 0; + cudaTextureDesc texDesc; + + memset(&texDesc, 0, sizeof(cudaTextureDesc)); + texDesc.addressMode[0] = cudaAddressModeClamp; + texDesc.filterMode = cudaFilterModePoint; + texDesc.readMode = cudaReadModeElementType; + texDesc.normalizedCoords = 0; thrust::host_vector tt(numr); @@ -107,11 +106,13 @@ void SLGridSph::initialize_cuda(std::vector& cuArray, cuda_safe_call(cudaMemcpyToArray(cuArray[0], 0, 0, &tt[0], tsize, cudaMemcpyHostToDevice), __FILE__, __LINE__, "copy texture to array"); // Specify texture - memset(&resDesc[0], 0, sizeof(cudaResourceDesc)); - resDesc[0].resType = cudaResourceTypeArray; - resDesc[0].res.array.array = cuArray[0]; + cudaResourceDesc resDesc; + + memset(&resDesc, 0, sizeof(cudaResourceDesc)); + resDesc.resType = cudaResourceTypeArray; + resDesc.res.array.array = cuArray[0]; - cuda_safe_call(cudaCreateTextureObject(&tex[0], &resDesc[0], &texDesc[0], NULL), __FILE__, __LINE__, "create texture object"); + cuda_safe_call(cudaCreateTextureObject(&tex[0], &resDesc, &texDesc, NULL), __FILE__, __LINE__, "create texture object"); for (int l=0; l<=lmax; l++) { for (int n=0; n& cuArray, cuda_safe_call(cudaMemcpyToArray(cuArray[i], 0, 0, &tt[0], tsize, cudaMemcpyHostToDevice), __FILE__, __LINE__, "copy texture to array"); // Specify texture - memset(&resDesc[i], 0, sizeof(cudaResourceDesc)); - resDesc[i].resType = cudaResourceTypeArray; - resDesc[i].res.array.array = cuArray[i]; + cudaResourceDesc resDesc; + + memset(&resDesc, 0, sizeof(cudaResourceDesc)); + resDesc.resType = cudaResourceTypeArray; + resDesc.res.array.array = cuArray[i]; + + cudaTextureDesc texDesc; - memset(&texDesc[i], 0, sizeof(cudaTextureDesc)); - texDesc[i].addressMode[0] = cudaAddressModeClamp; - texDesc[i].filterMode = cudaFilterModePoint; - texDesc[i].readMode = cudaReadModeElementType; - texDesc[i].normalizedCoords = 0; + memset(&texDesc, 0, sizeof(cudaTextureDesc)); + texDesc.addressMode[0] = cudaAddressModeClamp; + texDesc.filterMode = cudaFilterModePoint; + texDesc.readMode = cudaReadModeElementType; + texDesc.normalizedCoords = 0; - cuda_safe_call(cudaCreateTextureObject(&tex[i], &resDesc[i], &texDesc[i], NULL), __FILE__, __LINE__, "create texture object"); + cuda_safe_call(cudaCreateTextureObject(&tex[i], &resDesc, &texDesc, NULL), __FILE__, __LINE__, "create texture object"); } } diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index e61c1a9ce..da3cd51f5 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -100,14 +100,6 @@ protected: //! Read the HDF5 cache virtual bool ReadH5Cache(); -#if HAVE_LIBCUDA==1 - //@{ - //! Texture sturctures - std::vector resDesc; - std::vector texDesc; - //@} -#endif - public: //! Flag for MPI enabled (default: 0=off) diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index b2184e822..2a04a4c44 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -238,14 +238,6 @@ private: //! Read HDF5 cache bool ReadH5Cache(); -#if HAVE_LIBCUDA==1 - //@{ - //! Texture sturctures - std::vector resDesc; - std::vector texDesc; - //@} -#endif - public: //! Flag for MPI enabled (default: 0=off) diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index fe91c4e0c..218e70177 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -83,9 +83,6 @@ thrust::host_vector returnTestBioCyl1 return f_d; } -static std::vector resDesc; -static std::vector texDesc; - void BiorthCyl::initialize_cuda (std::vector& cuArray, thrust::host_vector& tex @@ -97,14 +94,6 @@ void BiorthCyl::initialize_cuda size_t ndim2 = 2; size_t ndim = ndim1 + ndim2; - // Cuda resouce descriptions - // - resDesc.resize(ndim); - - // Specify texture object parameters - // - texDesc.resize(ndim); - // Interpolation data array // cuArray.resize(ndim); @@ -128,12 +117,14 @@ void BiorthCyl::initialize_cuda for (size_t n=0; n Date: Sun, 18 Feb 2024 11:24:27 -0500 Subject: [PATCH 47/65] Disable cuda debugging [no ci] --- src/cudaPolarBasis.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index fd79d5b52..98855711c 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -11,7 +11,7 @@ // #define BOUNDS_CHECK // #define VERBOSE_CTR // #define VERBOSE_DBG -#define TEMP_DEBUG +// #define TEMP_DEBUG // Global symbols for coordinate transformation // From 9557eb9f5e2ff50a3201daecb6e28c5034b880f2 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 18 Feb 2024 11:25:29 -0500 Subject: [PATCH 48/65] Remove debugging code [no ci] --- src/cudaPolarBasis.cu | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index 98855711c..933fa050b 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -11,7 +11,6 @@ // #define BOUNDS_CHECK // #define VERBOSE_CTR // #define VERBOSE_DBG -// #define TEMP_DEBUG // Global symbols for coordinate transformation // @@ -1201,12 +1200,6 @@ forceKernelPlr3(dArray P, dArray I, fr += a*dp0 + b*dp1; } -#ifdef TEMP_DEBUG - if (npart < 10) { - printf("background: %8d %12.4e [%12.4e %12.4e] [%12.4e %12.4e] [%5.4f %5.4f]\n", - i, R, pp, log(R), fr, -1.0/R, a, b); - } -#endif continue; } From 3c46a8375b17114b232ff27b5a0ff30ca923fede Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 18 Feb 2024 12:34:23 -0500 Subject: [PATCH 49/65] Version bump [no ci] --- CMakeLists.txt | 2 +- doc/exp.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 410d81683..ab0c66397 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/doc/exp.cfg b/doc/exp.cfg index 38468d75a..db58900a8 100644 --- a/doc/exp.cfg +++ b/doc/exp.cfg @@ -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 From 03d816a14b347f82bc5f8904c1d3b1167b40bd79 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 18 Feb 2024 16:27:10 -0500 Subject: [PATCH 50/65] Fixed missing m=0 sine/cosine update in kernel loop [no ci] --- src/cudaPolarBasis.cu | 114 +++++++++++++++++++++--------------------- 1 file changed, 57 insertions(+), 57 deletions(-) diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index 933fa050b..817e16880 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -1200,81 +1200,81 @@ forceKernelPlr3(dArray P, dArray I, fr += a*dp0 + b*dp1; } - continue; - } + } else { - if (mm) norm = norm1; - else norm = norm0; + if (mm) norm = norm1; + else norm = norm0; - for (int n=0; n(tex._v[k], indX, indY , 0) * c00 + - tex3D(tex._v[k], indX+1, indY , 0) * c10 + - tex3D(tex._v[k], indX, indY+1, 0) * c01 + - tex3D(tex._v[k], indX+1, indY+1, 0) * c11 + tex3D(tex._v[k], indX, indY , 0) * c00 + + tex3D(tex._v[k], indX+1, indY , 0) * c10 + + tex3D(tex._v[k], indX, indY+1, 0) * c01 + + tex3D(tex._v[k], indX+1, indY+1, 0) * c11 #else - int2_as_double(tex3D(tex._v[k], indX, indY , 0)) * c00 + - int2_as_double(tex3D(tex._v[k], indX+1, indY , 0)) * c10 + - int2_as_double(tex3D(tex._v[k], indX, indY+1, 0)) * c01 + - int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 0)) * c11 + int2_as_double(tex3D(tex._v[k], indX, indY , 0)) * c00 + + int2_as_double(tex3D(tex._v[k], indX+1, indY , 0)) * c10 + + int2_as_double(tex3D(tex._v[k], indX, indY+1, 0)) * c01 + + int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 0)) * c11 #endif - ); + ); - cuFP_t rfrc = - ( + cuFP_t rfrc = + ( #if cuREAL == 4 - tex3D(tex._v[k], indX, indY , 1) * c00 + - tex3D(tex._v[k], indX+1, indY , 1) * c10 + - tex3D(tex._v[k], indX, indY+1, 1) * c01 + - tex3D(tex._v[k], indX+1, indY+1, 1) * c11 + tex3D(tex._v[k], indX, indY , 1) * c00 + + tex3D(tex._v[k], indX+1, indY , 1) * c10 + + tex3D(tex._v[k], indX, indY+1, 1) * c01 + + tex3D(tex._v[k], indX+1, indY+1, 1) * c11 #else - int2_as_double(tex3D(tex._v[k], indX, indY , 1)) * c00 + - int2_as_double(tex3D(tex._v[k], indX+1, indY , 1)) * c10 + - int2_as_double(tex3D(tex._v[k], indX, indY+1, 1)) * c01 + - int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 1)) * c11 + int2_as_double(tex3D(tex._v[k], indX, indY , 1)) * c00 + + int2_as_double(tex3D(tex._v[k], indX+1, indY , 1)) * c10 + + int2_as_double(tex3D(tex._v[k], indX, indY+1, 1)) * c01 + + int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 1)) * c11 #endif - ); + ); - cuFP_t zfrc = - ( + cuFP_t zfrc = + ( #if cuREAL == 4 - tex3D(tex._v[k], indX, indY , 2) * c00 + - tex3D(tex._v[k], indX+1, indY , 2) * c10 + - tex3D(tex._v[k], indX, indY+1, 2) * c01 + - tex3D(tex._v[k], indX+1, indY+1, 2) * c11 + tex3D(tex._v[k], indX, indY , 2) * c00 + + tex3D(tex._v[k], indX+1, indY , 2) * c10 + + tex3D(tex._v[k], indX, indY+1, 2) * c01 + + tex3D(tex._v[k], indX+1, indY+1, 2) * c11 #else - int2_as_double(tex3D(tex._v[k], indX, indY , 2)) * c00 + - int2_as_double(tex3D(tex._v[k], indX+1, indY , 2)) * c10 + - int2_as_double(tex3D(tex._v[k], indX, indY+1, 2)) * c01 + - int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 2)) * c11 + int2_as_double(tex3D(tex._v[k], indX, indY , 2)) * c00 + + int2_as_double(tex3D(tex._v[k], indX+1, indY , 2)) * c10 + + int2_as_double(tex3D(tex._v[k], indX, indY+1, 2)) * c01 + + int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 2)) * c11 #endif - ); + ); - if (zz < 0.0) zfrc *= -1.0; + if (zz < 0.0) zfrc *= -1.0; - // The trigonometric norm with a minus sign for the tabled values: - // -1/sqrt(2*pi) for m==0 or -1/sqrt(pi) for m>0 - // - cuFP_t Sfac = -norm; - cuFP_t facC = coef._v[IImn(mm, 'c', n, nmax)]; - cuFP_t facS = 0.0; - if (mm>0) { - facS = coef._v[IImn(mm, 's', n, nmax)]; + // The trigonometric norm with a minus sign for the tabled values: + // -1/sqrt(2*pi) for m==0 or -1/sqrt(pi) for m>0 + // + cuFP_t Sfac = -norm; + cuFP_t facC = coef._v[IImn(mm, 'c', n, nmax)]; + cuFP_t facS = 0.0; + if (mm>0) { + facS = coef._v[IImn(mm, 's', n, nmax)]; + } + + pp += potl * ( facC * ccos + facS * ssin) * Sfac; + fr += rfrc * ( facC * ccos + facS * ssin) * Sfac; + fz += zfrc * ( facC * ccos + facS * ssin) * Sfac; + fp += potl * ( facC * ssin - facS * ccos) * Sfac * mm; } - - pp += potl * ( facC * ccos + facS * ssin) * Sfac; - fr += rfrc * ( facC * ccos + facS * ssin) * Sfac; - fz += zfrc * ( facC * ccos + facS * ssin) * Sfac; - fp += potl * ( facC * ssin - facS * ccos) * Sfac * mm; } - + // Trig recursion to squeeze avoid internal FP fct call // cuFP_t cosM = ccos; From 3ed5b7dee5d4037f4cf307c44e0c2ba1112909e6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 18 Feb 2024 16:58:41 -0500 Subject: [PATCH 51/65] Wrapped evaluation in an if block to implement other exclusions [no ci] --- src/cudaPolarBasis.cu | 174 ++++++++++++++++++++++-------------------- 1 file changed, 91 insertions(+), 83 deletions(-) diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index 817e16880..6dc0fd0dd 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -1161,119 +1161,127 @@ forceKernelPlr3(dArray P, dArray I, for (int mm=0; mm<=muse; mm++) { - if (plrM0only and mm>0 ) continue; - if (plrNO_M0 and mm==0 ) continue; - if (plrNO_M1 and mm==1 ) continue; - if (plrEVEN_M and (mm/2)*2 != mm) continue; + bool compute = true; - if (plrM0back and mm==0 ) { - int ndim1 = (mmax+1)*nmax; + if (plrM0only and mm>0 ) compute = false; + if (plrNO_M0 and mm==0 ) compute = false; + if (plrNO_M1 and mm==1 ) compute = false; + if (plrEVEN_M and (mm/2)*2 != mm) compute = false; - cuFP_t xx = X*plrDxi/plrDx0; - int ind = floor(xx); + if (compute) { - if (ind<0) ind = 0; - if (ind>plrNumr-2) ind = plrNumr - 2; + if (plrM0back and mm==0) { + int ndim1 = (mmax+1)*nmax; - cuFP_t a = (cuFP_t)(ind+1) - xx; - cuFP_t b = 1.0 - a; + cuFP_t xx = X*plrDxi/plrDx0; + int ind = floor(xx); - // Do the interpolation for the prefactor potential - // + if (ind<0) ind = 0; + if (ind>plrNumr-2) ind = plrNumr - 2; + + cuFP_t a = (cuFP_t)(ind+1) - xx; + cuFP_t b = 1.0 - a; + + // Do the interpolation for the prefactor potential + // #if cuREAL == 4 - cuFP_t pp0 = tex1D(tex._v[ndim1+0], ind ); - cuFP_t pp1 = tex1D(tex._v[ndim1+0], ind+1); - cuFP_t dp0 = -tex1D(tex._v[ndim1+1], ind ); - cuFP_t dp1 = -tex1D(tex._v[ndim1+1], ind+1); + cuFP_t pp0 = tex1D(tex._v[ndim1+0], ind ); + cuFP_t pp1 = tex1D(tex._v[ndim1+0], ind+1); + cuFP_t dp0 = -tex1D(tex._v[ndim1+1], ind ); + cuFP_t dp1 = -tex1D(tex._v[ndim1+1], ind+1); #else - cuFP_t pp0 = int2_as_double(tex1D(tex._v[ndim1+0], ind )); - cuFP_t pp1 = int2_as_double(tex1D(tex._v[ndim1+0], ind+1)); - cuFP_t dp0 = -int2_as_double(tex1D(tex._v[ndim1+1], ind )); - cuFP_t dp1 = -int2_as_double(tex1D(tex._v[ndim1+1], ind+1)); + cuFP_t pp0 = int2_as_double(tex1D(tex._v[ndim1+0], ind )); + cuFP_t pp1 = int2_as_double(tex1D(tex._v[ndim1+0], ind+1)); + cuFP_t dp0 = -int2_as_double(tex1D(tex._v[ndim1+1], ind )); + cuFP_t dp1 = -int2_as_double(tex1D(tex._v[ndim1+1], ind+1)); #endif - if (xx<=0.0) { - pp += pp0; - fr += dp0; + if (xx<=0.0) { + pp += pp0; + fr += dp0; + } else { + pp += a*pp0 + b*pp1; + fr += a*dp0 + b*dp1; + } + } else { - pp += a*pp0 + b*pp1; - fr += a*dp0 + b*dp1; - } - } else { + if (mm) norm = norm1; + else norm = norm0; - if (mm) norm = norm1; - else norm = norm0; - - for (int n=0; n(tex._v[k], indX, indY , 0) * c00 + - tex3D(tex._v[k], indX+1, indY , 0) * c10 + - tex3D(tex._v[k], indX, indY+1, 0) * c01 + - tex3D(tex._v[k], indX+1, indY+1, 0) * c11 + tex3D(tex._v[k], indX, indY , 0) * c00 + + tex3D(tex._v[k], indX+1, indY , 0) * c10 + + tex3D(tex._v[k], indX, indY+1, 0) * c01 + + tex3D(tex._v[k], indX+1, indY+1, 0) * c11 #else - int2_as_double(tex3D(tex._v[k], indX, indY , 0)) * c00 + - int2_as_double(tex3D(tex._v[k], indX+1, indY , 0)) * c10 + - int2_as_double(tex3D(tex._v[k], indX, indY+1, 0)) * c01 + - int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 0)) * c11 + int2_as_double(tex3D(tex._v[k], indX, indY , 0)) * c00 + + int2_as_double(tex3D(tex._v[k], indX+1, indY , 0)) * c10 + + int2_as_double(tex3D(tex._v[k], indX, indY+1, 0)) * c01 + + int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 0)) * c11 #endif - ); + ); - cuFP_t rfrc = - ( + cuFP_t rfrc = + ( #if cuREAL == 4 - tex3D(tex._v[k], indX, indY , 1) * c00 + - tex3D(tex._v[k], indX+1, indY , 1) * c10 + - tex3D(tex._v[k], indX, indY+1, 1) * c01 + - tex3D(tex._v[k], indX+1, indY+1, 1) * c11 + tex3D(tex._v[k], indX, indY , 1) * c00 + + tex3D(tex._v[k], indX+1, indY , 1) * c10 + + tex3D(tex._v[k], indX, indY+1, 1) * c01 + + tex3D(tex._v[k], indX+1, indY+1, 1) * c11 #else - int2_as_double(tex3D(tex._v[k], indX, indY , 1)) * c00 + - int2_as_double(tex3D(tex._v[k], indX+1, indY , 1)) * c10 + - int2_as_double(tex3D(tex._v[k], indX, indY+1, 1)) * c01 + - int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 1)) * c11 + int2_as_double(tex3D(tex._v[k], indX, indY , 1)) * c00 + + int2_as_double(tex3D(tex._v[k], indX+1, indY , 1)) * c10 + + int2_as_double(tex3D(tex._v[k], indX, indY+1, 1)) * c01 + + int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 1)) * c11 #endif - ); + ); - cuFP_t zfrc = - ( + cuFP_t zfrc = + ( #if cuREAL == 4 - tex3D(tex._v[k], indX, indY , 2) * c00 + - tex3D(tex._v[k], indX+1, indY , 2) * c10 + - tex3D(tex._v[k], indX, indY+1, 2) * c01 + - tex3D(tex._v[k], indX+1, indY+1, 2) * c11 + tex3D(tex._v[k], indX, indY , 2) * c00 + + tex3D(tex._v[k], indX+1, indY , 2) * c10 + + tex3D(tex._v[k], indX, indY+1, 2) * c01 + + tex3D(tex._v[k], indX+1, indY+1, 2) * c11 #else - int2_as_double(tex3D(tex._v[k], indX, indY , 2)) * c00 + - int2_as_double(tex3D(tex._v[k], indX+1, indY , 2)) * c10 + - int2_as_double(tex3D(tex._v[k], indX, indY+1, 2)) * c01 + - int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 2)) * c11 + int2_as_double(tex3D(tex._v[k], indX, indY , 2)) * c00 + + int2_as_double(tex3D(tex._v[k], indX+1, indY , 2)) * c10 + + int2_as_double(tex3D(tex._v[k], indX, indY+1, 2)) * c01 + + int2_as_double(tex3D(tex._v[k], indX+1, indY+1, 2)) * c11 #endif - ); + ); - if (zz < 0.0) zfrc *= -1.0; + if (zz < 0.0) zfrc *= -1.0; - // The trigonometric norm with a minus sign for the tabled values: - // -1/sqrt(2*pi) for m==0 or -1/sqrt(pi) for m>0 - // - cuFP_t Sfac = -norm; - cuFP_t facC = coef._v[IImn(mm, 'c', n, nmax)]; - cuFP_t facS = 0.0; - if (mm>0) { - facS = coef._v[IImn(mm, 's', n, nmax)]; + // The trigonometric norm with a minus sign for the tabled values: + // -1/sqrt(2*pi) for m==0 or -1/sqrt(pi) for m>0 + // + cuFP_t Sfac = -norm; + cuFP_t facC = coef._v[IImn(mm, 'c', n, nmax)]; + cuFP_t facS = 0.0; + if (mm>0) { + facS = coef._v[IImn(mm, 's', n, nmax)]; + } + + pp += potl * ( facC * ccos + facS * ssin) * Sfac; + fr += rfrc * ( facC * ccos + facS * ssin) * Sfac; + fz += zfrc * ( facC * ccos + facS * ssin) * Sfac; + fp += potl * ( facC * ssin - facS * ccos) * Sfac * mm; } - - pp += potl * ( facC * ccos + facS * ssin) * Sfac; - fr += rfrc * ( facC * ccos + facS * ssin) * Sfac; - fz += zfrc * ( facC * ccos + facS * ssin) * Sfac; - fp += potl * ( facC * ssin - facS * ccos) * Sfac * mm; + // END: radial-order loop } + // END: 2-d interpolation block } + // END: compute block // Trig recursion to squeeze avoid internal FP fct call // From 64de90174618d4574d5facb89cbdd052d238de9f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 18 Feb 2024 17:24:44 -0500 Subject: [PATCH 52/65] More comments in cuda kernel [no ci] --- src/cudaPolarBasis.cu | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index 6dc0fd0dd..a10bdcd36 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -1171,8 +1171,12 @@ forceKernelPlr3(dArray P, dArray I, if (compute) { if (plrM0back and mm==0) { - int ndim1 = (mmax+1)*nmax; + // Offset of 1-d arrays in texture array + // + int offst = (mmax+1)*nmax; + // Compute X value for 1-d spacing from 2-d spacing + // cuFP_t xx = X*plrDxi/plrDx0; int ind = floor(xx); @@ -1185,16 +1189,16 @@ forceKernelPlr3(dArray P, dArray I, // Do the interpolation for the prefactor potential // #if cuREAL == 4 - cuFP_t pp0 = tex1D(tex._v[ndim1+0], ind ); - cuFP_t pp1 = tex1D(tex._v[ndim1+0], ind+1); - cuFP_t dp0 = -tex1D(tex._v[ndim1+1], ind ); - cuFP_t dp1 = -tex1D(tex._v[ndim1+1], ind+1); + cuFP_t pp0 = tex1D(tex._v[offst+0], ind ); + cuFP_t pp1 = tex1D(tex._v[offst+0], ind+1); + cuFP_t dp0 = -tex1D(tex._v[offst+1], ind ); + cuFP_t dp1 = -tex1D(tex._v[offst+1], ind+1); #else - cuFP_t pp0 = int2_as_double(tex1D(tex._v[ndim1+0], ind )); - cuFP_t pp1 = int2_as_double(tex1D(tex._v[ndim1+0], ind+1)); - cuFP_t dp0 = -int2_as_double(tex1D(tex._v[ndim1+1], ind )); - cuFP_t dp1 = -int2_as_double(tex1D(tex._v[ndim1+1], ind+1)); + cuFP_t pp0 = int2_as_double(tex1D(tex._v[offst+0], ind )); + cuFP_t pp1 = int2_as_double(tex1D(tex._v[offst+0], ind+1)); + cuFP_t dp0 = -int2_as_double(tex1D(tex._v[offst+1], ind )); + cuFP_t dp1 = -int2_as_double(tex1D(tex._v[offst+1], ind+1)); #endif if (xx<=0.0) { pp += pp0; From b3aa1f0ad7bb90684346cbce60525e781da99cf7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 19 Feb 2024 12:03:44 -0500 Subject: [PATCH 53/65] Fix long standing mistake in destroy_cuda members [no ci] --- src/cudaCylinder.cu | 16 ++++++++-------- src/cudaPolarBasis.cu | 30 +++++++++++++++++++----------- src/cudaSphericalBasis.cu | 17 ++++++++--------- 3 files changed, 35 insertions(+), 28 deletions(-) diff --git a/src/cudaCylinder.cu b/src/cudaCylinder.cu index bea08a813..8c978ea6b 100644 --- a/src/cudaCylinder.cu +++ b/src/cudaCylinder.cu @@ -2077,17 +2077,17 @@ void Cylinder::multistep_update_cuda() void Cylinder::destroy_cuda() { + // Deallocate the texture objects + // for (size_t i=0; i P, dArray I, for (int mm=0; mm<=muse; mm++) { - bool compute = true; + bool compute = true; // Compute the force for this m value + // by default + // Check restrictions + // if (plrM0only and mm>0 ) compute = false; if (plrNO_M0 and mm==0 ) compute = false; if (plrNO_M1 and mm==1 ) compute = false; @@ -1170,6 +1173,8 @@ forceKernelPlr3(dArray P, dArray I, if (compute) { + // Use the unperturbed fixed background model for m=0 + // if (plrM0back and mm==0) { // Offset of 1-d arrays in texture array // @@ -1208,7 +1213,10 @@ forceKernelPlr3(dArray P, dArray I, fr += a*dp0 + b*dp1; } - } else { + } + // END: unperturbed fixed background model for m=0 + // BEGIN: use basis function expansion + else { if (mm) norm = norm1; else norm = norm0; @@ -1283,7 +1291,7 @@ forceKernelPlr3(dArray P, dArray I, } // END: radial-order loop } - // END: 2-d interpolation block + // END: basis-function contribution block } // END: compute block @@ -2627,18 +2635,18 @@ void PolarBasis::multistep_update_cuda() void PolarBasis::destroy_cuda() { + // Deallocate textures + // for (size_t i=0; i Date: Mon, 19 Feb 2024 12:23:35 -0500 Subject: [PATCH 54/65] Change to keep the compute-sanitzer happy [no ci] --- src/cudaMultistep.cu | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/cudaMultistep.cu b/src/cudaMultistep.cu index ed39a205d..2425fd97d 100644 --- a/src/cudaMultistep.cu +++ b/src/cudaMultistep.cu @@ -484,8 +484,21 @@ void cuda_compute_levels() // Make thrust do device copy // - float minT = *(thrust::min_element(c->minDT.begin(), c->minDT.end())); - float maxT = *(thrust::max_element(c->maxDT.begin(), c->maxDT.end())); + thrust::host_vector minDT = c->minDT; + thrust::host_vector maxDT = c->maxDT; + + // Note: could not manage to do this on the device. So resorted + // to copying the block results and doing the reduction on the + // host. I assumed that I would be able to do something like + // this: + // + // float minT = *(thrust::min_element(c->minDT.begin(), c->minDT.end())); + // float maxT = *(thrust::max_element(c->maxDT.begin(), c->maxDT.end())); + // + // It seems to work, but the compute-sanitizer is not happy. + + float minT = *(thrust::min_element(minDT.begin(), minDT.end())); + float maxT = *(thrust::max_element(maxDT.begin(), maxDT.end())); if (minTmaxdt) maxdt = maxT; From bf6dc5d40c1e99d1fbfaeec44d00e76c201dd2a4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 22 Feb 2024 16:35:13 -0500 Subject: [PATCH 55/65] Spelling mistake --- src/Sphere.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Sphere.cc b/src/Sphere.cc index bd81633c3..95d4783ec 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -113,7 +113,7 @@ void Sphere::initialize() if (conf["dfac"]) dfac = conf["dfac"].as(); if (conf["modelname"]) model_file = conf["modelname"].as(); if (conf["cachename"]) cache_file = conf["cachename"].as(); - else throw std::runtime_error("Sphere: you must specify a cachname"); + else throw std::runtime_error("Sphere: you must specify a cachename"); if (conf["dtime"]) dtime = conf["dtime"].as(); if (conf["logr"]) logr = conf["logr"].as(); if (conf["plummer"]) plummer = conf["plummer"].as(); From b591decf384413c93f0d444b650f9243125d646c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 22 Feb 2024 16:35:47 -0500 Subject: [PATCH 56/65] Allow a default cachename for a true default cache in SLGridSph [no ci] --- include/SLGridMP2.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 2a04a4c44..399ffea83 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -249,7 +249,7 @@ public: SLGridSph(std::shared_ptr mod, int lmax, int nmax, int numr, double rmin, double rmax, bool cache, int Cmap, double RMAP, - std::string cachename="", + std::string cachename=".slgrid_sph_cache", bool Verbose=false); //! Constructor (uses file *model_file_name* for file) @@ -257,7 +257,7 @@ public: int lmax, int nmax, int numr, double rmin, double rmax, bool cache, int Cmap, double RMAP, int DIVERGE, double DFAC, - std::string cachename="", + std::string cachename=".slgrid_sph_cache", bool Verbose=false); //! Destructor From 04fc211647ef1e80a4daf7d1785026f3336f8edf Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 22 Feb 2024 16:37:31 -0500 Subject: [PATCH 57/65] Spelling mistake [no ci] Spelling mistake From 0865d9ec87e8d14803dee9b4017692f00b3ce033 Mon Sep 17 00:00:00 2001 From: mdw Date: Tue, 27 Feb 2024 07:27:36 -0500 Subject: [PATCH 58/65] Added NQDHT flag to allowed options in FlatDisk for experimentation with Bessel order [no ci] --- src/FlatDisk.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/FlatDisk.cc b/src/FlatDisk.cc index 2d2dc7345..ba5ab0ec2 100644 --- a/src/FlatDisk.cc +++ b/src/FlatDisk.cc @@ -15,6 +15,7 @@ FlatDisk::valid_keys = { "mmax", "numx", "numy", + "NQDHT", "knots", "logr", "model", From 9a4472c9f943574706f4f588c334fb6616d48db3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 27 Feb 2024 08:04:39 -0500 Subject: [PATCH 59/65] Added additional allowed flags for FlatDisk To BasisFactory [no ci] --- coefs/BiorthBasis.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/coefs/BiorthBasis.cc b/coefs/BiorthBasis.cc index a9e49f98a..b6622f643 100644 --- a/coefs/BiorthBasis.cc +++ b/coefs/BiorthBasis.cc @@ -1394,6 +1394,7 @@ namespace BasisClasses "numx", "numy", "numr", + "NQDHT", "knots", "logr", "model", @@ -1416,6 +1417,7 @@ namespace BasisClasses "Mmax", "nmax", "mmax", + "mlim", "dof", "subsamp", "samplesz", From 4754b0aa6bc3570c96ad728d67062a4ea2bf4979 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 27 Feb 2024 11:47:04 -0500 Subject: [PATCH 60/65] Make zeropos and zerovel the default options [no ci] --- utils/ICs/ZangICs.cc | 37 +++++++++++++++++++++++++++---------- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/utils/ICs/ZangICs.cc b/utils/ICs/ZangICs.cc index bece0fce8..638c6c755 100644 --- a/utils/ICs/ZangICs.cc +++ b/utils/ICs/ZangICs.cc @@ -38,7 +38,8 @@ main(int ac, char **av) options.add_options() ("h,help", "Print this help message") - ("z,zerovel", "Zero the mean velocity") + ("V,nozerovel", "Do not zero the mean velocity") + ("P,nozeropos", "Do not zero the center of mass") ("d,debug", "Print debug grid") ("N,number", "Number of particles to generate", cxxopts::value(N)->default_value("100000")) @@ -182,9 +183,9 @@ main(int ac, char **av) // Track number of iteration overflows int over = 0; - // Velocity zeroing + // Position and velocity zeroing // - std::vector> zerovel(nomp); + std::vector> zeropos(nomp), zerovel(nomp); // Generation loop with OpenMP // @@ -233,9 +234,12 @@ main(int ac, char **av) vel[n][1] = vr*sin(phi) + vt*cos(phi); vel[n][2] = 0.0; - // Accumulate mean velocity + // Accumulate mean position and velocity // - for (int k=0; k<3; k++) zerovel[tid][k] += vel[n][k]; + for (int k=0; k<3; k++) { + zeropos[tid][k] += pos[n][k]; + zerovel[tid][k] += vel[n][k]; + } // Print progress bar if (tid==0) ++(*progress); @@ -246,16 +250,29 @@ main(int ac, char **av) // double mass = (model->get_mass(Rmax) - model->get_mass(Rmin))/N; - // Reduce the mean velocity + // Reduce the mean position and velocity // for (int n=1; n Date: Tue, 27 Feb 2024 16:07:33 -0500 Subject: [PATCH 61/65] Fix for EVEN_M in BasisFactory for BirothCyl [no ci] --- coefs/BiorthBasis.cc | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/coefs/BiorthBasis.cc b/coefs/BiorthBasis.cc index b6622f643..5e32d51bc 100644 --- a/coefs/BiorthBasis.cc +++ b/coefs/BiorthBasis.cc @@ -1729,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(0, N1); n<=std::min(nmax-1, N2); n++) { From e3b663e6475f692fc8c890f33065445ee3296a73 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 Mar 2024 09:49:20 -0500 Subject: [PATCH 62/65] Updates to enforce the cache name requirement [no ci] --- exputil/EmpCylSL.cc | 2 +- include/EmpCylSL.H | 30 +++++++++++++++--------------- utils/ICs/eof_basis.cc | 6 ++++-- 3 files changed, 20 insertions(+), 18 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 90840b972..1db3164f2 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -3096,7 +3096,7 @@ void EmpCylSL::make_eof(void) // Cache table for restarts // - if (myid==0) cache_grid(1); + if (myid==0) cache_grid(1, cachefile); // Basis complete but still need to compute coefficients // diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 55fa787df..f9da0d718 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -235,7 +235,7 @@ protected: std::string cachefile; // 1=write, 0=read // return: 0=failure - int cache_grid(int, string file=""); + int cache_grid(int, std::string file); double integral(int, int, int, int); void get_pot(Eigen::MatrixXd&, Eigen::MatrixXd&, double, double); double massR(double R); @@ -495,8 +495,8 @@ public: */ EmpCylSL(int numr, int lmax, int mmax, int nord, - double ascale, double hscale, int Nodd=-1, - std::string cachename=""); + double ascale, double hscale, int Nodd, + std::string cachename); //! Destructor ~EmpCylSL(void); @@ -504,16 +504,16 @@ public: //! Reconstruct basis with new parameters void reset(int numr, int lmax, int mmax, int nord, double ascale, double hscale, int Nodd, - std::string cachename=""); + std::string cachename); //! Read EOF basis header from saved file - int read_eof_header(const string& eof_file); + int read_eof_header(const std::string& eof_file); //! Read EOF basis from saved file - int read_eof_file(const string& eof_file); + int read_eof_file(const std::string& eof_file); //! Dump the EOF basis in ascii format - void dump_eof_file(const string& eof_file, const string& dump_file); + void dump_eof_file(const std::string& eof_file, const std::string& dump_file); //! Read basis from cache file int read_cache(void); @@ -694,26 +694,26 @@ public: void dump_coefs_binary(std::ostream& out, double time); //! Plot basis - void dump_basis(const string& name, int step, double Rmax=-1.0); + void dump_basis(const std::string& name, int step, double Rmax=-1.0); //! Plot full fields for debugging - void dump_images(const string& OUTFILE, + void dump_images(const std::string& OUTFILE, double XYOUT, double ZOUT, int OUTR, int OUTZ, bool logscale); //! Plot basis images for debugging - void dump_images_basis(const string& OUTFILE, + void dump_images_basis(const std::string& OUTFILE, double XYOUT, double ZOUT, int OUTR, int OUTZ, bool logscale, int M1, int M2, int N1, int N2); //! Plot PCA basis images for debugging - void dump_images_basis_pca(const string& runtag, + void dump_images_basis_pca(const std::string& runtag, double XYOUT, double ZOUT, int OUTR, int OUTZ, int M, int N, int cnt); //! Plot EOF basis images for debugging - void dump_images_basis_eof(const string& runtag, + void dump_images_basis_eof(const std::string& runtag, double XYOUT, double ZOUT, int OUTR, int OUTZ, int M, int N, int cnt, const Eigen::VectorXd& tp); @@ -760,7 +760,7 @@ public: //! Set file name for EOF analysis and sample size for subsample //! computation - inline void setHall(string file, unsigned tot) + inline void setHall(std::string file, unsigned tot) { hallfile = file; nbodstot = tot; @@ -771,14 +771,14 @@ public: if (PCAVAR) { - const string types[] = + const std::string types[] = { "Hall", "Truncate", "None" }; - const string desc[] = + const std::string desc[] = { "Compute the S/N but do not modify coefficients", "Tapered signal-to-noise power defined by Hall", diff --git a/utils/ICs/eof_basis.cc b/utils/ICs/eof_basis.cc index 685b21b66..de38052c9 100644 --- a/utils/ICs/eof_basis.cc +++ b/utils/ICs/eof_basis.cc @@ -27,7 +27,7 @@ main(int argc, char **argv) std::string eof, tag; double rmax, zmax; - int nout, mmax, norder; + int nout, mmax, norder, nodd; cxxopts::Options options(argv[0], "Dump the entire disk orthgonal function file in ascii"); @@ -45,6 +45,8 @@ main(int argc, char **argv) cxxopts::value(mmax)->default_value("6")) ("n,nmax", "maximum radial order", cxxopts::value(norder)->default_value("18")) + ("nodd", "number of vertically antisymmtric functions", + cxxopts::value(nodd)->default_value("6")) ("N,nout", "number of grid points in each dimension", cxxopts::value(nout)->default_value("40")) ; @@ -83,7 +85,7 @@ main(int argc, char **argv) double acyl = 0.01; double hcyl = 0.001; - EmpCylSL test(nmax, lmax, mmax, nord, acyl, hcyl); + EmpCylSL test(nmax, lmax, mmax, nord, acyl, hcyl, nodd, eof); test.read_eof_file(eof); From cd04ea1d8cd61a993cd8cf61e3f300974406c91e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 Mar 2024 11:26:44 -0500 Subject: [PATCH 63/65] Change default nmaxfid to 128 in EmpCyl2d [no ci] --- exputil/BiorthCyl.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index f4f2abd1b..859dc13d9 100644 --- a/exputil/BiorthCyl.cc +++ b/exputil/BiorthCyl.cc @@ -40,7 +40,7 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) mmax = conf["Lmax"].as(); if (conf["nmaxfid"]) nmaxfid = conf["nmaxfid"].as(); - else nmaxfid = 40; + else nmaxfid = 256; if (conf["nmax"]) nmax = conf["nmax"].as(); else nmax = nmaxfid; From 614e98f4312e309b40f424b7056f285828070632 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 1 Mar 2024 11:43:49 -0500 Subject: [PATCH 64/65] Throw exception if yaml-cpp detects an error in parse.cc [no ci] --- src/parse.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/parse.cc b/src/parse.cc index 675a84572..160ac12a3 100644 --- a/src/parse.cc +++ b/src/parse.cc @@ -609,7 +609,9 @@ void YAML_parse_args(int argc, char** argv) if (done) { MPI_Finalize(); - exit(EXIT_SUCCESS); + throw EXPException("YAML configuration error", + "parsing failure, check your YAML config for syntax?", + __FILE__, __LINE__); } } From 9daee1613c36fbb0f4d45486c81b8e66e8327beb Mon Sep 17 00:00:00 2001 From: mdw Date: Fri, 1 Mar 2024 12:22:52 -0500 Subject: [PATCH 65/65] Remove MPI_Finalize() + exit() in favor of std::runtime_error exceptions in parse.cc [no ci] --- src/parse.cc | 47 +++++++++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/src/parse.cc b/src/parse.cc index 160ac12a3..3419a68d0 100644 --- a/src/parse.cc +++ b/src/parse.cc @@ -71,8 +71,8 @@ void initialize(void) << std::string(60, '-') << std::endl << parse << std::endl << std::string(60, '-') << std::endl; - MPI_Finalize(); - exit(-1); + throw EXPException("YAML configuration error", "Error parsing Global stanza", + __FILE__, __LINE__); } if (_G) { @@ -290,8 +290,8 @@ void initialize(void) MPI_Bcast(&nOK, 1, MPI_INT, 0, MPI_COMM_WORLD); if (nOK) { - MPI_Finalize(); - exit(10); + throw EXPException("Configuration error", "error opening outdir", + __FILE__, __LINE__); } } @@ -327,8 +327,8 @@ void initialize(void) MPI_Bcast(&iok, 1, MPI_INT, 0, MPI_COMM_WORLD); if (!iok) { - MPI_Finalize(); - exit(11); + throw EXPException("Configuration error", "error creating files in outdir", + __FILE__, __LINE__); } } @@ -382,8 +382,8 @@ void update_parm() catch (YAML::Exception & error) { if (myid==0) std::cout << "Error updating parameters in parse.update_parm: " << error.what() << std::endl; - MPI_Finalize(); - exit(-1); + throw EXPException("Configuration error", "error updating YAML config", + __FILE__, __LINE__); } } @@ -404,8 +404,8 @@ void write_parm(void) MPI_Bcast(&nOK, 1, MPI_INT, 0, MPI_COMM_WORLD); if (nOK) { - MPI_Finalize(); - exit(102); + throw EXPException("Configuration error", "error writing YAML config", + __FILE__, __LINE__); } if (myid!=0) return; @@ -502,8 +502,8 @@ void YAML_parse_args(int argc, char** argv) vm = options.parse(argc, argv); } catch (cxxopts::OptionException& e) { std::cout << "Option error: " << e.what() << std::endl; - MPI_Finalize(); - exit(-1); + throw EXPException("cxxopts error", "error parsing options", + __FILE__, __LINE__); } if (vm.count("help")) { @@ -530,8 +530,7 @@ void YAML_parse_args(int argc, char** argv) MPI_Bcast(&done, 1, MPI_INT, 0, MPI_COMM_WORLD); if (done) { - MPI_Finalize(); - exit(EXIT_SUCCESS); + throw EXPException("Clean exit", "done", __FILE__, __LINE__); } // If config file is not specified, use first trailing, unmatched @@ -573,8 +572,9 @@ void YAML_parse_args(int argc, char** argv) MPI_Bcast(&done, 1, MPI_INT, 0, MPI_COMM_WORLD); if (done) { - MPI_Finalize(); - exit(EXIT_SUCCESS); + throw EXPException("Configuration error", + "error opening YAML configuration file", + __FILE__, __LINE__); } try { @@ -590,6 +590,14 @@ void YAML_parse_args(int argc, char** argv) // MPI_Bcast(&done, 1, MPI_INT, 0, MPI_COMM_WORLD); + // Exit if done + // + if (done) { + throw EXPException("YAML configuration error", + "parsing failure, check your YAML config for syntax?", + __FILE__, __LINE__); + } + // The current YAML structure will now be serialized and broadcast // to all processes // @@ -608,13 +616,12 @@ void YAML_parse_args(int argc, char** argv) MPI_Bcast(&done, 1, MPI_INT, 0, MPI_COMM_WORLD); if (done) { - MPI_Finalize(); - throw EXPException("YAML configuration error", - "parsing failure, check your YAML config for syntax?", + throw EXPException("Configuration error in main process", + "exiting", __FILE__, __LINE__); } } - + // Receive the YAML serialized length // int line_size;