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/coefs/BiorthBasis.cc b/coefs/BiorthBasis.cc index 18427eb08..5e32d51bc 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; @@ -1394,6 +1394,7 @@ namespace BasisClasses "numx", "numy", "numr", + "NQDHT", "knots", "logr", "model", @@ -1405,7 +1406,10 @@ namespace BasisClasses "NO_M0", "NO_M1", "EVEN_M", + "M0_BACK", "M0_ONLY", + "NO_MONO", + "diskconf", "ssfrac", "playback", "coefMaster", @@ -1413,6 +1417,7 @@ namespace BasisClasses "Mmax", "nmax", "mmax", + "mlim", "dof", "subsamp", "samplesz", @@ -1456,6 +1461,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(); @@ -1722,18 +1729,20 @@ namespace BasisClasses } // Get the basis fields + // ortho->get_dens (dend[tid], R, z); ortho->get_pot (potd[tid], R, z); ortho->get_rforce (potR[tid], R, z); ortho->get_zforce (potZ[tid], R, z); // m loop + // for (int m=0, moffset=0; m<=mmax; m++) { - if (m==0 and NO_M0) continue; - if (m==1 and NO_M1) continue; - if (EVEN_M and m/2*2 != m) continue; - if (m>0 and M0_only) break; + if (m==0 and NO_M0) { moffset++; continue; } + if (m==1 and NO_M1) { moffset += 2; continue; } + if (EVEN_M and m/2*2 != m) { moffset += 2; continue; } + if (m>0 and M0_only) break; if (m==0) { for (int n=std::max(0, N1); n<=std::min(nmax-1, N2); n++) { @@ -1850,7 +1859,7 @@ namespace BasisClasses } // Radial grid spacing - double dx = (logxmax - logxmin)/numgrid; + double dx = (logxmax - logxmin)/(numgrid-1); // Basis storage Eigen::MatrixXd tabpot, tabden, tabrfc; 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 diff --git a/exputil/BiorthCyl.cc b/exputil/BiorthCyl.cc index 24e90f634..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; @@ -75,11 +75,8 @@ 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["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; @@ -101,6 +98,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: " @@ -116,6 +116,11 @@ BiorthCyl::BiorthCyl(const YAML::Node& conf) : conf(conf) geometry = "cylinder"; forceID = "BiorthCyl"; + + // Hardwired for testing. Could be parameters. + // + logr = false; + biorth = "bess"; initialize(); @@ -158,15 +163,9 @@ void BiorthCyl::initialize() void BiorthCyl::create_tables() { - // Hardwired for testing. Could be variables. - // - bool logr = false; - std::string target("expon"); - std::string biorth("bess"); - - emp = EmpCyl2d(mmax, nmaxfid, knots, numr, - rcylmin*scale, rcylmax*scale, acyltbl*scale, scale, - cmapR, logr, target, biorth); + emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, + rcylmin*scale, rcylmax*scale, + scale, cmapR, logr, diskconf, biorth); if (conf["basis"]) emp.basisTest(true); @@ -476,7 +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 ("cmapR", HighFive::DataSpace::From(cmapR)).write(cmapR); file.createAttribute ("cmapZ", HighFive::DataSpace::From(cmapZ)).write(cmapZ); } @@ -597,7 +595,7 @@ bool BiorthCyl::ReadH5Cache() double v; HighFive::Attribute vv = h5file.getAttribute(name); vv.read(v); if (fabs(value - v) < 1.0e-16) return true; else { - if (myid==0) std::cout << "ReadH5Cache(): wanted " + if (myid==0) std::cout << "---- BiortyCyl::ReadH5Cache: wanted " << name << "=" << value << " but found " << name << "=" << v << std::endl; @@ -635,7 +633,6 @@ bool BiorthCyl::ReadH5Cache() if (not checkDbl(rcylmin, "rcylmin")) return false; if (not checkDbl(rcylmax, "rcylmax")) return false; if (not checkDbl(scale, "scale")) return false; - if (not checkDbl(acyltbl, "acyltbl")) return false; if (not checkInt(cmapR, "cmapR")) return false; if (not checkInt(cmapZ, "cmapZ")) return false; @@ -649,6 +646,14 @@ bool BiorthCyl::ReadH5Cache() std::cerr << "---- BiorthCyl::ReadH5Cache: " << "read <" << cachename << ">" << std::endl; + // Create the basis instance + // + emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, + rcylmin*scale, rcylmax*scale, + scale, cmapR, logr, diskconf, biorth); + + if (conf["basis"]) emp.basisTest(true); + return true; } catch (HighFive::Exception& err) { @@ -755,7 +760,6 @@ YAML::Node BiorthCyl::getHeader(const std::string& cachefile) node["rcylmin"] = getDbl("rcylmin"); node["rcylmax"] = getDbl("rcylmax"); node["scale"] = getDbl("scale"); - node["acyltbl"] = getDbl("acyltbl"); node["cmapR"] = getInt("cmapR"); node["cmapZ"] = getInt("cmapZ"); } @@ -797,9 +801,9 @@ 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, - cmapR, false, "expon", "bess"); + emp = EmpCyl2d(mmax, nmaxfid, nmax, knots, numr, + rcylmin*scale, rcylmax*scale, scale, + cmapR, logr, diskconf, biorth); } return emp.orthoCheck(); 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/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index 3c2152076..9247bd0d3 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -336,12 +336,40 @@ class EmpCyl2d::ExponCyl : public EmpCyl2d::ModelCyl private: - double sigma0; + 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(double scl) - { acyl = scl; 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; @@ -365,9 +393,41 @@ class EmpCyl2d::ExponCyl : public EmpCyl2d::ModelCyl class EmpCyl2d::KuzminCyl : public EmpCyl2d::ModelCyl { +private: + + 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(double scl) { acyl = scl; id = "kuzmin"; } + KuzminCyl(const YAML::Node& par) + { + parse(par); + id = "kuzmin"; + } double pot(double R) { double a2 = acyl * acyl; @@ -392,63 +452,247 @@ class EmpCyl2d::KuzminCyl : public EmpCyl2d::ModelCyl class EmpCyl2d::MestelCyl : public EmpCyl2d::ModelCyl { +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(double scl) { acyl = scl; id = "mestel"; } + MestelCyl(const YAML::Node& par) + { + parse(par); + 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"); } }; -std::shared_ptr -EmpCyl2d::createModel(const std::string type, double acyl) +class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl { - // 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); }); + +private: + //! Parameters + double mu, nu, ri, ro; - 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); + //! Softening factor (not currently used) + double asoft = 1.0e-8; + + //! Ignore inner cut-off for N<0.05 + bool Inner = true; + + //! Taper factors + double Tifac, Tofac; + + //! Inner taper function + double Tinner(double Jp) + { + double fac = pow(Jp, nu); + return fac/(Tifac + fac); } - 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); + //! Outer taper function + double Touter(double Jp) + { + return 1.0/(1.0 + pow(Jp/Tofac, mu)); } - 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); + //! Deriv of inner taper function + double dTinner(double Jp) + { + double fac = pow(Jp, nu); + double fac2 = Tifac + fac; + return nu*fac/Jp/(fac2*fac2); + } + + //! Deriv of outer taper function + double dTouter(double Jp) + { + double fac = pow(Jp/Tofac, mu); + double fac2 = 1.0 + fac; + return -mu*fac/Jp/(fac2*fac2); + } + +protected: + + //! Assign values from YAML + void parse(const YAML::Node& conf) + { + try { + 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; + + if (conf["Ro"]) + ro = conf["Ro"].as(); + else + ro = 10.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 YAML::Node& par) : MestelCyl(par) + { + // Parse the YAML + parse(par); + + // Assign the id + id = "zang"; + + // Cache taper factors + Tifac = pow(ri*vrot, nu); + Tofac = ro*vrot; + + if (nu<0.05) { + // Exponent is now for mapping only + Inner = false; + } + } + + //! Surface density + double dens(double R) + { + double ret = MestelCyl::dens(R) * Touter(R); + if (Inner) ret *= Tinner(R); + return ret; + } + +}; + +std::shared_ptr EmpCyl2d::createModel() +{ + 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'"); + + // 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 (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 (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 (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]" << std::endl; - return std::make_shared(acyl); + std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk [Default]" + << 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); + } } @@ -516,47 +760,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, - 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), - model(model), biorth(biorth), cache_name_2d(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), 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, acyl); - basis = Basis2d::createBasis(mmax, nmax, rmax, biorth); + disk = createModel(); + 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, - 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), +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), 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 +819,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 +835,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 +861,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 +875,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)); @@ -759,10 +1007,15 @@ 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); - 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 +1023,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(params)).write(params); file.createAttribute("model", HighFive::DataSpace::From(model)). write(model); file.createAttribute("biorth", HighFive::DataSpace::From(biorth)).write(biorth); @@ -830,13 +1083,19 @@ bool EmpCyl2d::ReadH5Cache() 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; 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 +1104,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 checkStr(params, "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 +1278,7 @@ void EmpCyl2d::checkCoefs() if (myid) return; Mapping map(scale, cmap); - auto disk = createModel(model, acyl); + auto disk = createModel(); LegeQuad lw(knots); Eigen::VectorXd coefs(nmax), coef0(nmax); diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 55fa4a0c0..1db3164f2 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) { @@ -3094,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/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/exputil/cudaSLGridMP2.cu b/exputil/cudaSLGridMP2.cu index 58202f19f..8d51ac338 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; @@ -86,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); @@ -110,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; - cuda_safe_call(cudaCreateTextureObject(&tex[0], &resDesc[0], &texDesc[0], NULL), __FILE__, __LINE__, "create texture object"); + memset(&resDesc, 0, sizeof(cudaResourceDesc)); + resDesc.resType = cudaResourceTypeArray; + resDesc.res.array.array = cuArray[0]; + + 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/exputil/mestel.cc b/exputil/mestel.cc new file mode 100644 index 000000000..926edbc24 --- /dev/null +++ b/exputil/mestel.cc @@ -0,0 +1,205 @@ +#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/Tofac, mu)); +} + +double TaperedMestelDisk::dTinner(double Jp) +{ + double fac = pow(Jp, nu); + double fac2 = Tifac + fac; + return nu*fac/Jp/(fac2*fac2); +} + +double TaperedMestelDisk::dTouter(double Jp) +{ + double fac = pow(Jp/Tofac, mu); + double fac2 = 1.0 + fac; + return -mu*fac/Jp/(fac2*fac2); +} + +TaperedMestelDisk::TaperedMestelDisk(double nu, double mu, double Ri, double Ro, + double VROT, double RMIN, double RMAX) + : nu(nu), mu(mu), Ri(Ri), Ro(Ro), MestelDisk(VROT, RMIN, RMAX) +{ + Tifac = pow(Ri*vrot, nu); + Tofac = Ro*vrot; + 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) * 2.0*M_PI*rr; + 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; +} diff --git a/exputil/orbit.cc b/exputil/orbit.cc index 1d08a6618..7e5637308 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; } @@ -49,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(); @@ -73,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(); @@ -120,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) @@ -156,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; @@ -168,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/exputil/orbit_trans.cc b/exputil/orbit_trans.cc index db777cfff..35edbd32f 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,125 @@ 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 +192,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 +214,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 +338,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 +354,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 +381,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 +392,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 +436,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 +493,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/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; 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()); diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index 790201aa4..da3cd51f5 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -41,14 +41,14 @@ class BiorthCyl protected: - YAML::Node conf; + 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; + double rcylmin, rcylmax, scale, acyltbl, acylcut, Ninner, Mouter; - bool EVEN_M, verbose; + bool EVEN_M, verbose, logr; //@{ //! Grid parameters @@ -79,8 +79,8 @@ 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; //! Active cache file name std::string cachename; @@ -198,6 +198,13 @@ 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 std::tuple background(double r, double z) + { + auto [p, dr, d] = emp.background(r); + return {p, dr, 0.0}; + } + //! Get the table bounds double getRtable() { return rcylmax*scale; } diff --git a/include/EmpCyl2d.H b/include/EmpCyl2d.H index 8a12c6892..256c070cf 100644 --- a/include/EmpCyl2d.H +++ b/include/EmpCyl2d.H @@ -5,12 +5,14 @@ #include #include #include +#include #include #include #include #include +#include // Parameters parsing #include // For config macros // For reading and writing cache file @@ -26,13 +28,14 @@ class EmpCyl2d { public: + //! A two-dimensional disk model for computing the EOF class ModelCyl { protected: - double acyl; std::string id; + virtual void parse(const YAML::Node&) = 0; public: @@ -45,9 +48,12 @@ public: protected: - double acyl, rmin, rmax, scale; + //! Contains parameter database + YAML::Node Params; + + 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 +71,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 @@ -82,12 +88,12 @@ protected: std::shared_ptr disk; std::shared_ptr basis; - static std::shared_ptr - createModel(const std::string type, double A); + std::shared_ptr createModel(); 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, 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 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, 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. @@ -198,6 +204,20 @@ 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); + } + + std::tuple background(double r) + { + return {disk->pot(r), disk->dpot(r), disk->dens(r)}; + } + //@} }; diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 9c33136ff..f9da0d718 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(); @@ -238,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); @@ -498,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); @@ -507,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); @@ -697,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); @@ -763,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; @@ -774,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/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 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 77db06360..fa5541a22 100644 --- a/include/mestel.H +++ b/include/mestel.H @@ -1,10 +1,14 @@ #ifndef _mestel_H #define _mestel_H +#include + +//! Infiinte Mestel Disk class MestelDisk : public AxiSymModel { -private: +protected: + double vrot; double rot; double rmin; double rmax; @@ -14,88 +18,101 @@ private: public: - 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"; } - + //! Constructor + MestelDisk(double VROT = 1.0, double RMIN = 1.0e-6, double RMAX = 1.0e6); + + //@{ + //! Required member functions + virtual double get_mass(const double R); + + virtual double get_density(const double R); - // 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) { - 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; - } + 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); + //@} - // 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) { - 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) { - L = fabs(L); - if (L==0.0) return 0.0; - return F*pow(L, q) * exp(-E/sig2); - } - - 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) { - 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) { - 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); + //@} }; +//! Doubly tapered Mestel Disk +class TaperedMestelDisk : public MestelDisk +{ +private: + + std::shared_ptr interp; + +protected: + //! Taper exponents + double nu, mu; + + //! Inner taper radius + double Ri; + + //! Outer taper radius + double Ro; + + //! Taper factors + double Tifac, Tofac; + + //! Inner taper function + double Tinner(double Jp); + + //! Outer taper function + double Touter(double Jp); + + //! Deriv of inner taper function + double dTinner(double Jp); + + //! Deriv of outer taper function + double dTouter(double Jp); + +public: + + //! Constructor + TaperedMestelDisk(double nu, double mu, double Ri, double Ro, + double VROT = 1.0, + double RMIN = 1.0e-6, double RMAX = 1.0e6); + + //! Overloaded density function + 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); + + double dfde(double E, double L); + + double dfdl(double E, double L); + + double d2fde2(double E, double L); + //@} + +}; #endif 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] 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 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.H b/src/FlatDisk.H index ad8ae58df..64be763a4 100644 --- a/src/FlatDisk.H +++ b/src/FlatDisk.H @@ -74,6 +74,14 @@ private: } #endif + //! Background evaluation + 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); void get_potl_dens(double r, double z, diff --git a/src/FlatDisk.cc b/src/FlatDisk.cc index 6875ef853..ba5ab0ec2 100644 --- a/src/FlatDisk.cc +++ b/src/FlatDisk.cc @@ -12,14 +12,15 @@ FlatDisk::valid_keys = { "numr", "rcylmin", "rcylmax", - "acyltbl", "mmax", "numx", "numy", + "NQDHT", "knots", "logr", "model", "biorth", + "diskconf", "cachename" }; @@ -30,7 +31,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; @@ -46,7 +46,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 @@ -72,7 +71,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/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 // diff --git a/src/PolarBasis.H b/src/PolarBasis.H index aecffe886..caab18f0d 100644 --- a/src/PolarBasis.H +++ b/src/PolarBasis.H @@ -32,9 +32,11 @@ 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) @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 @@ -265,9 +267,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) @@ -376,6 +378,10 @@ protected: #endif + //! Background evaluation + 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); @@ -394,9 +400,12 @@ protected: //! Flag self_consitency bool self_consistent; - //! Flag fixed monopole + //! Omit monopole bool NO_M0; + //! Omit off-grid force algorithm + bool NO_MONO; + //! Flag drop m=1 term bool NO_M1; @@ -406,6 +415,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..4e7926c4f 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -35,11 +35,14 @@ PolarBasis::valid_keys = { "scale", "rmin", "rmax", + "Mmax", "self_consistent", "NO_M0", "NO_M1", "EVEN_M", "M0_ONLY", + "M0_BACK", + "NO_MONO", "mlim", "ssfrac", "playback", @@ -69,6 +72,8 @@ PolarBasis::PolarBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : NO_M1 = false; EVEN_M = false; M0_only = false; + M0_back = false; + NO_MONO = false; ssfrac = 0.0; subset = false; coefMaster = true; @@ -104,13 +109,18 @@ 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["NO_MONO"]) NO_MONO = conf["NO_MONO"].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; @@ -718,7 +728,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 @@ -1396,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) { @@ -1406,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) { @@ -1419,11 +1434,21 @@ 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) { + auto [p, drc, dzc] = get_pot_background(r, zz); + + 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]); + + potl = mfac*norm0 * p; + potr = mfac*norm0 * drc; + potz = mfac*norm0 * dzc; + } } // Asymmetric terms? @@ -1432,7 +1457,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..95d4783ec 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 cachename"); 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 diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index 9be032270..218e70177 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -53,6 +53,36 @@ 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; +} + void BiorthCyl::initialize_cuda (std::vector& cuArray, thrust::host_vector& tex @@ -60,7 +90,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 // @@ -71,15 +103,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; @@ -94,6 +117,15 @@ void BiorthCyl::initialize_cuda for (size_t n=0; n> 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 + + // Size of interpolation array + // + size_t tsize = numr*sizeof(cuFP_t); + + // Make the textures + // + for (int n=0; n::epsilon(); struct Element { @@ -221,4 +322,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; + } + + } 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; iminDT.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; diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index c3299d658..226252aaf 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -17,14 +17,17 @@ __device__ __constant__ cuFP_t plrRscale, plrHscale, plrXmin, plrXmax, plrYmin, plrYmax, plrDxi, plrDyi; +__device__ __constant__ +cuFP_t plrDx0; + __device__ __constant__ cuFP_t plrCen[3], plrBody[9], plrOrig[9]; __device__ __constant__ -int plrNumx, plrNumy, plrCmapR, plrCmapZ, plrOrient; +int plrNumx, plrNumy, plrCmapR, plrCmapZ, plrOrient, plrNumr; __device__ __constant__ -bool plrAcov, plrNO_M0, plrNO_M1, plrEVEN_M, plrM0only; +bool plrAcov, plrNO_M0, plrNO_M1, plrEVEN_M, plrM0only, plrM0back, plrNoMono; // Index function for sine and cosine coefficients // @@ -75,8 +78,10 @@ void testConstantsPlr() printf(" Ymax = %f\n", plrYmax ); printf(" Dxi = %f\n", plrDxi ); printf(" Dyi = %f\n", plrDyi ); + 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); @@ -84,6 +89,8 @@ void testConstantsPlr() printf(" NO_M1 = %d\n", plrNO_M1 ); printf(" EVEN_M = %d\n", plrEVEN_M); printf(" M0only = %d\n", plrM0only); + printf(" M0back = %d\n", plrM0back); + printf(" NoMono = %d\n", plrNoMono); printf("-------------------------\n"); } @@ -241,8 +248,22 @@ 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(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), + __FILE__, __LINE__, "Error copying plrM0back"); + cuda_safe_call(cudaMemcpyToSymbol(plrAcov, &subsamp, sizeof(bool), size_t(0), cudaMemcpyHostToDevice), __FILE__, __LINE__, "Error copying plrAcov"); + + cuda_safe_call(cudaMemcpyToSymbol(plrNumr, &f.numr, sizeof(int), size_t(0), cudaMemcpyHostToDevice), + __FILE__, __LINE__, "Error copying plrNumr"); + + cuFP_t Dx0 = (f.xmax - f.xmin)/(f.numr - 1); + cuda_safe_call(cudaMemcpyToSymbol(plrDx0, &Dx0, sizeof(cuFP_t), size_t(0), cudaMemcpyHostToDevice), + __FILE__, __LINE__, "Error copying plrDx0"); + } @@ -596,7 +617,9 @@ forceKernelPlr6(dArray 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) { @@ -1071,7 +1094,9 @@ forceKernelPlr3(dArray 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) { @@ -1136,83 +1161,140 @@ 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; // Compute the force for this m value + // by default - if (mm) norm = norm1; - else norm = norm0; + // 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; + if (plrEVEN_M and (mm/2)*2 != mm) compute = false; - for (int n=0; nplrNumr-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 - 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 + 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 - 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 + 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; + fr += dp0; + } else { + pp += a*pp0 + b*pp1; + fr += a*dp0 + b*dp1; + } + + } + // END: unperturbed fixed background model for m=0 + // BEGIN: use basis function expansion + else { + + 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 +#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 +#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; + } + // END: radial-order loop } - - 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: basis-function contribution block } - + // END: compute block + // Trig recursion to squeeze avoid internal FP fct call // cuFP_t cosM = ccos; @@ -2553,18 +2635,18 @@ void PolarBasis::multistep_update_cuda() void PolarBasis::destroy_cuda() { + // Deallocate textures + // for (size_t i=0; i > 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++) { diff --git a/src/parse.cc b/src/parse.cc index 675a84572..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,11 +616,12 @@ 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 in main process", + "exiting", + __FILE__, __LINE__); } } - + // Receive the YAML serialized length // int line_size; 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..638c6c755 --- /dev/null +++ b/utils/ICs/ZangICs.cc @@ -0,0 +1,318 @@ +/* + A tapered Mestel disk IC generator +*/ + // C++/STL headers +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include // Progress bar +#include // Option parsing + +int +main(int ac, char **av) +{ + //===================== + // Begin option parsing + //===================== + + int N; // Number of particles + double mu, nu, Ri, Ro; // 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") + ("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")) + ("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("1.0")) + ("o,Ro", "Outer radius for taper", + cxxopts::value(Ro)->default_value("20.0")) + ("r,Rmin", "Inner radius for model", + cxxopts::value(Rmin)->default_value("0.001")) + ("R,Rmax", "Outer radius for model", + 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", + cxxopts::value(seed)) + ("f,file", "Output body file", + cxxopts::value(bodyfile)->default_value("zang.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); + } + + SphericalOrbit::ZFRAC=0.3; // TEST + + // Create the model + // + auto model = std::make_shared(nu, mu, Ri, Ro, + 1.0, Rmin, Rmax); + model->setup_df(sigma); + + // 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 + // + 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(Rmin); + // 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 + // + 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; + for (int i=0; i<=numE; i++) { + double E = Emin + dE*i; + 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); + + 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); + + // Maximum number rejection-method iterations + int itmax = 100000; + + // Track number of iteration overflows + int over = 0; + + // Position and velocity zeroing + // + std::vector> zeropos(nomp), zerovel(nomp); + + // Generation loop with OpenMP + // +#pragma omp parallel for reduction(+:over) + for (int n=0; nnew_orbit(E, K); + double F = model->distf(E, orb[tid]->get_action(1)) / orb[tid]->get_freq(0); + if (F/peak > uniform(gen)) break; + } + + if (j==itmax) over++; + + 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))); + // ^ + // | + // For sanity-----+ + + 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; + + // Accumulate mean position and velocity + // + 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); + } + std::cout << std::endl << "** Main loop complete" << std::endl; + + // Compute the particle mass + // + double mass = (model->get_mass(Rmax) - model->get_mass(Rmin))/N; + + // Reduce the mean position and velocity + // + for (int n=1; n(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); diff --git a/utils/SL/EOF2d.cc b/utils/SL/EOF2d.cc index 32b207c01..9db5a1e35 100644 --- a/utils/SL/EOF2d.cc +++ b/utils/SL/EOF2d.cc @@ -6,24 +6,26 @@ #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, nmax, knots, M, N, nradial, nout; - double A, rmin, rmax; - std::string filename, type, biorth; + int numr, mmax, nmaxfid, nmax, knots, M, N, nradial, nout; + double A, rmin, rmax, rout; + std::string filename, config, 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() @@ -42,33 +44,38 @@ 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,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", + ("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)", - cxxopts::value(type)->default_value("expon")) - ("biorth", "Biorthogonal type (cb, bess)", + ("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", + ("o,filename", "Output filename", cxxopts::value(filename)->default_value("testeof")) ; @@ -101,10 +108,14 @@ int main(int argc, char** argv) // if (vm.count("cmap")) cmap = true; + // Parameter vector for the EmpCyl2d models + // + YAML::Node par = YAML::Load(config); + // Make the class instance // - EmpCyl2d emp(mmax, nmax, knots, numr, rmin, rmax, A, 1.0, 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; @@ -117,6 +128,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")) { @@ -134,20 +151,14 @@ 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 // 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); @@ -207,13 +218,22 @@ int main(int argc, char** argv) for (int i=0; i