diff --git a/CMakeLists.txt b/CMakeLists.txt index ffba0d1b6..9f5742909 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features project( EXP - VERSION "7.8.1" + VERSION "7.8.3" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index a9167396d..1598f3f7d 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -137,6 +137,7 @@ namespace BasisClasses Eigen::Vector3d currentAccel(double time); public: + //! The current pseudo acceleration Eigen::Vector3d pseudo {0, 0, 0}; @@ -249,8 +250,8 @@ namespace BasisClasses //@{ //! Initialize non-inertial forces - void setNonInertial(int N, Eigen::VectorXd& x, Eigen::MatrixXd& pos); - void setNonInertial(int N, const std::string& orient); + void setNonInertial(int Naccel, const Eigen::VectorXd& x, const Eigen::MatrixXd& pos); + void setNonInertial(int Naccel, const std::string& orient); //@} //! Set the current pseudo acceleration @@ -259,10 +260,22 @@ namespace BasisClasses if (Naccel > 0) pseudo = currentAccel(time); } + //! Returns true if non-inertial forces are active + bool usingNonInertial() { return Naccel > 0; } + + //! Reset to inertial coordinates + void setInertial() + { + Naccel = 0; + pseudo = {0, 0, 0}; + } + //! Get the field label vector std::vector getFieldLabels(void) { return getFieldLabels(coordinates); } + //! Get the basis expansion center + std::vector getCenter() { return coefctr; } }; using BasisPtr = std::shared_ptr; diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 483dcca46..b410c7861 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -283,8 +283,16 @@ namespace BasisClasses return makeFromArray(time); } - void Basis::setNonInertial(int N, Eigen::VectorXd& t, Eigen::MatrixXd& pos) + void Basis::setNonInertial(int N, const Eigen::VectorXd& t, const Eigen::MatrixXd& pos) { + // Sanity checks + if (t.size() < 1) + throw std::runtime_error("Basis: setNonInertial: no times in time array"); + + if (t.size() != pos.rows()) + throw std::runtime_error("Basis::setNonInertial: size mismatch in time and position arrays"); + + // Set the data Naccel = N; t_accel = t; p_accel = pos; @@ -351,16 +359,31 @@ namespace BasisClasses { Eigen::Vector3d ret; - if (time < t_accel(0) || time > 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; + auto n = t_accel.size(); + + // Allow a little bit of buffer in the allowable on-grid range but + // otherwise force termination + // + if ( time < t_accel(0 ) - 0.5*(t_accel(1 ) - t_accel(0 )) || + time > t_accel(n-1) + 0.5*(t_accel(n-1) - t_accel(n-2)) ) { + + std::ostringstream msg; + msg << "Basis::currentAccel: " << time + << " is outside the range of the non-inertial DB [" + << t_accel(0) << ", " << t_accel(n-1) << "]"; + + throw std::runtime_error(msg.str()); } + // Do the quadratic interpolation + // 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 imax = std::lower_bound(t_accel.data(), t_accel.data()+n, time) - t_accel.data(); + + // Get a range around the current time of approx size Naccel + // + imax = std::min(n-1, imax + Naccel/2); + imin = std::max(imax - Naccel, 0); int num = imax - imin + 1; Eigen::VectorXd t(num); diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 7705c7bc1..c0cd6e9e5 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -396,6 +396,10 @@ namespace BasisClasses CoefClasses::SphStruct* cf = dynamic_cast(coef.get()); + // Cache the current coefficient structure + // + coefret = coef; + // Assign internal coefficient table (doubles) from the complex struct // for (int l=0, L0=0, L1=0; l<=lmax; l++) { @@ -631,7 +635,7 @@ namespace BasisClasses Spherical::crt_eval (double x, double y, double z) { - double R = sqrt(x*x + y*y); + double R = sqrt(x*x + y*y) + 1.0e-18; double phi = atan2(y, x); auto v = cyl_eval(R, z, phi); @@ -1474,6 +1478,10 @@ namespace BasisClasses CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); + // Cache the current coefficient structure + // + coefret = coef; + for (int m=0; m<=mmax; m++) { // Set to zero on m=0 call only--------+ sl->set_coefs(m, (*cf->coefs).row(m).real(), (*cf->coefs).row(m).imag(), m==0); } @@ -1771,6 +1779,10 @@ namespace BasisClasses CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); auto & cc = *cf->coefs; + // Cache the current coefficient structure + // + coefret = coef; + // Assign internal coefficient table (doubles) from the complex struct // for (int m=0, m0=0; m<=mmax; m++) { @@ -2474,6 +2486,9 @@ namespace BasisClasses CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); auto & cc = *cf->coefs; + // Cache the current coefficient structure + coefret = coef; + // Assign internal coefficient table (doubles) from the complex struct // for (int m=0, m0=0; m<=mmax; m++) { @@ -2874,6 +2889,10 @@ namespace BasisClasses auto cf = dynamic_cast(coef.get()); expcoef = *cf->coefs; + // Cache the current coefficient structure + // + coefret = coef; + coefctr = {0.0, 0.0, 0.0}; } @@ -3306,6 +3325,10 @@ namespace BasisClasses auto cf = dynamic_cast(coef.get()); expcoef = *cf->coefs; + // Cache the cuurent coefficient structure + // + coefret = coef; + coefctr = {0.0, 0.0, 0.0}; } @@ -3685,15 +3708,46 @@ namespace BasisClasses // auto basis = std::get<0>(mod); + // Get expansion center + // + auto ctr = basis->getCenter(); + if (basis->usingNonInertial()) ctr = {0, 0, 0}; + // Get fields // int rows = accel.rows(); for (int n=0; ngetFields(ps(n, 0), ps(n, 1), ps(n, 2)); - // First 6 fields are density and potential, follewed by acceleration + auto v = basis->getFields(ps(n, 0) - ctr[0], + ps(n, 1) - ctr[1], + ps(n, 2) - ctr[2]); + // First 6 fields are density and potential, followed by acceleration for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k); } + // true for deep debugging + // | + // v + if (false and basis->usingNonInertial()) { + + auto coefs = basis->getCoefficients(); + auto time = coefs->time; + auto ctr = coefs->ctr; + + std::ofstream tmp; + if (time <= 0.0) tmp.open("pseudo.dat"); + else tmp.open("pseudo.dat", ios::app); + + if (tmp) + tmp << std::setw(16) << std::setprecision(5) << time + << std::setw(16) << std::setprecision(5) << ctr[0] + << std::setw(16) << std::setprecision(5) << ctr[1] + << std::setw(16) << std::setprecision(5) << ctr[2] + << std::setw(16) << std::setprecision(5) << basis->pseudo(0) + << std::setw(16) << std::setprecision(5) << basis->pseudo(1) + << std::setw(16) << std::setprecision(5) << basis->pseudo(2) + << std::endl; + } + return accel; } @@ -3720,13 +3774,12 @@ namespace BasisClasses throw std::runtime_error(sout.str()); } - auto it1 = std::lower_bound(times.begin(), times.end(), t); - auto it2 = it1 + 1; + auto it2 = std::lower_bound(times.begin(), times.end(), t); + auto it1 = it2; - if (it2 == times.end()) { - it2--; - it1 = it2 - 1; - } + if (it2 == times.end()) throw std::runtime_error("Basis::AllTimeAccel::evalcoefs: time t=" + std::to_string(t) + " out of bounds"); + else if (it2 == times.begin()) it2++; + else it1--; double a = (*it2 - t)/(*it2 - *it1); double b = (t - *it1)/(*it2 - *it1); @@ -3784,13 +3837,12 @@ namespace BasisClasses throw std::runtime_error(sout.str()); } - auto it1 = std::lower_bound(times.begin(), times.end(), t); - auto it2 = it1 + 1; - - if (it2 == times.end()) { - it2--; - it1 = it2 - 1; - } + auto it2 = std::lower_bound(times.begin(), times.end(), t); + auto it1 = it2; + + if (it2 == times.end()) throw std::runtime_error("Basis::AllTimeAccel::evalcoefs: time t=" + std::to_string(t) + " out of bounds"); + else if (it2 == times.begin()) it2++; + else it1--; double a = (*it2 - t)/(*it2 - *it1); double b = (t - *it1)/(*it2 - *it1); @@ -3837,25 +3889,93 @@ namespace BasisClasses { int rows = ps.rows(); - // Drift 1/2 - for (int n=0; n kf(4); + for (int i=0; i<4; i++) kf[i].resize(rows, 6); + + // Step 1 + // + accel.setZero(); + for (auto mod : bfe) F(t, ps, accel, mod); + for (int n=0; n(t+h, ps); } diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 4afca4cf9..26d011e28 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -27,7 +27,7 @@ namespace CoefClasses Eigen::VectorXcd store; //! Center data - std::vector ctr; + std::vector ctr= {0.0, 0.0, 0.0}; //! Destructor virtual ~CoefStruct() {} @@ -62,6 +62,26 @@ namespace CoefClasses //! Read-only access to coefficient data (no copy) Eigen::Ref getCoefs() { return store; } + //! Set new center data (no copy) + void setCenter(std::vector& STORE) + { + if (STORE.size() != ctr.size()) + throw std::invalid_argument("CoefStruct::setCenter: center vector size does not match"); + ctr = STORE; + } + + //! Read-only access to center (no copy) + std::vector getCenter() { return ctr; } + + //! Set coefficient time (no copy) + void setTime(double& STORE) + { + time = STORE; + } + + //! Read-only access to center (no copy) + double getTime() { return time; } + }; //! Specialization of CoefStruct for spheres diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index bb214cf20..fb1332a11 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -199,6 +199,18 @@ void BasisFactoryClasses(py::module &m) a fixed potential model. AccelFunc can be inherited by a native Python class and the evalcoefs() may be implemented in Python and passed to IntegrateOrbits in the same way as a native C++ class. + + Non-inertial frames of reference + -------------------------------- + Each component of a multiple component simulation may have its own expansion + center. Orbit integration in the frame of reference of the expansion is + accomplished by defining a moving frame of reference using the setNonInertial() + call with either an array of n times and center positions (as an nx3 array) + or by initializing with an EXP orient file. + + We provide a member function, setNonInertialAccel(t), to estimate the frame + acceleration at a given time. This may be useful for user-defined acceleration + routines. This is automatically called default C++ evalcoefs() routine. )"; using namespace BasisClasses; @@ -945,9 +957,104 @@ void BasisFactoryClasses(py::module &m) ------- list: str list of basis function labels - )" - ); + )" + ) + .def("setInertial", &BasisClasses::Basis::setInertial, + R"( + Reset to inertial coordinates + + Parameters + ---------- + None + + Returns + ------- + None + + See also + -------- + setNonInertial : set non-inertial data + setNonInertialAccel : set the non-inertial acceleration + )" + ) + .def("setNonInertial", + [](BasisClasses::Basis& A, + int N, const Eigen::VectorXd& times, const Eigen::MatrixXd& pos) { + A.setNonInertial(N, times, pos); + }, + R"( + Initialize for pseudo-force computation with a time series of positions + using (1) a time vector and (2) a center position matrix with rows of three + vectors + + Parameters + ---------- + N : int + number of previous positions to use for quadratic fit + times : list or numpy.ndarray + list of time points + pos : numpy.ndarray + an array with N rows and 3 columns of center positions + + Returns + ------- + None + + See also + -------- + setNonInertial : set non-inertial data from an Orient file + setNonInertialAccel : set the non-inertial acceration + )", + py::arg("N"), py::arg("times"), py::arg("pos") + ) + .def("setNonInertial", + [](BasisClasses::Basis& A, int N, const std::string orient) + { + A.setNonInertial(N, orient); + }, + R"( + Initialize for pseudo-force computation with a time series of positions + using a EXP orient file + + Parameters + ---------- + N : int + number of previous positions to use for quadratic fit + orient : str + name of the orient file + + Returns + ------- + None + + See also + -------- + setNonInertial : set non-inertial data from a time series of values + setNonInertialAccel : set the non-inertial acceration + )", + py::arg("N"), py::arg("orient") + ) + .def("setNonInertialAccel", &BasisClasses::Basis::setNonInertialAccel, + R"( + Set the pseudo acceleration for the non-inertial data at a given time + + Parameters + ---------- + time : float + evaluation time + Returns + ------- + None + + See also + -------- + setNonInertial : set non-inertial data from a time series of values + setNonInertial : set non-inertial data from an EXP orient file + )", + py::arg("time") + ); + py::class_, PyBiorthBasis, BasisClasses::Basis> (m, "BiorthBasis") .def(py::init(), diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 97c394db0..fc4fa42dd 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -712,11 +712,86 @@ void CoefficientClasses(py::module &m) { str geometry type )") - .def_readwrite("time", &CoefStruct::time, + .def_readonly("time", &CoefStruct::time, R"( float data's time stamp )") + .def_readonly("center", &CoefStruct::ctr, + R"( + float + data's center value + )") + .def("getCoefTime", &CoefStruct::getTime, + R"( + Read-only access to the coefficient time + + Returns + ------- + numpy.ndarray + vector of center data + + See also + -------- + setCoefTime : read-write access to the coefficient time + )") + .def("setCoefTime", + static_cast(&CoefStruct::setTime), + py::arg("tval"), + R"( + Set the coefficient time + + Parameters + ---------- + tval : float + time value + + Returns + ------- + None + + Notes + ----- + + See also + -------- + getCenter : read-only access to coefficient time + )") + .def("getCoefCenter", &CoefStruct::getCenter, + R"( + Read-only access to the center data + + Returns + ------- + numpy.ndarray + vector of center data + + See also + -------- + setCenter : read-write access to the center data + )") + .def("setCoefCenter", + static_cast&)>(&CoefStruct::setCenter), + py::arg("mat"), + R"( + Set the center vector + + Parameters + ---------- + mat : numpy.ndarray + center vector + + Returns + ------- + None + + Notes + ----- + + See also + -------- + getCenter : read-only access to center data + )") .def("getCoefs", &CoefStruct::getCoefs, R"( Read-only access to the underlying data store diff --git a/pyEXP/UtilWrappers.cc b/pyEXP/UtilWrappers.cc index 119005c84..04a80ae01 100644 --- a/pyEXP/UtilWrappers.cc +++ b/pyEXP/UtilWrappers.cc @@ -155,4 +155,22 @@ void UtilityClasses(py::module &m) { Report on the version and git commit. This is the same version information reported by the EXP N-body code.)"); + m.def("Version", + []() { + std::string version_str = VERSION; + std::istringstream iss(version_str); + std::string token; + int parts[3] = {0, 0, 0}; + int i = 0; + + while (std::getline(iss, token, '.') && i < 3) { + parts[i++] = std::stoi(token); + } + + return std::make_tuple(parts[0], parts[1], parts[2]); + }, + R"( + Return the version. + + This is the same version information reported by the EXP N-body code.)"); }