Permalink
Browse files

observables: Implementation of cylindrical velocity profile for parti…

…cles and LB fluid.
  • Loading branch information...
KaiSzuttor committed Jan 11, 2018
1 parent 6f0fadf commit 2a275a3b08f617603ad1f592aa06b29234286012
@@ -0,0 +1,104 @@
/*
Copyright (C) 2016,2017 The ESPResSo project
This file is part of ESPResSo.
ESPResSo is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ESPResSo is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "CylindricalLBVelocityProfileAtParticlePositions.hpp"
#include "lbgpu.hpp"
#include "utils.hpp"
#include "utils/Histogram.hpp"
namespace Observables {
std::vector<double> CylindricalLBVelocityProfileAtParticlePositions::
operator()(PartCfg &partCfg) const {
std::vector<size_t> n_bins{{static_cast<size_t>(n_r_bins),
static_cast<size_t>(n_phi_bins),
static_cast<size_t>(n_z_bins)}};
std::vector<std::pair<double, double>> limits{
{std::make_pair(min_r, max_r), std::make_pair(min_phi, max_phi),
std::make_pair(min_z, max_z)}};
Utils::CylindricalHistogram<double> histogram(n_bins, 3, limits);
#ifdef LB_GPU
// First collect all positions (since we want to call the LB function to
// get the fluid velocities only once).
std::vector<double> ppos(3 * ids().size());
std::vector<double> ppos_cyl(3 * ids().size());
std::vector<double> ppos_shifted(3 * ids().size());
for (int index = 0; index < ids().size(); ++index) {
auto const ppos_tmp =
::Vector<3, double>(folded_position(partCfg[ids()[index]]));
ppos[3 * index + 0] = ppos_tmp[0];
ppos[3 * index + 1] = ppos_tmp[1];
ppos[3 * index + 2] = ppos_tmp[2];
auto const ppos_shifted_tmp = ppos_tmp - center;
if (axis == "z") {
ppos_shifted[3 * index + 0] = ppos_shifted_tmp[0];
ppos_shifted[3 * index + 1] = ppos_shifted_tmp[1];
ppos_shifted[3 * index + 2] = ppos_shifted_tmp[2];
} else if (axis == "x") {
ppos_shifted[3 * index + 0] = -ppos_shifted_tmp[2];
ppos_shifted[3 * index + 1] = ppos_shifted_tmp[1];
ppos_shifted[3 * index + 2] = ppos_shifted_tmp[0];
} else if (axis == "y") {
ppos_shifted[3 * index + 0] = ppos_shifted_tmp[0];
ppos_shifted[3 * index + 1] = -ppos_shifted_tmp[2];
ppos_shifted[3 * index + 2] = ppos_shifted_tmp[1];
}
auto const ppos_cyl_tmp =
Utils::transform_to_cylinder_coordinates(ppos_shifted_tmp);
ppos_cyl[3 * index + 0] = ppos_cyl_tmp[0];
ppos_cyl[3 * index + 1] = ppos_cyl_tmp[1];
ppos_cyl[3 * index + 2] = ppos_cyl_tmp[2];
}
std::vector<double> velocities(3 * ids().size());
lb_lbfluid_get_interpolated_velocity_at_positions(
ppos.data(), velocities.data(), ids().size());
for (int index = 0; index < ids().size(); ++index) {
// 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[3 * index + 0] * velocities[3 * index + 0] +
ppos_shifted[3 * index + 1] * velocities[3 * index + 1]) /
std::sqrt(ppos_shifted[3 * index + 0] * ppos_shifted[3 * index + 0] +
ppos_shifted[3 * index + 1] * ppos_shifted[3 * index + 1]);
// v_phi = (x * v_y - y * v_x ) / (x^2 + y^2)
double v_phi = (ppos_shifted[3 * index + 0] * velocities[3 * index + 1] -
ppos_shifted[3 * index + 1] * velocities[3 * index + 0]) /
(ppos_shifted[3 * index + 0] * ppos_shifted[3 * index + 0] +
ppos_shifted[3 * index + 1] * ppos_shifted[3 * index + 1]);
// v_z = v_z
double v_z = velocities[3 * index + 2];
histogram.update(
std::vector<double>{{ppos_cyl[3 * index + 0], ppos_cyl[3 * index + 1],
ppos_cyl[3 * index + 2]}},
std::vector<double>{{v_r, v_phi, v_z}});
}
auto hist_tmp = histogram.get_histogram();
auto tot_count = histogram.get_tot_count();
for (size_t ind=0; ind < hist_tmp.size(); ++ind) {
if (hist_tmp[ind] > 0.0) {
hist_tmp[ind] /= tot_count[ind];
}
}
return hist_tmp;
#endif // LB_GPU
#ifndef LB_GPU
return histogram.get_histogram();
#endif
}
} // namespace Observables
@@ -0,0 +1,39 @@
/*
Copyright (C) 2016,2017 The ESPResSo project
This file is part of ESPResSo.
ESPResSo is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ESPResSo is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OBSERVABLES_CYLINDRICALLBVELOCITYPROFILEATPARTICLEPOSITIONS_HPP
#define OBSERVABLES_CYLINDRICALLBVELOCITYPROFILEATPARTICLEPOSITIONS_HPP
#include "CylindricalProfileObservable.hpp"
#include "partCfg_global.hpp"
#include "utils/Histogram.hpp"
namespace Observables {
class CylindricalLBVelocityProfileAtParticlePositions
: public CylindricalProfileObservable {
public:
virtual std::vector<double> operator()(PartCfg &partCfg) const override;
virtual int n_values() const override {
return 3 * n_r_bins * n_phi_bins * n_z_bins;
}
};
} // Namespace Observables
#endif
@@ -0,0 +1,73 @@
#ifndef OBSERVABLES_CYLINDRICALVELOCITYPROFILE_HPP
#define OBSERVABLES_CYLINDRICALVELOCITYPROFILE_HPP
#include "CylindricalProfileObservable.hpp"
#include "utils.hpp"
#include "utils/Histogram.hpp"
namespace Observables {
class CylindricalVelocityProfile : public CylindricalProfileObservable {
public:
virtual std::vector<double> operator()(PartCfg &partCfg) const override {
std::vector<size_t> n_bins{{static_cast<size_t>(n_r_bins),
static_cast<size_t>(n_phi_bins),
static_cast<size_t>(n_z_bins)}};
std::vector<std::pair<double, double>> limits{
{std::make_pair(min_r, max_r), std::make_pair(min_phi, max_phi),
std::make_pair(min_z, max_z)}};
Utils::CylindricalHistogram<double> histogram(n_bins, 3, limits);
for (int id : ids()) {
auto const ppos = ::Vector<3, double>(folded_position(partCfg[id]));
::Vector<3, double> ppos_shifted;
ppos_shifted = ppos - center;
::Vector<3, double> vel = {partCfg[id].m.v[0], partCfg[id].m.v[1],
partCfg[id].m.v[2]};
if (axis == "x") {
// x' = -z, y' = y, z'= x
ppos_shifted = ::Vector<3, double>{-ppos_shifted[2], ppos_shifted[1],
ppos_shifted[0]};
vel = {-vel[2], vel[1], vel[0]};
} else if (axis == "y") {
// x' = x, y' = -z, z' = y
ppos_shifted = ::Vector<3, double>{ppos_shifted[0], -ppos_shifted[2],
ppos_shifted[1]};
vel = {vel[0], -vel[2], vel[1]};
}
auto const ppos_cyl =
Utils::transform_to_cylinder_coordinates(ppos_shifted);
// 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] * vel[0] / time_step +
ppos_shifted[1] * vel[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] * vel[1] / time_step -
ppos_shifted[1] * vel[0] / time_step) /
(ppos_shifted[0] * ppos_shifted[0] +
ppos_shifted[1] * ppos_shifted[1]);
// v_z = v_z
double v_z = vel[2] / time_step;
// Write data to the histogram.
histogram.update(ppos_cyl, std::vector<double>{{v_r, v_phi, v_z}});
}
auto hist_tmp = histogram.get_histogram();
auto tot_count = histogram.get_tot_count();
for (size_t ind=0; ind < hist_tmp.size(); ++ind) {
if (hist_tmp[ind] > 0.0) {
hist_tmp[ind] /= tot_count[ind];
}
}
return hist_tmp;
}
virtual int n_values() const override {
return 3 * n_r_bins * n_phi_bins * n_z_bins;
}
};
} // Namespace Observables
#endif
@@ -110,7 +110,7 @@ template <typename T> class Histogram {
std::vector<std::pair<T, T>> limits);
std::vector<size_t> get_n_bins() const;
std::vector<T> get_histogram() const;
size_t get_tot_count() const;
std::vector<size_t> get_tot_count() const;
std::vector<std::pair<T, T>> get_limits() const;
std::vector<T> get_bin_sizes() const;
void update(std::vector<T> const &data);
@@ -131,8 +131,8 @@ template <typename T> class Histogram {
std::vector<T> m_hist;
// Number of dimensions for a single data point.
size_t m_n_dims_data;
// Track the number of total hits.
size_t m_tot_count;
// Track the number of total hits per bin entry.
std::vector<size_t> m_tot_count;
};
/**
@@ -146,7 +146,7 @@ template <typename T> class Histogram {
template <typename T>
Histogram<T>::Histogram(std::vector<size_t> n_bins, size_t n_dims_data,
std::vector<std::pair<T, T>> limits)
: m_n_bins(n_bins), m_limits(limits), m_n_dims_data(n_dims_data), m_tot_count(0) {
: m_n_bins(n_bins), m_limits(limits), m_n_dims_data(n_dims_data) {
if (n_bins.size() != limits.size()) {
throw std::invalid_argument("Argument for number of bins and limits do "
"not have same number of dimensions!");
@@ -156,6 +156,7 @@ Histogram<T>::Histogram(std::vector<size_t> n_bins, size_t n_dims_data,
m_n_dims_data * std::accumulate(std::begin(n_bins), std::end(n_bins), 1,
std::multiplies<size_t>());
m_hist = std::vector<T>(n_bins_total);
m_tot_count = std::vector<size_t>(n_bins_total);
}
/**
@@ -192,8 +193,8 @@ void Histogram<T>::update(std::vector<T> const &data,
throw std::invalid_argument("Wrong dimensions of given weights!");
for (size_t ind = 0; ind < m_n_dims_data; ++ind) {
m_hist[flat_index + ind] += weights[ind];
m_tot_count[flat_index + ind] += 1;
}
m_tot_count += 1;
}
}
@@ -229,7 +230,7 @@ template <typename T> std::vector<T> Histogram<T>::get_histogram() const {
/**
* \brief Get the histogram count data.
*/
template <typename T> size_t Histogram<T>::get_tot_count() const {
template <typename T> std::vector<size_t> Histogram<T>::get_tot_count() const {
return m_tot_count;
}
/**
@@ -258,7 +259,6 @@ class CylindricalHistogram : public Histogram<T> {
using Histogram<T>::get_limits;
using Histogram<T>::get_bin_sizes;
using Histogram<T>::m_hist;
using Histogram<T>::m_tot_count;
using Histogram<T>::m_n_dims_data;
private:
@@ -386,3 +386,73 @@ class CylindricalLBFluxDensityProfileAtParticlePositions(Observable):
"""
_so_name = "Observables::CylindricalLBFluxDensityProfileAtParticlePositions"
@script_interface_register
class CylindricalLBVelocityProfileAtParticlePositions(Observable):
"""Calculates the LB fluid velocity at the particle positions in polar coordinates.
Parameters
----------
ids : array_like of :obj:`int`
The ids of (existing) particles to take into account.
center : array_like of :obj:`float`
Position of the center of the polar coordinate system for the histogram.
axis : :obj:`str` (``x``, ``y``, or ``z``)
Orientation of the ``z``-axis of the polar coordinate system for the histogram.
n_r_bins : :obj:`int`
Number of bins in radial direction.
n_phi_bins : :obj:`int`
Number of bins for the azimuthal direction.
n_z_bins : :obj:`int`
Number of bins in ``z`` direction.
min_r : :obj:`float`
Minimum ``r`` to consider.
min_phi : :obj:`float`
Minimum ``phi`` to consider.
min_z : :obj:`float`
Minimum ``z`` to consider.
max_r : :obj:`float`
Maximum ``r`` to consider.
max_phi : :obj:`float`
Maximum ``phi`` to consider.
max_z : :obj:`float`
Maximum ``z`` to consider.
"""
_so_name = "Observables::CylindricalLBVelocityProfileAtParticlePositions"
@script_interface_register
class CylindricalVelocityProfile(Observable):
"""Calculates the LB particle velocity profile in polar coordinates.
Parameters
----------
ids : array_like of :obj:`int`
The ids of (existing) particles to take into account.
center : array_like of :obj:`float`
Position of the center of the polar coordinate system for the histogram.
axis : :obj:`str` (``x``, ``y``, or ``z``)
Orientation of the ``z``-axis of the polar coordinate system for the histogram.
n_r_bins : :obj:`int`
Number of bins in radial direction.
n_phi_bins : :obj:`int`
Number of bins for the azimuthal direction.
n_z_bins : :obj:`int`
Number of bins in ``z`` direction.
min_r : :obj:`float`
Minimum ``r`` to consider.
min_phi : :obj:`float`
Minimum ``phi`` to consider.
min_z : :obj:`float`
Minimum ``z`` to consider.
max_r : :obj:`float`
Maximum ``r`` to consider.
max_phi : :obj:`float`
Maximum ``phi`` to consider.
max_z : :obj:`float`
Maximum ``z`` to consider.
"""
_so_name = "Observables::CylindricalVelocityProfile"
@@ -28,8 +28,10 @@
#include "Observable.hpp"
#include "core/observables/CylindricalDensityProfile.hpp"
#include "core/observables/CylindricalVelocityProfile.hpp"
#include "core/observables/CylindricalFluxDensityProfile.hpp"
#include "core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.hpp"
#include "core/observables/CylindricalLBVelocityProfileAtParticlePositions.hpp"
#include "core/observables/CylindricalProfileObservable.hpp"
namespace ScriptInterface {
@@ -108,9 +110,12 @@ class CylindricalProfileObservable : public Observable {
};
NEW_CYLINDRICAL_PROFILE_OBSERVABLE(CylindricalDensityProfile)
NEW_CYLINDRICAL_PROFILE_OBSERVABLE(CylindricalVelocityProfile)
NEW_CYLINDRICAL_PROFILE_OBSERVABLE(CylindricalFluxDensityProfile)
NEW_CYLINDRICAL_PROFILE_OBSERVABLE(
CylindricalLBFluxDensityProfileAtParticlePositions)
NEW_CYLINDRICAL_PROFILE_OBSERVABLE(
CylindricalLBVelocityProfileAtParticlePositions)
} /* namespace Observables */
} /* namespace ScriptInterface */
@@ -77,8 +77,10 @@ void initialize() {
REGISTER(FluxDensityProfile);
REGISTER(LBVelocityProfile);
REGISTER(CylindricalDensityProfile);
REGISTER(CylindricalVelocityProfile);
REGISTER(CylindricalFluxDensityProfile);
REGISTER(CylindricalLBFluxDensityProfileAtParticlePositions);
REGISTER(CylindricalLBVelocityProfileAtParticlePositions);
#undef REGISTER
#undef REGISTER_PID_OBS

0 comments on commit 2a275a3

Please sign in to comment.