Permalink
Browse files

Merge pull request #1613 from KaiSzuttor/observable_lb

Observable lb
  • Loading branch information...
fweik committed Dec 6, 2017
2 parents 29189a1 + 962850e commit 44b3f95a80331004e03405866f6618677f4144ab
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 44b3f95

Please sign in to comment.