Permalink
Browse files

Merge branch 'python' of github.com:espressomd/espresso into accumulator

  • Loading branch information...
KaiSzuttor committed Dec 7, 2017
2 parents 98080ee + 2f66734 commit 6640648413aba47498d4c9f80e605b1dda8045b5
Showing with 1,298 additions and 445 deletions.
  1. +182 −114 doc/sphinx/run.rst
  2. +1 −1 samples/observables_correlators.py
  3. +2 −4 src/core/correlators/Correlator.cpp
  4. +46 −46 src/core/lbgpu_cuda.cu
  5. +48 −53 src/core/observables/CylindricalFluxDensityProfile.hpp
  6. +95 −0 src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp
  7. +64 −0 src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.hpp
  8. +1 −0 src/core/observables/CylindricalProfileObservable.hpp
  9. +13 −16 src/core/observables/DensityProfile.hpp
  10. +14 −16 src/core/observables/FluxDensityProfile.hpp
  11. +14 −17 src/core/observables/ForceDensityProfile.hpp
  12. +19 −19 src/core/observables/LbRadialVelocityProfile.hpp
  13. +20 −18 src/core/observables/LbVelocityProfile.hpp
  14. +1 −0 src/core/observables/Observable.hpp
  15. +10 −10 src/core/observables/ProfileObservable.hpp
  16. +18 −18 src/core/observables/profiles.hpp
  17. +0 −1 src/core/reaction_ensemble.cpp
  18. +1 −0 src/core/unit_tests/CMakeLists.txt
  19. +37 −0 src/core/unit_tests/histogram.cpp
  20. +273 −15 src/core/utils/Histogram.hpp
  21. +5 −9 src/python/espressomd/analyze.pyx
  22. +46 −22 src/python/espressomd/correlators.py
  23. +65 −0 src/python/espressomd/observables.py
  24. +4 −2 src/python/espressomd/reaction_ensemble.pyx
  25. +18 −9 src/script_interface/observables/CylindricalProfileObservable.hpp
  26. +27 −27 src/script_interface/observables/ProfileObservable.hpp
  27. +1 −0 src/script_interface/observables/initialize.cpp
  28. +2 −0 testsuite/CMakeLists.txt
  29. +49 −28 testsuite/observable_cylindricalfluxdensity.py
  30. +121 −0 testsuite/observable_cylindricallbfluxdensity.py
  31. +65 −0 testsuite/observable_profile.py
  32. +36 −0 testsuite/tests_common.py
View

Large diffs are not rendered by default.

Oops, something went wrong.
@@ -31,7 +31,7 @@
corr_operation="square_distance_componentwise", compress1="discard1")
# Instance a correlator calculating the FCS autocorrelation function from particle positions, using the symmetric focal spot with wx=wy=wz=10 (sigma)
fcs = Correlator(tau_lin=16, tau_max=10000, dt=0.1, obs1=p,
corr_operation="fcs_acf", args=[10, 10, 10], compress1="discard1")
corr_operation="fcs_acf", args=[10, 10, 10], compress1="discard2")
# Ask the correlator for its parameters
print c.get_params()
@@ -829,16 +829,14 @@ int componentwise_product ( double* A, unsigned int dim_A, double* B, unsigned i
}
int complex_conjugate_product ( double* A, unsigned int dim_A, double* B, unsigned int dim_B, double* C, unsigned int dim_corr, Vector3d ) {
unsigned int i,j;
unsigned int j;
if (!(dim_A == dim_B )) {
printf("Error in complex_conjugate product: The vector sizes do not match");
return 5;
}
j=0;
for ( i = 0; i < dim_A/2; i++ ) {
for ( j = 0; j < dim_A; j+=2 ) {
C[j] = A[j]*B[j] + A[j+1]*B[j+1];
C[j+1] = A[j+1]*B[j] - A[j]*B[j+1];
j=j+2;
}
return 0;
}
View
@@ -3804,34 +3804,34 @@ __global__ void fill_lb_radial_velocity_profile(LB_nodes_gpu n_a, radial_profile
unsigned int phibin=blockIdx.x;
unsigned int zbin=blockIdx.y;
float roffset=pdata->minr;
float r_incr=(pdata->maxr-pdata->minr)/(pdata->rbins-1);
float roffset=pdata->min_r;
float r_incr=(pdata->max_r-pdata->min_r)/(pdata->n_r_bins-1);
float r = roffset + rbin*r_incr;
unsigned int maxj;
float phioffset, phi_incr;
if ( pdata->phibins == 1 ) {
maxj = (int)floorf( 2*3.1415f*pdata->maxr/para.agrid ) ;
if ( pdata->n_phi_bins == 1 ) {
maxj = (int)floorf( 2*3.1415f*pdata->max_r/para.agrid ) ;
phioffset=0;
phi_incr=2*3.1415f/maxj;
} else {
maxj = pdata->phibins;
phioffset=pdata->minphi;
phi_incr=(pdata->maxphi-pdata->minphi)/(pdata->phibins);
maxj = pdata->n_phi_bins;
phioffset=pdata->min_phi;
phi_incr=(pdata->max_phi-pdata->min_phi)/(pdata->n_phi_bins);
}
float phi = phioffset + phibin*phi_incr;
unsigned int maxk;
float zoffset, z_incr;
if ( pdata->zbins == 1 ) {
if ( pdata->n_z_bins == 1 ) {
maxk = (int) para.dim_z;
zoffset=-pdata->center[2];
z_incr=para.agrid;
} else {
maxk = (int) pdata->zbins;
zoffset=pdata->minz;
z_incr=(pdata->maxz-pdata->minz)/(pdata->zbins-1);
maxk = (int) pdata->n_z_bins;
zoffset=pdata->min_z;
z_incr=(pdata->max_z-pdata->min_z)/(pdata->n_z_bins-1);
}
float z = zoffset + zbin*z_incr;
@@ -3873,34 +3873,34 @@ __global__ void fill_lb_velocity_profile(LB_nodes_gpu n_a, profile_data* pdata,
if ( pdata->xbins == 1 ) {
if ( pdata->n_x_bins == 1 ) {
/* maxi = (int) floor(gridDim.x/para.agrid); */
xoffset=0;
x_incr=para.agrid;
} else {
/* maxi = pdata->xbins; */
xoffset=pdata->minx;
x_incr=(pdata->maxx-pdata->minx)/(pdata->xbins-1);
/* maxi = pdata->n_x_bins; */
xoffset=pdata->min_x;
x_incr=(pdata->max_x-pdata->min_x)/(pdata->n_x_bins-1);
}
float x = xoffset + xbin*x_incr;
if ( pdata->ybins == 1 ) {
if ( pdata->n_y_bins == 1 ) {
maxj = (int) floorf(para.dim_y/para.agrid);
yoffset=0;
y_incr=para.agrid;
} else {
maxj = pdata->ybins;
yoffset=pdata->miny;
y_incr=(pdata->maxy-pdata->miny)/(pdata->ybins-1);
maxj = pdata->n_y_bins;
yoffset=pdata->min_y;
y_incr=(pdata->max_y-pdata->min_y)/(pdata->n_y_bins-1);
}
float y = yoffset + ybin*y_incr;
if ( pdata->zbins == 1 ) {
if ( pdata->n_z_bins == 1 ) {
maxk = (int) floorf(para.dim_z/para.agrid);
zoffset=0;
z_incr=para.agrid;
} else {
maxk = (int) pdata->zbins;
zoffset=pdata->minz;
z_incr=(pdata->maxz-pdata->minz)/(pdata->zbins-1);
maxk = (int) pdata->n_z_bins;
zoffset=pdata->min_z;
z_incr=(pdata->max_z-pdata->min_z)/(pdata->n_z_bins-1);
}
float z = zoffset + zbin*z_incr;
@@ -3928,23 +3928,23 @@ int statistics_observable_lbgpu_radial_velocity_profile(radial_profile_data* pda
unsigned int maxj, maxk;
float normalization_factor=1;
if ( pdata->rbins == 1 ) {
if ( pdata->n_r_bins == 1 ) {
return 1;
}
unsigned int maxi=pdata->rbins;
unsigned int maxi=pdata->n_r_bins;
if ( pdata->phibins == 1 ) {
maxj = (int)floorf( 2*3.1415f*pdata->maxr/lbpar_gpu.agrid ) ;
if ( pdata->n_phi_bins == 1 ) {
maxj = (int)floorf( 2*3.1415f*pdata->max_r/lbpar_gpu.agrid ) ;
normalization_factor/=maxj;
} else {
maxj = pdata->phibins;
maxj = pdata->n_phi_bins;
}
if ( pdata->zbins == 1 ) {
if ( pdata->n_z_bins == 1 ) {
maxk = (int) lbpar_gpu.dim_z;
normalization_factor/=maxk;
} else {
maxk = pdata->zbins;
maxk = pdata->n_z_bins;
}
for (int i = 0; i<n_A; i++) {
@@ -3980,11 +3980,11 @@ int statistics_observable_lbgpu_radial_velocity_profile(radial_profile_data* pda
for (int j =0; j<maxj; j++)
for (int k =0; k<maxk; k++) {
linear_index = 0;
if (pdata->rbins > 1)
linear_index += i*pdata->phibins*pdata->zbins;
if (pdata->phibins > 1)
linear_index += j*pdata->zbins;
if (pdata->zbins > 1)
if (pdata->n_r_bins > 1)
linear_index += i*pdata->n_phi_bins*pdata->n_z_bins;
if (pdata->n_phi_bins > 1)
linear_index += j*pdata->n_z_bins;
if (pdata->n_z_bins > 1)
linear_index +=k;
A[3*linear_index+0]+=host_data[3*(i*maxj*maxk + j*maxk + k)+0]*normalization_factor*lbpar_gpu.tau/lbpar_gpu.agrid;
A[3*linear_index+1]+=host_data[3*(i*maxj*maxk + j*maxk + k)+1]*normalization_factor*lbpar_gpu.tau/lbpar_gpu.agrid;
@@ -4011,23 +4011,23 @@ int statistics_observable_lbgpu_velocity_profile(profile_data* pdata, double* A,
float normalization_factor=1;
if ( pdata->xbins == 1 ) {
if ( pdata->n_x_bins == 1 ) {
maxi = (int) floor(lbpar_gpu.dim_x/lbpar_gpu.agrid);
normalization_factor/=maxi;
} else {
maxi = pdata->xbins;
maxi = pdata->n_x_bins;
}
if ( pdata->ybins == 1 ) {
if ( pdata->n_y_bins == 1 ) {
maxj = (int) floor(lbpar_gpu.dim_y/lbpar_gpu.agrid);
normalization_factor/=maxj;
} else {
maxj = pdata->ybins;
maxj = pdata->n_y_bins;
}
if ( pdata->zbins == 1 ) {
if ( pdata->n_z_bins == 1 ) {
maxk = (int) floor(lbpar_gpu.dim_z/lbpar_gpu.agrid);
normalization_factor/=maxk;
} else {
maxk = pdata->zbins;
maxk = pdata->n_z_bins;
}
for (int i = 0; i<n_A; i++) {
@@ -4065,11 +4065,11 @@ int statistics_observable_lbgpu_velocity_profile(profile_data* pdata, double* A,
for ( j = 0; j < maxj; j++ ) {
for ( k = 0; k < maxk; k++ ) {
linear_index = 0;
if (pdata->xbins > 1)
linear_index += i*pdata->ybins*pdata->zbins;
if (pdata->ybins > 1)
linear_index += j*pdata->zbins;
if (pdata->zbins > 1)
if (pdata->n_x_bins > 1)
linear_index += i*pdata->n_y_bins*pdata->n_z_bins;
if (pdata->n_y_bins > 1)
linear_index += j*pdata->n_z_bins;
if (pdata->n_z_bins > 1)
linear_index +=k;
A[3*linear_index+0]+=host_data[3*(i*maxj*maxk + j*maxk + k)+0]*normalization_factor*lbpar_gpu.tau/lbpar_gpu.agrid;
@@ -9,58 +9,53 @@ namespace Observables {
class CylindricalFluxDensityProfile : public CylindricalProfileObservable {
public:
virtual std::vector<double> operator()(PartCfg &partCfg) const override {
std::vector<double> res(n_values());
double bin_volume;
int r_bin, phi_bin, z_bin;
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 histogram(n_bins, 3, limits);
for (int id : ids()) {
auto const ppos = ::Vector<3, double>(folded_position(partCfg[id]));
auto const ppos_shifted = ppos - center;
::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);
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)) {
res[ind + 0] += v_r / bin_volume;
}
if (std::isfinite(v_phi)) {
res[ind + 1] += v_phi / bin_volume;
}
if (std::isfinite(v_z)) {
res[ind + 2] += v_z / bin_volume;
}
}
// 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}});
}
return res;
histogram.normalize();
return histogram.get_histogram();
}
virtual int n_values() const override {
return 3 * n_r_bins * n_phi_bins * n_z_bins;
@@ -72,19 +67,19 @@ class CylindricalFluxDensityProfile : public CylindricalProfileObservable {
// 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()}};
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];
std::vector<double> tmp = operator()(partCfg());
for (auto it = tmp.begin(); it != tmp.end(); it += 3) {
index = std::distance(tmp.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]}};
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";
}
Oops, something went wrong.

0 comments on commit 6640648

Please sign in to comment.