diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d22f3e97..cae3c8d66 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,12 @@ set(CMAKE_CXX_EXTENSIONS OFF) # Compiler flags. Not all tested thoroughly... if(CMAKE_CXX_COMPILER_ID MATCHES "Clang") # using Clang + message(STATUS "Clang compiler version is: ${CMAKE_CXX_COMPILER_VERSION}") + if(CMAKE_SYSTEM_NAME MATCHES "Linux") + if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 18.0.0) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libstdc++") + endif() + endif() elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU") # using GCC elseif(CMAKE_CXX_COMPILER_ID MATCHES "Intel") @@ -25,6 +31,9 @@ endif() # Required compiler features add_compile_options(-D_REENTRANT) +# Bake in library paths (esp. useful for HPC sites with modules) +set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) + # Check and enforce that we are a git repository. Necessary for # submodules to work correctly. if(EXISTS "${PROJECT_SOURCE_DIR}/.git") @@ -67,6 +76,7 @@ find_package(HDF5 COMPONENTS C CXX HL REQUIRED) find_package(TIRPC) # Check for alternative Sun rpc support find_package(Eigen3 REQUIRED) find_package(PNG) +find_package(Git) # Check for FE include(FEENABLE) @@ -183,7 +193,7 @@ add_compile_definitions(DATTRIB_CUDA=${CUDA_EXP_DATTRIB}) # Get the current working branch execute_process( - COMMAND git rev-parse --abbrev-ref HEAD + COMMAND ${GIT_EXECUTABLE} rev-parse --abbrev-ref HEAD WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} OUTPUT_VARIABLE GIT_BRANCH OUTPUT_STRIP_TRAILING_WHITESPACE @@ -191,7 +201,7 @@ execute_process( # Git submodule updates execute_process( - COMMAND git submodule update --init --recursive + COMMAND ${GIT_EXECUTABLE} submodule update --init --recursive WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} RESULT_VARIABLE GIT_SUBMOD_RESULT ) @@ -204,7 +214,7 @@ endif() # Get the latest abbreviated commit hash of the working branch execute_process( - COMMAND git rev-parse HEAD + COMMAND ${GIT_EXECUTABLE} rev-parse HEAD WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} OUTPUT_VARIABLE GIT_COMMIT OUTPUT_STRIP_TRAILING_WHITESPACE diff --git a/INSTALL.md b/INSTALL.md index 449daac4c..6bcb7e53c 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -1,16 +1,25 @@ # Configuring and building EXP -We are now using git submodules to provide `yaml-cpp`, which is not -common in the HPC environments. So, from the top-level directory, do -the following: +We are now using git submodules to provide a number of packages that +are not common in the HPC environments. These include + +- `HighFive`, a C++ wrapper +class to HDF5 +- `yaml-cpp`, a C++ class for reading and emitting YAML code +- `pybind11` which is used for Python bindings to the C++ classes +- `rapidxml`, used to write VTK files for rendering outside of EXP +- `png++`, C++ wrappers to the png library [Note: png support is optional] + +CMake will automatically download and configure these packages on the +first call. However, if you would to do this manually, from the +top-level directory, execute the following command: ``` git submodule update --init --recursive ``` -This will install `yaml-cpp` in the `extern` directory. The png++ C++ -wrappers to the png library are also installed in `extern`. Note: png -support is optional. +This will install the source packages in the `extern` directory. + ## EXP uses CMake for building a configuration diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 089887934..965c051c7 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -123,8 +124,22 @@ namespace BasisClasses //! Midplane escursion parameter double colh = 4.0; - public: + //@ + //! Pseudo-acceleration db + Eigen::VectorXd t_accel; + Eigen::MatrixXd p_accel; + //@} + //! Number of center points in acceleration estimator + int Naccel = 0; + + //! Get the current pseudo acceleration value + Eigen::Vector3d currentAccel(double time); + + public: + //! The current pseudo acceleration + Eigen::Vector3d pseudo {0, 0, 0}; + //! Constructor from YAML node Basis(const YAML::Node& conf, const std::string& name="Basis"); @@ -231,6 +246,19 @@ namespace BasisClasses //! Height above/below the plane for midplane search in disk scale //! lengths void setColumnHeight(double value) { colh = value; } + + //@{ + //! Initialize non-inertial forces + void setNonInertial(int N, Eigen::VectorXd& x, Eigen::MatrixXd& pos); + void setNonInertial(int N, const std::string& orient); + //@} + + //! Set the current pseudo acceleration + void setNonInertialAccel(double time) + { + if (Naccel > 0) pseudo = currentAccel(time); + } + }; using BasisPtr = std::shared_ptr; diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 5b574a8e2..483dcca46 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -184,6 +184,9 @@ namespace BasisClasses else if ( !name.compare("flatdisk") ) { basis = std::make_shared(conf); } + else if ( !name.compare("CBDisk") ) { + basis = std::make_shared(conf); + } else if ( !name.compare("slabSL") ) { basis = std::make_shared(conf); } @@ -280,5 +283,98 @@ namespace BasisClasses return makeFromArray(time); } + void Basis::setNonInertial(int N, Eigen::VectorXd& t, Eigen::MatrixXd& pos) + { + Naccel = N; + t_accel = t; + p_accel = pos; + } + + void Basis::setNonInertial(int N, const std::string& orient) + { + std::ifstream in(orient); + + if (not in) { + throw std::runtime_error("Cannot open Orient file with centering data: " + orient); + } + + const int cbufsiz = 16384; + std::unique_ptr cbuf(new char [cbufsiz]); + + // Look for data and write it while + // accumlating data for averaging + Eigen::Vector3d testread; + double time, dummy; + + std::vector times; + std::vector centers; + + while (in) { + + in.getline(cbuf.get(), cbufsiz); + if (in.rdstate() & (ios::failbit | ios::eofbit)) break; + + // Skip comment lines + // + if (cbuf[0] == '#') continue; + + std::istringstream line(cbuf.get()); + + // Read until current time is reached + line >> time; // + line >> dummy; + line >> dummy; + + bool allRead = true; + for (int i=0; i<8; i++) { + if (line.eof()) allRead = false; + for (int k; k<3; k++) line >> testread(k); + } + if (allRead) { + times.push_back(time); + centers.push_back(testread); + } + } + + // Repack data + Naccel = N; + t_accel.resize(times.size()); + p_accel.resize(times.size(), 3); + for (int i=0; i t_accel(t_accel.size()-1)) { + std::cout << "ERROR: " << time << " is outside of range of the non-inertial DB" + << std::endl; + ret.setZero(); + return ret; + } + else { + int imin = 0; + int imax = std::lower_bound(t_accel.data(), t_accel.data()+t_accel.size(), time) - t_accel.data(); + if (imax > Naccel) imin = imax - Naccel; + + int num = imax - imin + 1; + Eigen::VectorXd t(num); + Eigen::MatrixXd p(num, 3); + for (int i=imin; i<=imax; i++) { + t(i-imin) = t_accel(i); + for (int k=0; k<3; k++) p(i-imin, k) = p_accel(i, k); + } + + for (int k=0; k<3; k++) + ret(k) = 2.0*std::get<0>(QuadLS(t, p.col(k)).coefs()); + } + return ret; + } + } // END namespace BasisClasses diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index a4af27c1a..b5aa9ac89 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -529,7 +529,7 @@ namespace BasisClasses //! Compute the orthogonality of the basis by returning inner //! produce matrices - std::vector orthoCheck(); + std::vector orthoCheck(int knots=400); //! Biorthogonality sanity check bool orthoTest() @@ -542,6 +542,157 @@ namespace BasisClasses }; + /** + Uses the Clutton-Brock basis to evaluate expansion coeffients and + provide potential and density basis fields + */ + class CBDisk : public BiorthBasis + { + + public: + + using BasisMap = std::map; + using BasisArray = std::vector>; + + private: + + //! Helper for constructor + void initialize(); + + int mmax, nmax; + double scale; + + bool NO_M0, NO_M1, EVEN_M, M0_only; + + std::vector potd, potR, potZ, dend; + + Eigen::MatrixXd expcoef; + int N1, N2; + int used; + + using matT = std::vector; + using vecT = std::vector; + + double totalMass; + int npart; + + Eigen::VectorXd work; + + //! For coefficient writing + typedef Eigen::Matrix + EigenColMajor; + + //! Get potential + void get_pot(Eigen::MatrixXd& tab, double r); + + //! Get density + void get_dens(Eigen::MatrixXd& tab, double r); + + //! Get force + void get_force(Eigen::MatrixXd& tab, double r); + + //@{ + //! 2D Clutton-Brock basis + double phif(const int n, const int m, const double r); + double dphi(const int n, const int m, const double r); + + double potl(const int n, const int m, const double r); + double dpot(const int n, const int m, const double r); + double dens(const int n, const int m, const double r); + + double norm(const int n, const int m); + + void potl(const int m, const double r, Eigen::VectorXd& a); + void dpot(const int m, const double r, Eigen::VectorXd& a); + void dens(const int m, const double r, Eigen::VectorXd& a); + //@} + + //! Potential, force, and density scaling + double fac1, fac2; + + protected: + + //! Evaluate basis in cylindrical coordinates + virtual std::vector + cyl_eval(double R, double z, double phi); + + //! Evaluate basis in spherical coordinates. Conversion from the + //! cylindrical evaluation above. + virtual std::vector + sph_eval(double r, double costh, double phi); + + // Cartesian + virtual std::vector + crt_eval(double x, double y, double z); + + //! Load coefficients into the new CoefStruct + virtual void load_coefs(CoefClasses::CoefStrPtr coefs, double time); + + //! Set coefficients + virtual void set_coefs(CoefClasses::CoefStrPtr coefs); + + //! Valid keys for YAML configurations + static const std::set valid_keys; + + //! Return readable class name + virtual const std::string classname() { return "CBDisk";} + + //! Subspace index + virtual const std::string harmonic() { return "m";} + + public: + + //! Constructor from YAML node + CBDisk(const YAML::Node& conf); + + //! Constructor from YAML string + CBDisk(const std::string& confstr); + + //! Destructor + virtual ~CBDisk(void) {} + + //! Print and return the cache parameters + static std::map + cacheInfo(const std::string& cachefile) + { + return BiorthCyl::cacheInfo(cachefile); + } + + //! Zero out coefficients to prepare for a new expansion + void reset_coefs(void); + + //! Make coefficients after accumulation + void make_coefs(void); + + //! Accumulate new coefficients + virtual void accumulate(double x, double y, double z, double mass); + + //! Return current maximum harmonic order in expansion + int getMmax() { return mmax; } + + //! Return current maximum order in radial expansion + int getNmax() { return nmax; } + + //! Return a vector of a vector of 1d basis-function grids for + //! CBDisk, logarithmically spaced between [logxmin, logxmax] + //! (base 10). + BasisArray getBasis(double logxmin=-3.0, double logxmax=0.5, int numgrid=2000); + + //! Compute the orthogonality of the basis by returning inner + //! produce matrices + std::vector orthoCheck(int knots=4000); + + //! Biorthogonality sanity check + bool orthoTest() + { + auto [ret, worst, lworst] = orthoCompute(orthoCheck()); + // For the CTest log + std::cout << "CBDisk::orthoTest: worst=" << worst << std::endl; + return ret; + } + + }; + /** Uses EmpCylSL basis to evaluate expansion coeffients and provide potential and density basis fields @@ -676,7 +827,7 @@ namespace BasisClasses //! Compute the orthogonality of the basis by returning inner //! produce matrices - std::vector orthoCheck() + std::vector orthoCheck(int knots=40) { return sl->orthoCheck(); } @@ -813,7 +964,7 @@ namespace BasisClasses //! Compute the orthogonality of the basis by returning inner //! produce matrices - std::vector orthoCheck(); + std::vector orthoCheck(int knots=40); //! Biorthogonality sanity check bool orthoTest() @@ -937,12 +1088,12 @@ namespace BasisClasses //! Compute the orthogonality of the basis by returning inner //! produce matrices - Eigen::MatrixXcd orthoCheck(); + std::vector orthoCheck(int knots=40); //! Biorthogonality sanity check bool orthoTest() { - auto c = orthoCheck(); + auto c = orthoCheck()[0]; double worst = 0.0; for (int n1=0; n1 FlatDisk::orthoCheck() + std::vector FlatDisk::orthoCheck(int knots) { return ortho->orthoCheck(); } @@ -1996,6 +1997,677 @@ namespace BasisClasses return ret; } + const std::set + CBDisk::valid_keys = { + "self_consistent", + "NO_M0", + "NO_M1", + "EVEN_M", + "M0_BACK", + "M0_ONLY", + "NO_MONO", + "background", + "playback", + "coefMaster", + "scale", + "Lmax", + "Mmax", + "nmax", + "mmax", + "mlim", + "dof", + "subsamp", + "samplesz", + "vtkfreq", + "tksmooth", + "tkcum", + "tk_type" + }; + + CBDisk::CBDisk(const YAML::Node& CONF) : + BiorthBasis(CONF, "CBDisk") + { + initialize(); + } + + CBDisk::CBDisk(const std::string& confstr) : + BiorthBasis(confstr, "CBDisk") + { + initialize(); + } + + void CBDisk::initialize() + { + + // Assign some defaults + // + mmax = 6; + nmax = 18; + scale = 1.0; + + // Check for unmatched keys + // + auto unmatched = YamlCheck(conf, valid_keys); + if (unmatched.size()) + throw YamlConfigError("Basis::Basis::CBDisk", "parameter", unmatched, __FILE__, __LINE__); + + // Assign values from YAML + // + try { + if (conf["Lmax"]) mmax = conf["Lmax"].as(); // Proxy + if (conf["Mmax"]) mmax = conf["Mmax"].as(); // Proxy + if (conf["mmax"]) mmax = conf["mmax"].as(); + if (conf["nmax"]) nmax = conf["nmax"].as(); + if (conf["scale"]) scale = conf["scale"].as(); + + N1 = 0; + N2 = std::numeric_limits::max(); + NO_M0 = NO_M1 = EVEN_M = M0_only = false; + + if (conf["N1"] ) N1 = conf["N1"].as(); + if (conf["N2"] ) N2 = conf["N2"].as(); + if (conf["NO_M0"]) NO_M0 = conf["NO_M0"].as(); + if (conf["NO_M1"]) NO_M1 = conf["NO_M1"].as(); + if (conf["EVEN_M"]) EVEN_M = conf["EVEN_M"].as(); + if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); + } + catch (YAML::Exception & error) { + if (myid==0) std::cout << "Error parsing parameter stanza for <" + << name << ">: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + + throw std::runtime_error("CBDisk: error parsing YAML"); + } + + // Set characteristic radius defaults + // + if (not conf["scale"]) conf["scale"] = 1.0; + + // Potential, force, and density scaling + // + fac1 = pow(scale, -0.5); + fac2 = pow(scale, -1.5); + + // Orthogonality sanity check + // + orthoTest(); + + // Get max threads + // + int nthrds = omp_get_max_threads(); + + // Allocate memory + // + potd.resize(nthrds); + potR.resize(nthrds); + dend.resize(nthrds); + + for (auto & v : potd) v.resize(mmax+1, nmax); + for (auto & v : potR) v.resize(mmax+1, nmax); + for (auto & v : dend) v.resize(mmax+1, nmax); + + expcoef.resize(2*mmax+1, nmax); + expcoef.setZero(); + + work.resize(nmax); + + used = 0; + + // Set cylindrical coordindates + // + coordinates = Coord::Cylindrical; + } + + // Get potential + void CBDisk::get_pot(Eigen::MatrixXd& tab, double r) + { + tab.resize(mmax+1, nmax); + Eigen::VectorXd a(nmax); + for (int m=0; m<=mmax; m++) { + potl(m, r/scale, a); + tab.row(m) = a; + } + tab *= fac1; + } + + // Get density + void CBDisk::get_dens(Eigen::MatrixXd& tab, double r) + { + tab.resize(mmax+1, nmax); + Eigen::VectorXd a(nmax); + for (int m=0; m<=mmax; m++) { + dens(m, r/scale, a); + tab.row(m) = a; + } + tab *= fac2; + } + + // Get force + void CBDisk::get_force(Eigen::MatrixXd& tab, double r) + { + tab.resize(mmax+1, nmax); + Eigen::VectorXd a(nmax); + for (int m=0; m<=mmax; m++) { + dpot(m, r/scale, a); + tab.row(m) = a; + } + tab *= fac1/scale; + } + + // Routines for computing biorthonormal pairs based on + // Clutton-Brock's 2-dimensional series + // + double CBDisk::phif(const int n, const int m, const double r) + { + // By recurrance relation + // + double r2 = r*r; + double fac = 1.0/(1.0 + r2); + double cur = sqrt(fac); + + for (int mm=1; mm<=m; mm++) cur *= fac*(2*mm - 1); + + if (n==0) return cur; + + double curl1 = cur; + double curl2; + + fac *= r2 - 1.0; + cur *= fac*(2*m+1); + + if (n==1) return cur; + + for (int nn=2; nn<=n; nn++) { + curl2 = curl1; + curl1 = cur; + cur = (2.0 + (double)(2*m-1)/nn)*fac*curl1 - + (1.0 + (double)(2*m-1)/nn)*curl2; + } + + return cur; + } + + double CBDisk::potl(const int n, const int m, const double r) + { + return pow(r, m) * phif(n, m, r)/sqrt(norm(n, m)); + } + + double CBDisk::dpot(const int n, const int m, const double r) + { + double ret = dphi(n, m, r); + if (m) ret = (phif(n, m, r)*m/r + ret) * pow(r, m); + return ret/sqrt(norm(n, m)); + } + + double CBDisk::dphi(const int n, const int m, const double r) + { + double ret = phif(n, m+1, r); + if (n>0) ret -= 2.0*phif(n-1, m+1, r); + if (n>1) ret += phif(n-2, m+1, r); + return -r*ret; + } + + + // By recurrance relation + // + void CBDisk::potl(const int m, const double r, Eigen::VectorXd& a) + { + a.resize(nmax); + + double pfac = pow(r, m); + + double r2 = r*r; + double fac = 1.0/(1.0 + r2); + double cur = sqrt(fac); + + for (int mm=1; mm<=m; mm++) cur *= fac*(2*mm - 1); + + a(0) = pfac*cur; + + if (nmax>0) { + + double curl1 = cur; + double curl2; + + fac *= r2 - 1.0; + cur *= fac*(2*m+1); + + a(1) = pfac*cur; + + for (int nn=2; nn(m)/r; + + Eigen::VectorXd b(nmax); + potl(m+1, r, b); + + for (int n=0; n0) a(n) += 2.0*b(n-1); + if (n>1) a(n) -= b(n-2); + } + } + + double CBDisk::dens(int n, int m, double r) + { + double f = 0.5/sqrt(norm(n, m))/M_PI; + + if (n>=2) + return f*pow(r, m)*(phif(n, m+1, r) - phif(n-2, m+1, r)); + else + return f*pow(r, m)*phif(n, m+1, r); + } + + + void CBDisk::dens(int mm, double r, Eigen::VectorXd& a) + { + a.resize(nmax); + + double pfac = pow(r, (double)mm+1.0e-20); + + int m = mm + 1; + double r2 = r*r; + double fac = 1.0/(1.0 + r2); + double cur = sqrt(fac); + + for (int M=1; M<=m; M++) cur *= fac*(2*M - 1); + + a(0) = pfac*cur; + + if (nmax>0) { + + double curl1 = cur; + double curl2; + + fac *= r2 - 1.0; + cur *= fac*(2*m+1); + + a(1) = pfac*cur; + + for (int nn=2; nn1; nn--) + a(nn) -= a(nn-2); + } + + for (int n=0; n CBDisk::orthoCheck(int num) + { + const double tol = 1.0e-4; + LegeQuad lq(num); + + std::vector ret(mmax+1); + for (auto & v : ret) v.resize(nmax, nmax); + + double Rmax = scale*100.0; + + for (int m=0; m<=mmax; m++) { + + ret[m].setZero(); + + for (int i=0; i(coef.get()); + auto cf = cc->coefs; + int rows = cf->rows(); + int cols = cf->cols(); + if (rows != mmax+1 or cols != nmax) { + std::ostringstream sout; + sout << "CBDisk::set_coefs: the basis has (mmax+1, nmax)=(" + << mmax+1 << ", " << nmax + << "). The coef structure has (rows, cols)=(" + << rows << ", " << cols << ")"; + + throw std::runtime_error(sout.str()); + } + } + + CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); + auto & cc = *cf->coefs; + + // Assign internal coefficient table (doubles) from the complex struct + // + for (int m=0, m0=0; m<=mmax; m++) { + for (int n=0; nctr.size()) + coefctr = cf->ctr; + else + coefctr = {0.0, 0.0, 0.0}; + } + + void CBDisk::accumulate(double x, double y, double z, double mass) + { + // Normalization factors + // + constexpr double norm0 = 1.0; + constexpr double norm1 = M_SQRT2; + + //====================== + // Compute coefficients + //====================== + + double R2 = x*x + y*y; + double R = sqrt(R); + + // Get thread id + int tid = omp_get_thread_num(); + + used++; + totalMass += mass; + + double phi = atan2(y, x); + + get_pot(potd[tid], R); + + // M loop + for (int m=0, moffset=0; m<=mmax; m++) { + + if (m==0) { + for (int n=0; nCBDisk::cyl_eval(double R, double z, double phi) + { + // Get thread id + int tid = omp_get_thread_num(); + + // Fixed values + constexpr double norm0 = 1.0; + constexpr double norm1 = M_SQRT2; + + double den0=0, den1=0, pot0=0, pot1=0, rpot=0, zpot=0, ppot=0; + + // Get the basis fields + // + get_dens (dend[tid], R); + get_pot (potd[tid], R); + get_force (potR[tid], R); + + // m loop + // + for (int m=0, moffset=0; m<=mmax; m++) { + + if (m==0 and NO_M0) { moffset++; continue; } + if (m==1 and NO_M1) { moffset += 2; continue; } + if (EVEN_M and m/2*2 != m) { moffset += 2; continue; } + if (m>0 and M0_only) break; + + if (m==0) { + for (int n=std::max(0, N1); n<=std::min(nmax-1, N2); n++) { + den0 += expcoef(0, n) * dend[tid](0, n) * norm0; + pot0 += expcoef(0, n) * potd[tid](0, n) * norm0; + rpot += expcoef(0, n) * potR[tid](0, n) * norm0; + } + + moffset++; + } else { + double cosm = cos(phi*m), sinm = sin(phi*m); + double vc, vs; + + vc = vs = 0.0; + for (int n=std::max(0, N1); n<=std::min(nmax-1, N2); n++) { + vc += expcoef(moffset+0, n) * dend[tid](m, n); + vs += expcoef(moffset+1, n) * dend[tid](m, n); + } + + den1 += (vc*cosm + vs*sinm) * norm1; + + vc = vs = 0.0; + for (int n=std::max(0, N1); n<=std::min(nmax-1, N2); n++) { + vc += expcoef(moffset+0, n) * potd[tid](m, n); + vs += expcoef(moffset+1, n) * potd[tid](m, n); + } + + pot1 += ( vc*cosm + vs*sinm) * norm1; + ppot += (-vc*sinm + vs*cosm) * m * norm1; + + vc = vs = 0.0; + for (int n=std::max(0, N1); n<=std::min(nmax-1, N2); n++) { + vc += expcoef(moffset+0, n) * potR[tid](m, n); + vs += expcoef(moffset+1, n) * potR[tid](m, n); + } + + rpot += (vc*cosm + vs*sinm) * norm1; + + moffset +=2; + } + } + + den0 *= -1.0; + den1 *= -1.0; + pot0 *= -1.0; + pot1 *= -1.0; + rpot *= -1.0; + ppot *= -1.0; + + return {den0, den1, den0+den1, pot0, pot1, pot0+pot1, rpot, zpot, ppot}; + } + + + std::vector CBDisk::sph_eval(double r, double costh, double phi) + { + // Cylindrical coords + // + double sinth = sqrt(fabs(1.0 - costh*costh)); + double R = r*sinth, z = r*costh; + + auto v = cyl_eval(R, z, phi); + + // Spherical force element converstion + // + double potr = v[6]*sinth + v[7]*costh; + double pott = v[6]*costh - v[7]*sinth; + + return {v[0], v[1], v[2], v[3], v[4], v[5], potr, pott, v[8]}; + } + + std::vector CBDisk::crt_eval(double x, double y, double z) + { + // Cylindrical coords from Cartesian + // + double R = sqrt(x*x + y*y) + 1.0e-18; + double phi = atan2(y, x); + + auto v = cyl_eval(R, z, phi); + + double potx = v[4]*x/R - v[8]*y/R; + double poty = v[4]*y/R + v[8]*x/R; + + return {v[0], v[1], v[2], v[3], v[4], v[5], potx, poty, v[7]}; + } + + CBDisk::BasisArray CBDisk::getBasis + (double logxmin, double logxmax, int numgrid) + { + // Assing return storage + BasisArray ret(mmax+1); + for (auto & v : ret) { + v.resize(nmax); + for (auto & u : v) { + u["potential"].resize(numgrid); // Potential + u["density" ].resize(numgrid); // Density + u["rforce" ].resize(numgrid); // Radial force + } + } + + // Radial grid spacing + double dx = (logxmax - logxmin)/(numgrid-1); + + // Basis storage + Eigen::MatrixXd tabpot, tabden, tabrfc; + + // Evaluate on the plane + for (int i=0; i Slab::valid_keys = { "nmaxx", @@ -2442,7 +3114,7 @@ namespace BasisClasses return ret; } - std::vector Slab::orthoCheck() + std::vector Slab::orthoCheck(int knots) { return ortho->orthoCheck(); } @@ -2747,9 +3419,11 @@ namespace BasisClasses return {0, den1, den1, 0, pot1, pot1, potr, pott, potp}; } - Eigen::MatrixXcd Cube::orthoCheck() + std::vector Cube::orthoCheck(int knots) { - return ortho->orthoCheck(); + std::vector ret; + ret.push_back(ortho->orthoCheck().array().abs()); + return ret; } // Generate coeffients from a particle reader @@ -2979,7 +3653,7 @@ namespace BasisClasses for (int n=0; ngetFields(ps(n, 0), ps(n, 1), ps(n, 2)); // First 6 fields are density and potential, follewed by acceleration - for (int k=0; k<3; k++) accel(n, k) += v[6+k]; + for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k); } return accel; @@ -3048,6 +3722,10 @@ namespace BasisClasses // Install coefficients // basis->set_coefs(newcoef); + + // Set non-inertial force + basis->setNonInertialAccel(t); + } SingleTimeAccel::SingleTimeAccel(double t, std::vector mod) diff --git a/expui/BiorthBess.H b/expui/BiorthBess.H index ed0dcb57e..05a2cef41 100644 --- a/expui/BiorthBess.H +++ b/expui/BiorthBess.H @@ -1,6 +1,7 @@ #ifndef _BiorthBess_H_ #define _BiorthBess_H_ +#include #include #include #include diff --git a/expui/BiorthBess.cc b/expui/BiorthBess.cc index c73a2b7ec..0b16c8ce0 100644 --- a/expui/BiorthBess.cc +++ b/expui/BiorthBess.cc @@ -1,3 +1,5 @@ +#include + #include #include #include diff --git a/expui/Centering.cc b/expui/Centering.cc index 68515d7d2..8b00b1f07 100644 --- a/expui/Centering.cc +++ b/expui/Centering.cc @@ -21,8 +21,8 @@ namespace Utility std::vector ctr(3, 0.0); - typedef point point3; - typedef kdtree tree3; + using point3 = KDtree::point; + using tree3 = KDtree::kdtree; std::vector points; diff --git a/expui/CoefContainer.H b/expui/CoefContainer.H index 0c8bf39d0..fe10b64f3 100644 --- a/expui/CoefContainer.H +++ b/expui/CoefContainer.H @@ -129,6 +129,13 @@ namespace MSSA void pack_table(); void unpack_table(); //@} + + //@{ + //! Trajectories + void pack_traj(); + void unpack_traj(); + //@} + //@} //@{ @@ -168,6 +175,8 @@ namespace MSSA //! Copy data from backround back to the coefficient database void background(); + //! Access to coefficients + CoefClasses::CoefsPtr getCoefs() { return coefs; } }; /** @@ -287,7 +296,7 @@ components: std::cout << " ("; for (auto l : key) std::cout << l << "|"; std::cout << ")" << std::endl; - throw std::runtime_error("CoefDB::getData: desired key is not data"); + throw std::runtime_error("CoefContainer::getData: desired key is not data"); } return comps[c]->data[key]; } @@ -322,7 +331,7 @@ components: auto it = namemap.find(name); if (it == namemap.end()) { throw std::runtime_error - ("CoefDB::isComplex: cannot find component with given name"); + ("CoefContainer::isComplex: cannot find component with given name"); } return comps[it->second]->complexKey; @@ -344,6 +353,17 @@ components: return ret; } + //! Access to coefficients + CoefClasses::CoefsPtr getCoefs(const std::string& name) { + auto it = namemap.find(name); + if (it == namemap.end()) { + throw std::runtime_error + ("CoefContainer::getCoefs: cannot find component with given name"); + } + + return comps[it->second]->coefs; + } + }; // END CoefContainer diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index 0d48b17fe..9dda72b9d 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -107,6 +107,8 @@ namespace MSSA pack_cube(); else if (dynamic_cast(coefs.get())) pack_table(); + else if (dynamic_cast(coefs.get())) + pack_traj(); else if (dynamic_cast(coefs.get())) pack_sphfld(); else if (dynamic_cast(coefs.get())) @@ -788,6 +790,89 @@ namespace MSSA } + void CoefDB::pack_traj() + { + auto cur = dynamic_cast(coefs.get()); + + times = cur->Times(); + complexKey = false; + + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); + + int traj = cf->traj; + int rank = cf->rank; + int ntimes = times.size(); + + // Promote desired keys which are the data columns + // + keys.clear(); + for (auto v : keys0) { + if (v.size() != 1) { + std::ostringstream sout; + sout << "CoefDB::pack_traj: key vector should have rank 1; " + << "found rank " << v.size() << " instead"; + throw std::runtime_error(sout.str()); + } + else if (v[0] >= traj) { + std::ostringstream sout; + sout << "CoefDB::pack_traj: requested key=" << v[0] + << " exceeded the number of trajectories, " << traj; + throw std::runtime_error(sout.str()); + } + else keys.push_back(v); // OKAY + } + + // No bkeys for a trajectory + // + bkeys.clear(); + + for (unsigned m=0; m + ( cur->getCoefStruct(times[t]).get() ); + + for (int n=0; ncoefs)(m, n).real(); + } + + } + } + } + + void CoefDB::unpack_traj() + { + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[0]).get() ); + + int traj = cf->traj; + int rank = cf->rank; + int ntimes = times.size(); + + for (int t=0; t + ( coefs->getCoefStruct(times[t]).get() ); + + for (unsigned m=0; mcoefs)(m, n) = data[key][t*rank+n]; + } + // End field loop + } + // END time loop + + } + + CoefContainer::CoefContainer(const CoefContainer& p) { // Protected data diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index b50bb46f9..d438dc408 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -321,6 +321,53 @@ namespace CoefClasses }; + //! Specialization of CoefStruct for collection of trajectories + class TrajStruct : public CoefStruct + { + public: + //! Coefficient map type + using coefType = Eigen::Map; + + //! Coefficient map + std::shared_ptr coefs; + + //! Number of trajectories + int traj; + + //! Phase-space rank + int rank; + + //! Constructor + TrajStruct() : traj(0), rank(0) { geom = "trajectory"; } + + //! No need to read + bool read(std::istream& in, bool exp_type, bool verbose=false); + + //! Create an empty data storage + void create(); + + //! Copy all of data contents to a new instance + std::shared_ptr deepcopy(); + + //! Assign array + void assign(const Eigen::MatrixXd& arr) + { + traj = arr.rows(); + rank = arr.cols(); + + allocate(); + *coefs = arr; + } + + //! Allocate storage arrays + void allocate() + { + store.resize(traj*rank); + coefs = std::make_shared(store.data(), traj, rank); + } + + }; + //! Specialization of CoefStruct for spheres for general fields class SphFldStruct : public CoefStruct { @@ -441,6 +488,7 @@ namespace CoefClasses using SlabStrPtr = std::shared_ptr ; using CubeStrPtr = std::shared_ptr ; using TblStrPtr = std::shared_ptr ; + using TrajStrPtr = std::shared_ptr ; using SphFldPtr = std::shared_ptr ; using PlrFldPtr = std::shared_ptr ; diff --git a/expui/CoefStruct.cc b/expui/CoefStruct.cc index 20cb6b937..3537bde33 100644 --- a/expui/CoefStruct.cc +++ b/expui/CoefStruct.cc @@ -70,6 +70,14 @@ namespace CoefClasses throw std::runtime_error("TblStruct::create: cols must be >0"); } + void TrajStruct::create() + { + if (traj>0 and rank>0) + coefs = std::make_shared(store.data(), traj, rank); + else + throw std::runtime_error("TrajStruct::create: number of trajectories and phase-space size must be >0"); + } + void SphFldStruct::create() { if (nfld>0 and nmax>0) @@ -177,6 +185,23 @@ namespace CoefClasses return ret; } + std::shared_ptr TrajStruct::deepcopy() + { + auto ret = std::make_shared(); + + copyfields(ret); + + assert(("TrajStruct::deepcopy dimension mismatch", + traj*rank == store.size())); + + + ret->coefs = std::make_shared(ret->store.data(), traj, rank); + ret->traj = traj; + ret->rank = rank; + + return ret; + } + std::shared_ptr SphFldStruct::deepcopy() { auto ret = std::make_shared(); @@ -515,6 +540,44 @@ namespace CoefClasses } } + bool TrajStruct::read(std::istream& in, bool exp_type, bool verbose) + { + std::vector row; + std::string line; + + if (std::getline(in, line)) { + + std::istringstream sin(line); + + // Read first value as time + // + sin >> time; + + // If okay, read the rest of the line as columns + // + double val; + while (1) { + sin >> val; + if (sin) row.push_back(val); + else break; + } + + // Number of cols + int cols = row.size(); + + if (traj*rank != cols) + throw std::runtime_error("TrajStruct::read: data != traj*rank"); + + store.resize(cols); + coefs = std::make_shared(store.data(), traj, rank); + for (int i=0; i coefs; + + //! An alternate packing + std::vector data; + + //! Read the coefficients + virtual void readNativeCoefs(const std::string& file, + int stride, double tmin, double tmax); + + //! Get the YAML config + virtual std::string getYAML(); + + //! Write parameter attributes + virtual void WriteH5Params(HighFive::File& file); + + //! Write coefficient data in H5 + virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); + + public: + + //! Constructor + TrajectoryData(bool verbose=false) : Coefs("table", verbose) {} + + //! Constructor from vectors + TrajectoryData(const std::vector& time, + const std::vector& data, + bool verbose=false); + + //! Constructor from file + TrajectoryData(std::string& file, bool verbose=false); + + //! H5 constructor + TrajectoryData(HighFive::File& file, int stride=1, + double tmin=-std::numeric_limits::max(), + double tmax= std::numeric_limits::max(), + bool verbose=false); + + //! Copy constructor + TrajectoryData(TrajectoryData& p) : Coefs(p) + { + coefs = p.coefs; + data = p.data; + times = p.times; + } + + //! Clear coefficient container + virtual void clear() { coefs.clear(); } + + //! Add a coefficient structure to the container + virtual void add(CoefStrPtr coef); + + //! Get coefficient data at given time + virtual Eigen::VectorXcd& getData(double time); + + /** Set coefficient matrix at given time. This is for pybind11, + since the operator() will not allow lvalue assignment, it + seems. */ + virtual void setData(double time, Eigen::VectorXcd& arr); + + //! Get coefficient structure at a given time + virtual std::shared_ptr getCoefStruct(double time) + { return coefs[roundTime(time)]; } + + //! Get list of coefficient times + virtual std::vector Times() { return times; } + + //! Get all coefficients indexed in column, time + Eigen::Tensor getAllCoefs(); + + /** Get power for the series. Not implemented meaningfully for + trajctory data. */ + virtual Eigen::MatrixXd& Power + (int min=0, int max=std::numeric_limits::max()) + { + power.resize(0, 0); + return power; + } + + //! Make keys for the remaining indices in a subspace + virtual std::vector makeKeys(); + virtual std::vector makeKeys(Key k) { return makeKeys(); } + + //! Compare two collections of stanzas (this is for testing only) + virtual bool CompareStanzas(std::shared_ptr check); + + //! Copy all of the data; make new instances of shared pointer + //! objects. This will recursively call the deepcopy() functions + //! or all dependent objects + virtual std::shared_ptr deepcopy(); + + virtual void zerodata() { + for (auto v : coefs) v.second->zerodata(); + for (auto & v : data) v.setZero(); + } + + }; + /** Derived class for spherical field coefficients */ class SphFldCoefs : public Coefs { diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index fe55b0efe..a09362f91 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -259,6 +259,26 @@ namespace CoefClasses } + std::shared_ptr TrajectoryData::deepcopy() + { + auto ret = std::make_shared(); + + // Copy the base-class fields + copyfields(ret); + + // Copy the local structures from the map to the struct pointers + // by copyfing fields, not the pointer + for (auto v : coefs) + ret->coefs[v.first] = + std::dynamic_pointer_cast(v.second->deepcopy()); + + ret->data = data; + ret->times = times; + + return ret; + } + + SphFldCoefs::SphFldCoefs(HighFive::File& file, int stride, double Tmin, double Tmax, bool verbose) : Coefs("sphere", verbose) @@ -1879,6 +1899,271 @@ namespace CoefClasses } + TrajectoryData::TrajectoryData(const std::vector& Times, + const std::vector& data, + bool verbose) : + data(data), Coefs("table", verbose) + { + times = Times; + int traj = data.size(); + int ntim = data[0].rows(); + int rank = data[0].cols(); + + if (ntim != times.size()) { + std::ostringstream msg; + msg << "TrajectoryData ERROR (runtime) ntim [" << ntim + << "] != times.size() [" << times.size() << "]"; + throw std::runtime_error(msg.str()); + } + + for (int i=0; i(); + c->time = times[i]; + c->traj = traj; + c->rank = rank; + c->store.resize(c->traj*c->rank); + c->coefs = std::make_shared(c->store.data(), c->traj, c->rank); + for (int m=0; mcoefs)(m, n) = data[m](i, n); + } + } + coefs[roundTime(c->time)] = c; + } + } + + TrajectoryData::TrajectoryData(std::string& file, bool verbose) : + Coefs("trajectory", verbose) + { + std::ifstream in(file); + + if (not in) { + throw std::runtime_error("TrajectoryData ERROR (runtime) opening file <" + file + ">"); + } + + while (in) { + TrajStrPtr c = std::make_shared(); + if (not c->read(in, verbose)) break; + + coefs[roundTime(c->time)] = c; + } + + for (auto c : coefs) times.push_back(c.first); + } + + + TrajectoryData::TrajectoryData(HighFive::File& file, int stride, + double Tmin, double Tmax, bool verbose) : + Coefs("trajectory", verbose) + { + int traj, rank; + unsigned count; + std::string config; + + file.getAttribute("name" ).read(name); + file.getAttribute("traj" ).read(traj); + file.getAttribute("rank" ).read(rank); + file.getAttribute("config").read(config); + file.getDataSet ("count" ).read(count); + + auto snaps = file.getGroup("snapshots"); + + snaps.getDataSet("times").read(times); + snaps.getDataSet("datatable").read(data); + + for (unsigned n=0; n Tmax) continue; + + // Pack the data into the coefficient variable + // + auto coef = std::make_shared(); + coef->traj = traj; + coef->rank = rank; + coef->time = times[n]; + coef->store.resize(traj*rank); + coef->coefs = std::make_shared + (coef->store.data(), traj, rank); + *coef->coefs = data[n]; + + coefs[roundTime(times[n])] = coef; + } + } + + Eigen::VectorXcd& TrajectoryData::getData(double time) + { + auto it = coefs.find(roundTime(time)); + + if (it == coefs.end()) { + + arr.resize(0); + + } else { + + arr = it->second->store; + + } + + return arr; + } + + void TrajectoryData::setData(double time, Eigen::VectorXcd& dat) + { + auto it = coefs.find(roundTime(time)); + + if (it == coefs.end()) { + std::ostringstream str; + str << "TrajectoryData::setData: requested time=" << time << " not found"; + throw std::runtime_error(str.str()); + } else { + it->second->store = dat; + it->second->coefs = std::make_shared + (it->second->store.data(), it->second->traj, it->second->rank); + + } + } + + void TrajectoryData::readNativeCoefs(const std::string& file, int stride, + double tmin, double tmax) + { + std::ifstream in(file); + + if (not in) { + throw std::runtime_error("TrajectoryData ERROR (runtime) opening file <" + file + ">"); + } + + int count = 0; + while (in) { + TrajStrPtr c = std::make_shared(); + if (not c->read(in, verbose)) break; + + if (count++ % stride) continue; + if (c->time < tmin or c->time > tmax) continue; + + coefs[roundTime(c->time)] = c; + } + + times.clear(); + for (auto t : coefs) times.push_back(t.first); + + if (myid==0) + std::cerr << "---- Coefs::factory: " + << "read ascii and created TrajectoryData" + << std::endl; + } + + + void TrajectoryData::WriteH5Params(HighFive::File& file) + { + int traj = coefs.begin()->second->traj; + int rank = coefs.begin()->second->rank; + + std::string geometry = "trajectory"; + + file.createAttribute("traj", HighFive::DataSpace::From(traj)).write(traj); + file.createAttribute("rank", HighFive::DataSpace::From(rank)).write(rank); + } + + unsigned TrajectoryData::WriteH5Times(HighFive::Group& snaps, unsigned count) + { + snaps.createDataSet("times", times); + snaps.createDataSet("datatable", data); + + count = times.size(); + + return count; + } + + + std::string TrajectoryData::getYAML() + { + std::string ret; + if (coefs.size()) { + ret = coefs.begin()->second->buf; + } + return ret; + } + + bool TrajectoryData::CompareStanzas(std::shared_ptr check) + { + bool ret = true; + + auto other = dynamic_cast(check.get()); + + // Check that every time in this one is in the other + for (auto v : coefs) { + if (other->coefs.find(roundTime(v.first)) == other->coefs.end()) { + std::cout << "Can't find Time=" << v.first << std::endl; + ret = false; + } + } + + if (ret) { + std::cout << "Times are the same, now checking parameters at each time" + << std::endl; + for (auto v : coefs) { + auto it = other->coefs.find(v.first); + if (v.second->traj != it->second->traj) ret = false; + if (v.second->rank != it->second->rank) ret = false; + } + } + + if (ret) { + std::cout << "Parameters are the same, now checking coefficients time" + << std::endl; + for (auto v : coefs) { + auto it = other->coefs.find(v.first); + for (int m=0; mtraj; m++) { + for (int n=0; nrank; n++) { + if ((*v.second->coefs)(m, n) != (*it->second->coefs)(m, n)) { + ret = false; + } + } + } + } + } + + return ret; + } + + + Eigen::Tensor + TrajectoryData::getAllCoefs() + { + Eigen::Tensor ret; + + auto times = Times(); + + int traj = coefs.begin()->second->traj; + int rank = coefs.begin()->second->rank; + int ntim = times.size(); + + // Resize the tensor + ret.resize(traj, rank, ntim); + + for (int t=0; tcoefs)(m, n).real(); + } + } + } + + return ret; + } + + std::vector TrajectoryData::makeKeys() + { + std::vector ret; + if (coefs.size()==0) return ret; + + unsigned traj = coefs.begin()->second->traj; + for (unsigned m=0; m& Times, const std::vector>& data, bool verbose) : @@ -2408,6 +2693,14 @@ namespace CoefClasses coefs[roundTime(coef->time)] = p; } + void TrajectoryData::add(CoefStrPtr coef) + { + auto p = std::dynamic_pointer_cast(coef); + if (not p) throw std::runtime_error("TrajectoryData::add: Null coefficient structure, nothing added!"); + + coefs[roundTime(coef->time)] = p; + } + void SphFldCoefs::add(CoefStrPtr coef) { diff --git a/expui/FieldGenerator.H b/expui/FieldGenerator.H index 4ab573aef..d83cc4107 100644 --- a/expui/FieldGenerator.H +++ b/expui/FieldGenerator.H @@ -115,6 +115,11 @@ namespace Field histogram1d(PR::PRptr reader, double rmax, int nbins, std::string proj, std::vector center={0.0, 0.0, 0.0}); + //! Compute spherical log histogram from particles + std::tuple + histo1dlog(PR::PRptr reader, double rmin, double rmax, int nbins, + std::vector center={0.0, 0.0, 0.0}); + //! Write field slices to files. This will be VTK your build is //! compiled with VTK and ascii tables otherwise. void file_slices(BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs, diff --git a/expui/FieldGenerator.cc b/expui/FieldGenerator.cc index 73a1595f7..b12df7592 100644 --- a/expui/FieldGenerator.cc +++ b/expui/FieldGenerator.cc @@ -894,6 +894,58 @@ namespace Field } // END histogram1d + std::tuple + FieldGenerator::histo1dlog(PR::PRptr reader, double rmin, double rmax, + int nbins, std::vector ctr) + { + if (rmin <= 0.0) throw std::runtime_error("FieldGenerator::histo1dlog: rmin must be > 0.0"); + if (rmax <= rmin) throw std::runtime_error("FieldGenerator::histo1dlog: rmax must be > rmin"); + + const double pi = 3.14159265358979323846; + + Eigen::VectorXf rad = Eigen::VectorXf::Zero(nbins); + Eigen::VectorXf ret = Eigen::VectorXf::Zero(nbins); + + double lrmin = log(rmin), lrmax = log(rmax); + double del = (lrmax - lrmin)/nbins; + + // Make the histogram + // + for (auto p=reader->firstParticle(); p!=0; p=reader->nextParticle()) { + double rad = 0.0; + for (int k=0; k<3; k++) { + double pp = p->pos[k] - ctr[k]; + rad += pp*pp; + } + + int indx = floor((log(sqrt(rad)) - lrmin)/del); + if (indx>=0 and indxmass; + } + + // Accumulate between MPI nodes; return value to root node + // + if (use_mpi) { + if (myid==0) + MPI_Reduce(MPI_IN_PLACE, ret.data(), ret.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + else + MPI_Reduce(ret.data(), NULL, ret.size(), + MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + } + + // Compute density + // + double d3 = exp(3.0*del); + double rf = 4.0*pi/3.0*(d3 - 1.0); + for (int i=0; i> FieldGenerator::points(BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs) diff --git a/expui/expMSSA.H b/expui/expMSSA.H index 1200f0d97..a15cfeb4c 100644 --- a/expui/expMSSA.H +++ b/expui/expMSSA.H @@ -50,6 +50,13 @@ namespace MSSA //! Right singular vectors Eigen::MatrixXd U; + //@{ + //! Koopman modes with Hankel matrices + Eigen::VectorXcd L; + Eigen::MatrixXcd Phi; + int window; + //@} + //! Parameters //@{ bool flip, verbose, powerf; @@ -280,6 +287,13 @@ namespace MSSA return ret; } + //! Estimate Koopman modes from the trajectory eigenvectors + std::tuple + getKoopmanModes(const double tol, int window, bool debug); + + //! Return the reconstructed Koopman modes + std::map getReconstructedKoopman(int mode); + }; diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc index b136d2e7d..8940d98c5 100644 --- a/expui/expMSSA.cc +++ b/expui/expMSSA.cc @@ -1554,6 +1554,259 @@ namespace MSSA { } + std::tuple + expMSSA::getKoopmanModes(double tol, int D, bool debug) + { + bool use_fullKh = true; // Use the non-reduced computation of + // Koopman/eDMD + // Number of channels + // + nkeys = mean.size(); + + // Make sure parameters are sane + // + if (numW<=0) numW = numT/2; + if (numW > numT/2) numW = numT/2; + + numK = numT - numW + 1; + + Eigen::VectorXd S1; + Eigen::MatrixXd Y1; + Eigen::MatrixXd V1; + Eigen::MatrixXd VT1; + Eigen::MatrixXd VT2; + + // Make a new trajetory matrix with smoothing + // + Y1.resize(numK, numW*nkeys + D*(nkeys-1)); + + size_t n=0, offset=0; + for (auto k : mean) { + for (int i=0; i 0) { + // Back blending + for (int j=0; j(D-j)/D; + } + } + // Main series + for (int j=0; j(D-j)/D; + } + } + } + offset += numW + D; + n++; + } + + // double Scale = Y1.norm(); + // auto YY = Y1/Scale; + + auto YY = Y1; + + // Use one of the built-in Eigen3 algorithms + // + if (params["Jacobi"]) { + // -->Using Jacobi + Eigen::JacobiSVD + svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); + S1 = svd.singularValues(); + V1 = svd.matrixV(); + } + // -->Using BDC + else if (params["BDCSVD"]) { + Eigen::BDCSVD + svd(YY, Eigen::ComputeFullU | Eigen::ComputeFullV); + // svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); + S1 = svd.singularValues(); + V1 = svd.matrixV(); + } + // -->Use Random approximation algorithm from Halko, Martinsson, + // and Tropp + else { + int srank = std::min(YY.cols(), YY.rows()); + RedSVD::RedSVD svd(YY, srank); + S1 = svd.singularValues(); + V1 = svd.matrixV(); + } + + std::cout << "shape V1 = " << V1.rows() << " x " + << V1.cols() << std::endl; + + std::cout << "shape Y1 = " << Y1.rows() << " x " + << Y1.cols() << std::endl; + + int lags = V1.rows(); + int rank = V1.cols(); + + std::ofstream out; + if (debug) out.open("debug.txt"); + + if (out) out << "rank=" << rank << " lags=" << lags << std::endl; + + VT1.resize(rank, lags-1); + VT2.resize(rank, lags-1); + + for (int j=0; j uu; + for (int i=0; iUsing Jacobi + if (params["Jacobi"]) { + Eigen::JacobiSVD + // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); + svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); + SS = svd.singularValues(); + UU = svd.matrixU(); + VV = svd.matrixV(); + } + // -->Using BDC + else if (params["BDCSVD"]) { + Eigen::BDCSVD + // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); + svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); + SS = svd.singularValues(); + UU = svd.matrixU(); + VV = svd.matrixV(); + } + // -->Use Random approximation algorithm from Halko, Martinsson, + // and Tropp + else { + // RedSVD::RedSVD svd(VT1, std::min(rank, numK-1)); + RedSVD::RedSVD svd(VT1, std::max(VT1.rows(), VT2.cols())); + SS = svd.singularValues(); + UU = svd.matrixU(); + VV = svd.matrixV(); + } + + if (out) out << "Singular values" << std::endl << SS << std::endl; + + // Compute inverse + for (int i=0; i tol) SS(i) = 1.0/SS(i); + else SS(i) = 0.0; + } + + // Compute full Koopman operator + if (use_fullKh) { + + Eigen::MatrixXd DD(VV.cols(), UU.cols()); + DD.setZero(); + for (int i=0; i es(AT); + + L = es.eigenvalues(); + Phi = es.eigenvectors(); + + if (out) { + out << std::endl << "Eigenvalues" << std::endl << L << std::endl + << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; + } + + } + // Compute the reduced Koopman operator + else { + + Eigen::MatrixXd AT = UU.transpose() * (VT2 * VV) * SS.asDiagonal(); + + // Compute spectrum + Eigen::EigenSolver es(AT, true); + + L = es.eigenvalues(); + auto W = es.eigenvectors(); + + // Compute the EDMD modes + // + Eigen::VectorXcd Linv(L); + for (int i=0; i tol) Linv(i) = 1.0/Linv(i); + else Linv(i) = 0.0; + } + + Phi = VT2 * VV * SS.asDiagonal() * W * Linv.asDiagonal(); + + if (out) { + out << std::endl << "Eigenvalues" << std::endl << L << std::endl + << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; + } + } + + // Cache window size + // + window = D; + + return {L, Phi}; + } + + std::map + expMSSA::getReconstructedKoopman(int mode) + { + // Copy the original map for return + // + auto newdata = data; + + size_t n=0, offset=0; + + for (auto u : mean) { + + double disp = totVar; + if (type == TrendType::totPow) disp = totPow; + if (disp==0.0) disp = var[u.first]; + + std::complex phase = 1.0; + for (int i=0; i(); - else nmax = nmaxfid; + else nmax = 12; if (conf["numr"]) numr = conf["numr"].as(); else numr = 2000; @@ -836,3 +836,38 @@ std::vector BiorthCyl::orthoCheck() return emp.orthoCheck(); } + + +void BiorthCyl::dump_basis(const string& name) +{ + const int num = 1000; + double rmin = emp.getRmin(); + double rmax = emp.getRmax(); + double r, dr = rmax/(num-1); + + std::ofstream out; + + for (int m=0; m<=mmax; m++) { + + std::ostringstream ins; + ins << name << ".biorthcyl." << m; + out.open(ins.str().c_str()); + + // Ok, write data + + for (int j=0; j ${DEP_INC} ${EIGEN3_INCLUDE_DIR} ${HDF5_INCLUDE_DIRS} ${FFTW_INCLUDE_DIRS}) -set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX - yaml-cpp ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES} - ${FFTW_DOUBLE_LIB}) +set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp + ${VTK_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES} + ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB}) if(ENABLE_CUDA) list(APPEND common_LINKLIB CUDA::toolkit CUDA::cudart) diff --git a/exputil/EmpCyl2d.cc b/exputil/EmpCyl2d.cc index af0901e4e..efd29bcdf 100644 --- a/exputil/EmpCyl2d.cc +++ b/exputil/EmpCyl2d.cc @@ -9,9 +9,7 @@ #include #include -#include #include -#include // Set to true for orthogonality checking // @@ -331,293 +329,6 @@ EmpCyl2d::Basis2d::createBasis(int mmax, int nmax, double rmax, return std::make_shared(); } -class EmpCyl2d::ExponCyl : public EmpCyl2d::ModelCyl -{ - -private: - - double acyl, sigma0; - - // Assign values from YAML - // - void parse(const YAML::Node& conf) - { - try { - if (conf["acyl"]) - acyl = conf["acyl"].as(); - else - acyl = 1.0; - } - catch (YAML::Exception & error) { - if (myid==0) - std::cout << "Error parsing parameters in EmpCyl2d::ExponCyl: " - << error.what() << std::endl - << std::string(60, '-') << std::endl - << "Config node" << std::endl - << std::string(60, '-') << std::endl - << conf << std::endl - << std::string(60, '-') << std::endl; - throw std::runtime_error("EmpCyl2d::ExponCyl: error parsing YAML config"); - } - } - -public: - - ExponCyl(const YAML::Node& par) - { - parse(par); - sigma0 = 0.5/(M_PI*acyl*acyl); - id = "expon"; - } - - double pot(double r) { - double y = 0.5 * r / acyl; - return -M_PI*sigma0*r * - (EXPmath::cyl_bessel_i(0, y)*EXPmath::cyl_bessel_k(1, y) - - EXPmath::cyl_bessel_i(1, y)*EXPmath::cyl_bessel_k(0, y)); - } - - double dpot(double r) { - double y = 0.5 * r / acyl; - return 2.0*M_PI*sigma0*y* - (EXPmath::cyl_bessel_i(0, y)*EXPmath::cyl_bessel_k(0, y) - - EXPmath::cyl_bessel_i(1, y)*EXPmath::cyl_bessel_k(1, y)); - } - - double dens(double r) { - return sigma0*exp(-r/acyl); - } - -}; - -class EmpCyl2d::KuzminCyl : public EmpCyl2d::ModelCyl -{ -private: - - double acyl; - - // Assign values from YAML - // - void parse(const YAML::Node& conf) - { - try { - if (conf["acyl"]) - acyl = conf["acyl"].as(); - else - acyl = 1.0; - } - catch (YAML::Exception & error) { - if (myid==0) - std::cout << "Error parsing parameters in EmpCyl2d::KuzminCyl: " - << error.what() << std::endl - << std::string(60, '-') << std::endl - << "Config node" << std::endl - << std::string(60, '-') << std::endl - << conf << std::endl - << std::string(60, '-') << std::endl; - throw std::runtime_error("EmpCyl2d::KuzminCyl: error parsing YAML config"); - } - } - -public: - - KuzminCyl(const YAML::Node& par) - { - parse(par); - id = "kuzmin"; - } - - double pot(double R) { - double a2 = acyl * acyl; - return -1.0/sqrt(R*R + a2); - } - - double dpot(double R) { - double a2 = acyl * acyl; - return R/pow(R*R + a2, 1.5); - } - - double dens(double R) { - double a2 = acyl * acyl; - return 4.0*M_PI*acyl/pow(R*R + a2, 1.5)/(2.0*M_PI); - // ^ - // | - // This 4pi from Poisson's eqn - } - -}; - - -class EmpCyl2d::MestelCyl : public EmpCyl2d::ModelCyl -{ -protected: - - double vrot, rot; - - // Assign values from YAML - // - void parse(const YAML::Node& conf) - { - try { - if (conf["vrot"]) - vrot = conf["vrot"].as(); - else - vrot = 1.0; - } - catch (YAML::Exception & error) { - if (myid==0) - std::cout << "Error parsing parameters in EmpCyl2d::MestelCyl: " - << error.what() << std::endl - << std::string(60, '-') << std::endl - << "Config node" << std::endl - << std::string(60, '-') << std::endl - << conf << std::endl - << std::string(60, '-') << std::endl; - throw std::runtime_error("EmpCyl2d::MestelCyl: error parsing YAML config"); - } - } - -public: - - MestelCyl(const YAML::Node& par) - { - parse(par); - rot = vrot*vrot; - id = "mestel"; - } - - virtual double pot(double R) { - if (R>0.0) return rot*log(R); - else throw std::runtime_error("MestelCyl::pot: R<=0"); - } - - virtual double dpot(double R) { - if (R>0.0) return rot/R; - else throw std::runtime_error("MestelCyl::dpot: R<=0"); - } - - virtual double dens(double R) { - if (R>0.0) return rot/(2.0*M_PI*R); - else throw std::runtime_error("MestelCyl::dens: R<=0"); - } -}; - - -class EmpCyl2d::ZangCyl : public EmpCyl2d::MestelCyl -{ - -private: - //! Parameters - double mu, nu, ri, ro; - - //! Softening factor (not currently used) - double asoft = 1.0e-8; - - //! Ignore inner cut-off for N<0.05 - bool Inner = true; - - //! Taper factors - double Tifac, Tofac; - - //! Inner taper function - double Tinner(double Jp) - { - double fac = pow(Jp, nu); - return fac/(Tifac + fac); - } - - //! Outer taper function - double Touter(double Jp) - { - return 1.0/(1.0 + pow(Jp/Tofac, mu)); - } - - //! Deriv of inner taper function - double dTinner(double Jp) - { - double fac = pow(Jp, nu); - double fac2 = Tifac + fac; - return nu*fac/Jp/(fac2*fac2); - } - - //! Deriv of outer taper function - double dTouter(double Jp) - { - double fac = pow(Jp/Tofac, mu); - double fac2 = 1.0 + fac; - return -mu*fac/Jp/(fac2*fac2); - } - -protected: - - //! Assign values from YAML - void parse(const YAML::Node& conf) - { - try { - if (conf["Ninner"]) - nu = conf["Ninner"].as(); - else - nu = 2.0; - - if (conf["Mouter"]) - mu = conf["Mouter"].as(); - else - mu = 2.0; - - if (conf["Ri"]) - ri = conf["Ri"].as(); - else - ri = 1.0; - - if (conf["Ro"]) - ro = conf["Ro"].as(); - else - ro = 10.0; - } - catch (YAML::Exception & error) { - if (myid==0) - std::cout << "Error parsing parameters in EmpCyl2d::ZangCyl: " - << error.what() << std::endl - << std::string(60, '-') << std::endl - << "Config node" << std::endl - << std::string(60, '-') << std::endl - << conf << std::endl - << std::string(60, '-') << std::endl; - throw std::runtime_error("EmpCyl2d::ZangCyl: error parsing YAML config"); - } - } - -public: - - //! Constructor - ZangCyl(const YAML::Node& par) : MestelCyl(par) - { - // Parse the YAML - parse(par); - - // Assign the id - id = "zang"; - - // Cache taper factors - Tifac = pow(ri*vrot, nu); - Tofac = ro*vrot; - - if (nu<0.05) { - // Exponent is now for mapping only - Inner = false; - } - } - - //! Surface density - double dens(double R) - { - double ret = MestelCyl::dens(R) * Touter(R); - if (Inner) ret *= Tinner(R); - return ret; - } - -}; - std::shared_ptr EmpCyl2d::createModel() { try { @@ -638,7 +349,7 @@ std::shared_ptr EmpCyl2d::createModel() if (name.find("kuzmin") != std::string::npos) { if (myid==0) { std::cout << "---- EmpCyl2d::ModelCyl: Making a Kuzmin disk"; - if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + if (Params["parameters"]) std::cout << " with" << std::endl << Params["parameters"]; std::cout << std::endl; } return std::make_shared(Params["parameters"]); @@ -647,7 +358,7 @@ std::shared_ptr EmpCyl2d::createModel() if (name.find("mestel") != std::string::npos) { if (myid==0) { std::cout << "---- EmpCyl2d::ModelCyl: Making a finite Mestel disk"; - if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + if (Params["parameters"]) std::cout << " with" << std::endl << Params["parameters"]; std::cout << std::endl; } return std::make_shared(Params["parameters"]); @@ -656,7 +367,9 @@ std::shared_ptr EmpCyl2d::createModel() if (name.find("zang") != std::string::npos) { if (myid==0) { std::cout << "---- EmpCyl2d::ModelCyl: Making a double-tapered Zang"; - if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + YAML::Emitter out; + out << YAML::Flow << Params["parameters"]; + if (Params["parameters"]) std::cout << " with " << out.c_str(); std::cout << std::endl; } return std::make_shared(Params["parameters"]); @@ -665,7 +378,7 @@ std::shared_ptr EmpCyl2d::createModel() if (name.find("expon") != std::string::npos) { if (myid==0) { std::cout << "---- EmpCyl2d::ModelCyl: Making an Exponential disk"; - if (Params["parameters"]) std::cout << " with " << Params["parameters"]; + if (Params["parameters"]) std::cout << " with" << std::endl << Params["parameters"]; std::cout << std::endl; } return std::make_shared(Params["parameters"]); @@ -767,6 +480,7 @@ EmpCyl2d::EmpCyl2d if (cache_name_2d.size()==0) cache_name_2d = default_cache_name; disk = createModel(); + model = disk->ID(); basis = Basis2d::createBasis(mmax, nmaxfid, rmax, biorth); basis_test = false; diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 01b941501..17c9b51e3 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -1253,6 +1253,7 @@ YAML::Node EmpCylSL::getHeader_hdf5(const std::string& cachefile) node["mmax"] = getInt("mmax"); node["numx"] = getInt("numx"); node["numy"] = getInt("numy"); + node["lmaxfid"] = getInt("lmaxfid"); node["nmaxfid"] = getInt("nmaxfid"); node["nmax"] = getInt("nmax"); node["neven"] = getInt("neven"); diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index c0934aebb..f3aa05107 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -335,7 +335,8 @@ namespace PR { const H5std_string GROUP_NAME_what ("/Header"); H5::H5File file( FILE_NAME, H5F_ACC_RDONLY ); - H5::Group what(file.openGroup( GROUP_NAME_what )); + // H5::Group what(file.openGroup( GROUP_NAME_what )); + H5::Group what = file.openGroup( GROUP_NAME_what ); // Get time { @@ -415,7 +416,8 @@ namespace PR { const H5std_string GROUP_NAME_what ("/Header"); H5::H5File file( FILE_NAME, H5F_ACC_RDONLY ); - H5::Group what(file.openGroup( GROUP_NAME_what )); + // H5::Group what(file.openGroup( GROUP_NAME_what )); + H5::Group what = file.openGroup( GROUP_NAME_what ); // Get time { @@ -448,7 +450,8 @@ namespace PR { totalCount = npart[ptype]; std::string grpnam = "/" + sout.str(); - H5::Group grp(file.openGroup(grpnam)); + // H5::Group grp(file.openGroup(grpnam)); + H5::Group grp = file.openGroup(grpnam); H5::DataSet dataset = grp.openDataSet("Coordinates"); H5::DataSpace dataspace = dataset.getSpace(); diff --git a/exputil/VtkGrid.cc b/exputil/VtkGrid.cc index 6ba83799f..395e56317 100644 --- a/exputil/VtkGrid.cc +++ b/exputil/VtkGrid.cc @@ -1,5 +1,8 @@ +#include #include +#ifdef HAVE_VTK + VtkGrid::VtkGrid(int nx, int ny, int nz, double xmin, double xmax, double ymin, double ymax, @@ -126,3 +129,161 @@ void VtkGrid::Write(const std::string& name) // writer->SetDataModeToBinary(); writer->Write(); } + +#else + +VtkGrid::VtkGrid(int nx, int ny, int nz, + double xmin, double xmax, + double ymin, double ymax, + double zmin, double zmax) : + ThreeDGrid(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax) +{ + // Set knots + coord["X"].resize(nx); + coord["Y"].resize(ny); + coord["Z"].resize(nz); + + for (int i=0; i1) { + for (int i=0; i& data, const std::string& name) +{ + // Remove XML sensitive characters; paraview bug? + std::string newName(name); + replaceAll(newName, ">", ".gt."); + replaceAll(newName, "<", ".lt."); + + dataSet[newName].resize(nx*ny*nz); + + // Insert grid data + // + int I = 0; + + for (int i=0; i(data[(k*ny + j)*nx + i]); + } + } + } + +} + + +void VtkGrid::writeBeg(std::ofstream & fout) +{ + int zero = 0; + fout << "" << std::endl; + fout << "" << std::endl; + fout << " " << std::endl; + fout << " " << std::endl; +} + +void VtkGrid::writeFields(std::ofstream & fout) +{ + fout << " " << std::endl; + for (auto v : dataSet) { + // Get ranges + auto vmin = std::min_element(v.second.begin(), v.second.end()); + auto vmax = std::max_element(v.second.begin(), v.second.end()); + fout << " " << std::endl; + int cnt = 0; + for (auto & d : v.second) { + if (cnt % 6 == 0) fout << " "; + fout << std::scientific << std::setprecision(8) << d << " "; + if (++cnt % 6 == 0) fout << std::endl; + } + if (cnt % 6) fout << std::endl; + fout << " " << std::endl; + } + fout << " " << std::endl; + + // No cell data here so this stanza is empty, but it's part of the spec + fout << " " << std::endl; + fout << " " << std::endl; +} + +void VtkGrid::writeCoords(std::ofstream & fout) +{ + fout << " " << std::endl; + for (auto & c : coord) { + fout << " " + << std::endl; + + int cnt = 0; + for (auto & v : c.second) { + if (cnt % 6 == 0) fout << " "; + fout << std::scientific << std::setprecision(8) << v << " "; + if (++cnt % 6 == 0) fout << std::endl; + } + if (cnt % 6) fout << std::endl; + fout << " " << std::endl; + } + fout << " " << std::endl; +} + +void VtkGrid::writeEnd(std::ofstream & fout) +{ + fout << " " << std::endl; + fout << " " << std::endl; + fout << "" << std::endl; +} + +void VtkGrid::Write(const std::string& name) +{ + // Create the filename with the correct extension for Paraview + std::ostringstream filename; + filename << name << ".vtr"; + + std::ofstream fout(filename.str()); + if (fout) { + + } else { + throw std::runtime_error + (std::string("VtkGrid::Write: could not open file <") + filename.str() + ">"); + } + + writeBeg(fout); // VTK XML preamble; open stanzas + writeFields(fout); // Write each data field in the dataset map + writeCoords(fout); // Write the coordinate grids + writeEnd(fout); // VTK XML finish; close open stanzas +} + +#endif diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index b052c7979..5a4d25881 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -173,7 +173,10 @@ public: cacheInfo(const std::string& cachefile, bool verbose=true); #if HAVE_LIBCUDA==1 - void initialize_cuda(std::vector& cuArray, + //! Initialize CUDA arrays. 'disk' is the model object for the + //! background m=0 potential and radial force + void initialize_cuda(std::shared_ptr disk, + std::vector& cuArray, thrust::host_vector& tex); virtual cudaMappingConstants getCudaMappingConstants() @@ -201,13 +204,6 @@ public: //! Evaluate all orders in matrices; for n-body void get_pot(Eigen::MatrixXd& Vc, Eigen::MatrixXd& Vs, double r, double z); - //! Background evaluation - virtual std::tuple background(double r, double z) - { - auto [p, dr, d] = emp.background(r); - return {p, dr, 0.0}; - } - //! Get the table bounds double getRtable() { return rcylmax*scale; } @@ -220,6 +216,14 @@ public: double getYmin() { return ymin; } double getYmax() { return ymax; } + //! Dump the basis for plotting + void dump_basis(const string& name); + + //! Get model name + const std::string getModelName() { return emp.getModelName(); } + + //! Get biorthogonal function base name + const std::string getBiorthName() { return emp.getBiorthName(); } }; #endif diff --git a/include/DataGrid.H b/include/DataGrid.H index 99f457a84..befc11055 100644 --- a/include/DataGrid.H +++ b/include/DataGrid.H @@ -8,10 +8,7 @@ #include #include -#ifdef HAVE_VTK #include -#endif - /** This implementation of ThreeDGrid instantiates at VtkGrid if VTK is @@ -32,19 +29,22 @@ public: DataGrid(int nx, int ny, int nz, double xmin, double xmax, double ymin, double ymax, - double zmin, double zmax) + double zmin, double zmax, + std::string type="VTK") { -#ifdef HAVE_VTK - ptr = std::make_shared(nx, ny, nz, - xmin, xmax, - ymin, ymax, - zmin, zmax); -#else - ptr = std::make_shared(nx, ny, nz, + std::transform(type.begin(), type.end(), type.begin(), + [](unsigned char c){ return std::tolower(c); }); + + if (type.compare("vtk") == 0) + ptr = std::make_shared(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax); -#endif + else + ptr = std::make_shared(nx, ny, nz, + xmin, xmax, + ymin, ymax, + zmin, zmax); } //! Add data diff --git a/include/EmpCyl2d.H b/include/EmpCyl2d.H index 703e8e610..7a5ca6bb8 100644 --- a/include/EmpCyl2d.H +++ b/include/EmpCyl2d.H @@ -10,11 +10,14 @@ #include #include -#include #include // Parameters parsing #include // For config macros +#include +#include +#include + // For reading and writing cache file // #include @@ -69,6 +72,295 @@ public: double d_xi_to_r(double xi); }; + + class ExponCyl : public ModelCyl + { + + private: + + double acyl, sigma0; + + // Assign values from YAML + // + void parse(const YAML::Node& conf) + { + try { + if (conf["acyl"]) + acyl = conf["acyl"].as(); + else + acyl = 1.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::ExponCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("EmpCyl2d::ExponCyl: error parsing YAML config"); + } + } + + public: + + ExponCyl(const YAML::Node& par) + { + parse(par); + sigma0 = 0.5/(M_PI*acyl*acyl); + id = "expon"; + } + + double pot(double r) { + double y = 0.5 * r / acyl; + return -M_PI*sigma0*r * + (EXPmath::cyl_bessel_i(0, y)*EXPmath::cyl_bessel_k(1, y) - + EXPmath::cyl_bessel_i(1, y)*EXPmath::cyl_bessel_k(0, y)); + } + + double dpot(double r) { + double y = 0.5 * r / acyl; + return 2.0*M_PI*sigma0*y* + (EXPmath::cyl_bessel_i(0, y)*EXPmath::cyl_bessel_k(0, y) - + EXPmath::cyl_bessel_i(1, y)*EXPmath::cyl_bessel_k(1, y)); + } + + double dens(double r) { + return sigma0*exp(-r/acyl); + } + + }; + + class KuzminCyl : public ModelCyl + { + private: + + double acyl; + + // Assign values from YAML + // + void parse(const YAML::Node& conf) + { + try { + if (conf["acyl"]) + acyl = conf["acyl"].as(); + else + acyl = 1.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::KuzminCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("EmpCyl2d::KuzminCyl: error parsing YAML config"); + } + } + + public: + + KuzminCyl(const YAML::Node& par) + { + parse(par); + id = "kuzmin"; + } + + double pot(double R) { + double a2 = acyl * acyl; + return -1.0/sqrt(R*R + a2); + } + + double dpot(double R) { + double a2 = acyl * acyl; + return R/pow(R*R + a2, 1.5); + } + + double dens(double R) { + double a2 = acyl * acyl; + return 4.0*M_PI*acyl/pow(R*R + a2, 1.5)/(2.0*M_PI); + // ^ + // | + // This 4pi from Poisson's eqn + } + + }; + + + class MestelCyl : public ModelCyl + { + protected: + + double vrot, rot; + + // Assign values from YAML + // + void parse(const YAML::Node& conf) + { + try { + if (conf["vrot"]) + vrot = conf["vrot"].as(); + else + vrot = 1.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::MestelCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("EmpCyl2d::MestelCyl: error parsing YAML config"); + } + } + + public: + + MestelCyl(const YAML::Node& par) + { + parse(par); + rot = vrot*vrot; + id = "mestel"; + } + + virtual double pot(double R) { + if (R>0.0) return rot*log(R); + else throw std::runtime_error("MestelCyl::pot: R<=0"); + } + + virtual double dpot(double R) { + if (R>0.0) return rot/R; + else throw std::runtime_error("MestelCyl::dpot: R<=0"); + } + + virtual double dens(double R) { + if (R>0.0) return rot/(2.0*M_PI*R); + else throw std::runtime_error("MestelCyl::dens: R<=0"); + } + }; + + + class ZangCyl : public MestelCyl + { + + private: + //! Parameters + double mu, nu, ri, ro; + + //! Softening factor (not currently used) + double asoft = 1.0e-8; + + //! Ignore inner cut-off for N<0.05 + bool Inner = true; + + //! Taper factors + double Tifac, Tofac; + + //! Inner taper function + double Tinner(double Jp) + { + double fac = pow(Jp, nu); + return fac/(Tifac + fac); + } + + //! Outer taper function + double Touter(double Jp) + { + return 1.0/(1.0 + pow(Jp/Tofac, mu)); + } + + //! Deriv of inner taper function + double dTinner(double Jp) + { + double fac = pow(Jp, nu); + double fac2 = Tifac + fac; + return nu*fac/Jp/(fac2*fac2); + } + + //! Deriv of outer taper function + double dTouter(double Jp) + { + double fac = pow(Jp/Tofac, mu); + double fac2 = 1.0 + fac; + return -mu*fac/Jp/(fac2*fac2); + } + + protected: + + //! Assign values from YAML + void parse(const YAML::Node& conf) + { + try { + if (conf["Ninner"]) + nu = conf["Ninner"].as(); + else + nu = 2.0; + + if (conf["Mouter"]) + mu = conf["Mouter"].as(); + else + mu = 2.0; + + if (conf["Ri"]) + ri = conf["Ri"].as(); + else + ri = 1.0; + + if (conf["Ro"]) + ro = conf["Ro"].as(); + else + ro = 10.0; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in EmpCyl2d::ZangCyl: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("EmpCyl2d::ZangCyl: error parsing YAML config"); + } + } + + public: + + //! Constructor + ZangCyl(const YAML::Node& par) : MestelCyl(par) + { + // Parse the YAML + parse(par); + + // Assign the id + id = "zang"; + + // Cache taper factors + Tifac = pow(ri*vrot, nu); + Tofac = ro*vrot; + + if (nu<0.05) { + // Exponent is now for mapping only + Inner = false; + } + } + + //! Surface density + double dens(double R) + { + double ret = MestelCyl::dens(R) * Touter(R); + if (Inner) ret *= Tinner(R); + return ret; + } + + }; + + protected: //! Contains parameter database @@ -113,11 +405,6 @@ protected: std::shared_ptr createModel(); - class ExponCyl; - class MestelCyl; - class KuzminCyl; - class ZangCyl; - //! Mapping instance Mapping map; @@ -206,22 +493,15 @@ public: //! Check for configuration bool operator()() { return configured; } - //@{ - //! Background evaluation - void background(double r, double & p, double & dr) - { - p = disk->pot(r); - dr = disk->dpot(r); - } - - std::tuple background(double r) - { - return {disk->pot(r), disk->dpot(r), disk->dens(r)}; - } - //@} - //! Get coordinate mapping object Mapping getMapping() { return map; } + + //! Get model name + const std::string getModelName() { return model; } + + //! Get biorthogonal function base name + const std::string getBiorthName() { return biorth; } + }; diff --git a/include/KDtree.H b/include/KDtree.H index 40a845b32..1bce71225 100644 --- a/include/KDtree.H +++ b/include/KDtree.H @@ -6,381 +6,414 @@ #include #include -//! Wrapper class for std::map to keep a finite size -template -class Cache : public std::map +namespace KDtree { -public: - Cache(int N=1) : max_size(N) {} + //! Wrapper class for std::map to keep a finite size + template + class Cache : public std::map + { + public: + Cache(int N=1) : max_size(N) {} - void resize(int N) { max_size = N; } + void resize(int N) { max_size = N; } - void add( key_t key, value_t value ) - { - (*this)[key] = value; - if (this->size() > max_size) { - this->erase(--(*this).end()); + void add( key_t key, value_t value ) + { + (*this)[key] = value; + if (this->size() > max_size) { + this->erase(--(*this).end()); + } } - } - -private: - int max_size; -}; + + private: + int max_size; + }; -/** Class for representing a point - Coordinate_type must be a numeric type - Field (weight) is a double -*/ -template -class point -{ -public: - //! Null initializer - point() {} - - //! Constructor - point(std::array c, double m=1, - unsigned long indx=0) : coords_(c), mass_(m), indx_(indx) {} + /** Class for representing a point + Coordinate_type must be a numeric type + Field (weight) is a double + */ + template + class point + { + public: + //! Null initializer + point() {} + + //! Constructor + point(std::array c, double m=1, + unsigned long indx=0) : coords_(c), mass_(m), indx_(indx) {} - point(std::array c, - std::array v, - double m=1, unsigned long indx=0) : coords_(c), vels_(v), - mass_(m), indx_(indx) {} + point(std::array c, + std::array v, + double m=1, unsigned long indx=0) : coords_(c), vels_(v), + mass_(m), indx_(indx) {} - //! List constructor - point(std::initializer_list list, double m=1, - unsigned long indx=0) : mass_(m), indx_(indx) - { - size_t n = std::min(dimensions, list.size()); - std::copy_n(list.begin(), n, coords_.begin()); - } + //! List constructor + point(std::initializer_list list, double m=1, + unsigned long indx=0) : mass_(m), indx_(indx) + { + size_t n = std::min(dimensions, list.size()); + std::copy_n(list.begin(), n, coords_.begin()); + } + + point(std::initializer_list list, + std::initializer_list vlst, + double m=1, + unsigned long indx=0) : mass_(m), indx_(indx) + { + size_t n = std::min(dimensions, list.size()); + std::copy_n(list.begin(), n, coords_.begin()); + std::copy_n(vlst.begin(), n, vels_.begin()); + } - point(std::initializer_list list, - std::initializer_list vlst, - double m=1, - unsigned long indx=0) : mass_(m), indx_(indx) - { - size_t n = std::min(dimensions, list.size()); - std::copy_n(list.begin(), n, coords_.begin()); - std::copy_n(vlst.begin(), n, vels_.begin()); - } + //! Copy constructor + point(const point& p) + { + for (size_t n=0; n coords_; + std::array vels_; + double mass_; + unsigned long indx_; + }; - double mass() const + //! For iostream printing of points + template + std::ostream& operator<<(std::ostream& out, + const point& pt) { - return mass_; + out << '('; + for (size_t i = 0; i < dimensions; ++i) { + if (i > 0) out << ", "; + out << pt.get(i); + } + out << ')'; + return out; } - unsigned long indx() const + //! k-d tree implementation + template + class kdtree { - return indx_; - } + public: -private: - std::array coords_; - std::array vels_; - double mass_; - unsigned long indx_; -}; - -//! For iostream printing of points -template -std::ostream& operator<<(std::ostream& out, - const point& pt) -{ - out << '('; - for (size_t i = 0; i < dimensions; ++i) { - if (i > 0) out << ", "; - out << pt.get(i); - } - out << ')'; - return out; -} - -//! k-d tree implementation -template -class kdtree -{ -public: - typedef point point_type; + using point_type = point; -private: - struct node - { - node(const point_type& pt) : point_(pt), left_(nullptr), right_(nullptr) + struct node { - } + node(const point_type& pt) : point_(pt), left_(nullptr), right_(nullptr) + { + } - coordinate_type get(size_t index) const - { - return point_.get(index); - } + coordinate_type get(size_t index) const + { + return point_.get(index); + } + + double distance(const point_type& pt) const + { + return point_.distance(pt); + } + point_type point_; + node* left_; + node* right_; + }; - double distance(const point_type& pt) const - { - return point_.distance(pt); - } - point_type point_; - node* left_; - node* right_; - }; + using cache = Cache; - node* root_; + private: - Cache best_; + node* root_; - size_t visited_; - std::vector nodes_; + size_t visited_; + std::vector nodes_; - struct node_cmp - { - node_cmp(size_t index) : index_(index) + struct node_cmp { - } + node_cmp(size_t index) : index_(index) + { + } - bool operator()(const node& n1, const node& n2) const + bool operator()(const node& n1, const node& n2) const + { + return n1.point_.get(index_) < n2.point_.get(index_); + } + size_t index_; + }; + + node* make_tree(size_t begin, size_t end, size_t index) { - return n1.point_.get(index_) < n2.point_.get(index_); + if (end <= begin) return nullptr; + size_t n = begin + (end - begin)/2; + std::nth_element(&nodes_[begin], &nodes_[n], &nodes_[end], node_cmp(index)); + index = (index + 1) % dimensions; + nodes_[n].left_ = make_tree(begin, n, index); + nodes_[n].right_ = make_tree(n + 1, end, index); + return &nodes_[n]; } - size_t index_; - }; - node* make_tree(size_t begin, size_t end, size_t index) - { - if (end <= begin) return nullptr; - size_t n = begin + (end - begin)/2; - std::nth_element(&nodes_[begin], &nodes_[n], &nodes_[end], node_cmp(index)); - index = (index + 1) % dimensions; - nodes_[n].left_ = make_tree(begin, n, index); - nodes_[n].right_ = make_tree(n + 1, end, index); - return &nodes_[n]; - } - - void nearestN(node* root, const point_type& point, size_t index, int N) - { - if (root == nullptr) return; + void nearestN(node* root, const point_type& point, size_t index, int N, + cache& best_) + { + if (root == nullptr) return; - ++visited_; + ++visited_; // For diagnostic info only - double d = root->distance(point); - if (best_.size()first) { - best_.add(d, root); + double d = root->distance(point); + if (best_.size()first) { + best_.add(d, root); + } + + // This is only correct if the test point is never in the data set . . . + // if (best_.begin()->first == 0) return; + + double dx = root->get(index) - point.get(index); + index = (index + 1) % dimensions; + + nearestN(dx > 0 ? root->left_ : root->right_, point, index, N, best_); + + if (best_.size()>=N and dx * dx >= best_.rbegin()->first) return; + nearestN(dx > 0 ? root->right_ : root->left_, point, index, N, best_); + } + + void accum(node* root, std::vector>& ret, unsigned cur) + { + ret[cur].push_back(root->point_.indx()); + if (root->left_) accum(root->left_, ret, cur); + if (root->right_) accum(root->right_, ret, cur); } - // This is only correct if the test point is never in the data set . . . - // if (best_.begin()->first == 0) return; + void walk(node* root, std::vector>& ret, + unsigned int cur, unsigned int lev, unsigned int level) + { + if (lev == level) accum(root, ret, cur); + else { + ret[cur].push_back(root->point_.indx()); + walk(root->left_, ret, 2*cur + 0, lev + 1, level); + walk(root->right_, ret, 2*cur + 1, lev + 1, level); + } + } - double dx = root->get(index) - point.get(index); - index = (index + 1) % dimensions; - nearestN(dx > 0 ? root->left_ : root->right_, point, index, N); + public: - if (best_.size()>=N and dx * dx >= best_.rbegin()->first) return; - nearestN(dx > 0 ? root->right_ : root->left_, point, index, N); - } - -public: - //@{ - //! Copy constructor2 - kdtree(const kdtree&) = delete; - kdtree& operator=(const kdtree&) = delete; - //@} - - /** - * Constructor taking a pair of iterators. Adds each - * point in the range [begin, end) to the tree. - * - * @param begin start of range - * @param end end of range - */ - template - kdtree(iterator begin, iterator end) - { - best_.resize(1); - visited_ = 0; - nodes_.reserve(std::distance(begin, end)); - for (auto i = begin; i != end; ++i) - nodes_.emplace_back(*i); - root_ = make_tree(0, nodes_.size(), 0); - } - - /** - * Constructor taking a function object that generates - * points. The function object will be called n times - * to populate the tree. - * - * @param f function that returns a point - * @param n number of points to add - */ - template - kdtree(func&& f, size_t n) - { - best_.resize(1); - visited_ = 0; - nodes_.reserve(n); - for (size_t i = 0; i < n; ++i) - nodes_.emplace_back(f()); - root_ = make_tree(0, nodes_.size(), 0); - } - - /** - * Returns true if the tree is empty, false otherwise. - */ - bool empty() const - { - return nodes_.empty(); - } - - /** - * Returns the number of nodes visited by the last call - * to nearest(). - */ - size_t visited() const - { - return visited_; - } + //@{ + //! Copy constructor2 + kdtree(const kdtree&) = delete; + kdtree& operator=(const kdtree&) = delete; + //@} - /** - * Returns the distance between the input point and return value - * from the last call to nearest(). - */ - double distance() const - { - return std::sqrt(best_.begin()->first); - } + /** + * Constructor taking a pair of iterators. Adds each + * point in the range [begin, end) to the tree. + * + * @param begin start of range + * @param end end of range + */ + template + kdtree(iterator begin, iterator end) + { + visited_ = 0; + nodes_.reserve(std::distance(begin, end)); + for (auto i = begin; i != end; ++i) + nodes_.emplace_back(*i); + root_ = make_tree(0, nodes_.size(), 0); + } + + /** + * Constructor taking a function object that generates + * points. The function object will be called n times + * to populate the tree. + * + * @param f function that returns a point + * @param n number of points to add + */ + template + kdtree(func&& f, size_t n) + { + visited_ = 0; + nodes_.reserve(n); + for (size_t i = 0; i < n; ++i) + nodes_.emplace_back(f()); + root_ = make_tree(0, nodes_.size(), 0); + } + + /** + * Returns true if the tree is empty, false otherwise. + */ + bool empty() const + { + return nodes_.empty(); + } - /** - * Finds the nearest N points in the tree to the given point. It is - * not valid to call this function if the tree is empty. - * - * @param pt a point - * @param N is the number of nearest points - * - * Returns: tuple of the first points, summed weight, and the radius of the Nth - * point - */ - std::tuple - nearestN(const point_type& pt, int N) - { - if (root_ == nullptr) throw std::logic_error("tree is empty"); - best_.clear(); - best_.resize(N); - visited_ = 0; - nearestN(root_, pt, 0, N); + /** + * Returns the number of nodes visited by the last call + * to nearest(). + */ + size_t visited() const + { + return visited_; + } + + /** + * Returns the distance between the input point and return value + * from the last call to nearest(). + */ + double distance(cache& best_) const + { + return std::sqrt(best_.begin()->first); + } + + /** + * Finds the nearest N points in the tree to the given point. It is + * not valid to call this function if the tree is empty. + * + * @param pt a point + * @param N is the number of nearest points + * + * Returns: tuple of the first points, summed weight, and the radius of the Nth + * point + */ + std::tuple + nearestN(const point_type& pt, int N) + { + if (root_ == nullptr) throw std::logic_error("tree is empty"); + cache best_(N); + visited_ = 0; + nearestN(root_, pt, 0, N, best_); - double wgt = 0.0; // Sum weights - for (auto b : best_) wgt += b.second->point_.mass(); + double wgt = 0.0; // Sum weights + for (auto b : best_) wgt += b.second->point_.mass(); -#if __GNUC__ > 6 - return {best_.begin()->second->point_, wgt, std::sqrt(best_.rbegin()->first)}; +#if __GNUC__ > 6 or defined(__clang__) + return {best_.begin()->second->point_, wgt, std::sqrt(best_.rbegin()->first), best_}; #else - std::tuple ret; - std::get<0>(ret) = best_.begin()->second->point_; - std::get<1>(ret) = wgt; - std::get<2>(ret) = std::sqrt(best_.rbegin()->first); - return ret; + std::tuple ret; + std::get<0>(ret) = best_.begin()->second->point_; + std::get<1>(ret) = wgt; + std::get<2>(ret) = std::sqrt(best_.rbegin()->first); + std::get<3>(ret) = best_; + return ret; #endif - } - - /** - * Finds the nearest N points in the tree to the given point. It is - * not valid to call this function if the tree is empty. - * - * @param pt a point - * @param N is the number of nearest points + } - * Returns: tuple of the nearest point list and the radius of the - * Nth point - */ - std::tuple, double> - nearestList(const point_type& pt, int N) - { - if (root_ == nullptr) throw std::logic_error("tree is empty"); - best_.clear(); - best_.resize(N); - visited_ = 0; - nearestN(root_, pt, 0, N); + /** + * Finds the nearest N points in the tree to the given point. It is + * not valid to call this function if the tree is empty. + * + * @param pt a point + * @param N is the number of nearest points + + * Returns: tuple of the nearest point list and the radius of the + * Nth point + */ + std::tuple, double, cache> + nearestList(const point_type& pt, int N) + { + if (root_ == nullptr) throw std::logic_error("tree is empty"); + cache best_(N); + visited_ = 0; + nearestN(root_, pt, 0, N, best_); - std::vector pts; // The returned point list - for (auto b : best_) pts.push_back(b.second->point_); + std::vector pts; // The returned point list + for (auto b : best_) pts.push_back(b.second->point_); #if __GNUC__ > 6 - return {pts, std::sqrt(best_.rbegin()->first)}; + return {pts, std::sqrt(best_.rbegin()->first), best_}; #else - std::tuple, double> ret; - std::get<0>(ret) = pts; - std::get<1>(ret) = std::sqrt(best_.rbegin()->first); - return ret; + std::tuple, double> ret; + std::get<0>(ret) = pts; + std::get<1>(ret) = std::sqrt(best_.rbegin()->first); + std::get<2>(ret) = best_; + return ret; #endif - } + } - std::vector getDist() - { - std::vector ret; - for (auto v : best_) ret.push_back(v.first); - return ret; - } + std::vector getDist(cache& best_) + { + std::vector ret; + for (auto v : best_) ret.push_back(v.first); + return ret; + } + + std::vector> getPartition(unsigned level) + { + int partitions = (1 << level); + std::vector> ret(partitions); + walk(root_, ret, 0, 0, level); + return ret; + } -}; + }; +} +// END: namespace KDtree diff --git a/include/PseudoAccel.H b/include/PseudoAccel.H new file mode 100644 index 000000000..c7f885859 --- /dev/null +++ b/include/PseudoAccel.H @@ -0,0 +1,96 @@ +#ifndef _PseudoAccel_H_ +#define _PseudoAccel_H_ + +#include +#include + +#include + +#include + +class PseudoAccel +{ +private: + + //! Maximum number of time points + unsigned int Nsize; + + //! Type flags + bool CENTER=false, AXIS=false; + + //! Queue of elements + using Elem = std::array; + std::deque queue; + + //! Storage + Eigen::Vector3d accel = Eigen::Vector3d::Zero(); + Eigen::Vector3d omega = Eigen::Vector3d::Zero(); + Eigen::Vector3d domdt = Eigen::Vector3d::Zero(); + +public: + + //! Construct a pseudo acceleration estimator + PseudoAccel(unsigned int Nsize, bool center, bool axis) : + Nsize(Nsize), CENTER(center), AXIS(axis) {} + + //! Add a center element to the pseudo accelration estimator + void add(double t, const Eigen::Vector3d& c, const Eigen::Vector3d& a) { + queue.push_back({t, c(0), c(1), c(2), a(0), a(1), a(2)}); + if (queue.size() > Nsize) queue.pop_front(); + } + + std::tuple operator()() + { + if (CENTER or AXIS) { + + std::vector t, x, y, z, u, v, w; + for (auto &e : queue) { + t.push_back(e[0]); + if (CENTER) { + x.push_back(e[1]); + y.push_back(e[2]); + z.push_back(e[3]); + } + if (AXIS) { + u.push_back(e[4]); + v.push_back(e[5]); + w.push_back(e[6]); + } + } + + // Compute the acceleration only if there are enough elements + // + if (queue.size()==Nsize) { + + if (CENTER) { + accel << + 2.0*std::get<0>(QuadLS(t, x).coefs()), + 2.0*std::get<0>(QuadLS(t, y).coefs()), + 2.0*std::get<0>(QuadLS(t, z).coefs()); + } + + if (AXIS) { + auto [a, b, c] = QuadLS(t, u).coefs(); + auto [d, e, f] = QuadLS(t, v).coefs(); + auto [g, h, i] = QuadLS(t, w).coefs(); + + double T = t.back(); // Last eval time + + Eigen::Vector3d n, dndt, d2ndt2; + n << a*T*T+b*T+c, d*T*T+e*T+f, g*T*T+h*T+i; + dndt << 2.0*a*T+b, 2.0*d*T+e, 2.0*g*T+h; + d2ndt2 << 2.0*a, 2.0*d, 2.0*g; + + // Get angular acceleration and its derivative + omega = n.cross(dndt); + domdt = n.cross(d2ndt2); + } + } + } + + // Return center acceleration, angular velocity and its derivative + return {accel, omega, domdt}; + } +}; + +#endif diff --git a/include/QuadLS.H b/include/QuadLS.H new file mode 100644 index 000000000..1f5d4d368 --- /dev/null +++ b/include/QuadLS.H @@ -0,0 +1,74 @@ +#ifndef _QUADLS_H_ +#define _QUADLS_H_ + +#include +#include +#include + +// Take a pair of data arrays (x, y) and perform a quadratic regression +template +class QuadLS +{ +private: + double _a, _b, _c; + + void fit( const T &x, const T &y ) + { + auto n = x.size(); + + if (n==y.size() and n>0) { + + double sumx=0, sumy=0, sumxy=0, sumx2y=0, sumx2=0, sumx3=0, sumx4=0; + + for (int i=0; i0.0) { + _a = (Sx2y*Sxx - Sxy*Sxx2 ) / denom; + _b = (Sxy*Sx2x2 - Sx2y*Sxx2) / denom; + _c = (sumy - sumx2*_a - sumx*_b) / n; + } else _a = _b = _c = 0.0; + + } + else { + _a = _b = _c = 0.0; + } + } + +public: + + //! Null constructor + QuadLS() { } + + //! Constructor + QuadLS( const T &x, const T &y ) + { fit(x, y); } + + //! Destructor + virtual ~QuadLS() { } + + //! Access coefficients + std::tuple coefs() const + { return std::make_tuple( _a, _b, _c ); } + + //! Evaluate the quadratic function + double operator()( const double x ) const + { return _a * x * x + _b * x + _c; } +}; + +#endif diff --git a/include/VtkGrid.H b/include/VtkGrid.H index 0d5e5fb24..cbd01d922 100644 --- a/include/VtkGrid.H +++ b/include/VtkGrid.H @@ -1,10 +1,15 @@ #ifndef _VtkGrid_H #define _VtkGrid_H -#include +#include +#include #include +#include #include #include +#include + +#ifdef HAVE_VTK // // VTK stuff @@ -55,6 +60,56 @@ public: void Write(const std::string& name); }; +#else + +#include + +/** + A ThreeDGrid implementation of the rectangular grid database + implement in VTK and designed for consumption by PyVTK or ParaView. + */ +class VtkGrid : public ThreeDGrid +{ +private: + // Grid coordinates + std::map> coord; + + // Dataset + std::map> dataSet; + + // Write header + void writeBeg(std::ofstream& file); + + // Write footer + void writeEnd(std::ofstream& file); + + // Write scalar fields + void writeFields(std::ofstream& file); + + // Write coordinates + void writeCoords(std::ofstream& file); + + // String replacement utility + void replaceAll(std::string& str, + const std::string& from, + const std::string& to); + +public: + //! Constructor + VtkGrid(int nx, int ny, int nz, + double xmin, double xmax, + double ymin, double ymax, + double zmin, double zmax); + + //! Add data for two-dimensional cylindrical basis + void Add(const std::vector& data, const std::string& name); + + //! Write output file + void Write(const std::string& name); +}; + +#endif + typedef std::shared_ptr VtkGridPtr; #endif diff --git a/include/fpetrap.h b/include/fpetrap.h index 797512fb5..fc511d041 100644 --- a/include/fpetrap.h +++ b/include/fpetrap.h @@ -16,7 +16,7 @@ #endif /// Global to hold error code -volatile int my_err; +static int my_err; /// Very lean error handler. Only good for setting a debugger breakpoint. static void my_fpu_handler(int err) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index e0b72ab25..8393b1318 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -33,19 +33,21 @@ void BasisFactoryClasses(py::module &m) 3. Cylindrical, created created by computing empirical orthogonal functions over a densely sampled SphericalSL basis; - 4. FlatDisk, an EOF rotation of the finite Bessel basis; and + 4. FlatDisk, an EOF rotation of the finite Bessel basis; - 5. Slab, a biorthogonal basis for a slab geometry with a finite + 5. CBDisk, the Clutton-Brock disk basis for testing; + + 6. Slab, a biorthogonal basis for a slab geometry with a finite finite vertical extent. The basis is constructed from direct solution of the Sturm-Liouville equation. - 6. Cube, a periodic cube basis whose functions are the Cartesian + 7. Cube, a periodic cube basis whose functions are the Cartesian eigenfunctions of the Cartesian Laplacian: sines and cosines. - 7. FieldBasis, for computing user-provided quantities from a + 8. FieldBasis, for computing user-provided quantities from a phase-space snapshot. - 8. VelocityBasis, for computing the mean field velocity fields from + 9. VelocityBasis, for computing the mean field velocity fields from a phase-space snapshot. This is a specialized version of FieldBasis. Each of these bases take a YAML configuration file as input. These parameter @@ -159,11 +161,13 @@ void BasisFactoryClasses(py::module &m) 3. FlatDisk uses cylindrical coordinates - 4. Slab uses Cartesian coordinates + 4. CBDisk uses cylindrical coordinates + + 5. Slab uses Cartesian coordinates - 5. Cube uses Cartesian coordinates + 6. Cube uses Cartesian coordinates - 6. FieldBasis and VelocityBasis provides two natural geometries for + 7. FieldBasis and VelocityBasis provides two natural geometries for field evaluation: a two-dimensional (dof=2) polar disk and a three-dimensional (dof=3) spherical geometry that are chosen using the 'dof' parameter. These use cylindrical and spherical @@ -582,6 +586,71 @@ void BasisFactoryClasses(py::module &m) }; + class PyCBDisk : public CBDisk + { + protected: + + std::vector sph_eval(double r, double costh, double phi) override + { + PYBIND11_OVERRIDE(std::vector, CBDisk, sph_eval, r, costh, phi); + } + + std::vector cyl_eval(double R, double z, double phi) override + { + PYBIND11_OVERRIDE(std::vector, CBDisk, cyl_eval, R, z, phi); + } + + std::vector crt_eval(double x, double y, double z) override + { + PYBIND11_OVERRIDE(std::vector, CBDisk, crt_eval, x, y, z); + } + + void load_coefs(CoefClasses::CoefStrPtr coefs, double time) override + { + PYBIND11_OVERRIDE(void, CBDisk, load_coefs, coefs, time); + } + + void set_coefs(CoefClasses::CoefStrPtr coefs) override + { + PYBIND11_OVERRIDE(void, CBDisk, set_coefs, coefs); + } + + const std::string classname() override { + PYBIND11_OVERRIDE(std::string, CBDisk, classname); + } + + const std::string harmonic() override { + PYBIND11_OVERRIDE(std::string, CBDisk, harmonic); + } + + public: + + // Inherit the constructors + using CBDisk::CBDisk; + + std::vector getFields(double x, double y, double z) override + { + PYBIND11_OVERRIDE(std::vector, CBDisk, getFields, x, y, z); + } + + void accumulate(double x, double y, double z, double mass) override + { + PYBIND11_OVERRIDE(void, CBDisk, accumulate, x, y, z, mass); + } + + void reset_coefs(void) override + { + PYBIND11_OVERRIDE(void, CBDisk, reset_coefs,); + } + + void make_coefs(void) override + { + PYBIND11_OVERRIDE(void, CBDisk, make_coefs,); + } + + }; + + class PySlab : public Slab { protected: @@ -1064,9 +1133,9 @@ void BasisFactoryClasses(py::module &m) Set the coordinate system for force evaluations. The natural coordinates for the basis class are the default; spherical coordinates for SphericalSL and Bessel, cylindrical coordinates for - Cylindrical and FlatDisk, and Cartesian coordinates for the Slab - and Cube. This member function can be used to override the default. - The available coorindates are: 'spherical', 'cylindrical', + Cylindrical, FlatDisk, and CBDisk, and Cartesian coordinates for the + Slab and Cube. This member function can be used to override the + default. The available coorindates are: 'spherical', 'cylindrical', 'cartesian'. Parameters @@ -1578,6 +1647,94 @@ void BasisFactoryClasses(py::module &m) )", py::arg("cachefile")); + py::class_, PyCBDisk, BasisClasses::BiorthBasis>(m, "CBDisk") + .def(py::init(), + R"( + Create a Clutton-Brock 2d disk basis + + Parameters + ---------- + YAMLstring : str + The YAML configuration for the razor-thin EOF basis. The default + parameters will give an exponential disk with scale length of + 0.01 units. Set the disk scale length using the 'scale' parameter. + + Returns + ------- + CBDisk + the new instance + )", py::arg("YAMLstring")) + .def("getBasis", &BasisClasses::CBDisk::getBasis, + R"( + Evaluate the potential-density basis functions + + Returned functions will linearly spaced 2d-grid for inspection. The + min/max radii are given in log_10 units. The structure is a two-grid of + dimension mmax by nmax each pointing to a dictionary of 1-d arrays + ('potential', 'density', 'rforce') of dimension numr. + + Parameters + ---------- + logxmin : float, default=-4.0 + the minimum radius in log10 scaled units + logxmax : float, default=-1.0 + the maximum radius in log10 scaled units + numr : int + the number of output evaluations + + Returns + ------- + list(dict{str: numpy.ndarray}) + list of lists of dictionaries in harmonic and radial pointing to + density and potential basis functions + )", + py::arg("logxmin")=-4.0, + py::arg("logxmax")=-1.0, + py::arg("numr")=400) + // The following member needs to be a lambda capture because + // orthoCheck is not in the base class and needs to have different + // parameters depending on the basis type. Here, the quadrature + // is determined by the scale of the meridional grid. + .def("orthoCheck", [](BasisClasses::CBDisk& A) + { + return A.orthoCheck(); + }, + R"( + Check orthgonality of basis functions by quadrature + + Inner-product matrix of Sturm-Liouville solutions indexed by + harmonic order used to assess fidelity. + + Parameters + ---------- + knots : int, default=40 + Number of quadrature knots + + Returns + ------- + list(numpy.ndarray) + list of numpy.ndarrays from [0, ... , Mmax] + )" + ) + .def_static("cacheInfo", [](std::string cachefile) + { + return BasisClasses::CBDisk::cacheInfo(cachefile); + }, + R"( + Report the parameters in a basis cache file and return a dictionary + + Parameters + ---------- + cachefile : str + name of cache file + + Returns + ------- + out : dict({tag: value}) + cache parameters + )", + py::arg("cachefile")); + py::class_, PySlab, BasisClasses::BiorthBasis>(m, "Slab") .def(py::init(), R"( diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index aa58603b9..21408b3ba 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -597,6 +597,84 @@ void CoefficientClasses(py::module &m) { }; + class PyTrajectoryData : public TrajectoryData + { + protected: + void readNativeCoefs(const std::string& file, int stride, double tmin, double tmax) override { + PYBIND11_OVERRIDE(void, TrajectoryData, readNativeCoefs, file, stride, tmin, tmax); + } + + std::string getYAML() override { + PYBIND11_OVERRIDE(std::string, TrajectoryData, getYAML,); + } + + void WriteH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(void, TrajectoryData, WriteH5Params, file); + } + + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { + PYBIND11_OVERRIDE(unsigned, TrajectoryData, WriteH5Times, group, count); + } + + public: + // Inherit the constructors + using TrajectoryData::TrajectoryData; + + Eigen::VectorXcd& getData(double time) override { + PYBIND11_OVERRIDE(Eigen::VectorXcd&, TrajectoryData, getData, time); + } + + void setData(double time, Eigen::VectorXcd& array) override { + PYBIND11_OVERRIDE(void, TrajectoryData, setData, time, array); + } + + std::shared_ptr getCoefStruct(double time) override { + PYBIND11_OVERRIDE(std::shared_ptr, TrajectoryData, getCoefStruct, + time); + } + + std::vector Times() override { + PYBIND11_OVERRIDE(std::vector, TrajectoryData, Times,); + } + + void WriteH5Coefs(const std::string& prefix) override { + PYBIND11_OVERRIDE(void, TrajectoryData, WriteH5Coefs, prefix); + } + + void ExtendH5Coefs(const std::string& prefix) override { + PYBIND11_OVERRIDE(void, TrajectoryData, ExtendH5Coefs, prefix); + } + + Eigen::MatrixXd& Power(int min, int max) override { + PYBIND11_OVERRIDE(Eigen::MatrixXd&, TrajectoryData, Power, min, max); + } + + bool CompareStanzas(std::shared_ptr check) override { + PYBIND11_OVERRIDE(bool, TrajectoryData, CompareStanzas, check); + } + + void clear() override { + PYBIND11_OVERRIDE(void, TrajectoryData, clear,); + } + + void add(CoefStrPtr coef) override { + PYBIND11_OVERRIDE(void, TrajectoryData, add, coef); + } + + std::vector makeKeys(Key k) override { + PYBIND11_OVERRIDE(std::vector, TrajectoryData, makeKeys, k); + } + + std::shared_ptr deepcopy() override { + PYBIND11_OVERRIDE(std::shared_ptr, TrajectoryData, deepcopy,); + } + + void zerodata() override { + PYBIND11_OVERRIDE(void, TrajectoryData, zerodata,); + } + + }; + py::class_, PyCoefStruct> (m, "CoefStruct") .def(py::init<>(), @@ -1720,6 +1798,76 @@ void CoefficientClasses(py::module &m) { R"( Return a 2-dimensional ndarray indexed by column and time + Returns + ------- + numpy.ndarray + 2-dimensional numpy array containing the data table + )"); + + py::class_, PyTrajectoryData, CoefClasses::Coefs> + (m, "TrajectoryData", "Container for trajectory/orbit data") + .def(py::init(), + R"( + Construct a null TrajectoryData object + + Parameters + ---------- + verbose : bool + display verbose information. + + Returns + ------- + TrajectoryData instance + )", py::arg("verbose")=true) + .def(py::init(), + R"( + Construct a TrajectoryData object from a data file + + Parameters + ---------- + type : str + ascii table data file + + Returns + ------- + TrajectoryData instance + )") + .def(py::init(), + R"( + Construct a TrajectoryData object from a data file + + Parameters + ---------- + type : str + ascii table data file + verbose : bool + display verbose information. + + Returns + ------- + TrajectoryData instance + )", py::arg("filename"), py::arg("verbose")=true) + .def(py::init&, std::vector&, bool>(), + R"( + Construct a TrajectoryData object from data arrays + + Parameters + ---------- + time : ndarray + time data + array : ndarray + data columns + verbose : bool + display verbose information. + + Returns + ------- + TrajectoryData instance + )", py::arg("time"), py::arg("array"), py::arg("verbose")=true) + .def("getAllCoefs", &CoefClasses::TrajectoryData::getAllCoefs, + R"( + Return a 2-dimensional ndarray indexed by column and time + Returns ------- numpy.ndarray diff --git a/pyEXP/EDMDWrappers.cc b/pyEXP/EDMDWrappers.cc index 6a7f25539..8e79dc5d1 100644 --- a/pyEXP/EDMDWrappers.cc +++ b/pyEXP/EDMDWrappers.cc @@ -129,7 +129,8 @@ void EDMDtoolkitClasses(py::module &m) { using namespace MSSA; - py::class_> f(m, "Koopman"); + py::class_> + f(m, "Koopman"); f.def(py::init(), R"( @@ -292,4 +293,5 @@ void EDMDtoolkitClasses(py::module &m) { -------- Koopman )"); + } diff --git a/pyEXP/FieldWrappers.cc b/pyEXP/FieldWrappers.cc index 165e58daf..9cfca7e9d 100644 --- a/pyEXP/FieldWrappers.cc +++ b/pyEXP/FieldWrappers.cc @@ -304,6 +304,32 @@ void FieldGeneratorClasses(py::module &m) { py::arg("projection"), py::arg("center") = std::vector(3, 0.0)); + f.def("histo1dlog", &Field::FieldGenerator::histo1dlog, + R"( + Make a 1d density histogram (array) in spherical projection with + logarithmic scaling + + Parameters + ---------- + reader : ParticleReader + particle reader instance + rmin : float + minimum radius of the histogram (>0) + rmax : float + linear extent of the histogram + nbins : int + number of bins + center : list(float, float, float), default=[0, 0, 0] + origin for computing the histogram + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray) + the evaluation grid and computed 1d histogram + )", + py::arg("reader"), py::arg("rmin"), py::arg("rmax"), py::arg("nbins"), + py::arg("center") = std::vector(3, 0.0)); + f.def("file_lines", &Field::FieldGenerator::file_lines, R"( Write field arrays to files using the supplied string prefix. diff --git a/pyEXP/MSSAWrappers.cc b/pyEXP/MSSAWrappers.cc index 2b158ff2c..aa5ef5c57 100644 --- a/pyEXP/MSSAWrappers.cc +++ b/pyEXP/MSSAWrappers.cc @@ -541,6 +541,49 @@ void MSSAtoolkitClasses(py::module &m) { power value )"); + + f.def("getKoopmanModes", &expMSSA::getKoopmanModes, + R"( + Compute the Koopman mode estimate from the right-singular vectors + + Uses eDMD to estimate the modes + + Parameters + ---------- + tol : double + singular value truncation level + window: int + Smoothing between serialized channels (0 for no smoothing) + debug : bool + flag indicating whether to print debug information + + Notes + ----- + Use getReconstructedKoopman() to copy the reconstruction for a + particular mode back to the coefficient db + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray) + vector of eigenvalues and modes + )", py::arg("tol")=1.0e-12, py::arg("window")=0, py::arg("debug")=false); + + f.def("getReconstructedKoopman", &expMSSA::getReconstructedKoopman, + R"( + Reconstruct the coefficients for a particular Koopman mode + + Parameters + ---------- + mode: int + The index of the mode to be reconstructed + + Returns + ------- + dict({id: Coefs},...) + reconstructed time series in the original coefficient form + + )", py::arg("mode")); + f.def("getRC", &expMSSA::getRC, R"( Access the detrended reconstructed channel series by internal key diff --git a/src/CBDisk.H b/src/CBDisk.H new file mode 100644 index 000000000..dcf68c260 --- /dev/null +++ b/src/CBDisk.H @@ -0,0 +1,108 @@ +#ifndef _CBDisk_H +#define _CBDisk_H + +#include +#include + +#include + +/** Reference implemenation of the Clutton-Brock disk basis + + Parameters + + \param mmax is the maximum azimuthal order + + \param Lmax alias for mmax + + \param Mmax alias for mmax + + \param scale radius for coordinate scaling. By default, the core radius of the Clutton-Brock density is 1. +*/ +class CBDisk : public PolarBasis +{ + +private: + + //@{ + //! 2D Clutton-Brock basis + double phif(const int n, const int m, const double r); + double dphi(const int n, const int m, const double r); + + double potl(const int n, const int m, const double r); + double dpot(const int n, const int m, const double r); + double dens(const int n, const int m, const double r); + + double norm(const int n, const int m); + + void potl(const int m, const double r, Eigen::VectorXd& a); + void dens(const int m, const double r, Eigen::VectorXd& a); + //@} + + //! Required maximum table radious (no limit here) + virtual double getRtable() { return std::numeric_limits::max(); } + + //! Sanity check + void orthoCheck(); + + void initialize(void); + + void get_dpotl(double r, double z, + Eigen::MatrixXd& p, Eigen::MatrixXd& dpr, Eigen::MatrixXd& dpz, int tid); + + void get_potl(double r, double z, Eigen::MatrixXd& p, int tid); + + //! Background instance + using Disk2d = std::shared_ptr; + Disk2d disk; + + //! Set background from YAML + void setBackground(); + + //! Background evaluation + virtual std::tuple + background(double r, double z) + { return {disk->pot(r), -disk->dpot(r), 0.0}; } + + void get_dens(double r, double z, Eigen::MatrixXd& p, int tid); + + void get_potl_dens(double r, double z, + Eigen::MatrixXd& p, Eigen::MatrixXd& d, int tid); + + //@{ + //! Parameters + double scale; + int mmax; + string model; + //@} + + //! Potential, force, and density scale factors + double fac1, fac2; + + //! Valid keys for YAML configurations + static const std::set valid_keys; + +protected: + + //! Clutton-Brock potential + void get_pot(Eigen::MatrixXd& Vc, Eigen::MatrixXd& Vs, double r, double z) + { + get_potl(r, z, Vc, 0); + Vs = Vc; + } + +public: + // Global parameters + /** Constructor + @param c0 is the instantiating caller (Component) + @param conf passes any parameters to basis instance + @param scale is the radius for coordinate scaling + */ + CBDisk(Component* c0, const YAML::Node& conf, MixtureBasis* m=0); + + //! Destructor + virtual ~CBDisk(); +}; + +#endif + + diff --git a/src/CBDisk.cc b/src/CBDisk.cc new file mode 100644 index 000000000..193b4ab79 --- /dev/null +++ b/src/CBDisk.cc @@ -0,0 +1,475 @@ +#include + +#include + +#include "expand.H" + +#include +#include +#include + +const std::set +CBDisk::valid_keys = { + "mmax", + "scale", + "diskconf", + "background" +}; + +CBDisk::CBDisk(Component* c0, const YAML::Node& conf, MixtureBasis* m) : + PolarBasis(c0, conf, m) +{ + id = "CBDisk"; + is_flat = true; + + // Radial scale factor/softening radius + // + scale = 1.0; + + // Get initialization info + // + initialize(); + + if (myid==0) { + std::string sep("---- "); + std::cout << "---- CBDisk parameters: " + << std::endl << sep << "scale=" << scale + << std::endl << sep << "lmax=" << Lmax + << std::endl << sep << "nmax=" << nmax + << std::endl << sep << "NO_L0=" << std::boolalpha << NO_M0 + << std::endl << sep << "NO_L1=" << std::boolalpha << NO_M1 + << std::endl << sep << "EVEN_M=" << std::boolalpha << EVEN_M + << std::endl << sep << "M0_ONLY=" << std::boolalpha << M0_only + << std::endl << sep << "selfgrav=" << std::boolalpha << self_consistent + << std::endl; + } + + setup(); + + // Potential, force, and density scaling + // + fac1 = pow(scale, -0.5); + fac2 = pow(scale, -1.5); + + if (myid==0) orthoCheck(); +} + + +void CBDisk::initialize() +{ + // Remove matched keys + // + for (auto v : valid_keys) current_keys.erase(v); + + // Assign values from YAML + // + try { + if (conf["mmax"]) mmax = conf["mmax"].as(); + if (conf["Lmax"]) mmax = conf["Lmax"].as(); + if (conf["Mmax"]) mmax = conf["Mmax"].as(); + if (conf["scale"]) scale = conf["scale"].as(); + + Lmax = Mmax = mmax; // Override base-class values + } + catch (YAML::Exception & error) { + if (myid==0) std::cout << "Error parsing parameters in CBDisk: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("CBDisk::initialize: error in parsing YAML"); + } + + // Set background model + if (conf["background"]) setBackground(); +} + +CBDisk::~CBDisk(void) +{ + // NADA +} + + +void CBDisk::setBackground() +{ + try { + + YAML::Node Params = conf["background"]; + + std::string name = Params["name"].as(); + auto params = Params["parameters"]; + + // Convert ID string to lower case + // + std::transform(name.begin(), name.end(), name.begin(), + [](unsigned char c){ return std::tolower(c); }); + + if (name.find("kuzmin") != std::string::npos) { + if (myid==0) { + std::cout << "---- CBDisk: Making a Kuzmin disk"; + if (params) std::cout << " with" << std::endl << Params["parameters"]; + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + if (name.find("mestel") != std::string::npos) { + if (myid==0) { + std::cout << "---- CBDisk: Making a finite Mestel disk"; + if (params) std::cout << " with" << std::endl << Params["parameters"]; + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + if (name.find("zang") != std::string::npos) { + if (myid==0) { + std::cout << "---- CBDisk: Making a double-tapered Zang"; + YAML::Emitter out; + out << YAML::Flow << params; + if (params) std::cout << " with " << out.c_str(); + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + if (name.find("expon") != std::string::npos) { + if (myid==0) { + std::cout << "---- CBDisk: Making an Exponential disk"; + if (params) std::cout << " with" << std::endl << params; + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + // Default if nothing else matches + if (myid==0) + std::cout << "---- CBDisk: Making an Exponential disk [Default]" + << std::endl; + disk = std::make_shared(params); + return; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in CBDisk::setBackground: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf["background"] << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("CBDisk::setBackground: error parsing YAML config"); + } +} + + +void CBDisk::get_dpotl(double r, double z, + Eigen::MatrixXd& pot, + Eigen::MatrixXd& dpr, + Eigen::MatrixXd& dpz, int tid) +{ + pot.resize(mmax+1, nmax); // Resize the return arrays + dpr.resize(mmax+1, nmax); + dpz.resize(mmax+1, nmax); + + dpz.setZero(); + + double R = r/scale; + + for (int m=0; m<=mmax; m++) { // Pack the array + for (int j=0; j0) ret -= 2.0*phif(n-1, m+1, r); + if (n>1) ret += phif(n-2, m+1, r); + return -r*ret; +} + + +// By recurrance relation +// +void CBDisk::potl(const int m, const double r, Eigen::VectorXd& a) +{ + a.resize(nmax); + + double pfac = pow(r, m); + + double r2 = r*r; + double fac = 1.0/(1.0 + r2); + double cur = sqrt(fac); + + for (int mm=1; mm<=m; mm++) cur *= fac*(2*mm - 1); + + a(0) = pfac*cur; + + if (nmax>0) { + + double curl1 = cur; + double curl2; + + fac *= r2 - 1.0; + cur *= fac*(2*m+1); + + a(1) = pfac*cur; + + for (int nn=2; nn=2) + return f*pow(r, m)*(phif(n, m+1, r) - phif(n-2, m+1, r)); + else + return f*pow(r, m)*phif(n, m+1, r); +} + + +void CBDisk::dens(int mm, double r, Eigen::VectorXd& a) +{ + a.resize(nmax); + + double pfac = pow(r, (double)mm+1.0e-20); + + int m = mm + 1; + double r2 = r*r; + double fac = 1.0/(1.0 + r2); + double cur = sqrt(fac); + + for (int M=1; M<=m; M++) cur *= fac*(2*M - 1); + + a(0) = pfac*cur; + + if (nmax>0) { + + double curl1 = cur; + double curl2; + + fac *= r2 - 1.0; + cur *= fac*(2*m+1); + + a(1) = pfac*cur; + + for (int nn=2; nn1; nn--) + a(nn) -= a(nn-2); + } + + for (int n=0; n ortho(mmax+1); + Eigen::VectorXd worst(mmax+1); + Eigen::MatrixXd vpot(mmax+1, nmax), vden(mmax+1, nmax); + + double Rmax = scale*100.0; + + worst.setZero(); + + // Allocate result matrices + // + for (auto & v : ortho) { + v.resize(nmax, nmax); + v.setZero(); + } + + for (int i=0; i worst(m)) worst(m) = fabs(test); + } + } + + if (worst(m) > tol) { + std::ostringstream sout; + sout << "CBDisk.ortho." << m ; + std::ofstream out(sout.str()); + if (out) out << ortho[m] << std::endl; + } + + } + + std::cout << hh << std::endl + << hh << "Orthogonality check" << std::endl + << hh << std::endl + << hh << std::setw(6) << "m" << std::setw(16) << "worst" << std::endl + << hh << std::setw(6) << "-" << std::setw(16) << "-----" << std::endl; + + for (int m=0; m<=mmax; m++) + std::cout << hh << std::setw(6) << m << std::setw(16) << worst(m) << std::endl; + std::cout << std::endl; + + const int numr = 10; + double dR = 3.0*scale/numr, dH = 0.01*scale; + + std::ofstream out("CBDisk.ortho.gradient"); + for (int m=0; m<=mmax; m++) { + for (int i=0; isecond->pos, tp->second->vel); + tp->second->acc[j] += val - acc[j]; + } + + //! Add to accerlation (by component) + inline void AddAccExt(int i, int j, double val) { PartMap::iterator tp = particles.find(i); if (tp == particles.end()) { throw BadIndexException(i, particles.size(), __FILE__, __LINE__); @@ -868,6 +890,17 @@ public: //! Add to acceleration (by array) inline void AddAcc(int i, double *val) { + PartMap::iterator tp = particles.find(i); + if (tp == particles.end()) { + throw BadIndexException(i, particles.size(), __FILE__, __LINE__); + } + auto acc = getPseudoAccel(tp->second->pos, tp->second->vel); + for (int k=0; k<3; k++) + tp->second->acc[k] += val[k] - acc[k]; + } + + //! Add to acceleration (by array) + inline void AddAccExt(int i, double *val) { PartMap::iterator tp = particles.find(i); if (tp == particles.end()) { throw BadIndexException(i, particles.size(), __FILE__, __LINE__); @@ -877,6 +910,17 @@ public: //! Add to accerlation (by vector) inline void AddAcc(int i, vector& val) { + PartMap::iterator tp = particles.find(i); + if (tp == particles.end()) { + throw BadIndexException(i, particles.size(), __FILE__, __LINE__); + } + auto acc = getPseudoAccel(tp->second->pos, tp->second->vel); + for (int k=0; k<3; k++) + tp->second->acc[k] += val[k] - acc[k]; + } + + //! Add to accerlation (by vector) + inline void AddAccExt(int i, vector& val) { PartMap::iterator tp = particles.find(i); if (tp == particles.end()) { throw BadIndexException(i, particles.size(), __FILE__, __LINE__); diff --git a/src/Component.cc b/src/Component.cc index d696a43da..b52f30c30 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -46,6 +47,7 @@ const std::set Component::valid_keys_parm = "EJ", "nEJkeep", "nEJwant", + "nEJaccel", "EJkinE", "EJext", "EJdiag", @@ -183,6 +185,7 @@ Component::Component(YAML::Node& CONF) EJ = 0; nEJkeep = 100; nEJwant = 500; + nEJaccel = 0; EJkinE = true; EJext = false; EJdiag = false; @@ -319,6 +322,7 @@ void Component::set_default_values() if (!cconf["EJ"]) cconf["EJ"] = EJ; if (!cconf["nEJkeep"]) cconf["nEJkeep"] = nEJkeep; if (!cconf["nEJwant"]) cconf["nEJwant"] = nEJwant; + if (!cconf["nEJaccel"]) cconf["nEJaccel"] = nEJaccel; if (!cconf["EJkinE"]) cconf["EJkinE"] = EJkinE; if (!cconf["EJext"]) cconf["EJext"] = EJext; if (!cconf["EJdiag"]) cconf["EJdiag"] = EJdiag; @@ -682,6 +686,7 @@ Component::Component(YAML::Node& CONF, istream *in, bool SPL) : conf(CONF) EJ = 0; nEJkeep = 100; nEJwant = 500; + nEJaccel = 0; EJkinE = true; EJext = false; EJdiag = false; @@ -794,6 +799,7 @@ void Component::configure(void) std::cout << "Component: eEJ0 is no longer used, Ecurr is computed from the bodies using the expansion directly" << std::endl; if (cconf["nEJkeep" ]) nEJkeep = cconf["nEJkeep" ].as(); if (cconf["nEJwant" ]) nEJwant = cconf["nEJwant" ].as(); + if (cconf["nEJaccel"]) nEJaccel = cconf["nEJaccel"].as(); if (cconf["EJx0" ]) EJx0 = cconf["EJx0" ].as(); if (cconf["EJy0" ]) EJy0 = cconf["EJy0" ].as(); if (cconf["EJz0" ]) EJz0 = cconf["EJz0" ].as(); @@ -883,6 +889,9 @@ void Component::configure(void) else if ( !id.compare("flatdisk") ) { force = new FlatDisk(this, fconf); } + else if ( !id.compare("CBDisk") ) { + force = new CBDisk(this, fconf); + } else if ( !id.compare("direct") ) { force = new Direct(this, fconf); } @@ -925,6 +934,7 @@ void Component::initialize(void) for (int k=0; k<3; k++) { com[k] = center[k] = cov[k] = coa[k] = 0.0; com0[k] = cov0[k] = acc0[k] = angmom[k] = 0.0; + accel[k] = 0.0; } if (com_system) { @@ -1109,6 +1119,7 @@ void Component::initialize(void) if (EJdiag) cout << "Process " << myid << ": about to create Orient with" << " nkeep=" << nEJkeep << " nwant=" << nEJwant + << " naccel=" << nEJaccel << " EJkinE=" << EJkinE << " EJext=" << EJext; @@ -1135,10 +1146,12 @@ void Component::initialize(void) if (EJkinE) EJctl |= Orient::KE; if (EJext) EJctl |= Orient::EXTERNAL; - orient = new Orient(nEJkeep, nEJwant, EJ, EJctl, EJlogfile, EJdT, EJdamp); + orient = new Orient(nEJkeep, nEJwant, nEJaccel, + EJ, EJctl, EJlogfile, EJdT, EJdamp); if (restart && (EJ & Orient::CENTER)) { Eigen::VectorXd::Map(¢er[0], 3) = orient->currentCenter(); + std::tie(accel, omega, domdt) = orient->currentAccel(); } else { if (EJlinear) orient -> set_linear(); if (not com_system) { @@ -3107,6 +3120,7 @@ void Component::fix_positions_cpu(unsigned mlevel) if ((EJ & Orient::CENTER) && !EJdryrun) { auto ctr = orient->currentCenter(); + std::tie(accel, omega, domdt) = orient->currentAccel(); bool ok = true; for (int i=0; i<3; i++) { if (std::isnan(ctr[i])) ok = false; @@ -3942,3 +3956,23 @@ void Component::AddPart(PartPtr p) nbodies = particles.size(); } +Eigen::Vector3d& Component::getPseudoAccel(double* pos, double* vel) +{ + pseudo.setZero(); + + if (EJdryrun) return pseudo; + + if (EJ & Orient::CENTER) { + pseudo += accel; + } + + if (EJ & Orient::AXIS) { + Eigen::Map P(pos); + Eigen::Map V(vel); + + // Coriolis + Euler + Centrifugal forces + pseudo += 2.0*omega.cross(V) + domdt.cross(P) + omega.cross(omega.cross(P)); + } + + return pseudo; +} diff --git a/src/ComponentContainer.cc b/src/ComponentContainer.cc index 72030cddf..f399a59d0 100644 --- a/src/ComponentContainer.cc +++ b/src/ComponentContainer.cc @@ -771,7 +771,7 @@ void ComponentContainer::compute_potential(unsigned mlevel) if (timing) itmr->second.start(); ext->set_multistep_level(mlevel); - if (use_cuda and not c->force->cudaAware()) { + if (use_cuda and not ext->cudaAware()) { #if HAVE_LIBCUDA==1 c->CudaToParticles(); #endif diff --git a/src/Cylinder.cc b/src/Cylinder.cc index e355bcdf9..80eecc839 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -106,7 +106,8 @@ Cylinder::valid_keys = { "self_consistent", "playback", "coefCompute", - "coefMaster" + "coefMaster", + "dumpbasis" }; Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : @@ -552,6 +553,8 @@ void Cylinder::initialize() } } + if (conf["dumpbasis"]) dump_basis = true; + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing Cylinder parameters: " diff --git a/src/FlatDisk.H b/src/FlatDisk.H index 64be763a4..81c7b4eb7 100644 --- a/src/FlatDisk.H +++ b/src/FlatDisk.H @@ -16,8 +16,7 @@ typedef std::shared_ptr CylPtr; /** Computes the potential, acceleration and density using an empirical function analysis for a flat Bessel basis and a target density - \param acyltbl is the scale length of target density (default is - 0.6) + Parameters \param rcylmin is the minimum value in the table for the radial basis functions (default is 0.0) @@ -28,22 +27,33 @@ typedef std::shared_ptr CylPtr; \param scale is the expansion factor to get the physical scale from the internal scale (default is 0.01) + \param mmax is the maximum azimuthal order + \param numx is the number of grid points in the scaled radial direction \param numy is the number of grid points in the scaled vertical direction - \param numr is the number radial coordinate knots in the table - (default is 1000) - - \param nfid is the number of basis elements used in the EOF computation - \param knots is the number of quadrature points for EmpCyl2d \param logr scales the EmpCyl2d grid logarithmically - \param model name for EmpCyl2d (default: expon) + \param model name for conditioning the new basis using EmpCyl2d (default: expon) + + \param biorth is the biorthogonal basis set used by EmpCyl2d for conditioning (default: bess) + + \param background sets a fixed monopole potential and force for the model. Used for stability analysis of disks with active and inactive tapers, following Zang + + \param nmaxfid is the maximum radial order for the basis used by EmpCyl2d for conditioning + + \param numr is the number of grid points in the radial direction for the numerical basis in EmpCyl2d + + \param NQDHT is the number of grid points for the numerical Hankel transform - \param biorth set for EmpCyl2d (default: bess) + \param diskconf is a YAML configuration for the target disk model in EmpCyl2d + + \param cachename is the name of the cache file for the newly constructed basis + + \param dumpbasis provides the user with a ascii table of the basis potential-density pairs */ class FlatDisk : public PolarBasis { @@ -65,7 +75,7 @@ private: virtual void initialize_cuda() { sampT = floor(sqrt(component->CurTotal())); - ortho->initialize_cuda(cuInterpArray, tex); + ortho->initialize_cuda(disk, cuInterpArray, tex); } virtual cudaMappingConstants getCudaMappingConstants() @@ -74,14 +84,18 @@ private: } #endif + //! Background instance + using Disk2d = std::shared_ptr; + Disk2d disk; + + //! Set background from YAML + void setBackground(); + //! Background evaluation virtual std::tuple - get_pot_background(double r, double z) - { - auto [p, dr, dz] = ortho->background(r, z); - return {p, -dr, -dz}; - } - + background(double r, double z) + { return {disk->pot(r), -disk->dpot(r), 0.0}; } + void get_dens(double r, double z, Eigen::MatrixXd& p, int tid); void get_potl_dens(double r, double z, @@ -99,6 +113,8 @@ private: //! Valid keys for YAML configurations static const std::set valid_keys; + bool dump_basis; + protected: //! Use BiorthCyl to evaluate the basis at r, z diff --git a/src/FlatDisk.cc b/src/FlatDisk.cc index faf626fa8..706bd95ba 100644 --- a/src/FlatDisk.cc +++ b/src/FlatDisk.cc @@ -8,38 +8,46 @@ const std::set FlatDisk::valid_keys = { - "nmaxfid", - "numr", "rcylmin", "rcylmax", "mmax", "numx", "numy", - "NQDHT", "knots", "logr", "model", "biorth", + "background", + "nmaxfid", + "numr", + "NQDHT", "diskconf", - "cachename" + "cachename", + "dumpbasis" }; FlatDisk::FlatDisk(Component* c0, const YAML::Node& conf, MixtureBasis* m) : PolarBasis(c0, conf, m) { id = "FlatDisk EOF"; - // Defaults + + // Defaults + // knots = 40; numr = 2000; rcylmin = 0.0; rcylmax = 10.0; logr = false; is_flat = true; + dump_basis = false; - // Get initialization info + // Get initialization info + // initialize(); - id += "[" + model + "][" + biorth + "]"; + // Append info to basis ID + // + id += "[" + ortho->getModelName() + "][" + ortho->getBiorthName() + "]"; if (myid==0) { std::string sep("---- "); @@ -81,6 +89,8 @@ void FlatDisk::initialize() if (conf["biorth"]) biorth = conf["biorth"].as(); if (conf["mmax"]) Lmax = Mmax = mmax; // Override base-class values + + if (conf["dumpbasis"]) dump_basis = conf["dumpbasis"].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in FlatDisk: " @@ -96,6 +106,10 @@ void FlatDisk::initialize() // Create the BiorthCyl instance ortho = std::make_shared(conf); + if (dump_basis) ortho->dump_basis(runtag); + + // Set background model + if (conf["background"]) setBackground(); } FlatDisk::~FlatDisk(void) @@ -104,6 +118,83 @@ FlatDisk::~FlatDisk(void) } +void FlatDisk::setBackground() +{ + try { + + YAML::Node Params = conf["background"]; + + std::string name = Params["name"].as(); + auto params = Params["parameters"]; + + // Convert ID string to lower case + // + std::transform(name.begin(), name.end(), name.begin(), + [](unsigned char c){ return std::tolower(c); }); + + if (name.find("kuzmin") != std::string::npos) { + if (myid==0) { + std::cout << "---- FlatDisk: Making a Kuzmin disk"; + if (params) std::cout << " with" << std::endl << Params["parameters"]; + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + if (name.find("mestel") != std::string::npos) { + if (myid==0) { + std::cout << "---- FlatDisk: Making a finite Mestel disk"; + if (params) std::cout << " with" << std::endl << Params["parameters"]; + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + if (name.find("zang") != std::string::npos) { + if (myid==0) { + std::cout << "---- FlatDisk: Making a double-tapered Zang"; + YAML::Emitter out; + out << YAML::Flow << params; + if (params) std::cout << " with " << out.c_str(); + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + if (name.find("expon") != std::string::npos) { + if (myid==0) { + std::cout << "---- FlatDisk: Making an Exponential disk"; + if (params) std::cout << " with" << std::endl << params; + std::cout << std::endl; + } + disk = std::make_shared(params); + return; + } + + // Default if nothing else matches + if (myid==0) + std::cout << "---- FlatDisk: Making an Exponential disk [Default]" + << std::endl; + disk = std::make_shared(params); + return; + } + catch (YAML::Exception & error) { + if (myid==0) + std::cout << "Error parsing parameters in FlatDisk::setBackground: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf["background"] << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("FlatDisk::setBackground: error parsing YAML config"); + } +} + + void FlatDisk::get_dpotl(double r, double z, Eigen::MatrixXd& p, Eigen::MatrixXd& dpr, diff --git a/src/Orient.H b/src/Orient.H index b2bb1fa91..bb106bd51 100644 --- a/src/Orient.H +++ b/src/Orient.H @@ -8,6 +8,7 @@ #include #include +#include #include @@ -97,6 +98,9 @@ private: double damp; Eigen::Matrix3d body, orig; Eigen::Vector3d axis, center, axis1, center1, center0, cenvel0; + Eigen::Vector3d pseudo = Eigen::Vector3d::Zero(); + Eigen::Vector3d omega = Eigen::Vector3d::Zero(); + Eigen::Vector3d domdt = Eigen::Vector3d::Zero(); double sumX, sumX2; Eigen::Vector3d sumY, sumY2, sumXY, slope, intercept; double sigA, sigC, sigCz; @@ -110,8 +114,10 @@ private: string logfile; bool linear; - vector pos, psa, vel; + std::vector pos, psa, vel; + std::shared_ptr accel; + void accumulate_cpu(double time, Component* c); #if HAVE_LIBCUDA==1 void accumulate_gpu(double time, Component* c); @@ -125,7 +131,7 @@ public: enum ControlFlags {DIAG=1, KE=2, EXTERNAL=4}; //! Constructor - Orient(int number_to_keep, int target, + Orient(int number_to_keep, int target, int Naccel, unsigned orient_flags, unsigned control_flags, string logfile, double dt=0.0, double damping=1.0); @@ -167,6 +173,16 @@ public: //! Return current center Eigen::Vector3d& currentCenter(void) {return center;}; + //! Return current pseudo force, angular velocity, and angular + //! velocity derivative + const std::tuple + currentAccel(void) + { + if (accel) return (*accel)(); + pseudo.setZero(); // No acceleration or angular velocity + return {pseudo, pseudo, pseudo}; + } + //! Return variance of axis determination (linear least squares solution) double currentAxisVar(void) {return sigA;} diff --git a/src/Orient.cc b/src/Orient.cc index 676e40792..99597dc24 100644 --- a/src/Orient.cc +++ b/src/Orient.cc @@ -36,7 +36,7 @@ void EL3::debug() const } } -Orient::Orient(int n, int nwant, unsigned Oflg, unsigned Cflg, +Orient::Orient(int n, int nwant, int naccel, unsigned Oflg, unsigned Cflg, string Logfile, double dt, double damping) { keep = n; @@ -59,6 +59,13 @@ Orient::Orient(int n, int nwant, unsigned Oflg, unsigned Cflg, cenvel0.setZero(); axis .setZero(); + // Instantiate the pseudo-force helper + // class if naccel>0 + if (naccel) { + accel = std::make_shared + (naccel, oflags & CENTER, oflags & AXIS); + } + // Initialize last time to something // in the near infinite past lasttime = -std::numeric_limits::max(); @@ -70,8 +77,7 @@ Orient::Orient(int n, int nwant, unsigned Oflg, unsigned Cflg, // Set up identity body.setIdentity(); orig = body; - - // Check for previous state on + // check for previous state on // a restart int in_ok; std::vector in1(4), in2(4); @@ -169,6 +175,14 @@ Orient::Orient(int n, int nwant, unsigned Oflg, unsigned Cflg, if (static_cast(sumsC.size()) > keep) sumsC.pop_front(); } + bool allRead = true; + for (int i=0; i<3; i++) { + if (line.eof()) allRead = false; + for (int k; k<3; k++) line >> pseudo(k); + } + if (allRead) { + if (accel) accel->add(time, pseudo, axis1); + } } cout << " -- Orient: current log=" << logfile << " backup=" << backupfile << endl; @@ -244,12 +258,21 @@ Orient::Orient(int n, int nwant, unsigned Oflg, unsigned Cflg, << setw(15) << "| X-com(dif)" // 22 << setw(15) << "| Y-com(dif)" // 23 << setw(15) << "| Z-com(dif)" // 24 + << setw(15) << "| X-accel" // 25 + << setw(15) << "| Y-accel" // 26 + << setw(15) << "| Z-accel" // 27 + << setw(15) << "| Omega_X" // 28 + << setw(15) << "| Omega_Y" // 29 + << setw(15) << "| Omega_Z" // 30 + << setw(15) << "| dOmega/dt_X" // 31 + << setw(15) << "| dOmega/dt_Y" // 32 + << setw(15) << "| dOmega/dt_Z" // 33 << endl; out.fill('-'); int icnt = 1; out << "# " << setw(13) << icnt++; - for (int i=0; i<23; i++) out << "| " << setw(13) << icnt++; + for (int i=0; i<32; i++) out << "| " << setw(13) << icnt++; out << endl; out.close(); @@ -528,12 +551,16 @@ void Orient::accumulate(double time, Component *c) } if ((cflags & DIAG) && myid==0) { - cout << " Orient info [" << time << ", " << c->name << "]: " - << used << " particles used, Ecurr=" << Ecurr - << " Center=" - << center1[0] << ", " - << center1[1] << ", " - << center1[2] << endl; + std::cout << " Orient info [" << time << ", " << c->name << "]: " + << used << " particles used, Ecurr=" << Ecurr + << " Center=" + << center1[0] << ", " + << center1[1] << ", " + << center1[2] + << " Axis=" + << axis1[0] << ", " + << axis1[1] << ", " + << axis1[2] << std::endl; } if (static_cast(sumsA.size()) > keep + 1) { @@ -660,6 +687,12 @@ void Orient::accumulate(double time, Component *c) } else center = center1; + // Add to pseudo-acceleration estimator + // + if (accel) { + accel->add(time, center1, axis1); + } + // Increment initial center according // to user specified velocity center0 += cenvel0*dtime; @@ -732,6 +765,17 @@ void Orient::logEntry(double time, Component *c) // Columns 22 - 24 for (int k=0; k<3; k++) outl << setw(15) << c->com0[k]; + if (accel) std::tie(pseudo, omega, domdt) = (*accel)(); + + // Columns 25 - 27 + for (int k=0; k<3; k++) outl << setw(15) << pseudo[k]; + + // Columns 28 - 30 + for (int k=0; k<3; k++) outl << setw(15) << omega[k]; + + // Columns 31 - 33 + for (int k=0; k<3; k++) outl << setw(15) << domdt[k]; + outl << endl; } } diff --git a/src/OutputContainer.H b/src/OutputContainer.H index 9d512b32a..520611a79 100644 --- a/src/OutputContainer.H +++ b/src/OutputContainer.H @@ -32,7 +32,7 @@ public: void initialize(); //! Execute the all methods in the container - void Run(int nstep, int mstep=std::numeric_limits::max(), bool last=false); + void Run(int nstep, int mstep=std::numeric_limits::max(), bool final=false); }; #endif diff --git a/src/ParticleFerry.cc b/src/ParticleFerry.cc index 142ca2305..30aa3d7c0 100644 --- a/src/ParticleFerry.cc +++ b/src/ParticleFerry.cc @@ -126,13 +126,17 @@ void ParticleFerry::particlePack(PartPtr in, char* buffer) // std::vector iattrib // - memcpy(&buffer[pos], &in->iattrib[0], nimax*sizeof(int)); - pos += sizeof(int)*nimax; + if (nimax) { + memcpy(&buffer[pos], &in->iattrib[0], nimax*sizeof(int)); + pos += sizeof(int)*nimax; + } // std::vector dattrib // - memcpy(&buffer[pos], &in->dattrib[0], ndmax*sizeof(double)); - pos += sizeof(double)*ndmax; + if (ndmax) { + memcpy(&buffer[pos], &in->dattrib[0], ndmax*sizeof(double)); + pos += sizeof(double)*ndmax; + } // unsigned level // @@ -210,15 +214,19 @@ void ParticleFerry::particleUnpack(PartPtr out, char* buffer) // std::vector iattrib // - out->iattrib.resize(nimax); - memcpy(&out->iattrib[0], &buffer[pos], nimax*sizeof(int)); - pos += sizeof(int)*nimax; + if (nimax) { + out->iattrib.resize(nimax); + memcpy(&out->iattrib[0], &buffer[pos], nimax*sizeof(int)); + pos += sizeof(int)*nimax; + } // std::vector dattrib // - out->dattrib.resize(ndmax); - memcpy(&out->dattrib[0], &buffer[pos], ndmax*sizeof(double)); - pos += sizeof(double)*ndmax; + if (ndmax) { + out->dattrib.resize(ndmax); + memcpy(&out->dattrib[0], &buffer[pos], ndmax*sizeof(double)); + pos += sizeof(double)*ndmax; + } // unsigned level // diff --git a/src/PolarBasis.H b/src/PolarBasis.H index caab18f0d..bc835b4a0 100644 --- a/src/PolarBasis.H +++ b/src/PolarBasis.H @@ -380,7 +380,7 @@ protected: //! Background evaluation virtual std::tuple - get_pot_background(double r, double z) { return {0.0, 0.0, 0.0}; } + background(double r, double z) { return {0.0, 0.0, 0.0}; } //! Thread method for accerlation compuation virtual void * determine_acceleration_and_potential_thread(void * arg); diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index 67f2e454f..aae196aeb 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -401,8 +401,20 @@ void * PolarBasis::determine_coefficients_thread(void * arg) { // For biorthogonal density component and normalization // - constexpr double norm0 = 2.0*M_PI * 0.5*M_2_SQRTPI/M_SQRT2; - constexpr double norm1 = 2.0*M_PI * 0.5*M_2_SQRTPI; + constexpr double norm0_3d = 2.0*M_PI * 0.5*M_2_SQRTPI/M_SQRT2; + constexpr double norm1_3d = 2.0*M_PI * 0.5*M_2_SQRTPI; + + constexpr double norm0_2d = 1.0; + constexpr double norm1_2d = M_SQRT2; + + double norm0, norm1; + if (is_flat) { + norm0 = norm0_2d; + norm1 = norm1_2d; + } else { + norm0 = norm0_3d; + norm1 = norm1_3d; + } double r, r2, facL=1.0, fac1, fac2, phi, mass; double xx, yy, zz; @@ -1044,8 +1056,20 @@ void PolarBasis::multistep_update(int from, int to, Component *c, int i, int id) // For biorthogonal density component and normalization // - constexpr double norm0 = 2.0*M_PI * 0.5*M_2_SQRTPI/M_SQRT2; - constexpr double norm1 = 2.0*M_PI * 0.5*M_2_SQRTPI; + constexpr double norm0_3d = 2.0*M_PI * 0.5*M_2_SQRTPI/M_SQRT2; + constexpr double norm1_3d = 2.0*M_PI * 0.5*M_2_SQRTPI; + + constexpr double norm0_2d = 1.0; + constexpr double norm1_2d = M_SQRT2; + + double norm0, norm1; + if (is_flat) { + norm0 = norm0_2d; + norm1 = norm1_2d; + } else { + norm0 = norm0_3d; + norm1 = norm1_3d; + } double mass = c->Mass(i) * component->Adiabatic(); @@ -1323,8 +1347,20 @@ void PolarBasis::compute_multistep_coefficients() void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) { - constexpr double norm0 = 0.5*M_2_SQRTPI/M_SQRT2; - constexpr double norm1 = 0.5*M_2_SQRTPI; + constexpr double norm0_3d = 0.5*M_2_SQRTPI/M_SQRT2; + constexpr double norm1_3d = 0.5*M_2_SQRTPI; + + constexpr double norm0_2d = 1.0; + constexpr double norm1_2d = M_SQRT2; + + double norm0, norm1; + if (is_flat) { + norm0 = norm0_2d; + norm1 = norm1_2d; + } else { + norm0 = norm0_3d; + norm1 = norm1_3d; + } double r, r0=0.0, phi; double potr, potz, potl, potp, p, pc, drc, drs, dzc, dzs, ps, dfacp, facdp; @@ -1432,7 +1468,7 @@ void * PolarBasis::determine_acceleration_and_potential_thread(void * arg) if (not NO_M0) { if (M0_back) { - auto [p, drc, dzc] = get_pot_background(r, zz); + auto [p, drc, dzc] = background(r, zz); potl = mfac * p; potr = mfac * drc; diff --git a/src/cudaBiorthCyl.cu b/src/cudaBiorthCyl.cu index 218e70177..5a5790202 100644 --- a/src/cudaBiorthCyl.cu +++ b/src/cudaBiorthCyl.cu @@ -84,7 +84,8 @@ thrust::host_vector returnTestBioCyl1 } void BiorthCyl::initialize_cuda -(std::vector& cuArray, +(std::shared_ptr disk, + std::vector& cuArray, thrust::host_vector& tex ) { @@ -186,15 +187,19 @@ void BiorthCyl::initialize_cuda // Add background arrays // std::vector> tt(ndim2); - for (auto & v : tt) v.resize(numr); + for (auto & v : tt) { + v.resize(numr); + thrust::fill(v.begin(), v.end(), 0.0); + } - double dx0 = (xmax - xmin)/(numr - 1); + if (disk) { + double dx0 = (xmax - xmin)/(numr - 1); - for (int i=0; ipot(r); + tt[1][i] = disk->dpot(r); + } } // Allocate CUDA array in device memory (a one-dimension 'channel') diff --git a/src/cudaPolarBasis.cu b/src/cudaPolarBasis.cu index 0e51e34ce..ad6bb5349 100644 --- a/src/cudaPolarBasis.cu +++ b/src/cudaPolarBasis.cu @@ -884,13 +884,11 @@ __global__ void coefKernelPlr3 const int tid = blockDim.x * blockIdx.x + threadIdx.x; const int N = lohi.second - lohi.first; - // The inner product of the potential-density pair is 1/(2*pi). So - // the biorthgonal norm is 2*pi*density. The inner product of the - // trig functions is 2*pi. So the norm is 1/sqrt(2*pi). The total - // norm is therefore 2*pi/sqrt(2*pi) = sqrt(2*pi) = 2.5066... + // 2d basis are assumed to be normed on input upto the sqrt(2) for + // the sine and cosine angular terms. // - cuFP_t norm = 2.5066282746310007; - if (m) norm = 3.5449077018110322; + cuFP_t norm = 1.0; + if (m) norm = 1.4142135623730951; for (int n=0; n P, dArray I, // Normalization constants // - cuFP_t norm0 = 0.39894228040143270286; - cuFP_t norm1 = 0.56418958354775627928; + cuFP_t norm0 = 1.0; + cuFP_t norm1 = 1.4142135623730951; cuFP_t norm; int muse = mmax > mlim ? mlim : mmax; diff --git a/src/end.cc b/src/end.cc index 2f7c0afb9..008a03e74 100644 --- a/src/end.cc +++ b/src/end.cc @@ -15,9 +15,9 @@ void clean_up(void) { // Call for final output to files - output->Run(this_step, 0, true); + if (output) output->Run(this_step, 0, true); // Cache for restart - external->finish(); + if (external) external->finish(); MPI_Barrier(MPI_COMM_WORLD); diff --git a/src/user/CMakeLists.txt b/src/user/CMakeLists.txt index e340cce7e..1f7b4c85a 100644 --- a/src/user/CMakeLists.txt +++ b/src/user/CMakeLists.txt @@ -32,8 +32,8 @@ set(halo_SRC UserHalo.cc) set(mndisk_SRC UserMNdisk.cc) set(mwgala_SRC UserMW.cc) -if (ENABLE_CUDA) - set(cudatest_Src UserTestCuda.cc cudaUserTest.cu) +if(ENABLE_CUDA) + set(cudatest_SRC UserTestCuda.cc cudaUserTest.cu) endif() foreach(mlib ${USER_MODULES}) diff --git a/src/user/UserBar.cc b/src/user/UserBar.cc index 527732640..db3f6eeed 100644 --- a/src/user/UserBar.cc +++ b/src/user/UserBar.cc @@ -478,14 +478,14 @@ void * UserBar::determine_acceleration_and_potential_thread(void * arg) nn = pp * pow(rr/b5, 3.0)/(b5*b5); } - cC->AddAcc(i, 0, - ffac*( 2.0*( xx*cos2p + yy*sin2p)*fac - 5.0*nn*xx ) ); + cC->AddAccExt(i, 0, + ffac*( 2.0*( xx*cos2p + yy*sin2p)*fac - 5.0*nn*xx ) ); - cC->AddAcc(i, 1, - ffac*( 2.0*(-yy*cos2p + xx*sin2p)*fac - 5.0*nn*yy ) ); + cC->AddAccExt(i, 1, + ffac*( 2.0*(-yy*cos2p + xx*sin2p)*fac - 5.0*nn*yy ) ); - cC->AddAcc(i, 2, - ffac*( -5.0*nn*zz ) ); + cC->AddAccExt(i, 2, + ffac*( -5.0*nn*zz ) ); cC->AddPotExt(i, -ffac*pp*fac ); diff --git a/src/user/UserDisk.cc b/src/user/UserDisk.cc index e80dd789a..1bee560ad 100644 --- a/src/user/UserDisk.cc +++ b/src/user/UserDisk.cc @@ -395,9 +395,9 @@ void * UserDisk::determine_acceleration_and_potential_thread(void * arg) getTable(rr, zz, pot, fr, fz); // Add acceleration by disk - cC->AddAcc(i, 0, amp * fr*xx/(rr+1.0e-10) ); - cC->AddAcc(i, 1, amp * fr*yy/(rr+1.0e-10) ); - cC->AddAcc(i, 2, amp * fz ); + cC->AddAccExt(i, 0, amp * fr*xx/(rr+1.0e-10) ); + cC->AddAccExt(i, 1, amp * fr*yy/(rr+1.0e-10) ); + cC->AddAccExt(i, 2, amp * fz ); // Add external potential cC->AddPotExt(i, pot); diff --git a/src/user/UserHalo.cc b/src/user/UserHalo.cc index 09d774e3a..22f670c57 100644 --- a/src/user/UserHalo.cc +++ b/src/user/UserHalo.cc @@ -128,7 +128,7 @@ void UserHalo::determine_acceleration_and_potential(void) exp_thread_fork(false); - print_timings("UserHalo: accleration timings"); + print_timings("UserHalo: acceleration timings"); } @@ -214,9 +214,9 @@ void * UserHalo::determine_acceleration_and_potential_thread(void * arg) */ // END DEBUG - // Add external accerlation + // Add external acceleration for (int k=0; k<3; k++) - cC->AddAcc(i, k, -dpot*pos[k]/(qq[k]*r) ); + cC->AddAccExt(i, k, -dpot*pos[k]/(qq[k]*r) ); // Add external potential cC->AddPotExt(i, pot ); diff --git a/src/user/UserLogPot.cc b/src/user/UserLogPot.cc index 7f34c3af2..cf17d6e55 100644 --- a/src/user/UserLogPot.cc +++ b/src/user/UserLogPot.cc @@ -105,11 +105,11 @@ void * UserLogPot::determine_acceleration_and_potential_thread(void * arg) zz = cC->Pos(i, 2); rr = R*R + xx*xx + yy*yy/(b*b) + zz*zz/(c*c); - cC->AddAcc(i, 0,-v2*xx/rr ); + cC->AddAccExt(i, 0,-v2*xx/rr ); - cC->AddAcc(i, 1, -v2*yy/(rr*b*b) ); + cC->AddAccExt(i, 1, -v2*yy/(rr*b*b) ); - cC->AddAcc(i, 2, -v2*zz/(rr*c*c) ); + cC->AddAccExt(i, 2, -v2*zz/(rr*c*c) ); cC->AddPotExt(i, 0.5*v2*log(rr) ); } diff --git a/src/user/UserMNdisk.cc b/src/user/UserMNdisk.cc index 7aa8b935f..f7ceca8ca 100644 --- a/src/user/UserMNdisk.cc +++ b/src/user/UserMNdisk.cc @@ -174,9 +174,9 @@ void * UserMNdisk::determine_acceleration_and_potential_thread(void * arg) double fz = -mass*zz*ab/(zb*dn*dn*dn); // Add acceleration by disk - cC->AddAcc(i, 0, amp * fr*xx/(rr+1.0e-10) ); - cC->AddAcc(i, 1, amp * fr*yy/(rr+1.0e-10) ); - cC->AddAcc(i, 2, amp * fz ); + cC->AddAccExt(i, 0, amp * fr*xx/(rr+1.0e-10) ); + cC->AddAccExt(i, 1, amp * fr*yy/(rr+1.0e-10) ); + cC->AddAccExt(i, 2, amp * fz ); // Add external potential cC->AddPotExt(i, pot); diff --git a/src/user/UserMW.cc b/src/user/UserMW.cc index 011c4b174..0c4c36dac 100644 --- a/src/user/UserMW.cc +++ b/src/user/UserMW.cc @@ -207,9 +207,9 @@ void * UserMW::determine_acceleration_and_potential_thread(void * arg) double zacc = az1 + az2 + az3 + az4; // Add acceleration to particle - cC->AddAcc(i, 0, xacc ); - cC->AddAcc(i, 1, yacc ); - cC->AddAcc(i, 2, zacc ); + cC->AddAccExt(i, 0, xacc ); + cC->AddAccExt(i, 1, yacc ); + cC->AddAccExt(i, 2, zacc ); // Add external potential cC->AddPotExt(i, pot); diff --git a/utils/Analysis/KDcyltest.cc b/utils/Analysis/KDcyltest.cc index b272ccdbb..edf29c8bb 100644 --- a/utils/Analysis/KDcyltest.cc +++ b/utils/Analysis/KDcyltest.cc @@ -343,8 +343,8 @@ main(int argc, char **argv) // This is the kd- NN density estimate; skipped by default for Ndens=0 // if (Ndens<=0) Ndens = 2; - typedef point point3; - typedef kdtree tree3; + using point3 = KDtree::point; + using tree3 = KDtree::kdtree; std::vector points; diff --git a/utils/Analysis/KL_cyl.cc b/utils/Analysis/KL_cyl.cc index c77a9416b..ec7f4755d 100644 --- a/utils/Analysis/KL_cyl.cc +++ b/utils/Analysis/KL_cyl.cc @@ -436,11 +436,12 @@ main(int argc, char **argv) // This is the kd- NN density estimate; skipped by default for Ndens=0 // if (Ndens) { + if (myid==0) std::cout << "Computing KD density estimate for " << nbod << " points" << std::endl; - typedef point point3; - typedef kdtree tree3; + using point3 = KDtree::point; + using tree3 = KDtree::kdtree; std::vector points; diff --git a/utils/Analysis/KL_sph.cc b/utils/Analysis/KL_sph.cc index 7c1a385cf..a9d205839 100644 --- a/utils/Analysis/KL_sph.cc +++ b/utils/Analysis/KL_sph.cc @@ -302,8 +302,8 @@ main(int argc, char **argv) if (myid==0) std::cout << "Computing KD density estimate for " << nbod << " points" << std::endl; - typedef point point3; - typedef kdtree tree3; + using point3 = KDtree::point ; + using tree3 = KDtree::kdtree; std::vector points; diff --git a/utils/Analysis/haloprof.cc b/utils/Analysis/haloprof.cc index 8e3c69b98..5fcddab97 100644 --- a/utils/Analysis/haloprof.cc +++ b/utils/Analysis/haloprof.cc @@ -811,8 +811,8 @@ main(int argc, char **argv) // Density weight COM // if (KD) { - typedef point point3; - typedef kdtree tree3; + using point3 = KDtree::point ; + using tree3 = KDtree::kdtree; std::vector points; diff --git a/utils/ICs/initial2d.cc b/utils/ICs/initial2d.cc index b892d7172..106d4e7b5 100644 --- a/utils/ICs/initial2d.cc +++ b/utils/ICs/initial2d.cc @@ -208,6 +208,7 @@ main(int ac, char **av) int nhalo, ndisk; std::string hbods, dbods, suffix, centerfile, halofile1, halofile2; std::string cachefile, config, gentype, dtype, dmodel, mtype, ctype; + std::string diskconf; const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n"); @@ -363,13 +364,15 @@ main(int ac, char **av) ("centerfile", "File containing phase-space center", cxxopts::value(centerfile)) ("suffix", "Suffix for output files", - cxxopts::value(suffix)->default_value("diag")) + cxxopts::value(suffix)) ("threads", "Number of threads to run", cxxopts::value(nthrds)->default_value("1")) ("M,MP", "For testing an m>0 harmonic distribution using test2d", cxxopts::value(M)->default_value("2")) ("pert", "For testing a quadrupole distribution using test2d", cxxopts::value(pert)->default_value("0.0")) + ("diskconf", "Custom disk configuration file for basis construction", + cxxopts::value(diskconf)) ("allow", "Allow multimass algorithm to generature negative masses for testing") ("nomono", "Allow non-monotonic mass interpolation") ("report", "Print out progress in BiorthCyl table evaluation") @@ -582,16 +585,30 @@ main(int ac, char **av) yml << YAML::Key << "verbose" << YAML::Value << true; // Build the diskconf stanza - yml << YAML::Key << "diskconf" << YAML::BeginMap - << YAML::Key << "name" << YAML::Value << "expon" - << YAML::Key << "parameters" - << YAML::BeginMap - << YAML::Key << "aexp" << YAML::Value << ACYL - << YAML::EndMap - << YAML::EndMap; - - // End the configuration map - yml << YAML::EndMap; + yml << YAML::Key << "diskconf"; + if (vm.count("diskconf")) { + YAML::Node node = YAML::Load(vm["diskconf"].as()); + yml << YAML::Value << node; + } + else { + // Begin the map + yml << YAML::BeginMap + << YAML::Key << "name" << YAML::Value << "expon" + << YAML::Key << "parameters" + << YAML::BeginMap + << YAML::Key << "aexp" << YAML::Value << ACYL + << YAML::EndMap + << YAML::EndMap; + // End the configuration map + yml << YAML::EndMap; + } + + // TEST + std::cout << std::string(80, '-') << std::endl + << "Test output" << std::endl + << std::string(80, '-') << std::endl + << yml.c_str() << std::endl + << std::string(80, '-') << std::endl; // Create expansion only if needed . . . std::shared_ptr expandd; diff --git a/utils/PhaseSpace/psp2vtu.cc b/utils/PhaseSpace/psp2vtu.cc index 1d7800334..45a8f94de 100644 --- a/utils/PhaseSpace/psp2vtu.cc +++ b/utils/PhaseSpace/psp2vtu.cc @@ -215,8 +215,8 @@ main(int ac, char **av) dens->SetName("density"); if (Ndens) { - typedef point point3; - typedef kdtree tree3; + using point3 = KDtree::point ; + using tree3 = KDtree::kdtree; std::vector mass; std::vector points; diff --git a/utils/SL/CMakeLists.txt b/utils/SL/CMakeLists.txt index 59afb08fa..dfbbb1e1f 100644 --- a/utils/SL/CMakeLists.txt +++ b/utils/SL/CMakeLists.txt @@ -1,6 +1,6 @@ set(bin_PROGRAMS slcheck slshift orthochk diskpot qtest eoftest - oftest slabchk) + oftest slabchk expontst) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp exputil ${VTK_LIBRARIES}) @@ -35,6 +35,7 @@ add_executable(diskpot diskpot.cc CylindricalDisk.cc SLSphere.cc) add_executable(qtest qtest.cc) add_executable(eoftest EOF2d.cc) add_executable(oftest oftest.cc) +add_executable(expontst expontest.cc) foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) diff --git a/utils/SL/EOF2d.cc b/utils/SL/EOF2d.cc index 9db5a1e35..3cb8087b6 100644 --- a/utils/SL/EOF2d.cc +++ b/utils/SL/EOF2d.cc @@ -17,7 +17,7 @@ int main(int argc, char** argv) { bool logr = false, cmap = false, ortho = false, plane = false; int numr, mmax, nmaxfid, nmax, knots, M, N, nradial, nout; - double A, rmin, rmax, rout; + double scale, rmin, rmax, rout; std::string filename, config, biorth; // Parse command line @@ -51,8 +51,8 @@ int main(int argc, char** argv) cxxopts::value(N)->default_value("256")) ("n,nradial", "Radial order for vertical potential output", cxxopts::value(nradial)->default_value("0")) - ("A,length", "characteristic disk scale length", - cxxopts::value(A)->default_value("1.0")) + ("s,length", "characteristic disk scale length for mapping", + cxxopts::value(scale)->default_value("1.0")) ("mmax", "maximum number of angular harmonics in the expansion", cxxopts::value(mmax)->default_value("4")) ("nmaxfid", "maximum number of radial basis harmonics for EOF construction", @@ -60,13 +60,13 @@ int main(int argc, char** argv) ("nmax", "maximum number of radial harmonics in the expansion", cxxopts::value(nmax)->default_value("10")) ("numr", "radial knots for the SL grid", - cxxopts::value(numr)->default_value("4000")) + cxxopts::value(numr)->default_value("8000")) ("r,rmin", "minimum radius for the SL grid", cxxopts::value(rmin)->default_value("0.0001")) ("R,rmax", "maximum radius for the SL grid", cxxopts::value(rmax)->default_value("20.0")) ("knots", "Number of Legendre integration knots", - cxxopts::value(knots)->default_value("200")) + cxxopts::value(knots)->default_value("1000")) ("rout", "Outer radius for evaluation", cxxopts::value(rout)->default_value("10.0")) ("nout", "number of points in the output grid per side", @@ -114,7 +114,7 @@ int main(int argc, char** argv) // Make the class instance // - EmpCyl2d emp(mmax, nmaxfid, nmax, knots, numr, rmin, rmax, A, cmap, logr, + EmpCyl2d emp(mmax, nmaxfid, nmax, knots, numr, rmin, rmax, scale, cmap, logr, par, biorth); if (vm.count("basis")) emp.basisTest(true); @@ -253,9 +253,8 @@ int main(int argc, char** argv) for (int n=0; n +#include +#include +#include +#include +#include +#include + +#include + +#include // Hankel computation for potential +#include // Gauss-Legendre quadrature +#include + +int main(int argc, char** argv) +{ + int M=0, N, nout; + double A, rmax, rout; + std::string filename; + + // Parse command line + // + cxxopts::Options options + (argv[0], + "Check PotRZ and QDHT using the exact exponential disk solution\n"); + + options.add_options() + ("h,help", "Print this help message") + ("N,nsize", "Default radial grid size", + cxxopts::value(N)->default_value("256")) + ("A,length", "characteristic disk scale length", + cxxopts::value(A)->default_value("1.0")) + ("rmax", "Outer radius for transform", + cxxopts::value(rmax)->default_value("10.0")) + ("rout", "Outer radius for evaluation", + cxxopts::value(rout)->default_value("10.0")) + ("nout", "number of points in the output grid per side", + cxxopts::value(nout)->default_value("40")) + ("o,filename", "Output filename", + cxxopts::value(filename)->default_value("test.potrz")) + ; + + + //=================== + // Parse options + //=================== + + cxxopts::ParseResult vm; + + try { + vm = options.parse(argc, argv); + } catch (cxxopts::OptionException& e) { + std::cout << "Option error: " << e.what() << std::endl; + return 2; + } + + // Print help message and exit + // + if (vm.count("help")) { + std::cout << options.help() << std::endl << std::endl; + return 1; + } + + // Open output file + // + std::ofstream out(filename); + if (not out) throw std::runtime_error("Could not open output file: " + + filename); + + // Define some representative limits + // + double Rmax = rout; + + // Grid spacing + // + double dR = Rmax/(nout - 1); + + // Potential instance with radially sensitive convergence parameters + // + PotRZ pot(rmax, N, M); + + // Set the functor using a lambda + // + auto dens = [A](double R) { + return -exp(-R/A); + }; + + auto potl = [A](double R) { + double x = 0.5*R/A; + return M_PI*R*(std::cyl_bessel_i(1, x) * std::cyl_bessel_k(0, x) - + std::cyl_bessel_i(0, x) * std::cyl_bessel_k(1, x) ); + }; + + // Write the results + // + for (int i=0; i