Skip to content

Commit

Permalink
Implement pseudo-forces for 'Orient' moving 'CENTER' non-inertial frames
Browse files Browse the repository at this point in the history
  • Loading branch information
Martin D. Weinberg committed Dec 19, 2024
1 parent febd3b7 commit 46e7b2e
Show file tree
Hide file tree
Showing 10 changed files with 229 additions and 22 deletions.
30 changes: 29 additions & 1 deletion expui/BasisFactory.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include <ParticleReader.H>
#include <OrthoFunction.H>
#include <Coefficients.H>
#include <PseudoAccel.H>
#include <YamlCheck.H>
#include <localmpi.H>
#include <exputils.H>
Expand Down Expand Up @@ -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");

Expand Down Expand Up @@ -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<Basis>;
Expand Down
93 changes: 93 additions & 0 deletions expui/BasisFactory.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<char[]> cbuf(new char [cbufsiz]);

// Look for data and write it while
// accumlating data for averaging
Eigen::Vector3d testread;
double time, dummy;

std::vector<double> times;
std::vector<Eigen::Vector3d> 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<times.size(); i++) {
t_accel(i) = times[i];
for (int k=0; k<3; k++) p_accel(i, k) = centers[i](k);
}
}


Eigen::Vector3d Basis::currentAccel(double time)
{
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;
}
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<Eigen::VectorXd>(t, p.col(k)).coefs());
}
return ret;
}

}
// END namespace BasisClasses
8 changes: 6 additions & 2 deletions expui/BiorthBasis.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -3653,7 +3653,7 @@ namespace BasisClasses
for (int n=0; n<rows; n++) {
auto v = basis->getFields(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;
Expand Down Expand Up @@ -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<BasisCoef> mod)
Expand Down
37 changes: 36 additions & 1 deletion src/Component.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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__);
Expand All @@ -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__);
Expand All @@ -877,6 +903,15 @@ public:

//! Add to accerlation (by vector)
inline void AddAcc(int i, vector<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 accerlation (by vector)
inline void AddAccExt(int i, vector<double>& val) {
PartMap::iterator tp = particles.find(i);
if (tp == particles.end()) {
throw BadIndexException(i, particles.size(), __FILE__, __LINE__);
Expand Down
16 changes: 15 additions & 1 deletion src/Component.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ const std::set<std::string> Component::valid_keys_parm =
"EJ",
"nEJkeep",
"nEJwant",
"nEJaccel",
"EJkinE",
"EJext",
"EJdiag",
Expand Down Expand Up @@ -184,6 +185,7 @@ Component::Component(YAML::Node& CONF)
EJ = 0;
nEJkeep = 100;
nEJwant = 500;
nEJaccel = 0;
EJkinE = true;
EJext = false;
EJdiag = false;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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<int>();
if (cconf["nEJwant" ]) nEJwant = cconf["nEJwant" ].as<int>();
if (cconf["nEJaccel"]) nEJaccel = cconf["nEJaccel"].as<int>();
if (cconf["EJx0" ]) EJx0 = cconf["EJx0" ].as<double>();
if (cconf["EJy0" ]) EJy0 = cconf["EJy0" ].as<double>();
if (cconf["EJz0" ]) EJz0 = cconf["EJz0" ].as<double>();
Expand Down Expand Up @@ -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];
Expand All @@ -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) {
Expand Down Expand Up @@ -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;

Expand All @@ -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(&center[0], 3) = orient->currentCenter();
Eigen::VectorXd::Map(&accel [0], 3) = orient->currentAccel();
} else {
if (EJlinear) orient -> set_linear();
if (not com_system) {
Expand Down Expand Up @@ -1238,6 +1248,7 @@ Component::~Component(void)

delete [] com;
delete [] center;
delete [] accel;
delete [] cov;
delete [] coa;
delete [] angmom;
Expand Down Expand Up @@ -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;
Expand Down
Loading

0 comments on commit 46e7b2e

Please sign in to comment.