diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 6b8bc09eb..98d429c5f 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1206,9 +1206,15 @@ namespace BasisClasses //! Number of integration knots for orthogonal check int knots; + //! Even/odd exclusion flags + bool no_even = false, no_odd = false; + //! SLGridSlab mesh size int ngrid = 1000; + //! Coordinate map type for SLGridSlab + std::string cmap = "linear"; + //! Target model type for SLGridSlab std::string type = "isothermal"; @@ -1270,6 +1276,20 @@ namespace BasisClasses void computeAccel(double x, double y, double z, Eigen::Ref vstore); + //@{ + //! Covariance structures. First index is T, second is the + //! flattened 3-d k vector + std::vector meanV; + std::vector covrV; + std::vector dvarV; + + void init_covariance(); + void zero_covariance(); + void writeCovarH5Params(HighFive::File& file); + int sampT = 100; + bool diagcov = true; + //@} + public: //! Constructor from YAML node @@ -1326,6 +1346,64 @@ namespace BasisClasses return true; } + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) + { + if (covarStore) { + SubsampleCovariance::CovarData elem; + std::get<0>(elem) = sampleCounts; + std::get<1>(elem) = sampleMasses; + + int sampT = meanV.size(); + int Nsize = meanV[0].size(); + + std::get<2>(elem) = SubsampleCovariance::CoefType(sampT, 1, Nsize); + std::get<3>(elem) = SubsampleCovariance::CovrType(sampT, 1, Nsize, Nsize); + + for (int T=0; T(elem)(T, 0, n1) = meanV[T](n1); + for (int n2=0; n2 cov = 0.0; + if (diagcov) { + if (n1 == n2) cov = dvarV[T](n1); + } else if (covar && T < covrV.size()) { + cov = covrV[T](n1, n2); + } + std::get<3>(elem)(T, 0, n1, n2) = cov; + } + } + } + covarStore->writeCoefCovariance(compname, runtag, elem, time); + } else { + throw std::runtime_error("Slab::writeCoefCovariance: covariance storage not initialized"); + } + } + + + //! Enable coefficient covariance computation + virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, + bool ftype=false, + bool covr_tot=false, + bool covar_in=false) + { + pcavar = pcavar_in; + sampT = sampT_in; + scovr = covr_tot; + covar = covar_in; + if (pcavar) { + init_covariance(); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, scovr, covar); + } + } + + //! Return wave-number indices from flattened index + std::tuple index3D(unsigned k); + + //! Get flattened index from wave-number indicies + unsigned index1D(int kx, int ky, int kz); + }; /** diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 9c6d21a30..7546cafcb 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3833,21 +3833,36 @@ namespace BasisClasses "hslab", "zmax", "ngrid", + "cmap", "type", "knots", "verbose", - "check", - "method" + "method", + "no_even", + "no_odd", + "self_consistent", + "pcavar", + "samplesz", + "orthochk", + "cachename" }; Slab::Slab(const YAML::Node& CONF) : BiorthBasis(CONF, "slab") { initialize(); + + // Initialize covariance + // + if (pcavar) enableCoefCovariance(true, sampT); } Slab::Slab(const std::string& confstr) : BiorthBasis(confstr, "slab") { initialize(); + + // Initialize covariance + // + if (pcavar) enableCoefCovariance(true, sampT); } void Slab::initialize() @@ -3893,11 +3908,20 @@ namespace BasisClasses if (conf["hslab"]) hslab = conf["hslab"].as(); if (conf["zmax" ]) zmax = conf["zmax" ].as(); if (conf["ngrid"]) ngrid = conf["ngrid"].as(); + if (conf["cmap" ]) cmap = conf["cmap" ].as(); if (conf["type" ]) type = conf["type" ].as(); + if (conf["no_even"]) no_even = conf["no_even"].as(); + if (conf["no_odd"] ) no_odd = conf["no_odd" ].as(); + if (conf["knots"]) knots = conf["knots"].as(); - if (conf["check"]) check = conf["check"].as(); + if (conf["orthochk"]) check = conf["orthochk"].as(); + + if (conf["pcavar"]) pcavar = conf["pcavar"].as(); + if (conf["subsamp"]) sampT = conf["samplesz"].as(); + + if (conf["cachename"]) cachename = conf["cachename"].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -3910,6 +3934,19 @@ namespace BasisClasses throw std::runtime_error("Slab: error parsing YAML"); } + // Check for non-null cache file name. This must be specified + // to prevent recomputation and unexpected behavior. + // + if (cachename.size() == 0) { + throw std::runtime_error + ("SlabSL requires a specified cachename in your YAML config\n" + "for consistency with previous invocations and existing coefficient\n" + "sets. Please add explicitly add 'cachename: name' to your config\n" + "with new 'name' for creating a basis or an existing 'name' for\n" + "reading a previously generated basis cache\n"); + } + + // Finally, make the basis // SLGridSlab::mpi = 0; @@ -3919,7 +3956,8 @@ namespace BasisClasses int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy; - ortho = std::make_shared(nnmax, nmaxz, ngrid, zmax, type); + ortho = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, + cmap, type); // Orthogonality sanity check // @@ -4029,9 +4067,18 @@ namespace BasisClasses // Storage for basis evaluation Eigen::VectorXd zpot(nmaxz); + Eigen::VectorXcd g; + if (pcavar) { + g.resize(imx*imy*imz); + g.setZero(); + } + // Loop indices int ix, iy; + // Leading constants for biorthogonality + const double norm = -4.0*M_PI; + // Recursion multipliers std::complex stepx = exp(-kfac*x), facx; std::complex stepy = exp(-kfac*y), facy; @@ -4066,13 +4113,25 @@ namespace BasisClasses ortho->get_pot(zpot, z, iiy, iix); for (int iz=0; izorthoCheck(); } + void Slab::init_covariance() + { + if (pcavar) { + + int jmax = imx * imy * imz; + + meanV.resize(sampT); + for (auto& v : meanV) { + v.resize(jmax); + } + + if (covar) { + if (diagcov) { + dvarV.resize(sampT); + for (auto& v : dvarV) { + v.resize(jmax); + } + } else { + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(jmax, jmax); + } + } + } else { + covrV.clear(); + dvarV.clear(); + } + + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + zero_covariance(); + } + } + + + void Slab::zero_covariance() + { + for (int T=0; T 2*nmaxx) { + std::ostringstream sout; + sout << "Slab::index1D: x index [" << kx << "] must be in [0, " + << 2*nmaxx << "]"; + throw std::runtime_error(sout.str()); + } + + if (ky < 0 or ky > 2*nmaxy) { + std::ostringstream sout; + sout << "Slab::index1D: y index [" << ky << "] must be in [0, " + << 2*nmaxy << "]"; + throw std::runtime_error(sout.str()); + } + + if (kz < 0 or kx >= nmaxz) { + std::ostringstream sout; + sout << "Slab::index1D: z index [" << kz << "] must be in [0, " + << nmaxz-1 << "]"; + throw std::runtime_error(sout.str()); + } + + return + kx*(2*nmaxy+1)*nmaxz + + ky*nmaxz + + kz; + } + + std::tuple Slab::index3D(unsigned indx) + { + // Sanity check + // + if (indx >= imx*imy*imz) { + std::ostringstream sout; + sout << "Slab::index3D: index [" << indx << "] must be in 0 <= indx < " << imx*imy*imz; + throw std::runtime_error(sout.str()); + } + + // Compute the 3d index + // + int ix = indx/(imy*imz); + int iy = (indx - ix*imy*imz)/imz; + int iz = indx - ix*imy*imz - iy*imz; + + return {ix, iy, iz}; + } + const std::set Cube::valid_keys = { "nminx", @@ -4390,7 +4552,7 @@ namespace BasisClasses "nmaxz", "knots", "verbose", - "check", + "orthochk", "method", "pcavar", "subsamp", @@ -4463,7 +4625,7 @@ namespace BasisClasses if (conf["knots"]) knots = conf["knots" ].as(); - if (conf["check"]) check = conf["check" ].as(); + if (conf["orthochk"]) check = conf["orthochk"].as(); if (conf["pcavar"]) pcavar = conf["pcavar" ].as(); if (conf["subsamp"]) sampT = conf["subsamp" ].as(); @@ -4948,7 +5110,7 @@ namespace BasisClasses // int ix = indx/((2*nmaxy+1)*(2*nmaxz+1)); int iy = (indx - ix*(2*nmaxy+1)*(2*nmaxz+1))/(2*nmaxz+1); - int iz = indx - ix*(2*nmaxy+1)*(2*nmaxz+1)/(2*nmaxz+1) - iy*(2*nmaxz+1); + int iz = indx - ix*(2*nmaxy+1)*(2*nmaxz+1) - iy*(2*nmaxz+1); return {ix, iy, iz}; } @@ -4966,6 +5128,8 @@ namespace BasisClasses coef = std::make_shared(); else if (name.compare("flatdisk") == 0) coef = std::make_shared(); + else if (name.compare("slab") == 0) + coef = std::make_shared(); else if (name.compare("cube") == 0) coef = std::make_shared(); else { @@ -5031,6 +5195,10 @@ namespace BasisClasses coefret = std::make_shared(); else if (name.compare("flatdisk") == 0) coefret = std::make_shared(); + else if (name.compare("slab") == 0) + coefret = std::make_shared(); + else if (name.compare("cube") == 0) + coefret = std::make_shared(); else { std::ostringstream sout; sout << "Basis::createCoefficients: basis <" << name << "> not recognized" @@ -5658,6 +5826,16 @@ namespace BasisClasses file.createAttribute("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax); } + void Slab::writeCovarH5Params(HighFive::File& file) + { + file.createAttribute("nminx", HighFive::DataSpace::From(nminx)).write(nminx); + file.createAttribute("nminy", HighFive::DataSpace::From(nminy)).write(nminy); + file.createAttribute("nminz", HighFive::DataSpace::From(nminz)).write(nminz); + file.createAttribute("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx); + file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); + file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); + } + void Cube::writeCovarH5Params(HighFive::File& file) { file.createAttribute("nminx", HighFive::DataSpace::From(nminx)).write(nminx); diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index da76c6f28..b1ca806c7 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -3024,6 +3024,8 @@ namespace CoefClasses ret = std::make_shared(ret.get()); } else if (dynamic_cast(coef.get())) { ret = std::make_shared(ret.get()); + } else if (dynamic_cast(coef.get())) { + ret = std::make_shared(ret.get()); } else if (dynamic_cast(coef.get())) { ret = std::make_shared(ret.get()); } else if (dynamic_cast(coef.get())) { diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index e13c8cd46..d04d63316 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -21,6 +21,7 @@ #include "SLGridMP2.H" #include "massmodel.H" #include "EXPmath.H" +#include "libvars.H" // For parsed version info #ifdef USE_DMALLOC #include @@ -87,6 +88,14 @@ static int sl_dim; //====================================================================== +// Default model file names +const std::string SLGridSph::default_model = "SLGridSph.model"; + +// Default cache file names +const std::string SLGridSph::default_cache = ".slgrid_sph_cache"; +const std::string SLGridSlab::default_cache = ".slgrid_slab_cache"; + + int SLGridSph::mpi = 0; // initially off extern "C" { @@ -114,7 +123,7 @@ double sphdens(double r) { // This 4pi from Poisson's eqn // | - // | /-- This begins the true density profile + // | +--- This begins the true density profile // | | // v v return 4.0*M_PI * model->get_density(r); @@ -135,8 +144,8 @@ void SLGridSph::bomb(string oops) throw std::runtime_error(sout.str()); } - // Constructors - +// Constructors +// SLGridSph::SLGridSph(std::string modelname, int LMAX, int NMAX, int NUMR, double RMIN, double RMAX, @@ -1844,49 +1853,84 @@ static double poffset=0.0; class IsothermalSlab : public SlabModel { +private: + + std::string psa = + "---- IsothermalSlab NOW uses the traditional density profile proportional to sech^2(z/2H).\n" + "---- If you are using the old profile proportional to sech^2(z/H), please update your model\n" + "---- to use the new profile and set the scale height H to be half of the old value. This\n" + "---- will ensure that your model has the same density profile and potential as before, but\n" + "---- with a more standard functional form. If you have any questions or concerns about\n" + "---- this change, please contact the developers on GitHub."; + public: - IsothermalSlab() { id = "iso"; } + IsothermalSlab() { + id = "iso"; + if (myid==0 and (exp_build.major < 7 or + (exp_build.major == 7 and exp_build.minor < 11))) + std::cout << "---- SLGridSlab: IMPORTANT UPDATE for EXP " + << VERSION << '\n' << psa << std::endl; + } double pot(double z) { - return 2.0*M_PI*SLGridSlab::H*log(cosh(z/SLGridSlab::H)) - poffset; + return 4.0*M_PI*SLGridSlab::H*log(cosh(0.5*z/SLGridSlab::H)) - poffset; } double dpot(double z) { - return 2.0*M_PI*tanh(z/SLGridSlab::H); + return 2.0*M_PI*tanh(0.5*z/SLGridSlab::H); } double dens(double z) { - double tmp = 1.0/cosh(z/SLGridSlab::H); - return 4.0*M_PI * 0.5/SLGridSlab::H * tmp*tmp; + double tmp = 1.0/cosh(0.5*z/SLGridSlab::H); + return 4.0*M_PI * 0.25/SLGridSlab::H * tmp*tmp; + // ^ + // | + // +--- The 4*pi factor simplifies the SL solution } }; -//! Constant density slab with G = M = 1 -class ConstantSlab : public SlabModel +//! Uniform density slab with G = M = 1 +class UniformSlab : public SlabModel { public: - ConstantSlab() { id = "const"; } + UniformSlab() { id = "uniform"; } double pot(double z) { - return z*z/(4.0*SLGridSlab::H) - poffset; + // Offset defined so that potential is zero at |z|=H, and negative + // for |z| SLGridSlab::H) + return 2.0*M_PI*(zz - SLGridSlab::H); + else + return M_PI*(zz*zz/SLGridSlab::H - SLGridSlab::H); } double dpot(double z) { - return z/(2.0*SLGridSlab::H); + double zz = fabs(z); + if (zz > SLGridSlab::H) + return 2.0*M_PI*z/zz; + else + return 2.0*M_PI*z/SLGridSlab::H; } double dens(double z) { - return 4.0*M_PI / (2.0 * SLGridSlab::H); + if (fabs(z) > SLGridSlab::H) + return 0.0; + else + return 4.0*M_PI/(2.0 * SLGridSlab::H); + // ^ + // | + // +--- The 4*pi factor simplifies the SL solution } }; @@ -1919,6 +1963,49 @@ class ParabolicSlab : public SlabModel double h = SLGridSlab::H; double h2 = h*h; return 4.0*M_PI * 3.0*(1.0 - z*z/h2)/(4.0*h); + // ^ + // | + // +--- The 4*pi factor simplifies the SL solution + } +}; + +//! Cosine bell model with G = M = 1 +class CosineSlab : public SlabModel +{ +public: + + CosineSlab() { id = "cosine"; } + + double pot(double z) + { + double h = SLGridSlab::H; + if (fabs(z) > h) + return 2.0*M_PI*(fabs(z) - h); + else + return 2.0*M_PI/h*(0.5*z*z + h*h/(M_PI*M_PI)*(1.0 - cos(M_PI*z/h))) - M_PI*h*(4.0/(M_PI*M_PI) + 1.0); + } + + double dpot(double z) + { + double h = SLGridSlab::H; + if (fabs(z) > h) + return 2.0*M_PI*z/fabs(z); + else + return 2.0*M_PI/h*(z + h/M_PI*sin(M_PI*z/h)); + } + + double dens(double z) + { + double h = SLGridSlab::H; + if (fabs(z) > h) + return 0.0; + else { + double cosfac = cos(0.5*M_PI*z/h); + return 4.0*M_PI/h * cosfac*cosfac; + // ^ + // | + // +--- The 4*pi factor simplifies the SL solution + } } }; @@ -1937,8 +2024,12 @@ std::shared_ptr SlabModel::createModel(const std::string type) return std::make_shared(); } - if (data.find("const") != std::string::npos) { - return std::make_shared(); + if (data.find("unif") != std::string::npos) { + return std::make_shared(); + } + + if (data.find("cos") != std::string::npos) { + return std::make_shared(); } // Default @@ -1956,26 +2047,43 @@ void SLGridSlab::bomb(string oops) // Constructors SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, - const std::string TYPE, bool VERBOSE) + const std::string cachename, + const std::string CMAP, + const std::string TYPE, + bool VERBOSE) { + if (cachename.size()) slab_cache_name = cachename; + else throw std::runtime_error("SLGridSlab: you must specify a cachename"); + int kx, ky; numk = NUMK; nmax = NMAX; numz = NUMZ; + cmap = CMAP; type = TYPE; - zmax = ZMAX; slab = SlabModel::createModel(type); + if (myid==0) { + std::cout << "---- SLGridSlab::SLGridSlab: using slab model <" + << slab->ID() << ">" << std::endl; + } + poffset = slab->pot((1.0+ZEND)*zmax); tbdbg = VERBOSE; - // This could be controlled by a parameter...but at this point, this - // is a fixed tuning. - mM = CoordMap::factory(CoordMapTypes::Sech, H); + std::transform(cmap.begin(), cmap.end(), cmap.begin(), + [](unsigned char c){ return std::tolower(c); }); + + cmtype = + std::find_if(CoordMapNames.begin(), CoordMapNames.end(), + [&](const std::pair& pair) + { return pair.second == cmap; })->first; + + mM = CoordMap::factory(cmtype, H); init_table(); @@ -2193,9 +2301,6 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, } -const string slab_cache_name = ".slgrid_slab_cache"; - - bool SLGridSlab::ReadH5Cache(void) { if (!cache) return false; @@ -2255,8 +2360,21 @@ bool SLGridSlab::ReadH5Cache(void) if (not checkStr(geometry, "geometry")) return false; if (not checkStr(forceID, "forceID")) return false; + // Version check + // + if (h5file.hasAttribute("Version")) { + if (not checkStr(Version, "Version")) return false; + } else { + if (myid==0) + std::cout << "---- SLGridSlab::ReadH5Cache: " + << "recomputing cache for HighFive API change" + << std::endl; + return false; + } + // Parameter check // + if (not checkStr(cmap, "cmap")) return false; if (not checkStr(type, "type")) return false; if (not checkInt(numk, "numk")) return false; if (not checkInt(nmax, "nmax")) return false; @@ -2273,7 +2391,6 @@ bool SLGridSlab::ReadH5Cache(void) // Create table instances // - table = table_ptr_2D(new table_ptr_1D [numk+1]); for (int kx=0; kx<=numk; kx++) table[kx] = table_ptr_1D(new TableSlab [kx+1]); @@ -2348,9 +2465,11 @@ void SLGridSlab::WriteH5Cache(void) file.createAttribute("geometry", HighFive::DataSpace::From(geometry)).write(geometry); file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); + file.createAttribute("Version", HighFive::DataSpace::From(Version)).write(Version); // Write parameters // + file.createAttribute ("cmap", HighFive::DataSpace::From(cmap)).write(cmap); file.createAttribute ("type", HighFive::DataSpace::From(type)).write(type); file.createAttribute ("numk", HighFive::DataSpace::From(numk)).write(numk); file.createAttribute ("nmax", HighFive::DataSpace::From(nmax)).write(nmax); @@ -2413,18 +2532,19 @@ double SLGridSlab::LinearMap::d_xi_to_z(double xi) { return 1.0; } double SLGridSlab::get_pot(double x, int kx, int ky, int n, int which) { - int hold; - - // Flip sign for antisymmetric basis functions + // Flip sign for antisymmetric basis functions + // int sign=1; if (x<0 && 2*(n/2)!=n) sign=-1; x = fabs(x); - if (which) // Convert from z to x + // Convert from z to x + // + if (which) x = mM->z_to_xi(x); if (ky > kx) { - hold = ky; + int hold = ky; ky = kx; kx = hold; } @@ -2450,17 +2570,17 @@ double SLGridSlab::get_pot(double x, int kx, int ky, int n, int which) double SLGridSlab::get_dens(double x, int kx, int ky, int n, int which) { - int hold; - int sign=1; if (x<0 && 2*(n/2)!=n) sign=-1; x = fabs(x); - if (which) // Convert from z to x + // Convert from z to x + // + if (which) x = mM->z_to_xi(x); if (ky > kx) { - hold = ky; + int hold = ky; ky = kx; kx = hold; } @@ -2485,17 +2605,17 @@ double SLGridSlab::get_dens(double x, int kx, int ky, int n, int which) double SLGridSlab::get_force(double x, int kx, int ky, int n, int which) { - int hold; - int sign=1; if (x<0 && 2*(n/2)==n) sign = -1; x = fabs(x); - if (which) // Convert from z to x + // Convert from z to x + // + if (which) x = mM->z_to_xi(x); if (ky > kx) { - hold = ky; + int hold = ky; ky = kx; kx = hold; } @@ -2507,17 +2627,17 @@ double SLGridSlab::get_force(double x, int kx, int ky, int n, int which) double p = (x - xi[indx])/dxi; - // Use three point formula + // Use three point formula - // Point -1: indx-1 - // Point 0: indx - // Point 1: indx+1 + // Point -1: indx-1 + // Point 0: indx + // Point 1: indx+1 return mM->d_xi_to_z(x)/dxi * ( - (p - 0.5)*table[kx][ky].ef(n, indx-1)*p0[indx-1] - -2.0*p*table[kx][ky].ef(n, indx)*p0[indx] - + (p + 0.5)*table[kx][ky].ef(n, indx+1)*p0[indx+1] - ) / sqrt(table[kx][ky].ev[n]) * sign; + (p - 0.5)*table[kx][ky].ef(n, indx-1)*p0[indx-1] + -2.0*p*table[kx][ky].ef(n, indx)*p0[indx] + + (p + 0.5)*table[kx][ky].ef(n, indx+1)*p0[indx+1] + ) / sqrt(table[kx][ky].ev[n]) * sign; } @@ -2638,7 +2758,6 @@ void SLGridSlab::get_force(Eigen::MatrixXd& mat, double x, int which) l++; } } - } @@ -2661,13 +2780,20 @@ void SLGridSlab::get_pot(Eigen::VectorXd& vec, double x, int kx, int ky, int whi vec.resize(nmax); - int indx = (int)( (x-xmin)/dxi ); + int indx = 0; + if (x < xmin) { + indx = 0; + } else if (x > xmax) { + indx = numz - 2; + } else { + indx = (int)( (x-xmin)/dxi ); + } + if (indx<0) indx = 0; if (indx>numz-2) indx = numz - 2; double x1 = (xi[indx+1] - x)/dxi; double x2 = (x - xi[indx])/dxi; - sign2 = 1; for (int n=0; n xmax) { + indx = numz - 2; + } else { + indx = (int)( (x-xmin)/dxi ); + } + if (indx<0) indx = 0; if (indx>numz-2) indx = numz - 2; @@ -2737,11 +2871,18 @@ void SLGridSlab::get_force(Eigen::VectorXd& vec, double x, int kx, int ky, int w vec.resize(nmax); - int indx = (int)( (x-xmin)/dxi ); + int indx = 0; + if (x < xmin) { + indx = 0; + } else if (x > xmax) { + indx = numz - 2; + } else { + indx = (int)( (x-xmin)/dxi ); + } + if (indx<1) indx = 1; if (indx>numz-2) indx = numz - 2; - double p = (x - xi[indx])/dxi; double fac = mM->d_xi_to_z(x)/dxi; @@ -2800,7 +2941,6 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) cons[4] = (df + KKZ*f)*f; } cons[5] = 1.0; - // cons[5] = 1.0/(f*f); // // Initialize the vector INVEC(*): @@ -2908,7 +3048,7 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) << std::endl; - if (VERBOSE && iflag[i] != 0) { + if (tbdbg or (VERBOSE && iflag[i] != 0)) { if (iflag[i] > -10) { @@ -3034,7 +3174,7 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) << std::setw( 5) << iflag[i] << std::endl; - if (VERBOSE) { + if (tbdbg or (VERBOSE && iflag[i] != 0)) { if (iflag[i] > -10) { cout << std::setw(14) << "x" diff --git a/exputil/libvars.cc b/exputil/libvars.cc index 14ff71690..11ee23a56 100644 --- a/exputil/libvars.cc +++ b/exputil/libvars.cc @@ -3,6 +3,7 @@ standalone utilities */ +#include "config_exp.h" #include "libvars.H" namespace __EXP__ @@ -36,6 +37,7 @@ namespace __EXP__ //! Sanity tolerance for orthogonality double orthoTol = 1.0e-2; + }; diff --git a/include/EXPversion.H b/include/EXPversion.H new file mode 100644 index 000000000..2c17d3b5e --- /dev/null +++ b/include/EXPversion.H @@ -0,0 +1,48 @@ +#pragma once + +namespace __EXP__ +{ + // Compile-time parser for "X.Y.Z" format + + /* Example usage: + + #define VERSION_STR "2.4.12" + static constexpr EXPversion current_v = EXP_parse_version(VERSION_STR); + + int main() { + static_assert(current_v.major == 2, "Major version mismatch"); + std::cout << "Parsed: " << current_v.major << "." << current_v.minor << "\n"; + return 0;} + */ + + + struct EXPversion + { + int major, minor, patch; + }; + + // Compile-time parser for "X.Y.Z" format using c++-17 features + constexpr EXPversion EXP_parse_version(const char* str) + { + EXPversion v = {0, 0, 0}; + int* target = &v.major; + int current = 0; + + // Scan the version string until we hit a dot, then move to the next target + for (int i = 0; str[i] != '\0'; ++i) { + if (str[i] == '.') { + *target = current; + if (target == &v.major) target = &v.minor; + else if (target == &v.minor) target = &v.patch; + current = 0; + } else { + // Append the current decimal digit to the current version + current = current * 10 + (str[i] - '0'); + } + } + // Store the int of this version string + *target = current; + return v; + } +} + diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index dd46e28cf..7b135ea67 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -74,10 +74,10 @@ private: bool tbdbg; //! Default model file name - const std::string default_model = "SLGridSph.model"; + static const std::string default_model; //! Default cache file name - const std::string default_cache = ".slgrid_sph_cache"; + static const std::string default_cache; //! Model file name std::string model_file_name; @@ -292,12 +292,24 @@ private: void compute_table_worker(void); - // Local MPI stuff + //@{ + //! Local MPI stuff void mpi_setup(void); void mpi_unpack_table(void); int mpi_pack_table(TableSlab* table, int kx, int ky); + //@} + + //@{ + //! Cache reading and writing bool ReadH5Cache(void); void WriteH5Cache(void); + //@} + + //! Cache file name + std::string slab_cache_name; + + //! Cache versioning + inline static const std::string Version = "1.0"; int mpi_myid, mpi_numprocs; int mpi_bufsz; @@ -315,6 +327,13 @@ private: //! to finite interval enum CoordMapTypes {Tanh, Sech, Linear}; + //! Coordinate map chooser + const std::map CoordMapNames = { + {Tanh, "tanh" }, + {Sech, "sech" }, + {Linear, "linear"} + }; + //! The base class for all coordinate maps class CoordMap { @@ -374,12 +393,21 @@ private: virtual double d_xi_to_z(double z); }; + //! Coordinate map type (enum) + CoordMapTypes cmtype; + + //! Coordinate map type (string) + std::string cmap = "linear"; + //! Model type std::string type; //! The current map created by the SLGridSlab constructor std::unique_ptr mM; + //! Default cache file name + static const std::string default_cache; + public: //! Global MPI flag, default: 0=off @@ -403,7 +431,17 @@ public: //! Constructor SLGridSlab(int kmax, int nmax, int numz, double zmax, - const std::string type="isothermal", bool Verbose=false); + const std::string cachename, + const std::string coordmap="linear", + const std::string type="isothermal", + bool Verbose=false); + + //! Backward-compatible constructor for legacy call sites + SLGridSlab(int kmax, int nmax, int numz, double zmax, + const std::string type, + bool Verbose=false) + : SLGridSlab(kmax, nmax, numz, zmax, + default_cache, "linear", type, Verbose) {} //! Destructor ~SLGridSlab(); diff --git a/include/libvars.H b/include/libvars.H index 12dfd45d3..bbbcd6e11 100644 --- a/include/libvars.H +++ b/include/libvars.H @@ -7,6 +7,9 @@ #include +#include "config_exp.h" // For version string +#include "EXPversion.H" // For compile-time version parsing + namespace __EXP__ { //! @name Theading variables @@ -42,6 +45,9 @@ namespace __EXP__ //! Sanity tolerance for orthogonality extern double orthoTol; + //! Compile-time version parsing + static constexpr EXPversion exp_build = EXP_parse_version(VERSION); + }; #endif // END _LIBVARS_H diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index cd4f650cb..7c675e872 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -885,6 +885,11 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Slab, make_coefs,); } + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covr_tot, bool covar) override + { + PYBIND11_OVERRIDE(void, Slab, enableCoefCovariance, pcavar, nsamples, ftype, covr_tot, covar); + } + }; diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 1164ebf6b..8ed67fe1f 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -1000,8 +1000,13 @@ void CoefficientClasses(py::module &m) { py::class_, CoefStruct> (m, "SlabStruct") .def(py::init<>(), "Slab coefficient data structure object") - .def("assign", &SlabStruct::assign, - R"( + .def("assign", + [](CoefClasses::SlabStruct& A, py::array_t> mat) + { + auto M = make_tensor3>(mat); + A.assign(M); + }, + R"( Assign a coefficient matrix to CoefStruct. Parameters @@ -1022,8 +1027,13 @@ void CoefficientClasses(py::module &m) { py::class_, CoefStruct> (m, "CubeStruct") .def(py::init<>(), "Cube coefficient data structure object") - .def("assign", &CubeStruct::assign, - R"( + .def("assign", + [](CoefClasses::SlabStruct& A, py::array_t> mat) + { + auto M = make_tensor3>(mat); + A.assign(M); + }, + R"( Assign a coefficient matrix to CoefStruct. Parameters @@ -1045,13 +1055,13 @@ void CoefficientClasses(py::module &m) { (m, "TblStruct") .def(py::init<>(), "Multicolumn table data structure object") .def("assign", &TblStruct::assign, - R"( + R"( Assign a coefficient matrix to CoefStruct. Parameters ---------- mat : numpy.ndarray, complex - complex-valued NumPy array of table values + complex-valued NumPy array of table values Returns ------- @@ -1585,9 +1595,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1660,9 +1670,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1772,9 +1782,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1856,9 +1866,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient matrix - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray + the new coefficient array. Returns ------- @@ -1919,16 +1929,21 @@ void CoefficientClasses(py::module &m) { )", py::arg("time")) .def("setTensor", - &CoefClasses::SlabCoefs::setTensor, + [](CoefClasses::SlabCoefs& A, double time, + py::array_t> mat) + { + auto M = make_tensor3>(mat); + A.setTensor(time, M); + }, R"( Enter and/or rewrite the coefficient tensor at the provided time Parameters ---------- time : float - snapshot time corresponding to the the coefficient tensor - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient matrix + mat : numpy.ndarray of complex values + the new coefficient array. Returns ------- @@ -2020,9 +2035,9 @@ void CoefficientClasses(py::module &m) { Parameters ---------- time : float - snapshot time corresponding to the the coefficient tensor - mat : numpy.ndarray - the new coefficient array. + snapshot time corresponding to the the coefficient tensor + mat : numpy.ndarray + the new coefficient array. Returns ------- diff --git a/src/Component.cc b/src/Component.cc index c50884ce6..af9109217 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1393,7 +1393,7 @@ void Component::initialize(void) cout << "---- Component <" << name << ">: "; - if (nlevel<0) + if (nlevel<=0) std::cout << "no multistep level reporting"; else std::cout << "multistep level reporting every " << nlevel << " steps"; diff --git a/src/Cube.H b/src/Cube.H index 80b3ad5a8..03d8ee62f 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -10,7 +10,7 @@ #include #include "Coefficients.H" -#include "PotAccel.H" +#include "Basis.H" #if HAVE_LIBCUDA==1 #include @@ -66,7 +66,7 @@ subsampling is statistical assessment, this should not be a limitation. */ -class Cube : public PotAccel +class Cube : public Basis { private: @@ -323,6 +323,33 @@ public: //! Indicate whether or not subsampling is supported virtual bool hasSubsample() { return true; } + + //! Return density, potential, and potential gradient at a point in + //! Cartesian coordinates + void determine_fields_at_point(double x, double y, double z, + double *tdens0, double *tpotl0, + double *tdens, double *tpotl, + double *tpotx, double *tpoty, + double *tpotz); + + //! Returns the potential, its derivatives & density at a point in + //! spherical coordinates + void determine_fields_at_point_sph + (double r, double theta, double phi, + double *tdens0, double *dpotl0, double *tdens, double *tpotl, + double *tpotr, double *tpott, double *tpotp) { + throw std::runtime_error("Cube::determine_fields_at_point_sph: Not implemented for this basis"); + } + + //! Returns the potential, its derivatives & density at a point in + //! cylindrical coordinates + void determine_fields_at_point_cyl + (double r, double z, double phi, + double *tdens0, double *dpotl0, double *tdens, double *tpotl, + double *tpotr, double *tpotz, double *tpotp) { + throw std::runtime_error("Cube::determine_fields_at_point_cyl: Not implemented for this basis"); + } + }; diff --git a/src/Cube.cc b/src/Cube.cc index f4e184dca..07f173bbb 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -29,7 +29,7 @@ static bool deepDebug = false; static bool coefDebug = false; //@} -Cube::Cube(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) +Cube::Cube(Component* c0, const YAML::Node& conf) : Basis(c0, conf) { // ID parameters // @@ -1038,29 +1038,85 @@ PotAccel::CovarData Cube::getSubsample() std::get<2>(elem) = Eigen::Tensor, 3>(sampT, 1, osize); std::get<3>(elem) = Eigen::Tensor, 4>(sampT, 1, osize, osize); + if (not fullCovar) std::get<3>(elem).setZero(); + // Fill the covariance structure with subsamples + // for (int T=0; T(elem)(T, 0, n1) = meanV[T](n1); - for (int n2=0; n2(elem)(T, 0, n1, n2) = covrV[T](n1, n2); - else - std::get<3>(elem)(T, 0, n1, n2) = 0.0; - } - */ - std::get<2>(elem).chip(T, 0).chip(0, 0) = Eigen::TensorMap, 1>> (meanV[T].data(), osize); - std::get<3>(elem).chip(T, 0).chip(0, 0) = - Eigen::TensorMap, 2>> - (covrV[T].data(), osize, osize); - + // Only fill the covariance if requested + // + if (fullCovar) { + std::get<3>(elem).chip(T, 0).chip(0, 0) = + Eigen::TensorMap, 2>> + (covrV[T].data(), osize, osize); + } } return elem; } + +// Return density, potential, and potential gradient at a point +void Cube::determine_fields_at_point +(double x, double y, double z, + double *tdens0, double *tpotl0, double *tdens, double *tpotl, + double *tpotx, double *tpoty, double *tpotz) +{ + // Zero the fields + // + *tdens0 = *tpotl0 = 0.0; // k==0 contribution + *tdens = *tpotl = 0.0; // k!=0 contribution + *tpotx = *tpoty = *tpotz = 0.0; + + // Recursion multipliers + // + auto stepx = std::exp(kfac*x); + auto stepy = std::exp(kfac*y); + auto stepz = std::exp(kfac*z); + + // Initial values (note sign change from coefficient accumulation) + // + auto startx = std::exp(-kfac*(x*nmaxx)); + auto starty = std::exp(-kfac*(y*nmaxy)); + auto startz = std::exp(-kfac*(z*nmaxz)); + + std::complex facx, facy, facz; + int ix, iy, iz; + + for (facx=startx, ix=0; ix fac = facx*facy*facz*expcoef(ix, iy, iz); + + // Compute wavenumber; recall that the coefficients are + // stored as follows: -nmax,-nmax+1,...,0,...,nmax-1,nmax + // + int ii = ix-nmaxx; + int jj = iy-nmaxy; + int kk = iz-nmaxz; + + // No contribution to acceleration and potential ("swindle") + // for zero wavenumber + if (ii==0 && jj==0 && kk==0) continue; + + // Limit to minimum wave number + if (abs(ii)(0.0, dfac*ii)*fac*norm).real(); + *tpoty += (std::complex(0.0, dfac*jj)*fac*norm).real(); + *tpotz += (std::complex(0.0, dfac*kk)*fac*norm).real(); + } + } + } +} diff --git a/src/SlabSL.H b/src/SlabSL.H index 693ded977..792f1b985 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -9,7 +9,7 @@ #include "Coefficients.H" #include "SLGridMP2.H" #include "biorth1d.H" -#include "PotAccel.H" +#include "Basis.H" #if HAVE_LIBCUDA==1 #include @@ -17,10 +17,58 @@ #include #endif -/*! This routine computes the potential, acceleration and density - using expansion periodic in X & Y and outgoing vacuum boundary - condtions in Z */ -class SlabSL : public PotAccel +/** @class SlabSL + @brief This routine computes the potential, acceleration and density + using expansion periodic in X & Y and outgoing vacuum boundary + conditions in Z + + @details **YAML configuration** + + @param nmaxx is the maximum order of the expansion in x (default 6) + + @param nmaxy is the maximum order of the expansion in y (default 6) + + @param nmaxz is the maximum order of the expansion in z (default 6) + + @param nminx is the minimum order of the expansion in x (default 0) + + @param nminy is the minimum order of the expansion in y (default 0) + + @param hslab is the scale height of the slab (default 0.2) + + @param zmax is the maximum z for the slab (default 10.0) + + @param ngrid is the number of grid points in z for the + Sturm-Liouville solver (default 1000) + + @param cmap is the coordinate mapping for the SL solver (must be one of "tanh", "sech", or "linear", default "linear") + + @param type is the type of slab to solve for (default + "isothermal", must be "isothermal", "parabolic", or "constant") + + @param no_even is a flag to exclude even modes in the expansion or + K!=0 (default false) + + @param no_odd is a flag to exclude odd modes in the expansion for + K!=0 (default false) + + @param self_consistent set to true allows the particles to evolve + under the time-dependent basis expansion. For a basis fixed in time to + the initial time: set to false. + + @param nint integer number of steps between subsample computations + (default: 0, no subsampling) + + @param samplesz integer number of subsample bins (default: 0, + no subsampling). Must be greater than 0 when `nint>0`. + + @param orthochk boolean flag to compute orthogonality matrix for the + basis (default false). For debugging only; set to false for production. + + @param cachename is the name of the basis cache file. This is a + required parameter. +*/ +class SlabSL : public Basis { //! Pull in base-class overrides @@ -47,15 +95,21 @@ private: //! Current coefficient tensor std::vector expccof, expccofP; - int nminx, nminy; + //! Coefficient tensor for frozen potential (if self_consistent=false) + coefType expcofF; + + int nminx, nminy, sampT=0, nint=0; int nmaxx, nmaxy, nmaxz; double zmax, hslab; + bool no_even = false, no_odd = false; + int imx, imy, imz, jmax, nnmax; - double dfac; + double dfac, last=-std::numeric_limits::max(); std::complex kfac; + std::string cachename; - std::vector zfrc, zpot; + std::vector zfrc, zpot, zden; SlabSLCoefHeader coefheader; @@ -96,7 +150,11 @@ private: thrust::device_vector> dw_coef; thrust::device_vector> df_coef; - void resize_coefs(int N, int osize, int gridSize, int stride); + thrust::device_vector u_d; + + std::vector>> T_samp; + + void resize_coefs(int N, int osize, int gridSize, int sampT); }; //! A storage instance @@ -109,7 +167,7 @@ private: void cuda_initialize(); //! Initialize the cuda streams - void initialize_cuda() + void slab_initialize_cuda() { grid->initialize_cuda(cuInterpArray, tex); } @@ -120,15 +178,44 @@ private: #endif + //! Flag self_consistency + bool self_consistent = true; + + //! Flag whether coefficients have been initialized for the first time + bool firstime_coef = true; + //! Default number of grid points for SLGridSlab int ngrid = 1000; + //! Default coordinate mapping for SLGridSlab (must be "tanh", "sech", or "linear") + std::string cmap = "linear"; + //! Default slab type (must be "isothermal", "parabolic", or "constant") std::string type = "isothermal"; + //@{ + //! Covariance structures. First index is T, second is the + //! flattened 3-d k vector + std::vector meanV; + std::vector covrV; + + //@{ + //! Per thread covariance structures + std::vector> meanV1; + std::vector> covrV1; + std::vector workV1; + std::vector countV1; + std::vector massV1; + //@} + + void init_covariance(); + void zero_covariance(); + //@} + //@{ //! Usual evaluation interface void determine_coefficients(void); + void determine_acceleration_and_potential(void); void get_acceleration_and_potential(Component*); //@} @@ -188,6 +275,9 @@ private: // Biorth ID static const int ID=1; + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); + protected: //! Parse parameters and initialize on first call @@ -196,6 +286,9 @@ protected: //! Valid keys for YAML configurations static const std::set valid_keys; + //! Print orthogonality matrix for debugging + bool orthochk = false; + public: //! Id string @@ -209,6 +302,41 @@ public: //! Coefficient output void dump_coefs_h5(const std::string& file); + + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarData getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } + + //! Return density, potential, and potential gradient at a point in + //! Cartesian coordinates + void determine_fields_at_point(double x, double y, double z, + double *tdens0, double *tpotl0, + double *tdens, double *tpotl, + double *tpotx, double *tpoty, + double *tpotz); + + //! Returns the potential, its derivatives & density at a point in + //! spherical coordinates + void determine_fields_at_point_sph + (double r, double theta, double phi, + double *tdens0, double *dpotl0, double *tdens, double *tpotl, + double *tpotr, double *tpott, double *tpotp) { + throw std::runtime_error("SlabSL::determine_fields_at_point_sph: Not implemented for this basis"); + } + + //! Returns the potential, its derivatives & density at a point in + //! cylindrical coordinates + void determine_fields_at_point_cyl + (double r, double z, double phi, + double *tdens0, double *dpotl0, double *tdens, double *tpotl, + double *tpotr, double *tpotz, double *tpotp) { + throw std::runtime_error("SlabSL::determine_fields_at_point_cyl: Not implemented for this basis"); + } }; diff --git a/src/SlabSL.cc b/src/SlabSL.cc index 1de35675a..f6642bbf3 100644 --- a/src/SlabSL.cc +++ b/src/SlabSL.cc @@ -17,7 +17,16 @@ SlabSL::valid_keys = { "hslab", "zmax", "ngrid", - "type" + "cmap", + "type", + "nint", + "no_even", + "no_odd", + "samplesz", + "subsampleFloat", + "self_consistent", + "orthochk", + "cachename" }; //@{ @@ -26,21 +35,30 @@ static bool cudaAccumOverride = false; static bool cudaAccelOverride = false; //@} -SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) +SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : Basis(c0, conf) { id = "Slab (Sturm-Liouville)"; + + // Default values + // nminx = nminy = 0; nmaxx = nmaxy = nmaxz = 6; zmax = 10.0; hslab = 0.2; coef_dump = true; + cachename = ""; #if HAVE_LIBCUDA==1 cuda_aware = true; #endif + // Parse the YAML and initialize CUDA if needed + // initialize(); + if (cachename.size()==0) + throw std::runtime_error("SlabSL: you must specify a cachename"); + SLGridSlab::mpi = 1; SLGridSlab::ZBEG = 0.0; SLGridSlab::ZEND = 0.1; @@ -48,7 +66,10 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) int nnmax = (nmaxx > nmaxy) ? nmaxx : nmaxy; - grid = std::make_shared(nnmax, nmaxz, ngrid, zmax, type); + // Make the Sturm-Liouville grid and basis functions + // + grid = std::make_shared(nnmax, nmaxz, ngrid, zmax, cachename, + cmap, type); // Test for basis consistency (will generate an exception if maximum // error is out of tolerance) @@ -57,6 +78,8 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) int kxw = 0, kyw = 0; auto test = grid->orthoCheck(10000); for (int kx=0, indx=0; kx<=nnmax; kx++) { + // We only need the lower triangle of the kx-ky plane since the + // basis is symmetric in kx and ky for (int ky=0; ky<=kx; ky++, indx++) { for (int n1=0; n1 __EXP__::orthoTol) { if (myid==0) std::cout << "SlabSL: orthogonality failure, worst=" << worst @@ -93,22 +121,24 @@ SlabSL::SlabSL(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) << worst << std::endl; } - imx = 1 + 2*nmaxx; // Number of x wavenumber - imy = 1 + 2*nmaxy; // Number of y wavenumbers - imz = nmaxz; // Number of vertical functions - jmax = imx * imy * imz; // Total storage in tensor - + // Allocate temporary storage for coefficients + // expccof.resize(nthrds); for (auto & v : expccof) v.resize(imx, imy, imz); + // Lx and Ly are both 1.0, so the wavenumbers are 2*pi*n with n an + // integer dfac = 2.0*M_PI; kfac = std::complex(0.0, dfac); + // Arrays for multithreading force evaluation zpot.resize(nthrds); zfrc.resize(nthrds); + zden.resize(nthrds); for (auto & v : zpot) v.resize(nmaxz); for (auto & v : zfrc) v.resize(nmaxz); + for (auto & v : zden) v.resize(nmaxz); // Allocate coefficient tensor (one for each multistep level) and // zero-out contents @@ -163,7 +193,26 @@ void SlabSL::initialize() if (conf["ngrid"]) ngrid = conf["ngrid"].as(); if (conf["hslab"]) hslab = conf["hslab"].as(); if (conf["zmax" ]) zmax = conf["zmax" ].as(); + if (conf["cmap" ]) cmap = conf["cmap" ].as(); if (conf["type" ]) type = conf["type" ].as(); + if (conf["no_even"]) no_even = conf["no_even"].as(); + if (conf["no_odd"]) no_odd = conf["no_odd"].as(); + if (conf["orthochk"]) orthochk = conf["orthochk"].as(); + if (conf["cachename"]) cachename = conf["cachename"].as(); + + if (conf["self_consistent"]) { + self_consistent = conf["self_consistent"].as(); + } else + self_consistent = true; + + if (conf["nint"]) { + nint = conf["nint"].as(); + if (nint>0) computeSubsample = true; + } + + if (conf["samplesz"]) { + sampT = conf["samplesz"].as(); + } } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in SlabSL: " @@ -176,19 +225,109 @@ void SlabSL::initialize() throw std::runtime_error("SlabSL::initialze: error parsing YAML"); } + // Set dimensions of coefficient tensor and total storage in tensor + // + imx = 1 + 2*nmaxx; // Number of x wavenumber + imy = 1 + 2*nmaxy; // Number of y wavenumbers + imz = nmaxz; // Number of vertical functions + jmax = imx * imy * imz; // Total storage in tensor + + // Initialize covariance + // + init_covariance(); + + // Convert cmap to lower case (needed by cuda) + // + std::transform(cmap.begin(), cmap.end(), cmap.begin(), + [](unsigned char c){ return std::tolower(c); }); + #if HAVE_LIBCUDA==1 + // Cuda initialization (if needed) + // cuda_initialize(); #endif } +void SlabSL::init_covariance() +{ + if (computeSubsample) { + + meanV.resize(sampT); + for (auto& v : meanV) { + v.resize(jmax); + } + + workV1.resize(nthrds); + for (auto& v : workV1) v.resize(jmax); + + if (fullCovar) { + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(jmax, jmax); + } + } else { + covrV.clear(); + } + + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + meanV1.resize(nthrds); + covrV1.resize(nthrds); + countV1.resize(nthrds); + massV1.resize(nthrds); + + for (int n=0; nsetZero(); } + // Determine whether or not to compute a subsample + // + if (mstep==0 or mstep==std::numeric_limits::max()) { + if (nint>0 && this_step % nint == 0) { + if (tnow > last) { + requestSubsample = true; + last = tnow; + zero_covariance(); + } + } + } else { + subsampleComputed = false; + } + #if HAVE_LIBCUDA==1 (*barrier)("SlabSL::entering cuda coefficients", __FILE__, __LINE__); if (component->cudaDevice>=0 and use_cuda) { @@ -226,12 +379,11 @@ void SlabSL::determine_coefficients(void) exp_thread_fork(true); #endif - int used1 = 0, rank = expccof[0].size(); + int used1 = 0; // Reduce over threads used = 0; for (int i=1; i stepx, stepy; unsigned nbodies = component->levlist[mlevel].size(); - int id = *((int*)arg); - int nbeg = nbodies*id/nthrds; - int nend = nbodies*(id+1)/nthrds; + int id = *((int*)arg); + int nbeg = nbodies*id/nthrds; + int nend = nbodies*(id+1)/nthrds; double adb = cC->Adiabatic(); for (int q=nbeg; q(nmaxx)*kfac*cC->Pos(i, 0)); starty = exp(static_cast(nmaxy)*kfac*cC->Pos(i, 1)); - double zz = cC->Pos(i, 2), mm = -4.0*M_PI * cC->Mass(i) * adb; + double zz = cC->Pos(i, 2), ms = cC->Mass(i); + + // Biorthogonaity gives a -4𝛑 here + constexpr double nrm = - 4.0*M_PI; + double mm = nrm * ms * adb; for (facx=startx, ix=0; ixNumber(); // And compute number of bodies + // Call the parallel acceleration and potential + // MPL_start_timer(); + determine_acceleration_and_potential(); + MPL_stop_timer(); + // Clear external potential flag + use_external = false; +} + +void SlabSL::determine_acceleration_and_potential(void) +{ if (play_back) { swap_coefs(expccofP, expccof); } if (use_external == false) { - if (multistep && initializing) { + if (multistep && (self_consistent || initializing)) { compute_multistep_coefficients(); } @@ -454,11 +696,23 @@ void * SlabSL::determine_acceleration_and_potential_thread(void * arg) for (int iz=0; iz("nminx", HighFive::DataSpace::From(nminx)).write(nminx); + file.createAttribute("nminy", HighFive::DataSpace::From(nminy)).write(nminy); + file.createAttribute("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx); + file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); + file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); +} + +PotAccel::CovarData SlabSL::getSubsample() +{ + CovarData elem; + + std::get<0>(elem) = sampleCounts; + std::get<1>(elem) = sampleMasses; + std::get<2>(elem) = Eigen::Tensor, 3>(sampT, 1, jmax); + std::get<3>(elem) = Eigen::Tensor, 4>(sampT, 1, jmax, jmax); + + if (not fullCovar) std::get<3>(elem).setZero(); + + // Fill the covariance structure with subsamples + // + for (int T=0; T(elem).chip(T, 0).chip(0, 0) = + Eigen::TensorMap, 1>> + (meanV[T].data(), jmax); + + // Use the map to copy the covariance if requested + // + if (fullCovar) { + std::get<3>(elem).chip(T, 0).chip(0, 0) = + Eigen::TensorMap, 2>> + (covrV[T].data(), jmax, jmax); + } + } + + return elem; +} + +// Return density, potential, and potential gradient at a point +void SlabSL::determine_fields_at_point +(double x, double y, double z, + double *tdens0, double *tpotl0, double *tdens, double *tpotl, + double *tpotx, double *tpoty, double *tpotz) +{ + // Zero the fields + // + *tdens0 = *tpotl0 = 0.0; // k==0 contribution + *tdens = *tpotl = 0.0; // k!=0 contribution + *tpotx = *tpoty = *tpotz = 0.0; + + // Recursion multipliers + // + std::complex stepx = exp(kfac*x); + std::complex stepy = exp(kfac*y); + + // Initial values (note sign change) + // + std::complex startx = exp(-static_cast(nmaxx)*kfac*x); + std::complex starty = exp(-static_cast(nmaxy)*kfac*y); + + std::complex facx, facy, facz, facf, facd, fac; + int ix, iy; + + // Compute wavenumber; recall that the coefficients are stored + // as follows: -nmax,-nmax+1,...,0,...,nmax-1,nmax + // + for (facx=startx, ix=0; ix nmaxx) { + std::cerr << "Out of bounds: ii=" << ii << std::endl; + } + if (iiy > nmaxy) { + std::cerr << "Out of bounds: jj=" << jj << std::endl; + } + + if (iix>=iiy) { + grid->get_pot (zpot[0], z, iix, iiy); + grid->get_dens (zden[0], z, iix, iiy); + grid->get_force(zfrc[0], z, iix, iiy); + } + else { + grid->get_pot (zpot[0], z, iiy, iix); + grid->get_dens (zden[0], z, iiy, iix); + grid->get_force(zfrc[0], z, iiy, iix); + } + + for (int iz=0; iz(ii)*fac).real(); + *tpoty += (kfac*static_cast(jj)*fac).real(); + *tpotz += facf.real(); + } + } + } +} diff --git a/src/cudaCube.cu b/src/cudaCube.cu index 9f3098066..0d13a190e 100644 --- a/src/cudaCube.cu +++ b/src/cudaCube.cu @@ -1,5 +1,3 @@ -// -*- C++ -*- - #include #include @@ -1462,4 +1460,5 @@ void Cube::destroy_cuda() // Nothing } +// -*- C++ -*- diff --git a/src/cudaSlabSL.cu b/src/cudaSlabSL.cu index cd00f09d4..6ad2c80d8 100644 --- a/src/cudaSlabSL.cu +++ b/src/cudaSlabSL.cu @@ -1,5 +1,4 @@ -// -*- C++ -*- - +#include #include #include @@ -19,6 +18,9 @@ __device__ __constant__ int slabNumX, slabNumY, slabNumZ, slabNX, slabNY, slabNZ, slabNdim, slabNum; +__device__ __constant__ +bool slabNoEven, slabNoOdd; + __device__ __constant__ int slabCmap; @@ -81,6 +83,10 @@ void testConstantsSlab() printf(" Xmin = %e\n", slabXmin ); printf(" Xmax = %e\n", slabXmax ); printf(" Dxi = %e\n", slabDxi ); + if (slabNoOdd) printf(" NoOdd = true\n" ); + else printf(" NoOdd = false\n"); + if (slabNoEven) printf(" NoEven = true\n" ); + else printf(" NoEven = false\n"); printf("-------------------------\n"); } @@ -171,7 +177,9 @@ cuFP_t cu_d_xi_to_z(cuFP_t xi) // void SlabSL::cuda_initialize() { - // Nothing + // Turn off full covariance + // + fullCovar = false; } // Copy constants to device @@ -213,8 +221,30 @@ void SlabSL::initialize_constants() size_t(0), cudaMemcpyHostToDevice), __FILE__, __LINE__, "Error copying slabNum"); - // Sech map - int Cmap = 1; + cuda_safe_call(cudaMemcpyToSymbol(slabNoEven, &no_even, sizeof(bool), + size_t(0), cudaMemcpyHostToDevice), + __FILE__, __LINE__, "Error copying slabNoEven"); + + cuda_safe_call(cudaMemcpyToSymbol(slabNoOdd, &no_odd, sizeof(bool), + size_t(0), cudaMemcpyHostToDevice), + __FILE__, __LINE__, "Error copying slabNoOdd"); + + // Coordinate map + std::map CoordMap = + { + {"tanh" , 0}, + {"sech" , 1}, + {"linear", 2} + }; + + auto it = CoordMap.find(cmap); + if (it == CoordMap.end()) { + std::ostringstream sout; + sout << "cudaSlabSL: unknown coordinate map type '" << cmap + << "'. Valid values are: tanh, sech, linear"; + throw std::runtime_error(sout.str()); + } + int Cmap = it->second; cuda_safe_call(cudaMemcpyToSymbol(slabCmap, &Cmap, sizeof(int), size_t(0), cudaMemcpyHostToDevice), @@ -242,9 +272,47 @@ void SlabSL::initialize_constants() } +__global__ void wrapKernelSlab +(dArray P, dArray I, int stride, PII lohi) +{ + // Thread ID + // + const int tid = blockDim.x * blockIdx.x + threadIdx.x; + const int N = lohi.second - lohi.first; + + for (int n=0; n=P._s) printf("out of bounds: %s:%d\n", __FILE__, __LINE__); +#endif + cudaParticle & p = P._v[I._v[npart]]; + + // Truncate to box with sides in [0,1] + // + for (int k=0; k<2; k++) { + if (p.pos[k]<0.0) + p.pos[k] += floor(-p.pos[k]) + 1.0; + else + p.pos[k] += -floor(p.pos[k]); + } + } + // END: particle index limit + } + // END: stride loop +} + + __global__ void coefKernelSlab (dArray P, dArray I, dArray coef, - dArray tex, int stride, PII lohi) + dArray tex, dArray used, int stride, PII lohi) { // Thread ID // @@ -267,7 +335,10 @@ __global__ void coefKernelSlab cudaParticle & p = P._v[I._v[npart]]; cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; - cuFP_t mm = - p.mass * 2.0*slabDfac; + // Biorthogonality gives a minus here + cuFP_t mm = - 2.0 * slabDfac * p.mass; + + used._v[i] = p.mass; // Mark particle as used // Restore particles to the unit slab // @@ -322,15 +393,18 @@ __global__ void coefKernelSlab Y = cy; // Assign the min Y wavenumber conjugate - int kx = abs(ii); // X index into texture array - for (int jj=-slabNumY; jj<=slabNumY; jj++, Y*=sy) { - int ky = abs(jj); // Y index into texture array - + // Compute the offset into the texture array for the (kx, + // ky) pair. The storage is the upper triangular part of + // the (kx, ky) plane. + // + int kx = max(abs(ii), abs(jj)); + int ky = min(abs(ii), abs(jj)); int offset = 1 + (kx*(kx+1)/2 + ky)*slabNumZ; - // The vertical basis iteration + // The vertical basis iteration + // for (int n=0; n P, dArray I, for (int n=0; n P, dArray I, a*int2_as_double(tex1D(tex._v[k], in0 )) + b*int2_as_double(tex1D(tex._v[k], in0+1)) #endif - ) * p0 * sign; + ) * p0 * signV; cuFP_t f = ( #if cuREAL == 4 @@ -480,7 +565,7 @@ forceKernelSlab(dArray P, dArray I, -2.0*s*int2_as_double(tex1D(tex._v[k], jn0)) * int2_as_double(tex1D(tex._v[0], jn0)) + (s + 0.5)*int2_as_double(tex1D(tex._v[k], jn0+1)) * int2_as_double(tex1D(tex._v[0], jn0+1)) #endif - ) * sign * dx; + ) * signF * dx; fac = X * Y * v * coef._v[slabIndex(ii, jj, n)]; facf = X * Y * f * coef._v[slabIndex(ii, jj, n)]; @@ -516,7 +601,7 @@ public: } }; -void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride) +void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int sampT) { // Reserve space for coefficient reduction // @@ -531,6 +616,15 @@ void SlabSL::cudaStorage::resize_coefs(int N, int osize, int gridSize, int strid dN_coef.resize(osize*N); dc_coef.resize(osize*gridSize); dw_coef.resize(osize); // This will stay fixed + + u_d.resize(N); + + // Resize covariance storage, if sampT>0 + // Resize covariance storage, if sampT>0 + if (sampT) { + T_samp.resize(sampT); + for (int T=0; Tstream), cuS.df_coef.begin(), cuS.df_coef.end(), 0.0); + + if (requestSubsample) { + + cuS.T_samp.resize(sampT); + + for (int T=0; Tstream), + cuS.T_samp[T].begin(), cuS.T_samp[T].end(), 0.0); + } + } } void SlabSL::determine_coefficients_cuda() @@ -552,7 +657,7 @@ void SlabSL::determine_coefficients_cuda() // must be done every time // if (initialize_cuda_slab) { - initialize_cuda(); + slab_initialize_cuda(); initialize_cuda_slab = false; // Copy coordinate mapping @@ -621,6 +726,11 @@ void SlabSL::determine_coefficients_cuda() // cuda_zero_coefs(); + // Zero mass accumulation array + // + thrust::fill(cuS.u_d.begin(), cuS.u_d.end(), 0.0); + + // Get sorted particle range for mlevel // Get sorted particle range for mlevel // PII lohi = component->CudaGetLevelRange(mlevel, mlevel), cur; @@ -684,17 +794,30 @@ void SlabSL::determine_coefficients_cuda() // int sMemSize = BLOCK_SIZE * sizeof(CmplxT); + std::vector::iterator> bm; + if (requestSubsample) { + for (int T=0; Tstream>>> + (toKernel(cs->cuda_particles), toKernel(cs->indx1), stride, cur); + // Compute the coefficient contribution for each order // auto beg = cuS.df_coef.begin(); coefKernelSlab<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), toKernel(t_d), stride, cur); + toKernel(cuS.dN_coef), toKernel(t_d), toKernel(cuS.u_d), + stride, cur); // Begin the reduction by blocks [perhaps this should use a // stride?] @@ -726,6 +849,78 @@ void SlabSL::determine_coefficients_cuda() beg, beg, thrust::plus()); thrust::advance(beg, jmax); + + if (requestSubsample) { + + int sN = N/sampT; + int nT = sampT; + + if (sN==0) { // Fail-safe underrun + sN = 1; + nT = N; + } + + for (int T=0; T(mend - mbeg); + massV1 [0][T] += static_cast(thrust::reduce(mbeg, mend)); + + // Begin the reduction per grid block + // + /* A reminder to consider implementing strides in reduceSum */ + /* + unsigned int stride1 = s/BLOCK_SIZE/deviceProp.maxGridSize[0] + 1; + unsigned int gridSize1 = s/BLOCK_SIZE/stride1; + + if (s > gridSize1*BLOCK_SIZE*stride1) gridSize1++; + */ + + unsigned int gridSize1 = s/BLOCK_SIZE; + if (s > gridSize1*BLOCK_SIZE) gridSize1++; + + // Sum reduce into gridsize1 blocks taking advantage of + // GPU warp structure + // + + // Mean computation only + // + reduceSumS + <<stream>>> + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), jmax, N, k, k+s); + + // Finish the reduction for this order in parallel + // + thrust::counting_iterator index_begin(0); + thrust::counting_iterator index_end(gridSize1*jmax); + + thrust::reduce_by_key + ( + thrust::cuda::par.on(cs->stream), + thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), + thrust::make_transform_iterator(index_end, key_functor(gridSize1)), + cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() + ); + + + thrust::transform(thrust::cuda::par.on(cs->stream), + cuS.dw_coef.begin(), cuS.dw_coef.end(), + bm[T], bm[T], thrust::plus()); + + thrust::advance(bm[T], imz); + } + // END: subsample loop + } } // use1 += N; // Increment particle count @@ -735,6 +930,18 @@ void SlabSL::determine_coefficients_cuda() // host_coefs = cuS.df_coef; + if (requestSubsample) { + // T loop + // + for (int T=0; T retV = cuS.T_samp[T]; + + for (int n=0; n(retV[n]); + } + } + } + // Create a wavenumber tuple from a flattened index // auto indices = [&](int indx) @@ -982,7 +1189,7 @@ void SlabSL::determine_acceleration_cuda() // must be done every time // if (initialize_cuda_slab) { - initialize_cuda(); + slab_initialize_cuda(); initialize_cuda_slab = false; } @@ -1047,8 +1254,10 @@ void SlabSL::HtoD_coefs() host_coefs.resize(jmax); // Copy from Slab - for (int i=0; istream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), toKernel(t_d), stride, cur); + toKernel(cuS.dN_coef), toKernel(t_d), toKernel(cuS.u_d), + stride, cur); unsigned int gridSize1 = N/BLOCK_SIZE; if (N > gridSize1*BLOCK_SIZE) gridSize1++; @@ -1196,4 +1406,4 @@ void SlabSL::destroy_cuda() // Nothing } - +// -*- C++ -*- diff --git a/utils/ICs/bonnerebert.cc b/utils/ICs/bonnerebert.cc index a7589b6d0..8d7cceb00 100644 --- a/utils/ICs/bonnerebert.cc +++ b/utils/ICs/bonnerebert.cc @@ -28,8 +28,6 @@ */ -#define VERSION 0.1 - #include #include #include @@ -366,8 +364,8 @@ decode_switches (int argc, char **argv) "R:" /* runit */ "N:" /* number */ "S:" /* seed */ - "h" /* help */ - "V", /* version */ + "V" /* version */ + "h", /* help */ long_options, (int *) 0)) != EOF) { switch (c) @@ -399,13 +397,13 @@ decode_switches (int argc, char **argv) case 'S': /* --seed */ S = atoi(optarg); break; - case 'V': - cout << program_name << " " << VERSION << endl; - exit (0); - case 'h': usage (0); + case 'V': + cout << program_name << " version " << VERSION << endl; + exit (0); + default: usage (-1); } diff --git a/utils/SL/CMakeLists.txt b/utils/SL/CMakeLists.txt index dfbbb1e1f..5fbeb8977 100644 --- a/utils/SL/CMakeLists.txt +++ b/utils/SL/CMakeLists.txt @@ -28,7 +28,7 @@ set(common_INCLUDE ${CMAKE_CURRENT_SOURCE_DIR}/..) add_executable(slcheck slcheck.cc) -add_executable(slabchk slabchk.cc Model1d.cc) +add_executable(slabchk slabchk.cc Model1d.cc Orbit1d.cc) add_executable(slshift slshift.cc SLSphere.cc) add_executable(orthochk orthochk.cc) add_executable(diskpot diskpot.cc CylindricalDisk.cc SLSphere.cc) diff --git a/utils/SL/Model1d.H b/utils/SL/Model1d.H index 7af14197c..9fc747601 100644 --- a/utils/SL/Model1d.H +++ b/utils/SL/Model1d.H @@ -5,9 +5,23 @@ #include "interp.H" -class OneDModel +class OneDModel : public std::enable_shared_from_this { +protected: + + int Nfgrid = 800; // Number of points for tabulating + // distribution function for + // interpolation in getE(M) + + Linear1d dfgrid; // Interpolation object for + // distribution function grid (for + // getE(M)) + + std::vector Egrid, Mgrid; + + void compute_dfgrid(void); + public: bool dist_defined; @@ -38,6 +52,65 @@ public: virtual double distf(const double, const double V=0.0) = 0; virtual double dfde (const double, const double V=0.0) = 0; virtual double dfdv (const double, const double V=0.0) = 0; + + //! Model types + enum class ModelType {table, lowiso, sech2, sech2u, sech2halo, sech2mu, + uniform, cosine}; + + + //! Model names for sysout display (avoid temporaries and dynamic + //! allocations on every class instantiation) + inline static const std::map model_names = { + {ModelType::table, "file"}, + {ModelType::lowiso, "LowIso"}, + {ModelType::sech2, "Sech2"}, + {ModelType::sech2u, "Sech2u"}, + {ModelType::sech2halo, "Sech2Halo"}, + {ModelType::sech2mu, "Sech2(mu)"}, + {ModelType::uniform, "Uniform"}, + {ModelType::cosine, "Cosine"} + }; + + //! Model types for command line parsing + //! 'inline static const' allows initialization inside the class definition. + //! 'std::string_view' keys completely prevent temporary string creations. + inline static const std::map model_map = { + {"table", ModelType::table}, + {"lowiso", ModelType::lowiso}, + {"sech2", ModelType::sech2}, + {"sech2u", ModelType::sech2u}, + {"sech2halo", ModelType::sech2halo}, + {"sech2mu", ModelType::sech2mu}, + {"uniform", ModelType::uniform}, + {"cosine", ModelType::cosine} + }; + + //! Model type for sysout display + static std::string get_model_name(ModelType type) { + return std::string(model_names.at(type)); + } + + //! Command line parsing for model types + static ModelType parse_model_type(const std::string& str) { + std::string MODEL = str; + std::transform(MODEL.begin(), MODEL.end(), MODEL.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // A std::string implicitly and cheaply converts to a std::string_view, + // so finding it in our map works cleanly and creates zero temporaries. + auto it = model_map.find(MODEL); + if (it != model_map.end()) { + return it->second; + } else { + std::ostringstream oss; + oss << "Invalid model type: " << str << ". Valid types are: "; + for (const auto& pair : model_map) { + oss << pair.first << " "; + } + throw std::invalid_argument(oss.str()); + } + } + }; class OneDModelTable : public OneDModel @@ -427,6 +500,382 @@ public: if (!model_computed) reset(); return -exp(-E/dispz - 0.5*p*p/dispx)*p/dispx * norm; } +}; + +class Uniform : public OneDModel +{ +private: + // Half-thickness of the uniform slab + double H; + // Velocity widths in x and y directions (for distribution function) + double Vx, Vy; + // Vertical frequency for the uniform slab potential + double Omega; + // Normalization constant for the distribution function + double norm; + +public: + + Uniform(void) + { + Vx = 1.0; + Vy = 1.0; + // + // Units: G = Sigma_o = 1 + // + H = 1.0; + Omega = std::sqrt(2.0*M_PI/H); + norm = 1.0/(std::sqrt(2.0)*M_PI*H); + // ModelID = "Uniform"; + } + + Uniform(double h, double Vx) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vx; + // + // Units: G = Sigma_o = 1 + // + Omega = std::sqrt(2.0*M_PI/H); + norm = 1.0/(std::sqrt(2.0)*M_PI*H); + // ModelID = "Uniform"; + } + + Uniform(double h, double Vx, double Vy) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vy; + // + // Units: G = Sigma_o = 1 + // + Omega = std::sqrt(2.0*M_PI/H); + norm = 1.0/(std::sqrt(2.0)*M_PI*H); + // ModelID = "Uniform"; + } + + double get_mass(const double z) + { + if (z<-H) + return 0.0; + else if (z > H) + return 1.0; + else { + return (z + H)/(2.0*H); + } + } + + double get_density(const double z) + { + if (fabs(z) > H) + return 0.0; + else + return 1.0/(2.0*H); + } + + double get_pot(const double z) + { + // Offset defined so that potential is zero at |z|=H, and negative + // for |z| H) + return 2.0*M_PI*zz - 2.0*M_PI*H; + else + return M_PI*zz*zz/H - M_PI*H; + } + + double get_dpot(const double z) + { + double zz = fabs(z); + if (zz > H) + return 2.0*M_PI*z/zz; + else + return 2.0*M_PI*z/H; + } + + double get_dpot2(const double z) + { + double zz = fabs(z); + if (zz > H) + return 0.0; + else + return 2.0*M_PI/H; + } + + void get_pot_dpot(const double z, double& p, double& dp) + { + double zz = fabs(z); + if (zz > H) { + p = 2.0*M_PI*zz - 2.0*M_PI*H; + dp = 2.0*M_PI*z/zz; + } else { + p = M_PI*zz*zz/H - M_PI*H; + dp = 2.0*M_PI*z/H; + } + } + + double get_mass(const double x, const double y, const double z) + { + return get_mass(z); + } + + double get_density(const double x, const double y, const double z) + { + return get_density(z); + } + + double get_pot(const double x, const double y, const double z) + { + return get_pot(z); + } + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return H; } + double get_scale_height(void) { return H; } + + double distf(const double E, const double p=0.0) + { + // Convert energy to potential, so that the distribution function + // is nonzero only. That is, assume E = Omega*Ja, not zero at +/- + // H. + double OJ = fabs(E + M_PI*H); + return norm / std::sqrt(OJ); + } + + double dfde(const double E, const double p=0.0) + { + // Convert energy to potential, so that the distribution function + // is nonzero only. That is, assume E = Omega*Ja, not zero at +/- + // H. + double OJ = fabs(E + M_PI*H); + return -0.5*Omega*norm / std::pow(OJ, 1.5); + } + + double dfdv(const double E, const double p=0.0) + { + // The distribution function is independent of velocity, so this + // derivative is zero. There are delta functions at v=-Vx and + // v=+Vx, but those are not captured by this simple functional + // form. + return 0.0; + } + +}; + +//! A simple consine-bell slab model with finite thickness +class Cosine : public OneDModel +{ +private: + //! Half-thickness of the uniform slab + double H; + + //! Central density of the cosine slab + double rho0; + + //! Potential zero point is defined at |z|=H, so that potential is + //! negative inside the slab and zero at the edges + double Phi0; + + //! Velocity widths in x and y directions (for distribution function) + double Vx, Vy; + + //! Vertical frequency for the uniform slab potential + double Omega; + + //! Normalization constant for the distribution function + double norm; + + //! Regularization parameter for the distribution function (if needed) + double epsilon; + + //! Number of points in the tabulated distribution function + int numdf; + + //! Number of quadrature points for computing DF + int nint; + + //! Tabulated action and distribution function for interpolation + std::vector Etab, Ftab; + + //! DF table interpolation object + Linear1d df_interp; + + //! Precompute the distribution function table for interpolation + void compute_model(); + +public: + + //! Null constructor with default parameters + Cosine(void) + { + Vx = 1.0; + Vy = 1.0; + // + // Units: G = Sigma_o = 1 + // + H = 1.0; + rho0 = 1.0/H; + Phi0 = M_PI*H*(4.0/(M_PI*M_PI) + 1.0); + // ModelID = "Cosine"; + numdf = 800; + nint = 800; + dist_defined = true; + compute_model(); + } + + //! Constructor with specified half-thickness and velocity width + //! (same for x and y) + Cosine(double h, double Vx, double epsilon) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vx; + this->epsilon = epsilon; + // + // Units: G = Sigma_o = 1 + // + rho0 = 1.0/H; + Phi0 = M_PI*H*(4.0/(M_PI*M_PI) + 1.0); + // ModelID = "Cosine"; + numdf = 800; + nint = 800; + dist_defined = true; + compute_model(); + } + + //! Constructor with specified half-thickness and velocity widths in + //! x and y + Cosine(double h, double Vx, double Vy, double epsilon) + { + this->H = h; + this->Vx = Vx; + this->Vy = Vy; + this->epsilon = epsilon; + // + // Units: G = Sigma_o = 1 + // + rho0 = 1.0/H; + Phi0 = M_PI*H*(4.0/(M_PI*M_PI) + 1.0); + // ModelID = "Cosine"; + numdf = 800; + nint = 800; + dist_defined = true; + compute_model(); + } + + void setDF(int numdf, int nint=800) + { + this->numdf = numdf; + this->nint = nint; + compute_model(); + } + + double get_mass(const double z) + { + if (z<-H) + return 0.0; + else if (z > H) + return 1.0; + else { + return 0.5*rho0*(z + H + H/M_PI*sin(M_PI*z/H)); + } + } + + double get_density(const double z) + { + if (fabs(z) > H) + return 0.0; + else { + double cosfac = cos(M_PI*z/(2.0*H)); + return rho0*cosfac*cosfac; + } + } + + double get_pot(const double z) + { + // Offset defined so that potential is zero at |z|=H, and negative + // for |z| H) + return 2.0*M_PI*rho0*H*(zz - H); + else + return 2.0*M_PI*rho0*(0.5*z*z + H*H/(M_PI*M_PI)*(1 - cos(M_PI*z/H))) - Phi0; + } + + double get_dpot(const double z) + { + double zz = fabs(z); + if (zz > H) + return 2.0*M_PI*rho0*H*z/zz; + else + return 2.0*M_PI*rho0*(z + H/M_PI*sin(M_PI*z/H)); + } + + double get_dpot2(const double z) + { + double zz = fabs(z); + if (zz > H) + return 0.0; + else + return 2.0*M_PI*rho0*(1 + cos(M_PI*z/H)); + } + + void get_pot_dpot(const double z, double& p, double& dp) + { + double zz = fabs(z); + if (zz > H) { + p = 2.0*M_PI*rho0*H*(zz - H); + dp = 2.0*M_PI*rho0*H*z/zz; + } else { + p = 2.0*M_PI*rho0*(0.5*z*z + H*H/(M_PI*M_PI)*(1 - cos(M_PI*z/H))) - Phi0; + dp = 2.0*M_PI*rho0*(z + H/M_PI*sin(M_PI*z/H)); + } + } + + double get_mass(const double x, const double y, const double z) + { + return get_mass(z); + } + + double get_density(const double x, const double y, const double z) + { + return get_density(z); + } + + double get_pot(const double x, const double y, const double z) + { + return get_pot(z); + } + + double get_min_radius(void) { return 0.0; } + double get_max_radius(void) { return H; } + double get_scale_height(void) { return H; } + + double distf(const double E, const double p=0.0) + { + if (E>=0.0) + return 0.0; + else + return df_interp.eval(E); + } + + double dfde(const double E, const double p=0.0) + { + if (E>=0.0) + return 0.0; + else + return df_interp.deriv(E); + } + + double dfdv(const double E, const double p=0.0) + { + // The distribution function is independent of velocity, so this + // derivative is zero. There are delta functions at v=-Vx and + // v=+Vx, but those are not captured by this simple functional + // form. + return 0.0; + } }; diff --git a/utils/SL/Model1d.cc b/utils/SL/Model1d.cc index cea2faa16..1dfc409ef 100644 --- a/utils/SL/Model1d.cc +++ b/utils/SL/Model1d.cc @@ -2,11 +2,12 @@ #include #include #include +#include #include #include "gaussQ.H" #include "interp.H" - +#include "Orbit1d.H" #include "Model1d.H" #include "RK4.H" @@ -331,3 +332,233 @@ void Sech2Halo::reset() model_computed = true; dist_defined = true; } + +void Cosine::compute_model(void) +{ + double Emin = get_pot(0.0); + double Emax = get_pot(get_max_radius()); + + Etab.resize(numdf+1); + Ftab.resize(numdf+1); + + const int numz = 10000; + std::vector pot(numz+1), z(numz+1); + for (int i=0; i<=numz; i++) { + z[i] = H * i/numz; + pot[i] = get_pot(z[i]); + } + +#ifdef DEBUG + { + std::ofstream out("cosine_pot.dat"); + if (out) { + out << "#" + << std::setw(15) << std::right << "Z" + << std::setw(16) << std::right << "Potential" + << std::endl; + for (int i=0; i<=numz; i++) { + out << std::setw(16) << std::right << z[i] + << std::setw(16) << std::right << pot[i] + << std::endl; + } + } else { + std::cerr << "Could not open cosine_pot.dat for writing" << std::endl; + } + } +#endif + + Linear1d getZ(pot, z); + + const int numL = 400; + double fac = rho0/(2.0*sqrt(2.0)*H); + LegeQuad lege(numL); + + for (int i=0; i<=numdf; i++) { + Etab[i] = Emin + (Emax - Emin)*i/numdf; + Ftab[i] = 0.0; + double umax = sqrt(Emax - Etab[i]); + for (int j=0; j Z(numZ+1), Rho(numZ+1, 0.0); + const int nint = 800; + LegeQuad lg(nint); + + for (int i=0; i<=numZ; i++) { + Z[i] = get_max_radius() * i/numZ; + double P = get_pot(Z[i]); + double vmax = sqrt(2.0*(Emax - P)); + for (int j=0; j maxdif) { + maxdif = dif; + maxz = Z[i]; + } + } + std::cout << "Max density difference from DF integration: " << maxdif + << " at z=" << maxz << std::endl; + + std::ofstream out("cosine_df_check_density.dat"); + if (out) { + out << "#" + << std::setw(15) << std::right << "Z" + << std::setw(16) << std::right << "Density" + << std::setw(16) << std::right << "DF Density" + << std::setw(16) << std::right << "Difference" + << std::endl; + for (int i=0; i<=numZ; i++) { + out << std::setw(16) << std::right << Z[i] + << std::setw(16) << std::right << get_density(Z[i]) + << std::setw(16) << std::right << Rho[i] + << std::setw(16) << std::right << Rho[i] - get_density(Z[i]) + << std::endl; + } + + } else { + std::cerr << "Could not open cosine_df_check.dat for writing" << std::endl; + } + + out.close(); + out.open("cosine_df_check_harmonic.dat"); + if (out) { + out << "#" + << std::setw(15) << std::right << "Z" + << std::setw(16) << std::right << "E" + << std::setw(16) << std::right << "DF" + << std::setw(16) << std::right << "Harmonic DF" + << std::setw(16) << std::right << "Difference" + << std::endl; + + double Omega = sqrt(4.0*M_PI*rho0); + double Sigma2 = 8.0*rho0*H*H/M_PI; + + for (int i=0; i<=numdf; i++) { + double Z = getZ(Etab[i]); + double E = get_pot(Z) + Phi0; + double DF = df_interp(Etab[i]); + double DFh = rho0*Omega/(2.0*M_PI*Sigma2)*exp(-E/Sigma2); + + out << std::setw(16) << std::right << getZ(Z) + << std::setw(16) << std::right << Etab[i] + << std::setw(16) << std::right << DF + << std::setw(16) << std::right << DFh + << std::setw(16) << std::right << DF - DFh + << std::endl; + } + } else { + std::cerr << "Could not open cosine_df_check_harmonic.dat for writing" << std::endl; + } + } +#endif + // END DEBUG +} + +void OneDModel::compute_dfgrid(void) +{ + if (not dist_defined) + throw std::runtime_error("Model must have a DF defined in compute_dfgrid()"); + + double zmax = get_max_radius(); + double Emin = get_pot(0.0), Emax = get_pot(zmax); + double E, dE = (Emax - Emin)/double(Nfgrid+1); + + Egrid.resize(Nfgrid+1); + Mgrid.resize(Nfgrid+1); + + for (int i=0; i<=Nfgrid; i++) Egrid[i] = Emin + dE*i; + Mgrid[0] = 0.0; + + const int Nl = 32; + LegeQuad wkE(Nl); + + auto safeSelf = std::static_pointer_cast(shared_from_this()); + OneDOrbit orb(safeSelf); + + for (int i=1; i<=Nfgrid; i++) { + Mgrid[i] = Mgrid[i-1]; + for (int j=1; j +#include +#include + +#include + +#include "orbit.H" +// #include "massmodel1d.H" + +class OneDModel; +class OneDBiorth; + +class OneDOrbit : public RegularOrbit +{ + + using Vector = Eigen::VectorXd; + using Matrix = Eigen::MatrixXd; + using CVector = Eigen::VectorXcd; + using CMatrix = Eigen::MatrixXcd; + using Complex = std::complex; + + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + +private: + double energy; + double vp; // keep track orbit information in in-plane + // direction + + // computed in compute_freq + double zmax; + double freq; + double action; + + std::shared_ptr model; + std::shared_ptr biorth; + + void compute_action(void) { compute_freq(); } + void compute_freq(void); + void compute_angles(void); + void compute_biorth(void); + + struct OneD_Angle_Grid { + int size; + Vector time; + Matrix angl; + Matrix hght; + } angle_grid; + + constexpr static double TOL = 1.0e-10; + + // Number of points in angle grid + int NGRID; + // Number of knots for angle integrals + int NINT; + // Number of points for frequency integral + int NFREQ; + // Number of points for potential transform + int NREC; + + std::shared_ptr lkw; // Cache this for computing angles + + bool pot_cache = false; // Cache potential transforms for a + // given biorthogonal set + + std::vector pot_eval; // Cache potential evaluations for pot_trans + + +public: + + enum class Type { Angle, Height, Velocity }; + + //@{ + //! Constructors + OneDOrbit(void); + OneDOrbit(std::shared_ptr model); + OneDOrbit(std::shared_ptr model, double Energy, double vparallel=0.0); + //@} + + virtual ~OneDOrbit() = default; + + //! Copy constructor and assignment operator (defaulted) + OneDOrbit(OneDOrbit &orb); + // OneDOrbit &operator=(OneDOrbit &); + + // Required functions + + double get_action(const int i=0) { + if (action_defined == false) compute_action(); + return action; } + + double get_freq(const int i=0) { + if (freq_defined == false) compute_freq(); + return freq; } + + double get_angle(const int i, const double time); + + double get_angle(const Type f, const double time); + + // Specific functions + + void new_orbit(double Energy, double vparallel=0.0); + + void set_biorth(std::shared_ptr type) { + biorth = type; + biorth_defined = true; + } + + void set_numerical_params(int nrec, int ngrid=100, int nfreq=64, + int nint=20) + { NREC = nrec; NGRID = ngrid; NFREQ = nfreq; NINT = nint;} + + double pot_trans(const int k, double (*func)(double) ); + double pot_trans(const int k, const int n); + void pot_trans(const int k, Eigen::VectorXd& t); + + Complex pot_trans_complex(const int k, std::function ); + Complex pot_trans_complex(const int k, const int n); + Complex pot_trans_complex_check(const int k, const int n); + Complex pot_trans_complex_check(const int k, const int nx, const int n); + void pot_trans_complex(const int k, Eigen::VectorXcd& t); + void set_pot_cache(bool val, int n=0); + + double pot_trans(const int kx, const int ky, double (*func)(double) ); + double pot_trans(const int kx, const int ky, const int n); + void pot_trans(const int kx, const int ky, Eigen::VectorXd& t); + + Complex pot_trans_complex(const int kx, const int ky, const int kz, + double (*func)(double) ); + Complex pot_trans_complex(const int kx, const int ky, const int kz, + const int n); + void pot_trans_complex(const int kx, const int ky, const int kz, + Eigen::VectorXcd& t); + void pot_trans_complex_check(const int kx, const int ky, const int kz, + Eigen::VectorXcd& t); + + //@{ + // Safe access + const double Energy(void) { return energy; } + const double Vpara(void) { return vp; } + + const double apo(void) { + if (freq_defined == false) compute_freq(); + return zmax; + } + + OneDModel& modl(void) { return *model; } + + OneDBiorth& orth(void) { return *biorth; } + //@} +}; + +#endif diff --git a/utils/SL/Orbit1d.cc b/utils/SL/Orbit1d.cc new file mode 100644 index 000000000..785c23e46 --- /dev/null +++ b/utils/SL/Orbit1d.cc @@ -0,0 +1,595 @@ +// This routine computes orbit information for 1-d slab + +// #define NO_SPLINE 1 +#undef NO_SPLINE + +#include + +#include "numerical.H" +#include "gaussQ.H" +#include "interp.H" +#include "biorth1d.H" +// #include "massmodel1d.H" +#include "Orbit1d.H" +#include "Model1d.H" + +using namespace std::literals::complex_literals; + +class fct_zmax : public std::function +{ +private: + + std::shared_ptr mdl; + double E; + +public: + + fct_zmax(std::shared_ptr p, double e) : mdl(p), E(e) {} + + double operator()(double z) + { + return E - mdl->get_pot(z); + } +}; + + +OneDOrbit::OneDOrbit(void) +{ + model_defined = false; + action_defined = false; + freq_defined = false; + angle_defined = false; + biorth_defined = false; + + NGRID = 100; + NINT = 20; + NFREQ = 64; + NREC = 64; + + dof = 1; + OrbitID = "OneDimensionalOrbit"; +} + +OneDOrbit::OneDOrbit(std::shared_ptr MODEL) +{ + model = MODEL; + + model_defined = true; + action_defined = false; + freq_defined = false; + angle_defined = false; + biorth_defined = false; + + NGRID = 100; + NINT = 20; + NFREQ = 64; + NREC = 64; + + dof = 1; + OrbitID = "OneDimensionalOrbit"; +} + +OneDOrbit::OneDOrbit(std::shared_ptr MODEL, double ENERGY, double VPARA) +{ + model = MODEL; + energy = ENERGY; + vp = VPARA; + + model_defined = true; + action_defined = false; + freq_defined = false; + angle_defined = false; + biorth_defined = false; + + NGRID = 100; + NINT = 20; + NFREQ = 64; + NREC = 64; + + dof = 1; + OrbitID = "OneDimensionalOrbit"; +} + +OneDOrbit::OneDOrbit(OneDOrbit& p) +{ + model = p.model; + energy = p.energy; + vp = p.vp; + + zmax = p.zmax; + freq = p.freq; + action = p.action; + + model = p.model; + biorth = p.biorth; + + model_defined = p.model_defined; + action_defined = p.action_defined; + freq_defined = p.freq_defined; + angle_defined = p.angle_defined; + biorth_defined = p.biorth_defined; + + NGRID = p.NGRID; + NINT = p.NINT; + NFREQ = p.NFREQ; + NREC = p.NREC; + + lkw = p.lkw; + pot_cache = p.pot_cache; + pot_eval = p.pot_eval; + + dof = 1; + OrbitID = "OneDimensionalOrbit"; +} + + +void OneDOrbit::new_orbit(double ENERGY, double VPARA) +{ + energy = ENERGY; + vp = VPARA; + + freq_defined = false; + angle_defined = false; + if (pot_cache) pot_eval.clear(); +} + + +// Perform integral using centered rectangles + +void OneDOrbit::compute_freq(void) +{ + double E = energy; + fct_zmax fct(model, E); + + double mmin = model->get_min_radius(); + double mmax = model->get_max_radius(); + + zmax = zbrent(fct, mmin, 4.0*mmax, TOL); + + if (energy-model->get_pot(mmin) < TOL || zmax < TOL) { + freq = sqrt(4.0*M_PI*model->get_density(mmin)); + return; + } + + double dt = M_PI/(2.0*NFREQ); + + double ansf = 0.0; + double ansi = 0.0; + for (int j=0; jget_pot(z)) ); + ansi += cos(t) * sqrt( 2.0*(E-model->get_pot(z)) ); + } + + freq = 0.5*M_PI/(ansf*zmax*dt); + action = 4.0*ansi*zmax*dt; + + freq_defined = true; +} + + + +void OneDOrbit::compute_angles(void) +{ + if (!freq_defined) compute_freq(); + + double E = energy; + + // Gaussian integration (prevent recomputation) + // + if (lkw == nullptr or lkw->get_n() != NINT) + lkw = std::make_shared(NINT); + + angle_grid.size = NGRID; + angle_grid.time.resize(NGRID+1); + angle_grid.angl.resize(2, NGRID+1); + angle_grid.hght.resize(2, NGRID+1); + + double z, dz = zmax/NGRID; + + angle_grid.time(0) = 0.0; + angle_grid.angl(0, 0) = 0.0; + angle_grid.hght(0, 0) = 0.0; + + double ans, tmin, tmax, t; + + for (int i=1; i<=NGRID; i++) { + + z = angle_grid.hght(0, i) = dz*i; + + tmin = asin(dz*(i-1)/zmax); + if (dz*i/zmax>=1.0) + tmax = 0.5*M_PI; + else + tmax = asin(dz*i/zmax); + + ans = 0.0; + for (int j=0; jknot(j); + double denom = 2.0*(E-model->get_pot(zmax*sin(t))); + if (denom<=0.0) { + double dt = 0.5*M_PI - t; + denom = model->get_dpot(zmax) * zmax * dt*dt; + } + if (denom > 0.0) + ans += lkw->weight(j) * cos(t) / sqrt(denom); + else + std::cerr << "Warning: negative denominator in compute_angles: " + << "E=" << E << " pot=" << model->get_pot(zmax*sin(t)) + << " denom=" << denom << std::endl; + } + ans *= zmax*(tmax - tmin); + + angle_grid.time(i) = angle_grid.time(i-1) + ans; + angle_grid.angl(0, i) = freq * angle_grid.time(i); + } + +#ifndef NO_SPLINE + Eigen::VectorXd work(NGRID+1); + Spline(angle_grid.time, angle_grid.angl.row(0), -1.0e30, -1.0e30, work); + angle_grid.angl.row(1) = work; + + Spline(angle_grid.time, angle_grid.hght.row(0), -1.0e30, -1.0e30, work); + angle_grid.hght.row(1) = work; +#endif + + angle_defined = true; +} + + +double OneDOrbit::get_angle(const Type f, const double T) +{ + switch (f) { + case Type::Angle: + return get_angle(1, T); + case Type::Height: + return get_angle(2, T); + case Type::Velocity: + return get_angle(3, T); + default: + throw "OneDOrbit::get_angle: invalid angle type"; + } +} + +double OneDOrbit::get_angle(const int i, const double T) +{ + if (!angle_defined) compute_angles(); + + // Compute the time modulo the quarter period and identify the + // quadrant + // + double time, qperiod = 0.5*M_PI/freq; + + if (T<0.0) + time = T + 4.0*qperiod * (1.0 + (int)(0.25*fabs(T)/qperiod)); + else + time = T - 4.0*qperiod * (int)(0.25*fabs(T)/qperiod); + + int phase = (int)(time/qperiod); + double ptime = time - phase*qperiod; + + phase = phase % 4; + + double zsign = 1.0; + if (phase > 1) zsign = -1.0; + if (phase == 1 || phase == 3) ptime = qperiod - ptime; + + double ans, vel; + + switch (i) { + case 3: +#ifndef NO_SPLINE + Splint1(angle_grid.time, angle_grid.hght.row(0), angle_grid.hght.row(1), ptime, ans); +#else + ans = odd2(ptime, angle_grid.time, angle_grid.hght); +#endif + vel = std::sqrt(2.0*std::fabs((energy - model->get_pot(ans)))); + if (phase == 1 || phase == 2) vel = -vel; + /* + if (std::fabs(energy - model->get_pot(ans) - 0.5*vel*vel) > 1.0e-6) { + std::cout << "Velocity error: energy=" << energy + << " pot=" << model->get_pot(ans) + << " kinetic=" << 0.5*vel*vel + << " error=" << energy - model->get_pot(ans) - 0.5*vel*vel + << std::endl;return vel; + } + */ + return vel; + case 2: +#ifndef NO_SPLINE + Splint1(angle_grid.time, angle_grid.hght.row(0), angle_grid.hght.row(1), ptime, ans); +#else + ans = odd2(ptime, angle_grid.time, angle_grid.hght); +#endif + ans *= zsign; + break; + case 1: + default: +#ifndef NO_SPLINE + Splint1(angle_grid.time, angle_grid.angl.row(0), angle_grid.angl.row(1), ptime, ans); +#else + ans = odd2(ptime, angle_grid.time, angle_grid.angl); +#endif + switch (phase) { + case 1: + ans = M_PI - ans; + break; + case 2: + ans = M_PI + ans; + break; + case 3: + ans = 2.0*M_PI - ans; + break; + } + break; + } + + return ans; +} + +double OneDOrbit::pot_trans(const int k, double (*func)(double z)) +{ + double freq = get_freq(); + double qperiod = 0.5*M_PI/freq; + + // Use centered rectangles + double t, dt = 2.0*qperiod/NREC; + double ans = 0.0; + for (int i=0; ipotl(n, 0, get_angle(2, t)); + } + + ans *= dt*freq/M_PI; + + return ans; +} + +void OneDOrbit::pot_trans(const int k, Eigen::VectorXd& t) +{ + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling pot_trans(int, int)"); + + double freq = get_freq(); + double qperiod = 0.5*M_PI/freq; + + // Use centered rectangles + double T, dT = 2.0*qperiod/NREC; + Eigen::VectorXd p(t.size()); t.setZero(); + + for (int i=0; ipotl_vec(0, 0, get_angle(2, T), p); + biorth->potl(0, 0, get_angle(2, T), p); + t += cos(T*freq*k) * p; + } + + t *= dT*freq/M_PI; +} + +std::complex OneDOrbit::pot_trans_complex +(const int k, std::function func) +{ + double freq = get_freq(); + double hperiod = M_PI/freq; + + // Use centered rectangles + double t, dt = 2.0*hperiod/NREC; + Complex ans = 0.0; + + for (int i=0; i(k)) * func(get_angle(2, t)); + } + + ans *= 0.5*dt*freq/M_PI; + + return ans; +} + + +std::complex OneDOrbit::pot_trans_complex(const int k, const int n) +{ + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling " + "pot_trans_complex(int, int)"); + + double freq = get_freq(); + double hperiod = M_PI/freq; + + // Result + // + Complex ans = 0.0; + + // Use centered rectangles + // + double dt = 2.0*hperiod/NREC; + for (int i=0; i(k)) * + biorth->potl(n, 0, get_angle(2, t)); + } + + return ans * 0.5*dt*freq/M_PI; +} + +std::complex OneDOrbit::pot_trans_complex_check(const int k, const int n) +{ + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling " + "pot_trans_complex_check(int, int)"); + + double freq = get_freq(); + + // Result + Complex ans = 0.0; + + // Use centered rectangles + double dw = 2.0*M_PI/NREC; + + for (int i=0; i(k)) * + biorth->potl(n, 0, get_angle(2, w/freq)); + } + + return ans/double(NREC); +} + +std::complex OneDOrbit::pot_trans_complex_check +(const int k, const int nx, const int n) +{ + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling " + "pot_trans_complex_check(int, int)"); + + double freq = get_freq(); + + // Result + Complex ans = 0.0; + + // Use centered rectangles + double dw = 2.0*M_PI/NREC; + + for (int i=0; i(k)) * + biorth->potl(n, nx, get_angle(OneDOrbit::Type::Height, w/freq)); + } + + return ans/double(NREC); +} + +void OneDOrbit::set_pot_cache(bool val, int norder) +{ + pot_eval.clear(); + pot_cache = true; + if (val == false) { + pot_cache = false; + return; + } + + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling pot_trans(int, int)"); + + double freq = get_freq(); + double hperiod = M_PI/freq; + + // Use centered rectangles + double T, dT = 2.0*hperiod/NREC; + Eigen::VectorXd p(norder); + + pot_eval.resize(NREC); + for (int i=0; ipotl(0, 0, get_angle(2, T), p); + pot_eval[i] = p; + } +} + +void OneDOrbit::pot_trans_complex(const int k, Eigen::VectorXcd& t) +{ + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling pot_trans(int, int)"); + + double freq = get_freq(); + double hperiod = M_PI/freq; + + // Use centered rectangles + double dT = 2.0*hperiod/NREC; + + Eigen::VectorXd p(t.size()); // Evaluation vector for potential transform + t.setZero(); // Zero the return vector + + for (int i=0; i(k)); + if (pot_cache) { + t += factor * pot_eval[i].cast >(); + } else { + biorth->potl(0, 0, get_angle(2, T), p); + t += factor * p.cast >(); + } + } + + t *= 0.5*dT*freq/M_PI; +} + +void OneDOrbit::pot_trans_complex(const int nx, const int ny, const int nz, + Eigen::VectorXcd& t) +{ + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling pot_trans(int, int)"); + + double freq = get_freq(); + double hperiod = M_PI/freq; + // Use centered rectangles + double dT = 2.0*hperiod/NREC; + + t.setZero(); // Zero the return vector + + Eigen::VectorXd p(t.size()); + + for (int i=0; ipotl(nx, ny, get_angle(2, T), p); + t += exp(-1i*T*freq*static_cast(nz)) * + p.cast >(); + } + + t *= 0.5*dT*freq/M_PI; +} + +// This version is for checking the results of the above routine +// +void OneDOrbit::pot_trans_complex_check +(const int nx, const int ny, const int nz, Eigen::VectorXcd& t) +{ + if (!biorth_defined) + bomb ("Must define biorthogonal set before calling pot_trans(int, int)"); + + double freq = get_freq(); + double hperiod = M_PI/freq; + + // Zero the return vector + // + t.setZero(); + + // Use centered rectangles + // + Eigen::VectorXd p(t.size()); + double dw = 2.0*M_PI/NREC; + for (int i=0; ipotl(nx, ny, z, p); + t += exp(-1i*w*static_cast(nz)) * p.cast>(); + } + + t /= NREC; +} + diff --git a/utils/SL/orthochk.cc b/utils/SL/orthochk.cc index 3e6a86b0c..72295fffe 100644 --- a/utils/SL/orthochk.cc +++ b/utils/SL/orthochk.cc @@ -4,169 +4,90 @@ #include #include -#include - #include "biorth1d.H" #include "SLGridMP2.H" #include "gaussQ.H" #include "localmpi.H" +#include "cxxopts.H" -//=========================================================================== - -void usage(char *prog) -{ - cout << "Usage:\n\n" - << prog << " [options]\n\n" - << setw(15) << "Option" << setw(10) << "Argument" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Description" << endl << endl - << resetiosflags(ios::left) - << setw(15) << "-m or --mpi" << setw(10) << "No" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Turn on MPI for SL computation" << endl - << resetiosflags(ios::left) - << setw(15) << "-t or --Trig" << setw(10) << "No" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Use trigonometric basis" << endl - << resetiosflags(ios::left) - << setw(15) << "-s or --SL" << setw(10) << "No" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Use Sturm-Liouville basis" << endl - << setw(15) << "-T or --type" << setw(10) << "string" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Density target (isothermal, constant, parabolic)" << endl - << resetiosflags(ios::left) - << setw(15) << "-n " << setw(10) << "int" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Number of basis functions" << endl - << resetiosflags(ios::left) - << setw(15) << "-H " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Slab scale height" << endl - << resetiosflags(ios::left) - << setw(15) << "-k " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Wave number for Trig basis" << endl - << resetiosflags(ios::left) - << setw(15) << "-x " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Wave number in X for SL basis" << endl - << resetiosflags(ios::left) - << setw(15) << "-y " << setw(10) << "double" << setw(10) << " " - << setiosflags(ios::left) - << setw(40) << "Wave number in Y for SL basis" << endl - << resetiosflags(ios::left) - << "" << endl; - - exit(0); -} enum BioType1d {Trig, SL}; +std::map bioTypeMap = { + {"trig", Trig}, + {"sl", SL} +}; + +std::map bioStrMap = { + {Trig, "trig"}, + {SL, "sl"} +}; + int main(int argc, char** argv) { - bool use_mpi = false; double KX = 0.5; double H = 0.1; double ZMAX = 1.0; int NMAX = 10; int IKX = 1; int IKY = 3; - BioType1d Type = Trig; + int numz = 2000; + std::string bioTypeStr = "trig"; + std::string cachename = ".slab_sl_cache"; + std::string cmap = "linear"; std::string slabID = "iso"; + bool use_mpi = false; - int c; - while (1) { - int this_option_optind = optind ? optind : 1; - int option_index = 0; - static struct option long_options[] = { - {"mpi", 0, 0, 0}, - {"Trig", 0, 0, 0}, - {"SL", 0, 0, 0}, - {0, 0, 0, 0} - }; - - c = getopt_long (argc, argv, "msT:tx:y:k:n:z:H:h", - long_options, &option_index); - - if (c == -1) break; - - switch (c) { - case 0: - { - string optname(long_options[option_index].name); - - if (!optname.compare("mpi")) { - use_mpi = true; - } else if (!optname.compare("Trig")) { - Type = Trig; - } else if (!optname.compare("SL")) { - Type = SL; - } else if (!optname.compare("type")) { - slabID = optarg; - } else { - cout << "Option " << long_options[option_index].name; - if (optarg) cout << " with arg " << optarg; - cout << " is not defined " << endl; - exit(0); - } - } - break; - - case 'm': - use_mpi = true; - break; - - case 's': - Type = SL; - break; - - case 'T': - slabID = optarg; - break; - - case 't': - Type = Trig; - break; - - case 'x': - IKX = atoi(optarg); - break; - - case 'y': - IKY = atoi(optarg); - break; + cxxopts::Options options("orthochk", "Check orthogonality of 1D basis functions"); + + options.add_options() + ("m,mpi", "Use MPI") + ("T,type", "Slab type (trig or SL)", cxxopts::value(bioTypeStr)->default_value("trig")) + ("M,matrix", "Print orthogonality matrix for a particular kx, ky choice") + ("x,ikx", "IKX for SLGridSlab (default: 1)", cxxopts::value(IKX)) + ("y,iky", "IKY for SLGridSlab (default: 3)", cxxopts::value(IKY)) + ("k,kx", "KX for OneDTrig (default: 0.5)", cxxopts::value(KX)) + ("z,zmax", "ZMAX for OneDTrig and SLGridSlab (default: 1.0)", cxxopts::value(ZMAX)) + ("H,height", "Scale height H for SLGridSlab (default: 0.1)", cxxopts::value(H)) + ("n,nmax", "NMAX for SLGridSlab (default: 10)", cxxopts::value(NMAX)) + ("N,numz", "Number of z points for SLGridSlab (default: 2000)", cxxopts::value(numz)->default_value("2000")) + ("C,cmap", "Coordinate map for SLGridSlab (linear, tanh, or sech; default: linear)", cxxopts::value(cmap)->default_value("linear")) + ("s,slabid", "Slab model ID for SLGridSlab (iso, parabolic, or constant; default: iso)", cxxopts::value(slabID)->default_value("iso")) + ("c,cachename", "Cache file name for SLGridSlab (default: .slab_sl_cache)", cxxopts::value(cachename)->default_value(".slab_sl_cache")) + ("h,help", "Print usage"); + + auto result = options.parse(argc, argv); + + if (result.count("mpi")) { + local_init_mpi(argc, argv); + use_mpi = true; + } - case 'k': - KX = atof(optarg); - break; + if (result.count("help")) { + if (myid==0) + std::cout << options.help() << std::endl; + return 0; + } - case 'z': - ZMAX = atof(optarg); - break; - case 'H': - H = atof(optarg); - break; + // Determine biorthogonal function type + std::transform(bioTypeStr.begin(), bioTypeStr.end(), bioTypeStr.begin(), + [](unsigned char c){ return std::tolower(c); }); - case 'n': - NMAX = atoi(optarg); - break; + BioType1d Type = std::find_if(bioTypeMap.begin(), bioTypeMap.end(), + [bioTypeStr](const std::pair& pair) { return pair.first == bioTypeStr; }) != bioTypeMap.end() ? bioTypeMap[bioTypeStr] : Trig; - case 'h': - default: - usage(argv[0]); - } + // Check for valid coordinate map choice + std::transform(cmap.begin(), cmap.end(), cmap.begin(), + [](unsigned char c){ return std::tolower(c); }); - } - - //=================== - // MPI preliminaries - //=================== - - if (use_mpi) { - local_init_mpi(argc, argv); + std::vector valid_cmaps = { "tanh", "sech", "linear" }; + if (valid_cmaps.end() == + std::find_if(valid_cmaps.begin(), valid_cmaps.end(), + [cmap](const std::string& vmap) { return vmap == cmap; })) { + std::cerr << "Invalid coordinate map choice: " << cmap << std::endl; + return 1; } //=================== @@ -183,14 +104,14 @@ main(int argc, char** argv) case SL: { - const int NUMZ=800; int KMAX = max(IKX+1, IKY+1); SLGridSlab::ZBEG = 0.0; SLGridSlab::ZEND = 0.1; SLGridSlab::H = H; if (use_mpi) SLGridSlab::mpi = 1; - orthoSL = std::make_shared(KMAX, NMAX, NUMZ, ZMAX, slabID, true); + orthoSL = std::make_shared(KMAX, NMAX, numz, ZMAX, cachename, + cmap, slabID, true); } break; @@ -201,6 +122,8 @@ main(int argc, char** argv) } + int ikx = max(IKX, IKY), iky = min(IKX, IKY); + //=================== // Get info //=================== @@ -216,11 +139,12 @@ main(int argc, char** argv) cout << "1: Print out density, potential pairs" << endl; cout << "2: Check density" << endl; cout << "3: Check orthogonality" << endl; - cout << "4: Quit" << endl; + cout << "4: Internal ortho check" << endl; + cout << "5: Quit" << endl; cout << "?? "; cin >> iwhich; - if (iwhich < 1 || iwhich > 4) iwhich = 4; + if (iwhich < 1 || iwhich > 5) iwhich = 5; switch(iwhich) { case 1: @@ -249,16 +173,16 @@ main(int argc, char** argv) if (Type == Trig) { x = ortho->r_to_rb(z); out << setw(15) << z - << setw(15) << ortho->potl(N, i, z) + << setw(15) << ortho->potl (N, i, z) << setw(15) << ortho->force(N, i, z) - << setw(15) << ortho->dens(N, i, z) + << setw(15) << ortho->dens (N, i, z) << endl; } else { x = orthoSL->z_to_xi(z); out << setw(15) << z - << setw(15) << orthoSL->get_pot(x, IKX, IKY, N) - << setw(15) << orthoSL->get_force(x, IKX, IKY, N) - << setw(15) << orthoSL->get_dens(x, IKX, IKY, N) + << setw(15) << orthoSL->get_pot (x, ikx, iky, N) + << setw(15) << orthoSL->get_force(x, ikx, iky, N) + << setw(15) << orthoSL->get_dens (x, ikx, iky, N) << endl; } } @@ -281,64 +205,66 @@ main(int argc, char** argv) int num; cin >> num; - cout << "N? "; - int N; - cin >> N; - - cout << "dz? "; - double dz; - cin >> dz; + cout << "eps? "; + double eps; + cin >> eps; + double dz = 2.0*ZMAX/num; + double h = dz*eps; double x, z, d1, d2, d3; for (int i=0; ir_to_rb(z); - d1 = ( - ortho->potl(N, i, z+dz) - - ortho->potl(N, i, z)*2.0 + - ortho->potl(N, i, z-dz) - ) / (dz*dz); + for (int n=0; nforce(N, i, z-0.5*dz) - - ortho->force(N, i, z+0.5*dz) - ) / dz; + d1 = ( + ortho->potl(n, i, z+h) - + ortho->potl(n, i, z)*2.0 + + ortho->potl(n, i, z-h) + ) / (h*h); + + d2 = ( + ortho->force(n, i, z-0.5*h) - + ortho->force(n, i, z+0.5*h) + ) / h; - d3 = -KX*KX*ortho->potl(N, i, z); + d3 = -KX*KX*ortho->potl(n, i, z); - out << setw(15) << z - << setw(15) << d1+d3 - << setw(15) << d2+d3 - << setw(15) << -ortho->dens(N, i, z) - << endl; + out << setw(15) << d1+d3 + << setw(15) << d2+d3 + << setw(15) << -ortho->dens(n, i, z); + } + } else { - x = orthoSL->z_to_xi(z); - d1 = ( - orthoSL->get_pot(x+dz, IKX, IKY, N) - - orthoSL->get_pot(x, IKX, IKY, N)*2.0 + - orthoSL->get_pot(x-dz, IKX, IKY, N) - ) / (dz*dz); + Eigen::VectorXd pot0(NMAX), potN(NMAX), potP(NMAX), den0(NMAX); + Eigen::VectorXd frcP(NMAX), frcN(NMAX); - d2 = ( - orthoSL->get_force(x+0.5*dz, IKX, IKY, N) - - orthoSL->get_force(x-0.5*dz, IKX, IKY, N) - ) / dz; + orthoSL->get_pot (pot0, z , ikx, iky); + orthoSL->get_pot (potN, z-h, ikx, iky); + orthoSL->get_pot (potP, z+h, ikx, iky); + orthoSL->get_force(frcN, z-0.5*h, ikx, iky); + orthoSL->get_force(frcP, z+0.5*h, ikx, iky); + orthoSL->get_dens (den0, z, ikx, iky); - d3 = -4.0*M_PI*M_PI*(IKX*IKX+IKY*IKY)* - orthoSL->get_pot(x, IKX, IKY, N); + for (int n=0; nget_dens(x, IKX, IKY, N) - << endl; + d1 = (potP(n) - 2.0*pot0(n) + potN(n)) / (h*h); + d2 = (frcP(n) - frcN(n)) / h; + d3 = -4.0*M_PI*M_PI*(IKX*IKX+IKY*IKY) * pot0(n); + out << setw(15) << d1+d3 + << setw(15) << d2+d3 + << setw(15) << den0(n); + } } + out << endl; } } @@ -352,58 +278,101 @@ main(int argc, char** argv) LegeQuad lw(num); - std::cout << "N1, N2? "; - int N1, N2; - std::cin >> N1; - std::cin >> N2; - - double ximin, ximax; - switch (Type) { - case Trig: - ximin = ortho->r_to_rb(-ZMAX); - ximax = ortho->r_to_rb( ZMAX); - break; - case SL: - ximin = orthoSL->z_to_xi(-ZMAX); - ximax = orthoSL->z_to_xi( ZMAX); - break; + if (result.count("matrix")) { + + Eigen::VectorXd pot(NMAX), den(NMAX); + Eigen::MatrixXd ans(NMAX, NMAX); + ans.setZero(); + + for (int i=0; ipotl(i, i, z, pot); + ortho->dens(i, i, z, den); + } else { + orthoSL->get_pot (pot, z, ikx, iky); + orthoSL->get_dens(den, z, ikx, iky); + } + + for (int n1=0; n1> N1; + std::cin >> N2; - double x, r, ans=0.0; - for (int i=0; ipotl(N1, i, x); - double tmp2 = ortho->dens(N1, i, x); - double tmp3 = ortho->d_r_to_rb(x); - - ans += ortho->potl(N1, i, x)* - ortho->dens(N2, i, x) * - ortho->d_r_to_rb(x) * (ximax - ximin)*lw.weight(i); - } - + ximin = ortho->r_to_rb(-ZMAX); + ximax = ortho->r_to_rb( ZMAX); break; - case SL: + ximin = orthoSL->z_to_xi(-ZMAX); + ximax = orthoSL->z_to_xi( ZMAX); + break; + } + + double z, r, ans=0.0; + for (int i=0; ipotl(N1, i, z); + double tmp2 = ortho->dens(N1, i, z); + + ans += ortho->potl(N1, i, z)* + ortho->dens(N2, i, z) * 2.0*ZMAX*lw.weight(i); + } - ans += orthoSL->get_pot(x, IKX, IKY, N1)* - orthoSL->get_dens(x, IKX, IKY, N2) / - orthoSL->d_xi_to_z(x) * (ximax - ximin)*lw.weight(i); + break; - break; + case SL: + + ans += orthoSL->get_pot(z, IKX, IKY, N1)* + orthoSL->get_dens(z, IKX, IKY, N2) * 2.0*ZMAX*lw.weight(i); + + break; + } } + std::cout << "<" << N1 << "|" << N2 << "> = " << ans << std::endl; + } + break; + } + case 4: + { + switch (Type) { + case Trig: + std::cout << "No internal orthogonality check for OneDTrig" + << std::endl; + break; + case SL: + { + int knots; + std::cout << "Number of knots? "; + std::cin >> knots; + auto orthoMat = orthoSL->orthoCheck(knots); + std::cout << "Orthogonality check for SLGridSlab" << std::endl; + for (size_t i=0; i = " << ans << endl; } - - break; - default: done = true; break; @@ -412,7 +381,7 @@ main(int argc, char** argv) if (done) break; } } - + if (use_mpi) MPI_Finalize(); return 0; diff --git a/utils/SL/slabchk.cc b/utils/SL/slabchk.cc index 3baf1589f..f1bfbecbf 100644 --- a/utils/SL/slabchk.cc +++ b/utils/SL/slabchk.cc @@ -9,35 +9,44 @@ #include "gaussQ.H" #include "SLGridMP2.H" #include "Model1d.H" +#include "numerical.H" int main(int argc, char** argv) { - double H, zmax, zend; - int kmax, nmax, knots, numz; - std::string filename; + double H, zmax, zend, sigma, eps; + int kx, nmax, knots, numz, nfid; + std::string model, filename; // Parse command line // - cxxopts::Options options(argv[0], "Check the consistency a spherical SL basis"); + cxxopts::Options options(argv[0], "Check the consistency a slab SL basis"); options.add_options() ("h,help", "Print this help message") - ("ortho", "Compute orthogonality matrix") + ("o,ortho", "Compute orthogonality matrix") ("H,height", "Slab scale height", cxxopts::value(H)->default_value("1.0")) - ("K,kmax", "maximum order of in-plane harmonics", - cxxopts::value(kmax)->default_value("4")) + ("K,kx", "order of in-plane harmonics", + cxxopts::value(kx)->default_value("0")) ("N,nmax", "maximum number of vertical harmonics", cxxopts::value(nmax)->default_value("18")) ("n,numz", "size of vertical grid", cxxopts::value(numz)->default_value("1000")) ("Z,zmax", "maximum extent of vertical grid", cxxopts::value(zmax)->default_value("10.0")) - ("zend", "potential offset", + ("s,sigma","Velocity dispersion for model", + cxxopts::value(sigma)->default_value("1.0")) + ("d,zend", "potential offset", cxxopts::value(zend)->default_value("0.0")) - ("knots", "Number of Legendre integration knots", + ("e,eps", "edge offset for basis and sampling", + cxxopts::value(eps)->default_value("1.0e-3")) + ("k,knots", "Number of Legendre integration knots", cxxopts::value(knots)->default_value("40")) + ("j,ncheck", "Order of coefficient to check for orthogonality", + cxxopts::value(nfid)->default_value("2")) + ("m,model", "Model type (Sech2mu, Uniform)", + cxxopts::value(model)->default_value("Sech2mu")) ("p,prefix", "Output filename prefix", cxxopts::value(filename)->default_value("slabchk_test")) ; @@ -63,17 +72,42 @@ main(int argc, char** argv) return 1; } - // Generate Sech2 disk grid + // Define model type and instantiate model // - double dispz = 2.0*M_PI*H*H; - - Sech2 sech2(dispz); - double h = sech2.get_scale_height(); + OneDModel::ModelType mtype = OneDModel::parse_model_type(model); + std::shared_ptr mdl; + std::string mname = "isothermal"; + + if (mtype == OneDModel::ModelType::uniform) { + + std::cout << "Harmonic potential: setting scale height to " + << H << " and Zmax to " << H*(1.0 - eps) << std::endl; + mdl = std::make_shared(H, sigma, sigma); + zmax = H*(1.0 - eps); + mname = "uniform"; + } + else if (mtype == OneDModel::ModelType::cosine) { + + std::cout << "Cosine bell potential: setting scale height to " + << H << " and Zmax to " << H*0.999 << std::endl; + mdl = std::make_shared(H, sigma, sigma); + zmax = H*0.999; + mname = "cosine"; + } else if (mtype == OneDModel::ModelType::sech2mu) { + + std::cout << "Sech2 potential: using scale height " << H << std::endl; - SLGridSlab::H = h; - SLGridSlab::ZEND = zend; - std::cout << "Check...scale height is: " << SLGridSlab::H - << std::endl << std::endl; + auto tmp = std::make_shared(sigma*sigma, H); + tmp->set_hmax(10.0*H); + mdl = tmp; + + SLGridSlab::H = H; + SLGridSlab::ZEND = zend; + mname = "sech2mu"; + } else { + std::cerr << "model must be 'Sech2mu' or 'Harmonic'\n"; + return 1; + } // Particle position generator // @@ -81,15 +115,44 @@ main(int argc, char** argv) std::mt19937 gen(rd()); std::uniform_real_distribution<> dis(0.0, 1.0); - auto sample = [&dis, &gen, &h]() + auto sample = [&dis, &gen, &H, &eps, &mtype]() { - double m = dis(gen); - return 0.5*h*log(m/(1.0 - m)); + if (mtype == OneDModel::ModelType::uniform) { + return H*(2.0*dis(gen) - 1.0); + } + else if (mtype == OneDModel::ModelType::cosine) { + // The cumulative mass profile for the cosine bell model is: + // + auto M = [&H](double z) -> double { + return 0.5/H*(z + H + H/M_PI*sin(M_PI*z/H)); + }; + + // Generate a random mass and invert M(z) using Brent's method to get z + // + double del = 1.0 - eps; + double Mbot = M(-H*del), Mtop = M(H*del); + double m = Mbot + (Mtop - Mbot)*dis(gen); + + auto R = [&m, &H](double z) -> double { + return 0.5/H*(z + H + H/M_PI*sin(M_PI*z/H)) - m; + }; + + return zbrent(R, -H*del, H*del, 1e-10); + } + else if (mtype == OneDModel::ModelType::sech2mu) { + double m = dis(gen); + return 0.5*H*log(m/(1.0 - m)); + } + else { + throw std::runtime_error("Invalid model type for sampling"); + } }; // Generate Sturm-Liouville grid // - auto ortho = std::make_shared(kmax, nmax, numz, zmax, "isothermal"); + auto ortho = std::make_shared(kx+1, nmax, numz, zmax, + ".slgrid_slab_cache", + "linear", mname, false); LegeQuad lw(knots); @@ -105,12 +168,12 @@ main(int argc, char** argv) for (int n=0; nget_pot(x, 0, 0, n, 0)* - ortho->get_dens(x, 0, 0, 1, 0) / + ans1[n] += -ortho->get_pot(x, kx, 0, n, 0)* + ortho->get_dens(x, 0, 0, nfid, 0) / ortho->d_xi_to_z(x) * (ximax - ximin)*lw.weight(i); - ans2[n] += -ortho->get_pot(x, 0, 0, n, 0)* - 4.0*M_PI*sech2.get_density(z) / + ans2[n] += -ortho->get_pot(x, kx, 0, n, 0)* + 4.0*M_PI*mdl->get_density(z) / ortho->d_xi_to_z(x) * (ximax - ximin)*lw.weight(i); } } @@ -118,7 +181,7 @@ main(int argc, char** argv) // Monte Carlo version // int Number = 100000; - double fac = 4.0*M_PI*2.0*h/Number; + double fac = 4.0*M_PI/Number; for (int i=0; iz_to_xi(z); for (int n=0; nget_pot(x, 0, 0, n, 0) * fac; + ans3[n] += -ortho->get_pot(x, kx, 0, n, 0) * fac; } } + std::cout << "---- Coefficients\n"; + std::cout << std::setw(6) << "n" + << std::setw(18) << "Dens(1)" + << std::setw(18) << "Dens(model)" + << std::setw(18) << "Dens(MC)" + << std::setw(18) << "Model dens at z=0" + << std::endl; + for (int n=0; nget_dens(0.0, 0, 0, n)/(4.0*M_PI) + << std::setw(18) << ans2[n]*ortho->get_dens(0.0, kx, 0, n)/(4.0*M_PI) << std::endl; } std::cout << std::endl; int NUM = 20; - double dz = 3.0*H/(NUM-1); + double dz = zmax/(NUM-1); + + std::cout << "---- Comparing SL solution to model . . .\n"; + std::cout << std::setw(16) << "z" + << std::setw(16) << "rho (SL)" + << std::setw(16) << "rho (mod)" + << std::setw(16) << "phi (SL)" + << std::setw(16) << "phi (mod)" + << std::setw(16) << "dphi (SL)" + << std::setw(16) << "dphi (mod)" + << std::endl; for (int i=0; iget_dens (x, 0, 0, n, 0); - t += ans2[n]*ortho->get_pot (x, 0, 0, n, 0); - v += ans2[n]*ortho->get_force(x, 0, 0, n, 0); + s += ans2[n]*ortho->get_dens (x, kx, 0, n, 0); + t += ans2[n]*ortho->get_pot (x, kx, 0, n, 0); + v += ans2[n]*ortho->get_force(x, kx, 0, n, 0); } std::cout << std::setw(16) << z << std::setw(16) << s/(4.0*M_PI) - << std::setw(16) << sech2.get_density(z) + << std::setw(16) << mdl->get_density(z) << std::setw(16) << t - << std::setw(16) << sech2.get_pot(z) + << std::setw(16) << mdl->get_pot(z) << std::setw(16) << v - << std::setw(16) << sech2.get_dpot(z) + << std::setw(16) << mdl->get_dpot(z) << std::endl; } @@ -170,20 +251,20 @@ main(int argc, char** argv) if (out) { int NUM = 200; - double zmin = -3.0*H, zmax = 3.0*H; - double dz = (zmax - zmin)/(NUM-1); + double Zmin = -zmax, Zmax = zmax; + double dz = (Zmax - Zmin)/(NUM-1); for (int i=0; iz_to_xi(z); double s = 0.0; out << std::setw(16) << z; for (int n=0; nget_pot (z, 0, 0, n); - out << std::setw(16) << ortho->get_dens (z, 0, 0, n); - out << std::setw(16) << ortho->get_force(z, 0, 0, n); + out << std::setw(16) << ortho->get_pot (z, kx, 0, n); + out << std::setw(16) << ortho->get_dens (z, kx, 0, n); + out << std::setw(16) << ortho->get_force(z, kx, 0, n); } out << std::endl; } @@ -192,19 +273,25 @@ main(int argc, char** argv) } out.close(); - out.open(filename + ".ortho"); - if (out) { - auto test = ortho->orthoCheck(); - int cnt = 0; - for ( auto & v : test) { - out << "==== " << cnt++ << std::endl << v << std::endl; - } - - } else { - throw std::runtime_error("Error opening filename <" + filename + ".ortho>"); - } + if (vm.count("ortho")) { + std::cout << "---- Computing orthogonality matrix\n"; + + out.open(filename + ".ortho"); + if (out) { + auto test = ortho->orthoCheck(); + int cnt = 0; + for ( auto & v : test) { + out << "==== " << cnt++ << std::endl << v << std::endl; + } + + std::cout << "---- Orthogonality matrix written to <" << filename << ".ortho>\n"; + } else { + throw std::runtime_error("Error opening filename <" + filename + ".ortho>"); + } + } + return 0; } diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index cc735fe08..56dc4959f 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,6 +1,6 @@ set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj - testEmpDeproj testEmp testED testmd5) + testEmpDeproj testEmp testED testmd5 vtest) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil yaml-cpp ${VTK_LIBRARIES}) @@ -40,6 +40,7 @@ add_executable(testEmpDeproj testEmpDeproj.cc CubicSpline.cc add_executable(testEmp testEmp.cc EmpDeproj.cc) add_executable(testmd5 testmd5.cc) add_executable(testED testED.cc EmpDeproj.cc) +add_executable(vtest version_test.cc) foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) diff --git a/utils/Test/version_test.cc b/utils/Test/version_test.cc new file mode 100644 index 000000000..16f2b585b --- /dev/null +++ b/utils/Test/version_test.cc @@ -0,0 +1,15 @@ +// Example of using the version string and parsed integer triplet from libvars.H + +#include +#include "libvars.H" +using namespace __EXP__; + +int main() +{ + std::cout << "Version string: " << VERSION << '\n'; + std::cout << "Parsed into integer triplet:\n"; + std::cout << "-- Major=" << exp_build.major << '\n'; + std::cout << "-- Minor=" << exp_build.minor << '\n'; + std::cout << "-- Patch=" << exp_build.patch << '\n'; + return 0; +}