From 970474fc7afc980ceef5e37cd4f886b29be5a8bc Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 15 Apr 2025 19:15:13 -0400 Subject: [PATCH 01/22] Expose coefctr from Basis to client classes --- expui/BasisFactory.H | 2 ++ 1 file changed, 2 insertions(+) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index a9167396d..19b2fe6e6 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -263,6 +263,8 @@ namespace BasisClasses std::vector getFieldLabels(void) { return getFieldLabels(coordinates); } + //! Get the basis expansion center + std::vector getCenter() { return coefctr; } }; using BasisPtr = std::shared_ptr; From ec05fd98bff3452d512e5ecbc3c9203b458b09b6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 15 Apr 2025 19:15:20 -0400 Subject: [PATCH 02/22] Add omitted center subtraction for getFields evaluation in evalaccel() --- expui/BiorthBasis.cc | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index f9237bc65..e6eab3a2d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3685,11 +3685,17 @@ namespace BasisClasses // auto basis = std::get<0>(mod); + // Get expansion center + // + auto ctr = basis->getCenter(); + // Get fields // int rows = accel.rows(); for (int n=0; ngetFields(ps(n, 0), ps(n, 1), ps(n, 2)); + 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, follewed by acceleration for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k); } @@ -3899,10 +3905,14 @@ namespace BasisClasses // int numT = floor( (tfinal - tinit)/h ); - // Compute output step + // Number of output steps. Will attempt to find the best stride... + // + int stride = std::ceil(static_cast(numT)/static_cast(nout)); + if (stride>1) numT = nout * stride; + + // Compute the interval-matching step // - nout = std::min(numT, nout); - double H = (tfinal - tinit)/nout; + h = (tfinal - tinit)/(numT - 1); // Return data // @@ -3934,9 +3944,11 @@ namespace BasisClasses for (int k=0; k<6; k++) ret(n, k, 0) = ps(n, k); double tnow = tinit; - for (int s=1, cnt=1; s= H*cnt-h*1.0e-8) { + if (s++ % stride == 0) { times(cnt) = tnow; for (int n=0; n Date: Tue, 15 Apr 2025 19:23:14 -0400 Subject: [PATCH 03/22] Restore BiorthBasis.cc to pre tfinalFix condition for this branch --- expui/BiorthBasis.cc | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e6eab3a2d..5397b3f2b 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3905,14 +3905,10 @@ namespace BasisClasses // int numT = floor( (tfinal - tinit)/h ); - // Number of output steps. Will attempt to find the best stride... + // Compute output step // - int stride = std::ceil(static_cast(numT)/static_cast(nout)); - if (stride>1) numT = nout * stride; - - // Compute the interval-matching step - // - h = (tfinal - tinit)/(numT - 1); + nout = std::min(numT, nout); + double H = (tfinal - tinit)/nout; // Return data // @@ -3944,11 +3940,9 @@ namespace BasisClasses for (int k=0; k<6; k++) ret(n, k, 0) = ps(n, k); double tnow = tinit; - int s = 0, cnt = 0; - while (s < numT) { - if (tfinal - tnow < h) h = tfinal - tnow; + for (int s=1, cnt=1; s= H*cnt-h*1.0e-8) { times(cnt) = tnow; for (int n=0; n Date: Sat, 19 Apr 2025 21:47:51 +0100 Subject: [PATCH 04/22] Expose expansion center in pyEXP --- pyEXP/CoefWrappers.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 97c394db0..765db774b 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -717,6 +717,11 @@ void CoefficientClasses(py::module &m) { float data's time stamp )") + .def_readwrite("center", &CoefStruct::ctr, + R"( + array of floats + the expansion center + )") .def("getCoefs", &CoefStruct::getCoefs, R"( Read-only access to the underlying data store From 09e6c2955734b3f5df00e3bf2884725e435ee8ca Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sun, 20 Apr 2025 14:05:28 +0100 Subject: [PATCH 05/22] add getCenter/setCenter members --- expui/CoefStruct.H | 11 +++++++++++ pyEXP/CoefWrappers.cc | 43 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 50 insertions(+), 4 deletions(-) diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 4afca4cf9..3f662e8c3 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -62,6 +62,17 @@ 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; } + }; //! Specialization of CoefStruct for spheres diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 765db774b..8dd17b178 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -717,11 +717,46 @@ void CoefficientClasses(py::module &m) { float data's time stamp )") - .def_readwrite("center", &CoefStruct::ctr, + .def_readonly("center", &CoefStruct::ctr, R"( - array of floats - the expansion center - )") + array of floats + the expansion center + )") + .def("getCenter", &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("setCenter", + 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 From 093818c9487974c095bc93ba72e547909e287302 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 21 Apr 2025 10:56:13 +0100 Subject: [PATCH 06/22] change center member names; create time members --- expui/CoefStruct.H | 9 +++++++++ pyEXP/CoefWrappers.cc | 42 ++++++++++++++++++++++++++++++++++++------ 2 files changed, 45 insertions(+), 6 deletions(-) diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 3f662e8c3..4276cb95e 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -73,6 +73,15 @@ namespace CoefClasses //! 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/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 8dd17b178..d8e62127c 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -717,12 +717,42 @@ void CoefficientClasses(py::module &m) { float data's time stamp )") - .def_readonly("center", &CoefStruct::ctr, - R"( - array of floats - the expansion center + .def("getCoefTime", &CoefStruct::getTime, + R"( + Read-only access to the coefficient time + + Returns + ------- + numpy.ndarray + vector of center data + + See also + -------- + setCenter : read-write access to the coefficient time + )") + .def("setCoefTime", + static_cast(&CoefStruct::setTime), + py::arg("tval"), + R"( + Set the coefficien time + + Parameters + ---------- + tval : float + time value + + Returns + ------- + None + + Notes + ----- + + See also + -------- + getCenter : read-only access to coefficient time )") - .def("getCenter", &CoefStruct::getCenter, + .def("getCoefCenter", &CoefStruct::getCenter, R"( Read-only access to the center data @@ -735,7 +765,7 @@ void CoefficientClasses(py::module &m) { -------- setCenter : read-write access to the center data )") - .def("setCenter", + .def("setCoefCenter", static_cast&)>(&CoefStruct::setCenter), py::arg("mat"), R"( From 08957f38e4f2d3b5382b7373575d0a7a0025fcc7 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 21 Apr 2025 10:59:34 +0100 Subject: [PATCH 07/22] version bump only [no CI] --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ffba0d1b6..dbf34eaf7 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.2" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) From 9fe6efa40b841ae6d9874626aee5b44b5c09c81f Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 21 Apr 2025 19:57:51 +0100 Subject: [PATCH 08/22] add readonly attributes --- pyEXP/CoefWrappers.cc | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index d8e62127c..d01718371 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -712,11 +712,16 @@ 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 From fbc25a4d9f9e148f26f231cece55ca14f54ea8e2 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 22 Apr 2025 11:05:04 +0100 Subject: [PATCH 09/22] add default center --- expui/CoefStruct.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 4276cb95e..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() {} From aafec0f3920d20135de85213fb841f499ffbc8c1 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Tue, 22 Apr 2025 13:47:51 +0100 Subject: [PATCH 10/22] Apply suggestions from code review Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/CoefWrappers.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index d01718371..fc4fa42dd 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -733,13 +733,13 @@ void CoefficientClasses(py::module &m) { See also -------- - setCenter : read-write access to the coefficient time + setCoefTime : read-write access to the coefficient time )") .def("setCoefTime", static_cast(&CoefStruct::setTime), py::arg("tval"), R"( - Set the coefficien time + Set the coefficient time Parameters ---------- From ec11f5481c3c8e5c4271d6473e6a43ca6c93dbff Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 22 Apr 2025 14:11:22 -0400 Subject: [PATCH 11/22] Expose the pseudo-force infrastructure through pybind11 --- expui/BasisFactory.H | 7 +++- expui/BasisFactory.cc | 16 ++++++-- expui/BiorthBasis.cc | 3 +- pyEXP/BasisWrappers.cc | 89 +++++++++++++++++++++++++++++++++++++++++- 4 files changed, 106 insertions(+), 9 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 19b2fe6e6..2ee5ea584 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -249,8 +249,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(const Eigen::VectorXd& x, const Eigen::MatrixXd& pos); + void setNonInertial(const std::string& orient); //@} //! Set the current pseudo acceleration @@ -259,6 +259,9 @@ namespace BasisClasses if (Naccel > 0) pseudo = currentAccel(time); } + //! Returns true if non-intertial forces are active + bool usingNonInertial() { return Naccel > 0; } + //! Get the field label vector std::vector getFieldLabels(void) { return getFieldLabels(coordinates); } diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 483dcca46..7760eb038 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -283,14 +283,22 @@ namespace BasisClasses return makeFromArray(time); } - void Basis::setNonInertial(int N, Eigen::VectorXd& t, Eigen::MatrixXd& pos) + void Basis::setNonInertial(const Eigen::VectorXd& t, const Eigen::MatrixXd& pos) { - Naccel = N; + // 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 = t.size(); t_accel = t; p_accel = pos; } - void Basis::setNonInertial(int N, const std::string& orient) + void Basis::setNonInertial(const std::string& orient) { std::ifstream in(orient); @@ -337,7 +345,7 @@ namespace BasisClasses } // Repack data - Naccel = N; + Naccel = times.size(); t_accel.resize(times.size()); p_accel.resize(times.size(), 3); for (int i=0; igetCenter(); + std::vector ctr {0, 0, 0}; + if (basis->usingNonInertial()) ctr = basis->getCenter(); // Get fields // diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index bb214cf20..c713d10c4 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 times (N) 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,82 @@ void BasisFactoryClasses(py::module &m) ------- list: str list of basis function labels - )" - ); + )" + ) + .def("setNonInertial", + [](BasisClasses::Basis& A, + const Eigen::VectorXd& times, const Eigen::MatrixXd& pos) { + A.setNonInertial(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 + ---------- + 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("times"), py::arg("pos") + ) + .def("setNonInertial", + [](BasisClasses::Basis& A, const std::string orient) + { + A.setNonInertial(orient); + }, + R"( + Initialize for pseudo-force computation with a time series of positions + using a EXP orient file + + Parameters + ---------- + 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("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(), From dd1105cb6fb3d48483d83145eef0eebebb7b17b7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 22 Apr 2025 17:15:28 -0400 Subject: [PATCH 12/22] Fixed the flipped logic of inertial vs non-inertial center subtraction --- expui/BiorthBasis.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 54e08c063..fb0a0dcef 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3687,8 +3687,8 @@ namespace BasisClasses // Get expansion center // - std::vector ctr {0, 0, 0}; - if (basis->usingNonInertial()) ctr = basis->getCenter(); + auto ctr = basis->getCenter(); + if (basis->usingNonInertial()) ctr = {0, 0, 0}; // Get fields // From 8accceb37d5f08bd8137eaa750f2312bbf1550fd Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 22 Apr 2025 18:32:24 -0400 Subject: [PATCH 13/22] Restore Naccel setting, oops --- expui/BasisFactory.H | 4 ++-- expui/BasisFactory.cc | 8 ++++---- pyEXP/BasisWrappers.cc | 18 +++++++++++------- 3 files changed, 17 insertions(+), 13 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 2ee5ea584..511fc5669 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -249,8 +249,8 @@ namespace BasisClasses //@{ //! Initialize non-inertial forces - void setNonInertial(const Eigen::VectorXd& x, const Eigen::MatrixXd& pos); - void setNonInertial(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 diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 7760eb038..088fd8244 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -283,7 +283,7 @@ namespace BasisClasses return makeFromArray(time); } - void Basis::setNonInertial(const Eigen::VectorXd& t, const Eigen::MatrixXd& pos) + void Basis::setNonInertial(int N, const Eigen::VectorXd& t, const Eigen::MatrixXd& pos) { // Sanity checks if (t.size() < 1) @@ -293,12 +293,12 @@ namespace BasisClasses throw std::runtime_error("Basis::setNonInertial: size mismatch in time and position arrays"); // Set the data - Naccel = t.size(); + Naccel = N; t_accel = t; p_accel = pos; } - void Basis::setNonInertial(const std::string& orient) + void Basis::setNonInertial(int N, const std::string& orient) { std::ifstream in(orient); @@ -345,7 +345,7 @@ namespace BasisClasses } // Repack data - Naccel = times.size(); + Naccel = N; t_accel.resize(times.size()); p_accel.resize(times.size(), 3); for (int i=0; i Date: Tue, 22 Apr 2025 22:09:03 -0400 Subject: [PATCH 14/22] Put the evaluation point at the center of the interval when possible --- expui/BasisFactory.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 088fd8244..cedd05fdc 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -368,7 +368,10 @@ namespace BasisClasses 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; + + // Get a range around the current time of approx size Naccel + imax = std::min(t_accel.size()-1, imax + Naccel/2); + imin = std::max(imax - Naccel, 0); int num = imax - imin + 1; Eigen::VectorXd t(num); From 47edaf1234abc2d0d76e4426ae5cc84efda38f7f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 23 Apr 2025 13:06:49 -0400 Subject: [PATCH 15/22] Added a inertial coodindate reset member; fixed an interpolation fence-post error that did not affect the results. --- expui/BasisFactory.H | 8 ++++++++ expui/BiorthBasis.cc | 24 +++++++++++------------- pyEXP/BasisWrappers.cc | 18 ++++++++++++++++++ 3 files changed, 37 insertions(+), 13 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 511fc5669..73059a1e6 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}; @@ -262,6 +263,13 @@ namespace BasisClasses //! Returns true if non-intertial 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); } diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index fb0a0dcef..96164d328 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3727,13 +3727,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); @@ -3791,13 +3790,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); diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 27b3d1528..cef5ac170 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -959,6 +959,24 @@ void BasisFactoryClasses(py::module &m) 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 acceration + )" + ) .def("setNonInertial", [](BasisClasses::Basis& A, int N, const Eigen::VectorXd& times, const Eigen::MatrixXd& pos) { From 2cb664a8507b63469495c15ed5043106a9ae5506 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 23 Apr 2025 13:29:37 -0400 Subject: [PATCH 16/22] Comment typo only --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 96164d328..eaac49ab9 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3697,7 +3697,7 @@ namespace BasisClasses 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, follewed by acceleration + // 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); } From f792b9ca6133b2f6cb092b43c806269ec245fd49 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 23 Apr 2025 16:33:34 -0400 Subject: [PATCH 17/22] Add a bit of tolerance to the center array grid interpolation before flagging out of bounds --- expui/BasisFactory.cc | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index cedd05fdc..b410c7861 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -359,18 +359,30 @@ 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(); + 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(t_accel.size()-1, imax + Naccel/2); + // + imax = std::min(n-1, imax + Naccel/2); imin = std::max(imax - Naccel, 0); int num = imax - imin + 1; From c3a227589db6bf5c29dace9d9626b84f9814be22 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 23 Apr 2025 16:34:56 -0400 Subject: [PATCH 18/22] Assign passed coefficients as the coefret internal basis instance so this can be queried by clients; added pseudo accel debug output for testing --- expui/BiorthBasis.cc | 48 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index eaac49ab9..e87c172c0 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++) { @@ -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,10 @@ namespace BasisClasses CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); auto & cc = *cf->coefs; + // Cache the cuurent coefficient structure + // + coefret = coef; + // Assign internal coefficient table (doubles) from the complex struct // for (int m=0, m0=0; m<=mmax; m++) { @@ -2874,6 +2890,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 +3326,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}; } @@ -3701,6 +3725,30 @@ namespace BasisClasses for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k); } + // true for deep debugging + // | + // v + if (true 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; } From c74d5b3d580b4b68067c8842ab549fec51b983f0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 23 Apr 2025 20:47:40 -0400 Subject: [PATCH 19/22] Some cleanup; turn of RK4 test and pseudo force test output --- expui/BiorthBasis.cc | 104 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 86 insertions(+), 18 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e87c172c0..1972ad9bd 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -635,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); @@ -3728,7 +3728,7 @@ namespace BasisClasses // true for deep debugging // | // v - if (true and basis->usingNonInertial()) { + if (false and basis->usingNonInertial()) { auto coefs = basis->getCoefficients(); auto time = coefs->time; @@ -3890,25 +3890,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); } From e4a977c215ca93b435a02ce8008c95cfc9bdd7bc Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Thu, 24 Apr 2025 10:54:52 +0100 Subject: [PATCH 20/22] Apply suggestions from code review Typos in docstrings only. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BasisFactory.H | 2 +- expui/BiorthBasis.cc | 3 +-- pyEXP/BasisWrappers.cc | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 73059a1e6..1598f3f7d 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -260,7 +260,7 @@ namespace BasisClasses if (Naccel > 0) pseudo = currentAccel(time); } - //! Returns true if non-intertial forces are active + //! Returns true if non-inertial forces are active bool usingNonInertial() { return Naccel > 0; } //! Reset to inertial coordinates diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 1972ad9bd..c0cd6e9e5 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -2486,8 +2486,7 @@ namespace BasisClasses CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); auto & cc = *cf->coefs; - // Cache the cuurent coefficient structure - // + // Cache the current coefficient structure coefret = coef; // Assign internal coefficient table (doubles) from the complex struct diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index cef5ac170..fb1332a11 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -974,7 +974,7 @@ void BasisFactoryClasses(py::module &m) See also -------- setNonInertial : set non-inertial data - setNonInertialAccel : set the non-inertial acceration + setNonInertialAccel : set the non-inertial acceleration )" ) .def("setNonInertial", From 35d26a6ae8c78faff3e6290ddf907f5a4ebedd53 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Thu, 24 Apr 2025 21:13:36 +0100 Subject: [PATCH 21/22] Return Version tuple for compatibility checking --- pyEXP/UtilWrappers.cc | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) 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.)"); } From 1cb6c2de74bab092744a819762d17ebe3ccd547f Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Thu, 24 Apr 2025 21:19:47 +0100 Subject: [PATCH 22/22] version bump [no CI] --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index dbf34eaf7..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.2" + VERSION "7.8.3" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran)