Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Modernize core #3314

Merged
merged 6 commits into from
Nov 15, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions src/core/BoxGeometry.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,12 @@ class BoxGeometry {
* @param box_l Length that should be set.
*/
void set_length(Utils::Vector3d const &box_l) { m_length = box_l; }

/**
* @brief Box volume
* @return Return the volume of the box.
*/
double volume() const { return m_length[0] * m_length[1] * m_length[2]; }
};

#endif // CORE_BOX_GEOMETRY_HPP
100 changes: 45 additions & 55 deletions src/core/accumulators/Correlator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -190,18 +190,15 @@ constexpr const char init_errors[][64] = {

int Correlator::get_correlation_time(double *correlation_time) {
// We calculate the correlation time for each m_dim_corr by normalizing the
// correlation,
// integrating it and finding out where C(tau)=tau;
double C_tau;
int ok_flag;
// correlation, integrating it and finding out where C(tau)=tau
for (unsigned j = 0; j < m_dim_corr; j++) {
correlation_time[j] = 0.;
}

// here we still have to fix the stuff a bit!
for (unsigned j = 0; j < m_dim_corr; j++) {
C_tau = 1 * m_dt;
ok_flag = 0;
double C_tau = 1 * m_dt;
bool ok_flag = false;
for (unsigned k = 1; k < m_n_result - 1; k++) {
if (n_sweeps[k] == 0)
break;
Expand All @@ -215,7 +212,7 @@ int Correlator::get_correlation_time(double *correlation_time) {
2 * sqrt(tau[k - 1] * m_dt / n_data)) {
correlation_time[j] =
C_tau * (1 + (2 * (double)tau[k] + 1) / (double)n_data);
ok_flag = 1;
ok_flag = true;
break;
}
}
Expand All @@ -228,7 +225,6 @@ int Correlator::get_correlation_time(double *correlation_time) {
}

void Correlator::initialize() {
unsigned int i, j, k;
hierarchy_depth = 0;
// Class members are assigned via the initializer list

Expand Down Expand Up @@ -260,8 +256,9 @@ void Correlator::initialize() {
dim_A = 0;
dim_B = 0;

if (A_obs)
if (A_obs) {
dim_A = A_obs->n_values();
}
if (!B_obs) {
B_obs = A_obs;
}
Expand Down Expand Up @@ -352,24 +349,26 @@ void Correlator::initialize() {

result.resize(std::array<int, 2>{{m_n_result, m_dim_corr}});

for (i = 0; i < m_n_result; i++) {
for (j = 0; j < m_dim_corr; j++)
for (int i = 0; i < m_n_result; i++) {
for (int j = 0; j < m_dim_corr; j++) {
// and initialize the values
result[i][j] = 0;
}
}

newest = std::vector<unsigned int>(hierarchy_depth, m_tau_lin);

tau.resize(m_n_result);
for (i = 0; i < m_tau_lin + 1; i++) {
for (int i = 0; i < m_tau_lin + 1; i++) {
tau[i] = i;
}

for (j = 1; j < hierarchy_depth; j++)
for (k = 0; k < m_tau_lin / 2; k++) {
for (int j = 1; j < hierarchy_depth; j++) {
for (int k = 0; k < m_tau_lin / 2; k++) {
tau[m_tau_lin + 1 + (j - 1) * m_tau_lin / 2 + k] =
(k + (m_tau_lin / 2) + 1) * (1 << j);
}
}
}

void Correlator::update() {
Expand All @@ -378,20 +377,15 @@ void Correlator::update() {
"No data can be added after finalize() was called.");
}
// We must now go through the hierarchy and make sure there is space for the
// new
// datapoint. For every hierarchy level we have to decide if it necessary to
// move
// something
int i, j;
int highest_level_to_compress;
unsigned int index_new, index_old, index_res;
// new datapoint. For every hierarchy level we have to decide if it necessary
// to move something
int highest_level_to_compress = -1;

t++;

highest_level_to_compress = -1;
i = 0;
// Lets find out how far we have to go back in the hierarchy to make space for
// the new value
// Let's find out how far we have to go back in the hierarchy to make space
// for the new value
int i = 0;
while (true) {
if (((t - ((m_tau_lin + 1) * ((1 << (i + 1)) - 1) + 1)) % (1 << (i + 1)) ==
0)) {
Expand All @@ -405,9 +399,9 @@ void Correlator::update() {
}

// Now we know we must make space on the levels 0..highest_level_to_compress
// Now lets compress the data level by level.
// Now let's compress the data level by level.

for (i = highest_level_to_compress; i >= 0; i--) {
for (int i = highest_level_to_compress; i >= 0; i--) {
// We increase the index indicating the newest on level i+1 by one (plus
// folding)
newest[i + 1] = (newest[i + 1] + 1) % (m_tau_lin + 1);
Expand Down Expand Up @@ -441,9 +435,9 @@ void Correlator::update() {
}

// Now update the lowest level correlation estimates
for (j = 0; j < min(m_tau_lin + 1, n_vals[0]); j++) {
index_new = newest[0];
index_old = (newest[0] - j + m_tau_lin + 1) % (m_tau_lin + 1);
for (unsigned j = 0; j < min(m_tau_lin + 1, n_vals[0]); j++) {
auto const index_new = newest[0];
auto const index_old = (newest[0] - j + m_tau_lin + 1) % (m_tau_lin + 1);
auto const temp =
(corr_operation)(A[0][index_old], B[0][index_new], m_correlation_args);
assert(temp.size() == m_dim_corr);
Expand All @@ -457,9 +451,9 @@ void Correlator::update() {
for (int i = 1; i < highest_level_to_compress + 2; i++) {
for (unsigned j = (m_tau_lin + 1) / 2 + 1;
j < min(m_tau_lin + 1, n_vals[i]); j++) {
index_new = newest[i];
index_old = (newest[i] - j + m_tau_lin + 1) % (m_tau_lin + 1);
index_res =
auto const index_new = newest[i];
auto const index_old = (newest[i] - j + m_tau_lin + 1) % (m_tau_lin + 1);
auto const index_res =
m_tau_lin + (i - 1) * m_tau_lin / 2 + (j - m_tau_lin / 2 + 1) - 1;
auto const temp = (corr_operation)(A[i][index_old], B[i][index_new],
m_correlation_args);
Expand All @@ -480,20 +474,15 @@ int Correlator::finalize() {
throw std::runtime_error("Correlator::finalize() can only be called once.");
}
// We must now go through the hierarchy and make sure there is space for the
// new
// datapoint. For every hierarchy level we have to decide if it necessary to
// move
// something
int i, j;
int ll; // current lowest level
int vals_ll; // number of values remaining in the lowest level
// new datapoint. For every hierarchy level we have to decide if it necessary
// to move something
int highest_level_to_compress;
unsigned int index_new, index_old, index_res;

// make a flag that the correlation is finalized
finalized = 1;
// mark the correlation as finalized
finalized = true;

for (ll = 0; ll < hierarchy_depth - 1; ll++) {
for (int ll = 0; ll < hierarchy_depth - 1; ll++) {
int vals_ll; // number of values remaining in the lowest level
if (n_vals[ll] > m_tau_lin + 1)
vals_ll = m_tau_lin + n_vals[ll] % 2;
else
Expand All @@ -507,9 +496,9 @@ int Correlator::finalize() {
highest_level_to_compress = -1;
}

i = ll + 1; // lowest level, for which we have to check for compression
// Lets find out how far we have to go back in the hierarchy to make space
// for the new value
int i = ll + 1; // lowest level for which we have to check for compression
// Let's find out how far we have to go back in the hierarchy to make
// space for the new value
while (highest_level_to_compress > -1) {
if (n_vals[i] % 2) {
if (i < (hierarchy_depth - 1) && n_vals[i] > m_tau_lin) {
Expand All @@ -526,9 +515,9 @@ int Correlator::finalize() {

// Now we know we must make space on the levels
// 0..highest_level_to_compress
// Now lets compress the data level by level.
// Now let's compress the data level by level.

for (i = highest_level_to_compress; i >= ll; i--) {
for (int i = highest_level_to_compress; i >= ll; i--) {
// We increase the index indicating the newest on level i+1 by one (plus
// folding)
newest[i + 1] = (newest[i + 1] + 1) % (m_tau_lin + 1);
Expand All @@ -542,12 +531,13 @@ int Correlator::finalize() {
newest[ll] = (newest[ll] + 1) % (m_tau_lin + 1);

// We only need to update correlation estimates for the higher levels
for (i = ll + 1; i < highest_level_to_compress + 2; i++) {
for (j = (m_tau_lin + 1) / 2 + 1; j < min(m_tau_lin + 1, n_vals[i]);
for (int i = ll + 1; i < highest_level_to_compress + 2; i++) {
for (int j = (m_tau_lin + 1) / 2 + 1; j < min(m_tau_lin + 1, n_vals[i]);
j++) {
index_new = newest[i];
index_old = (newest[i] - j + m_tau_lin + 1) % (m_tau_lin + 1);
index_res =
auto const index_new = newest[i];
auto const index_old =
(newest[i] - j + m_tau_lin + 1) % (m_tau_lin + 1);
auto const index_res =
m_tau_lin + (i - 1) * m_tau_lin / 2 + (j - m_tau_lin / 2 + 1) - 1;

auto const temp = (corr_operation)(A[i][index_old], B[i][index_new],
Expand All @@ -569,7 +559,7 @@ std::vector<double> Correlator::get_correlation() {
std::vector<double> res;

// time + n_sweeps + corr_1...corr_n
int cols = 2 + m_dim_corr;
int const cols = 2 + m_dim_corr;
res.resize(m_n_result * cols);

for (int i = 0; i < m_n_result; i++) {
Expand Down
6 changes: 3 additions & 3 deletions src/core/accumulators/Correlator.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ class Correlator : public AccumulatorBase {
Correlator(int tau_lin, double tau_max, int delta_N, std::string compress1_,
std::string compress2_, std::string corr_operation, obs_ptr obs1,
obs_ptr obs2, Utils::Vector3d correlation_args_ = {})
: AccumulatorBase(delta_N), finalized(0), t(0),
: AccumulatorBase(delta_N), finalized(false), t(0),
m_correlation_args(correlation_args_), m_tau_lin(tau_lin),
m_dt(delta_N * time_step), m_tau_max(tau_max),
compressA_name(std::move(compress1_)),
Expand Down Expand Up @@ -233,8 +233,8 @@ class Correlator : public AccumulatorBase {
void set_internal_state(std::string const &);

private:
unsigned int finalized; ///< flag whether non-zero of correlation is finalized
unsigned int t; ///< global time in number of frames
bool finalized; ///< whether the correlation is finalized
unsigned int t; ///< global time in number of frames

Utils::Vector3d m_correlation_args; ///< additional arguments, which the
///< correlation may need (currently
Expand Down
2 changes: 1 addition & 1 deletion src/core/cells.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ CellPList ghost_cells = {nullptr, 0, 0};
/** Type of cell structure in use */
CellStructure cell_structure;

/** On of Cells::Resort, announces the level of resort needed.
/** One of @ref Cells::Resort, announces the level of resort needed.
*/
unsigned resort_particles = Cells::RESORT_NONE;
bool rebuild_verletlist = true;
Expand Down
53 changes: 19 additions & 34 deletions src/core/cells.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -145,19 +145,15 @@ struct CellStructure {

bool use_verlet_list = true;

/** Maximal pair range supported by current
* cell system.
*/
/** Maximal pair range supported by current cell system. */
Utils::Vector3d max_range = {};

/**
* Minimum range that has to be supported.
*/
/** Minimum range that has to be supported. */
double min_range;

/** returns the global local_cells */
/** Return the global local_cells */
CellPList local_cells() const;
/** returns the global ghost_cells */
/** Return the global ghost_cells */
CellPList ghost_cells() const;

/** Communicator to exchange ghost particles. */
Expand All @@ -168,7 +164,7 @@ struct CellStructure {
/** Cell system dependent function to find the right cell for a
* particle.
* \param p Particle.
* \return pointer to cell where to put the particle, nullptr
* \return pointer to cell where to put the particle, nullptr
* if the particle does not belong on this node.
*/
Cell *(*particle_to_cell)(const Particle &p) = nullptr;
Expand All @@ -184,18 +180,16 @@ struct CellStructure {
/** list of all cells. */
extern std::vector<Cell> cells;

/** list of all cells containing particles physically on the local
node */
/** list of all cells containing particles physically on the local node */
extern CellPList local_cells;
/** list of all cells containing ghosts */
extern CellPList ghost_cells;

/** Type of cell structure in use ( \ref Cell Structure ). */
/** Type of cell structure in use. */
extern CellStructure cell_structure;

/** If non-zero, cell systems should reset the position for checking
* the Verlet criterion. Moreover, the Verlet list has to be
* rebuilt.
* the Verlet criterion. Moreover, the Verlet list has to be rebuilt.
*/
extern bool rebuild_verletlist;

Expand All @@ -207,9 +201,9 @@ extern bool rebuild_verletlist;
/*@{*/

/** Reinitialize the cell structures.
* @param new_cs gives the new topology to use afterwards. May be set to
* \ref CELL_STRUCTURE_CURRENT for not changing it.
* @param range Desired interaction range
* @param new_cs The new topology to use afterwards. May be set to
* @ref CELL_STRUCTURE_CURRENT for not changing it.
* @param range Desired interaction range
*/
void cells_re_init(int new_cs, double range);

Expand All @@ -233,9 +227,6 @@ inline void realloc_cellplist(CellPList *cpl, int size) {

/** Sort the particles into the cells and initialize the ghost particle
* structures.
* @param global_flag if this is CELLS_GLOBAL_EXCHANGE, particle positions can
* have changed arbitrarily, otherwise the change should
* have been smaller then skin.
*/
void cells_resort_particles(int global_flag);

Expand All @@ -244,43 +235,37 @@ void cells_resort_particles(int global_flag);
* It calculates the maximal interaction range, and as said reinitializes
* the cells structure if something significant has changed.
*
* If bit CELL_FLAG_FAST is set, the routine should try to save time.
* If bit @ref CELL_FLAG_FAST is set, the routine should try to save time.
* Currently this means that if the maximal range decreased, it does
* not reorganize the particles. This is used in the NpT algorithm to
* avoid frequent reorganization of particles.
*
* If bit CELL_FLAG_GRIDCHANGED is set, it means the nodes' topology
* If bit @ref CELL_FLAG_GRIDCHANGED is set, it means the nodes' topology
* has changed, i. e. the grid or periodicity. In this case a full
* reorganization is due.
*
* @param flags a bitmask of CELL_FLAG_GRIDCHANGED,
* and/or CELL_FLAG_FAST, see above.
* @param flags a bitmask of @ref CELL_FLAG_GRIDCHANGED,
* and/or @ref CELL_FLAG_FAST, see above.
*/
void cells_on_geometry_change(int flags);

/** Update ghost information. If \ref resort_particles == 1,
* also a resorting of the particles takes place.
/** Update ghost information. If @ref resort_particles is not
* @ref Cells::RESORT_NONE, the particles are also resorted.
*/
void cells_update_ghosts();

/** Calculate and return the total number of particles on this node. */
int cells_get_n_particles();

/**
* @brief Get pairs closer than distance from the cells.
*
* This is mostly for testing purposes and uses link_cell
* to get pairs out of the cellsystem by a simple distance
* criterion.
* @brief Get pairs closer than @p distance from the cells.
*
* Pairs are sorted so that first.id < second.id
*/
std::vector<std::pair<int, int>> mpi_get_pairs(double distance);

/**
* @brief Increase the local resort level at least to level.
*
* The changed level has to be communicated via annouce_resort_particles.
* @brief Increase the local resort level at least to @p level.
*/
void set_resort_particles(Cells::Resort level);

Expand Down
Loading