diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 08988793..965c051c 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 d3c11304..483dcca4 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -283,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.cc b/expui/BiorthBasis.cc index d91d1142..513c2fd9 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -592,7 +592,7 @@ namespace BasisClasses double densfac = 1.0/(scale*scale*scale) * 0.25/M_PI; double potlfac = 1.0/scale; - + return {den0 * densfac, // 0 den1 * densfac, // 1 @@ -3653,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; @@ -3722,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/src/Component.H b/src/Component.H index 5dfc1ceb..170f1082 100644 --- a/src/Component.H +++ b/src/Component.H @@ -87,6 +87,9 @@ struct loadb_datum ensemble used to determine the paramters. The energy cutoff is dynamically adjusted to achieve this target. + @param nEJwant is the size of the past states used to estimate the + the acceleration of the expansion frame + @param EJx0 is the initial EJ center x-coordinate (default: 0) @param EJy0 is the initial EJ center y-coordinate (default: 0) @@ -399,10 +402,12 @@ public: double *com; //! Center of velocity double *cov; - //! Center of accelration + //! Center of acceleration double *coa; //! Center (e.g. set by Orient) double *center; + //! Pseudo-accleration (set by Orient) + double *accel; //! Angular momentum vector double *angmom; //! Phase space vector @@ -453,6 +458,9 @@ public: //! Target number of particles for average int nEJwant; + //! Target number of states for pseudo-acceleration estimation + int nEJaccel; + //! Initial EJ center //@{ //! x-coord @@ -859,6 +867,15 @@ public: //! Add to accerlation (by component) inline void AddAcc(int i, int j, double val) { + PartMap::iterator tp = particles.find(i); + if (tp == particles.end()) { + throw BadIndexException(i, particles.size(), __FILE__, __LINE__); + } + tp->second->acc[j] += val - accel[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 +885,15 @@ 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__); + } + for (int k=0; k<3; k++) tp->second->acc[k] += val[k] - accel[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 +903,15 @@ 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__); + } + for (int k=0; k<3; k++) tp->second->acc[k] += val[k] - accel[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 b1cab94b..c4fd1365 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -47,6 +47,7 @@ const std::set Component::valid_keys_parm = "EJ", "nEJkeep", "nEJwant", + "nEJaccel", "EJkinE", "EJext", "EJdiag", @@ -184,6 +185,7 @@ Component::Component(YAML::Node& CONF) EJ = 0; nEJkeep = 100; nEJwant = 500; + nEJaccel = 0; EJkinE = true; EJext = false; EJdiag = false; @@ -320,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; @@ -683,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; @@ -795,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(); @@ -916,6 +921,7 @@ void Component::initialize(void) { com = new double [3]; center = new double [3]; + accel = new double [3]; cov = new double [3]; coa = new double [3]; angmom = new double [3]; @@ -929,6 +935,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) { @@ -1113,6 +1120,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; @@ -1139,10 +1147,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(); + Eigen::VectorXd::Map(&accel [0], 3) = orient->currentAccel(); } else { if (EJlinear) orient -> set_linear(); if (not com_system) { @@ -1238,6 +1248,7 @@ Component::~Component(void) delete [] com; delete [] center; + delete [] accel; delete [] cov; delete [] coa; delete [] angmom; @@ -3111,12 +3122,15 @@ void Component::fix_positions_cpu(unsigned mlevel) if ((EJ & Orient::CENTER) && !EJdryrun) { auto ctr = orient->currentCenter(); + auto acc = orient->currentAccel(); bool ok = true; for (int i=0; i<3; i++) { if (std::isnan(ctr[i])) ok = false; + if (std::isnan(acc[i])) ok = false; } if (ok) { for (int i=0; i<3; i++) center[i] += ctr[i]; + for (int i=0; i<3; i++) accel [i] += acc[i]; } else if (myid==0) { cout << "Orient: center failure, T=" << tnow << ", adjustment skipped" << endl; diff --git a/src/Orient.H b/src/Orient.H index b2bb1fa9..110c9fff 100644 --- a/src/Orient.H +++ b/src/Orient.H @@ -8,6 +8,7 @@ #include #include +#include #include @@ -96,7 +97,7 @@ private: int keep, current; double damp; Eigen::Matrix3d body, orig; - Eigen::Vector3d axis, center, axis1, center1, center0, cenvel0; + Eigen::Vector3d axis, center, axis1, center1, center0, cenvel0, pseudo; double sumX, sumX2; Eigen::Vector3d sumY, sumY2, sumXY, slope, intercept; double sigA, sigC, sigCz; @@ -110,8 +111,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 +128,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 +170,14 @@ public: //! Return current center Eigen::Vector3d& currentCenter(void) {return center;}; + //! Return current pseudo force + Eigen::Vector3d& currentAccel(void) + { + pseudo.setZero(); + if (accel) return pseudo = accel->accel(); + return 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 676e4079..c865221a 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,8 @@ Orient::Orient(int n, int nwant, unsigned Oflg, unsigned Cflg, cenvel0.setZero(); axis .setZero(); + if (naccel) accel = std::make_shared(naccel); + // Initialize last time to something // in the near infinite past lasttime = -std::numeric_limits::max(); @@ -70,8 +72,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); @@ -124,6 +125,8 @@ Orient::Orient(int n, int nwant, unsigned Oflg, unsigned Cflg, // Look for data and write it while // accumlating data for averaging + Eigen::Vector3d testread; + while (in && restart) { in.getline(cbuffer, cbufsiz); @@ -169,6 +172,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 >> testread(k); + } + if (allRead) { + if (accel) accel->add(time, testread); + } } cout << " -- Orient: current log=" << logfile << " backup=" << backupfile << endl; @@ -244,12 +255,15 @@ 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 << 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<26; i++) out << "| " << setw(13) << icnt++; out << endl; out.close(); @@ -660,6 +674,10 @@ void Orient::accumulate(double time, Component *c) } else center = center1; + if (accel) { + accel->add(time, center); + } + // Increment initial center according // to user specified velocity center0 += cenvel0*dtime; @@ -732,6 +750,10 @@ void Orient::logEntry(double time, Component *c) // Columns 22 - 24 for (int k=0; k<3; k++) outl << setw(15) << c->com0[k]; + // Columns 25 - 27 + pseudo = accel->accel(); + for (int k=0; k<3; k++) outl << setw(15) << pseudo[k]; + outl << endl; } } diff --git a/src/user/UserBar.cc b/src/user/UserBar.cc index 52773264..db3f6eee 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 e80dd789..1bee560a 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 09d774e3..b3b1a2c7 100644 --- a/src/user/UserHalo.cc +++ b/src/user/UserHalo.cc @@ -216,7 +216,7 @@ void * UserHalo::determine_acceleration_and_potential_thread(void * arg) // Add external accerlation 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 );