Permalink
Browse files

Merge pull request #1594 from KaiSzuttor/observable

Implementation of cylindrical observable
  • Loading branch information...
KaiSzuttor committed Nov 14, 2017
2 parents b6fa6a9 + 9acc310 commit 3f55d4dbfb6f45073b339a26efb3a10fbe552762
@@ -0,0 +1,88 @@
#ifndef OBSERVABLES_CYLINDRICALFLUXDENSITYPROFILE_HPP
#define OBSERVABLES_CYLINDRICALFLUXDENSITYPROFILE_HPP
#include "CylindricalProfileObservable.hpp"
#include "utils.hpp"
#include "utils/Histogram.hpp"
namespace Observables {
class CylindricalFluxDensityProfile : public CylindricalProfileObservable {
public:
virtual int actual_calculate(PartCfg &partCfg) override {
double bin_volume;
int r_bin, phi_bin, z_bin;
for (int id : ids()) {
auto const ppos = ::Vector<3, double>(folded_position(partCfg[id]));
auto const ppos_shifted = ppos - center;
auto const ppos_cyl = Utils::transform_to_cylinder_coordinates(ppos_shifted);
r_bin = std::floor((ppos_cyl[0] - min_r) / r_bin_size());
phi_bin = std::floor((ppos_cyl[1] - min_phi) / phi_bin_size());
z_bin = std::floor((ppos_cyl[2] - min_z) / z_bin_size());
bin_volume =
PI *
((min_r + (r_bin + 1) * r_bin_size()) *
(min_r + (r_bin + 1) * r_bin_size()) -
(min_r + r_bin * r_bin_size()) * (min_r + r_bin * r_bin_size())) *
z_bin_size() * phi_bin_size() / (2 * PI);
if (r_bin >= 0 && r_bin < n_r_bins && phi_bin >= 0 &&
phi_bin < n_phi_bins && z_bin >= 0 && z_bin < n_z_bins) {
// Coordinate transform the velocities and divide core velocities by time_step to get MD
// units.
// v_r = (x * v_x + y * v_y) / sqrt(x^2 + y^2)
double v_r = (ppos_shifted[0] * partCfg[id].m.v[0] / time_step +
ppos_shifted[1] * partCfg[id].m.v[1] / time_step) /
std::sqrt(ppos_shifted[0] * ppos_shifted[0] +
ppos_shifted[1] * ppos_shifted[1]);
// v_phi = (x * v_y - y * v_x ) / (x^2 + y^2)
double v_phi = (ppos_shifted[0] * partCfg[id].m.v[1] / time_step -
ppos_shifted[1] * partCfg[id].m.v[0] / time_step) /
(ppos_shifted[0] * ppos_shifted[0] +
ppos_shifted[1] * ppos_shifted[1]);
// v_z = v_z
double v_z = partCfg[id].m.v[2] / time_step;
// Write a flat histogram.
// index calculation: using the following formula for N dimensions:
// ind = ind_{N-1} + sum_{j=0}^{N-2} (ind_j * prod_{k=j+1}^{N-1} n_k),
// where ind is the flattened index, ind_i is the ith unflattened index
// and n_i is the size of the ith dimension.
int ind =
3 * (r_bin * n_phi_bins * n_z_bins + phi_bin * n_z_bins + z_bin);
if (std::isfinite(v_r)) {last_value[ind + 0] += v_r / bin_volume;}
if (std::isfinite(v_phi)) {last_value[ind + 1] += v_phi / bin_volume;}
if (std::isfinite(v_z)) {last_value[ind + 2] += v_z / bin_volume;}
}
}
return 0;
}
virtual int n_values() const override {
return 3 * n_r_bins * n_phi_bins * n_z_bins;
}
private:
virtual void do_write() override {
// We override the implementation to actually write positions not plain
// indices.
static const int len_dims[4] = {n_r_bins, n_phi_bins, n_z_bins, 3};
static const int n_dims = 4;
static const std::array<double, 3> bin_sizes = {
r_bin_size(), phi_bin_size(), z_bin_size()};
std::array<double, 3> position;
int index;
int unravelled_index[4];
for (auto it = last_value.begin(); it != last_value.end(); it += 3) {
index = std::distance(last_value.begin(), it);
::Utils::unravel_index(len_dims, n_dims, index, unravelled_index);
position = {
(static_cast<double>(unravelled_index[0]) + 0.5) * bin_sizes[0],
(static_cast<double>(unravelled_index[1]) + 0.5) * bin_sizes[1],
(static_cast<double>(unravelled_index[2]) + 0.5) * bin_sizes[2]};
m_ofile << position[0] << " " << position[1] << " " << position[2] << " "
<< *it << " " << *(it + 1) << " " << *(it + 2) << "\n";
}
m_ofile << std::endl;
}
};
} // Namespace Observables
#endif
@@ -0,0 +1,30 @@
#ifndef OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP
#define OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP
#include <cmath>
#include "PidObservable.hpp"
#include "Vector.hpp"
#include "integrate.hpp"
#include "particle_data.hpp"
namespace Observables {
class CylindricalProfileObservable : public PidObservable {
public:
::Vector<3, double> center;
double min_r, max_r;
double min_phi, max_phi;
double min_z, max_z;
// Number of bins for each coordinate.
int n_r_bins, n_phi_bins, n_z_bins;
double r_bin_size() { return (max_r - min_r) / n_r_bins; }
double phi_bin_size() { return (max_phi - min_phi) / n_phi_bins; }
double z_bin_size() { return (max_z - min_z) / n_z_bins; }
virtual int n_values() const override {
return n_r_bins * n_phi_bins * n_z_bins;
};
};
} // Namespace Observables
#endif
@@ -30,6 +30,7 @@
namespace Observables {
class Observable {
public:
friend class CylindricalFluxDensityProfile;
Observable();
virtual ~Observable() = default;
int calculate();
@@ -33,7 +33,7 @@ class PidObservable : public Observable {
public:
std::vector<int> &ids() { return m_ids; }
std::vector<int> const &ids() const { return m_ids; }
int n_values() const override { return 3 * ids().size(); }
virtual int n_values() const override { return 3 * ids().size(); }
};
} // Namespace Observables
@@ -18,6 +18,7 @@
#include <iostream> //for std::cout
#include <fstream> //for std::ifstream, std::ofstream for input output into files
#include "utils.hpp" // for PI and random vectors
#include "utils/Histogram.hpp"
#include "partCfg_global.hpp"
namespace ReactionEnsemble{
@@ -1166,7 +1167,7 @@ int ReactionEnsemble::initialize_wang_landau(){
for(int flattened_index=0;flattened_index<m_current_wang_landau_system.len_histogram;flattened_index++){
//unravel index
int unraveled_index[m_current_wang_landau_system.nr_collective_variables];
unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index,unraveled_index);
Utils::unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index,unraveled_index);
//use unraveled index
double current_energy=unraveled_index[energy_collective_variable_index]*m_current_wang_landau_system.collective_variables[energy_collective_variable_index]->delta_CV+m_current_wang_landau_system.collective_variables[energy_collective_variable_index]->CV_minimum;
if(current_energy>max_boundaries_energies[get_flattened_index_wang_landau_without_energy_collective_variable(flattened_index,energy_collective_variable_index)] || current_energy<min_boundaries_energies[get_flattened_index_wang_landau_without_energy_collective_variable(flattened_index,energy_collective_variable_index)]-m_current_wang_landau_system.collective_variables[energy_collective_variable_index]->delta_CV ){
@@ -1414,18 +1415,6 @@ bool ReactionEnsemble::achieved_desired_number_of_refinements_one_over_t() {
}
/**
*Returns the unraveled index of the provided flattened index (needed for writing the Wang-Landau results to file)
*/
void ReactionEnsemble::unravel_index(int* len_dims, int ndims, int flattened_index, int* unraveled_index_out){
//idea taken from http://codinghighway.com/2014/02/22/c-multi-dimensional-arrays-part-2-flattened-to-unflattened-index/
int mul[ndims];
mul[ndims-1]=1;
for (int j = ndims-2; j >= 0; j--)
mul[j] = mul[j+1]*len_dims[j+1];
for (int j = 0; j < ndims; j++)
unraveled_index_out[j]=(flattened_index/mul[j])%len_dims[j];
}
/**
@@ -1443,7 +1432,7 @@ void ReactionEnsemble::write_wang_landau_results_to_file(char* full_path_to_outp
//unravel index
if(std::abs(m_current_wang_landau_system.wang_landau_potential[flattened_index]-m_current_wang_landau_system.double_fill_value)>1){ //only output data if they are not equal to m_current_reaction_system.double_fill_value. This if ensures that for the energy observable not allowed energies (energies in the interval [global_E_min, global_E_max]) in the multidimensional wang landau potential are printed out, since the range [E_min(nbar), E_max(nbar)] for each nbar may be a different one
int unraveled_index[m_current_wang_landau_system.nr_collective_variables];
unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index,unraveled_index);
Utils::unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index,unraveled_index);
//use unraveled index
for(int i=0;i<m_current_wang_landau_system.nr_collective_variables;i++){
fprintf(pFile, "%f ",unraveled_index[i]*m_current_wang_landau_system.collective_variables[i]->delta_CV+m_current_wang_landau_system.collective_variables[i]->CV_minimum);
@@ -1500,7 +1489,7 @@ void ReactionEnsemble::write_out_preliminary_energy_run_results (char* full_path
for(int flattened_index=0;flattened_index<m_current_wang_landau_system.len_histogram;flattened_index++){
//unravel index
int unraveled_index[m_current_wang_landau_system.nr_collective_variables];
unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index,unraveled_index);
Utils::unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index,unraveled_index);
//use unraveled index
for(int i=0;i<m_current_wang_landau_system.nr_collective_variables;i++){
fprintf(pFile, "%f ",unraveled_index[i]*m_current_wang_landau_system.collective_variables[i]->delta_CV+m_current_wang_landau_system.collective_variables[i]->CV_minimum);
@@ -1520,7 +1509,7 @@ int ReactionEnsemble::get_flattened_index_wang_landau_without_energy_collective_
int* nr_subindices_of_collective_variable=m_current_wang_landau_system.nr_subindices_of_collective_variable;
//unravel index
int unraveled_index[m_current_wang_landau_system.nr_collective_variables];
unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index_with_energy_collective_variable,unraveled_index);
Utils::unravel_index(nr_subindices_of_collective_variable,m_current_wang_landau_system.nr_collective_variables,flattened_index_with_energy_collective_variable,unraveled_index);
//use unraveled index
const int nr_collective_variables=m_current_wang_landau_system.nr_collective_variables-1; //forget the last collective variable (the energy collective variable)
double current_state[nr_collective_variables];
@@ -296,10 +296,6 @@ class ReactionEnsemble {
int index_of_current_collective_variable);
int *initialize_histogram();
double *initialize_wang_landau_potential();
void unravel_index(int *len_dims, int ndims, int flattened_index,
int *unraveled_index_out); // needed for writing results
// and energy collective
// variable
void reset_histogram();
int m_WL_accepted_moves = 0;
int m_WL_tries = 0;
View
@@ -42,6 +42,7 @@
#include "utils/List.hpp"
#include "utils/math/sqr.hpp"
#include "utils/memory.hpp"
#include "Vector.hpp"
/*************************************************************/
/** \name Mathematical, physical and chemical constants. */
@@ -105,7 +106,16 @@ template <unsigned n, typename T> inline T int_pow(T x) {
/** Calculate signum of val, if supported by T */
template <typename T> int sgn(T val) { return (T(0) < val) - (val < T(0)); }
/** \brief Transform the given 3D Vector to cylinder coordinates.
*/
inline ::Vector<3, double>
transform_to_cylinder_coordinates(::Vector<3, double> const &pos) {
double r = std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]);
double phi = std::atan2(pos[1], pos[0]);
return ::Vector<3, double>{r, phi, pos[2]};
}
} // Namespace Utils
/*************************************************************/
/** \name List operations . */
@@ -233,7 +243,7 @@ inline void sort_int_array(int *data, int size) {
}
}
/** permute an interger array field of size size about permute positions. */
/** permute an integer array field of size size about permute positions. */
inline void permute_ifield(int *field, int size, int permute) {
int i, tmp;
@@ -1091,7 +1101,8 @@ void vecsub(T const *const a, T const *const b, T *const c) {
template <typename T> int sign(T value) {
return (T(0) < value) - (value < T(0));
}
}
}// namespace utils
/*@}*/
@@ -0,0 +1,20 @@
namespace Utils {
/**
* \brief Returns the unravelled index of the provided flat index.
* Therefore is the inversion of flattening an ndims dimensional index.
* @param len_dims an int array of length ndims containing the lengths of the dimensions. (Input)
* @param ndims int denoting the number of dimensions. (Input)
* @flattened_index an int denoting the flat index. (Input)
* @unravelled_index_out an int array with length ndims where the unflat indices are written to. (Output)
*/
inline void unravel_index(const int* const len_dims, const int ndims, const int flattened_index, int* unravelled_index_out){
//idea taken from http://codinghighway.com/2014/02/22/c-multi-dimensional-arrays-part-2-flattened-to-unflattened-index/
std::vector<int> mul(ndims);
mul[ndims-1]=1;
for (int j = ndims-2; j >= 0; j--)
mul[j] = mul[j+1]*len_dims[j+1];
for (int j = 0; j < ndims; j++)
unravelled_index_out[j]=(flattened_index/mul[j])%len_dims[j];
}
} // Namespace Utils
@@ -120,3 +120,8 @@ class StressTensor(Observable):
@script_interface_register
class StressTensorAcf(Observable):
_so_name="Observables::StressTensorAcf"
@script_interface_register
class CylindricalFluxDensityProfile(Observable):
_so_name="Observables::CylindricalFluxDensityProfile"
@@ -41,6 +41,9 @@ cdef extern from "utils.hpp":
cdef void alloc_intlist(int_list * il, int size)
cdef void realloc_intlist(int_list * il, int size)
cdef extern from "utils/Histogram.hpp" namespace "Utils":
cdef void unravel_index(const int* const len_dims, const int ndims, const int flattened_index, int* unravelled_index_out)
cdef int_list * create_int_list_from_python_object(obj)
cdef np.ndarray create_nparray_from_int_list(int_list * il)
cdef np.ndarray create_nparray_from_double_list(double_list * dl)
@@ -18,8 +18,10 @@
#
from __future__ import print_function, absolute_import
cimport numpy as np
cimport cython
import numpy as np
from cpython.version cimport PY_MAJOR_VERSION
from libcpp.vector cimport vector
cdef extern from "stdlib.h":
void free(void * ptr)
@@ -191,3 +193,37 @@ cdef handle_errors(msg):
# Cast because cython does not support typed enums completely
if <int> err.level() == <int> ERROR:
raise Exception(msg)
def get_unravelled_index(len_dims, n_dims, flattened_index):
"""
Getting the unravelled index for a given flattened index in ``n_dims`` dimensions.
Parameters
----------
len_dims : array_like :obj:`int`
The length of each of the ``n_dims`` dimensions.
n_dims : :obj:`int`
The number of dimensions.
flattened_index : :obj:`int`
The flat index that should be converted back to an
``n_dims`` dimensional index.
Returns
-------
unravelled_index : array_like :obj:`int`
An array containing the index for each dimension.
"""
cdef vector[int] c_len_dims
for i in range(len(len_dims)):
c_len_dims.push_back(len_dims[i])
cdef int c_n_dims = n_dims
cdef int c_flattened_index = flattened_index
cdef vector[int] unravelled_index_out
unravelled_index_out.assign(n_dims, 0)
unravel_index(c_len_dims.data(), c_n_dims, c_flattened_index, unravelled_index_out.data())
out = np.empty(n_dims)
for i in range(n_dims):
out[i] = unravelled_index_out[i]
return out
@@ -73,7 +73,7 @@ void transform_vectors(Variant &v) {
return;
}
/* only double, tranform to vector<int> */
/* only double, tranform to vector<double> */
if (std::all_of(variant_vector.begin(), variant_vector.end(), is_double)) {
v = to_vector<double>(variant_vector);
return;
Oops, something went wrong.

0 comments on commit 3f55d4d

Please sign in to comment.