Cylinder transformation and observables

doc/sphinx/analysis.rst

@@ -411,10 +411,21 @@ or bin edges for the axes. Example::
                        density_profile.min_y, density_profile.max_y])
     plt.show()
 
+Observables based on cylindrical coordinates are also available.
+They require special parameters if the cylindrical coordinate system is non-standard, e.g. if you want the origin of the cylindrical coordinates to be at a special location of the box or if you want to make use of symmetries along an axis that is not parallel to the z-axis.
+For this purpose, use :class:`espressomd.math.CylindricalTransformationParameters` to create a consistent set of the parameters needed. Example::
+
+    import espressomd.math
+
+    # shifted and rotated cylindrical coordinates
+    cyl_transform_params = espressomd.math.CylindricalTransformationParameters(
+        center=[5.0, 5.0, 0.0], axis=[0, 1, 0], orientation=[0, 0, 1])
+
     # histogram in cylindrical coordinates
     density_profile = espressomd.observables.CylindricalDensityProfile(
-        ids=[0, 1], center=[5.0, 5.0, 0.0], axis=[0, 0, 1],
-        n_r_bins=8, min_r=0.0, max_r=4.0,
+        ids=[0, 1],
+        transform_params = cyl_transform_params,
+        n_r_bins=8, min_r=1.0, max_r=4.0,
         n_phi_bins=16, min_phi=-np.pi, max_phi=np.pi,
         n_z_bins=4, min_z=4.0, max_z=8.0)
     obs_data = density_profile.calculate()
@@ -779,4 +790,3 @@ Note that the cluster objects do not contain copies of the particles, but refer
 
 
 
-

doc/tutorials/charged_system/charged_system-1.ipynb

@@ -36,7 +36,7 @@
     "import espressomd\n",
     "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n",
     "\n",
-    "from espressomd import System, interactions, electrostatics, observables, accumulators\n",
+    "from espressomd import System, interactions, electrostatics, observables, accumulators, math\n",
     "\n",
     "import numpy as np\n",
     "from scipy import optimize\n",
@@ -438,20 +438,18 @@
     "```python\n",
     "def setup_profile_calculation(system, delta_N, ion_types, r_min, n_radial_bins):\n",
     "    radial_profile_accumulators = {}\n",
+    "    ctp = math.CylindricalTransformationParameters(center = np.array(system.box_l) / 2.,\n",
+    "                                                   axis = [0, 0, 1],\n",
+    "                                                   orientation = [1, 0, 0])\n",
     "    for ion_type in ion_types:\n",
     "        ion_ids = system.part.select(type=ion_type).id\n",
     "        radial_profile_obs = observables.CylindricalDensityProfile(\n",
     "            ids=ion_ids,\n",
-    "            center=np.array(system.box_l) / 2.,\n",
-    "            axis=[0, 0, 1, ],\n",
+    "            transform_params = ctp,\n",
     "            n_r_bins=n_radial_bins,\n",
-    "            n_phi_bins=1,\n",
-    "            n_z_bins=1,\n",
     "            min_r=r_min,\n",
-    "            min_phi=-np.pi,\n",
     "            min_z=-system.box_l[2] / 2.,\n",
     "            max_r=system.box_l[0] / 2.,\n",
-    "            max_phi=np.pi,\n",
     "            max_z=system.box_l[2] / 2.)\n",
     "\n",
     "        bin_edges = radial_profile_obs.bin_edges()\n",
@@ -945,7 +943,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.9"
+   "version": "3.8.5"
   }
  },
  "nbformat": 4,

samples/lb_profile.py

@@ -32,6 +32,7 @@
 import espressomd.shapes
 import espressomd.lbboundaries
 import espressomd.accumulators
+import espressomd.math
 
 system = espressomd.System(box_l=[10.0, 10.0, 5.0])
 system.time_step = 0.01
@@ -42,16 +43,14 @@
     agrid=1.0, dens=1.0, visc=1.0, tau=0.01, ext_force_density=[0, 0, 0.15], kT=1.0, seed=32)
 system.actors.add(lb_fluid)
 system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23)
-fluid_obs = espressomd.observables.CylindricalLBVelocityProfile(
+ctp = espressomd.math.CylindricalTransformationParameters(
     center=[5.0, 5.0, 0.0],
     axis=[0, 0, 1],
+    orientation=[1, 0, 0])
+fluid_obs = espressomd.observables.CylindricalLBVelocityProfile(
+    transform_params=ctp,
     n_r_bins=100,
-    n_phi_bins=1,
-    n_z_bins=1,
-    min_r=0.0,
     max_r=4.0,
-    min_phi=-np.pi,
-    max_phi=np.pi,
     min_z=0.0,
     max_z=10.0,
     sampling_density=0.1)

src/core/grid_based_algorithms/lb_collective_interface.cpp

@@ -81,6 +81,15 @@ mpi_lb_get_interpolated_velocity(Utils::Vector3d const &pos) {
 
 REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_interpolated_velocity)
 
+boost::optional<double>
+mpi_lb_get_interpolated_density(Utils::Vector3d const &pos) {
+  return detail::lb_calc_for_pos(pos, [&](auto pos) {
+    return lb_lbinterpolation_get_interpolated_density(pos);
+  });
+}
+
+REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_interpolated_density)
+
 auto mpi_lb_get_density(Utils::Vector3i const &index) {
   return detail::lb_calc_fluid_kernel(
       index, [&](auto const &modes, auto const &force_density) {

src/core/grid_based_algorithms/lb_collective_interface.hpp

@@ -27,6 +27,8 @@
 /* collective getter functions */
 boost::optional<Utils::Vector3d>
 mpi_lb_get_interpolated_velocity(Utils::Vector3d const &pos);
+boost::optional<double>
+mpi_lb_get_interpolated_density(Utils::Vector3d const &pos);
 boost::optional<double> mpi_lb_get_density(Utils::Vector3i const &index);
 boost::optional<Utils::Vector19d>
 mpi_lb_get_populations(Utils::Vector3i const &index);

src/core/grid_based_algorithms/lb_interface.cpp

@@ -1273,3 +1273,23 @@ lb_lbfluid_get_interpolated_velocity(const Utils::Vector3d &pos) {
   }
   throw NoLBActive();
 }
+
+double lb_lbfluid_get_interpolated_density(const Utils::Vector3d &pos) {
+  auto const folded_pos = folded_position(pos, box_geo);
+  auto const interpolation_order = lb_lbinterpolation_get_interpolation_order();
+  if (lattice_switch == ActiveLB::GPU) {
+    throw std::runtime_error(
+        "Density interpolation is not implemented for the GPU LB.");
+  }
+  if (lattice_switch == ActiveLB::CPU) {
+    switch (interpolation_order) {
+    case (InterpolationOrder::quadratic):
+      throw std::runtime_error("The non-linear interpolation scheme is not "
+                               "implemented for the CPU LB.");
+    case (InterpolationOrder::linear):
+      return mpi_call(::Communication::Result::one_rank,
+                      mpi_lb_get_interpolated_density, folded_pos);
+    }
+  }
+  throw NoLBActive();
+}

src/core/grid_based_algorithms/lb_interface.hpp

@@ -265,4 +265,11 @@ Utils::Vector3d lb_lbfluid_calc_fluid_momentum();
 const Utils::Vector3d
 lb_lbfluid_get_interpolated_velocity(const Utils::Vector3d &pos);
 
+/**
+ * @brief Calculates the interpolated fluid density on the master process.
+ * @param pos Position at which the density is to be calculated.
+ * @retval interpolated fluid density.
+ */
+double lb_lbfluid_get_interpolated_density(const Utils::Vector3d &pos);
+
 #endif

src/core/grid_based_algorithms/lb_interpolation.cpp

@@ -82,6 +82,16 @@ Utils::Vector3d node_u(Lattice::index_t index) {
   return Utils::Vector3d{modes[1], modes[2], modes[3]} / local_density;
 }
 
+double node_dens(Lattice::index_t index) {
+#ifdef LB_BOUNDARIES
+  if (lbfields[index].boundary) {
+    return lbpar.density;
+  }
+#endif // LB_BOUNDARIES
+  auto const modes = lb_calc_modes(index, lbfluid);
+  return lbpar.density + modes[0];
+}
+
 } // namespace
 
 const Utils::Vector3d
@@ -98,6 +108,19 @@ lb_lbinterpolation_get_interpolated_velocity(const Utils::Vector3d &pos) {
   return interpolated_u;
 }
 
+double lb_lbinterpolation_get_interpolated_density(const Utils::Vector3d &pos) {
+  double interpolated_dens = 0.;
+
+  /* Calculate fluid density at the position.
+     This is done by linear interpolation (eq. (11) @cite ahlrichs99a) */
+  lattice_interpolation(lblattice, pos,
+                        [&interpolated_dens](Lattice::index_t index, double w) {
+                          interpolated_dens += w * node_dens(index);
+                        });
+
+  return interpolated_dens;
+}
+
 void lb_lbinterpolation_add_force_density(
     const Utils::Vector3d &pos, const Utils::Vector3d &force_density) {
   switch (interpolation_order) {

src/core/grid_based_algorithms/lb_interpolation.hpp

@@ -41,6 +41,13 @@ InterpolationOrder lb_lbinterpolation_get_interpolation_order();
 const Utils::Vector3d
 lb_lbinterpolation_get_interpolated_velocity(const Utils::Vector3d &p);
 
+/**
+ * @brief Calculates the fluid density at a given position of the
+ * lattice.
+ * @note It can lead to undefined behaviour if the
+ * position is not within the local lattice. */
+double lb_lbinterpolation_get_interpolated_density(const Utils::Vector3d &p);
+
 /**
  * @brief Add a force density to the fluid at the given position.
  */

src/core/observables/CylindricalDensityProfile.hpp

@@ -41,7 +41,8 @@ class CylindricalDensityProfile : public CylindricalPidProfileObservable {
 
     for (auto p : particles) {
       histogram.update(Utils::transform_coordinate_cartesian_to_cylinder(
-          folded_position(traits.position(p), box_geo) - center, axis));
+          folded_position(traits.position(p), box_geo) - trafo_params->center(),
+          trafo_params->axis(), trafo_params->orientation()));
     }
 
     histogram.normalize();

src/core/observables/CylindricalFluxDensityProfile.hpp

@@ -43,11 +43,13 @@ class CylindricalFluxDensityProfile : public CylindricalPidProfileObservable {
 
     // Write data to the histogram
     for (auto p : particles) {
-      auto const pos = folded_position(traits.position(p), box_geo) - center;
+      auto const pos =
+          folded_position(traits.position(p), box_geo) - trafo_params->center();
       histogram.update(
-          Utils::transform_coordinate_cartesian_to_cylinder(pos, axis),
-          Utils::transform_vector_cartesian_to_cylinder(traits.velocity(p),
-                                                        axis, pos));
+          Utils::transform_coordinate_cartesian_to_cylinder(
+              pos, trafo_params->axis(), trafo_params->orientation()),
+          Utils::transform_vector_cartesian_to_cylinder(
+              traits.velocity(p), trafo_params->axis(), pos));
     }
     histogram.normalize();
     return histogram.get_histogram();

src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp

@@ -42,13 +42,23 @@ CylindricalLBFluxDensityProfileAtParticlePositions::evaluate(
     auto const pos = folded_position(traits.position(p), box_geo);
     auto const v = lb_lbfluid_get_interpolated_velocity(pos) *
                    lb_lbfluid_get_lattice_speed();
+    auto const flux_dens = lb_lbfluid_get_interpolated_density(pos) * v;
 
     histogram.update(
-        Utils::transform_coordinate_cartesian_to_cylinder(pos - center, axis),
-        Utils::transform_vector_cartesian_to_cylinder(v, axis, pos - center));
+        Utils::transform_coordinate_cartesian_to_cylinder(
+            pos - trafo_params->center(), trafo_params->axis(),
+            trafo_params->orientation()),
+        Utils::transform_vector_cartesian_to_cylinder(
+            flux_dens, trafo_params->axis(), pos - trafo_params->center()));
   }
 
-  histogram.normalize();
-  return histogram.get_histogram();
+  // normalize by number of hits per bin
+  auto hist_tmp = histogram.get_histogram();
+  auto tot_count = histogram.get_tot_count();
+  std::transform(hist_tmp.begin(), hist_tmp.end(), tot_count.begin(),
+                 hist_tmp.begin(), [](auto hi, auto ci) {
+                   return ci > 0 ? hi / static_cast<double>(ci) : 0.;
+                 });
+  return hist_tmp;
 }
 } // namespace Observables

src/core/observables/CylindricalLBProfileObservable.hpp

@@ -21,6 +21,7 @@
 
 #include "CylindricalProfileObservable.hpp"
 
+#include <utility>
 #include <utils/Vector.hpp>
 #include <utils/math/coordinate_transformation.hpp>
 #include <utils/math/vec_rotate.hpp>
@@ -30,15 +31,14 @@ namespace Observables {
 
 class CylindricalLBProfileObservable : public CylindricalProfileObservable {
 public:
-  CylindricalLBProfileObservable(Utils::Vector3d const &center,
-                                 Utils::Vector3d const &axis, int n_r_bins,
-                                 int n_phi_bins, int n_z_bins, double min_r,
-                                 double max_r, double min_phi, double max_phi,
-                                 double min_z, double max_z,
-                                 double sampling_density)
-      : CylindricalProfileObservable(center, axis, n_r_bins, n_phi_bins,
-                                     n_z_bins, min_r, max_r, min_phi, max_phi,
-                                     min_z, max_z),
+  CylindricalLBProfileObservable(
+      std::shared_ptr<Utils::CylindricalTransformationParameters> trafo_params,
+      int n_r_bins, int n_phi_bins, int n_z_bins, double min_r, double max_r,
+      double min_phi, double max_phi, double min_z, double max_z,
+      double sampling_density)
+      : CylindricalProfileObservable(std::move(trafo_params), n_r_bins,
+                                     n_phi_bins, n_z_bins, min_r, max_r,
+                                     min_phi, max_phi, min_z, max_z),
         sampling_density(sampling_density) {
     calculate_sampling_positions();
   }
@@ -47,17 +47,16 @@ class CylindricalLBProfileObservable : public CylindricalProfileObservable {
         limits[0], limits[1], limits[2], n_bins[0], n_bins[1], n_bins[2],
         sampling_density);
     for (auto &p : sampling_positions) {
-      double theta;
-      Utils::Vector3d rotation_axis;
-      auto p_cart = Utils::transform_coordinate_cylinder_to_cartesian(
-          p, Utils::Vector3d{{0.0, 0.0, 1.0}});
+      auto p_cart = Utils::transform_coordinate_cylinder_to_cartesian(p);
       // We have to rotate the coordinates since the utils function assumes
       // z-axis symmetry.
-      std::tie(theta, rotation_axis) =
-          Utils::rotation_params(Utils::Vector3d{{0.0, 0.0, 1.0}}, axis);
+      constexpr Utils::Vector3d z_axis{{0.0, 0.0, 1.0}};
+      auto const theta = Utils::angle_between(z_axis, trafo_params->axis());
+      auto const rot_axis =
+          Utils::vector_product(z_axis, trafo_params->axis()).normalize();
       if (theta > std::numeric_limits<double>::epsilon())
-        p_cart = Utils::vec_rotate(rotation_axis, theta, p_cart);
-      p = p_cart + center;
+        p_cart = Utils::vec_rotate(rot_axis, theta, p_cart);
+      p = p_cart + trafo_params->center();
     }
   }
   std::vector<Utils::Vector3d> sampling_positions;

src/core/observables/CylindricalLBVelocityProfile.cpp

@@ -35,11 +35,11 @@ std::vector<double> CylindricalLBVelocityProfile::operator()() const {
   for (auto const &p : sampling_positions) {
     auto const velocity = lb_lbfluid_get_interpolated_velocity(p) *
                           lb_lbfluid_get_lattice_speed();
-    auto const pos_shifted = p - center;
-    auto const pos_cyl =
-        Utils::transform_coordinate_cartesian_to_cylinder(pos_shifted, axis);
+    auto const pos_shifted = p - trafo_params->center();
+    auto const pos_cyl = Utils::transform_coordinate_cartesian_to_cylinder(
+        pos_shifted, trafo_params->axis(), trafo_params->orientation());
     histogram.update(pos_cyl, Utils::transform_vector_cartesian_to_cylinder(
-                                  velocity, axis, pos_shifted));
+                                  velocity, trafo_params->axis(), pos_shifted));
   }
   auto hist_data = histogram.get_histogram();
   auto const tot_count = histogram.get_tot_count();

src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp

@@ -41,17 +41,20 @@ std::vector<double> CylindricalLBVelocityProfileAtParticlePositions::evaluate(
                    lb_lbfluid_get_lattice_speed();
 
     histogram.update(
-        Utils::transform_coordinate_cartesian_to_cylinder(pos - center, axis),
-        Utils::transform_vector_cartesian_to_cylinder(v, axis, pos - center));
+        Utils::transform_coordinate_cartesian_to_cylinder(
+            pos - trafo_params->center(), trafo_params->axis(),
+            trafo_params->orientation()),
+        Utils::transform_vector_cartesian_to_cylinder(
+            v, trafo_params->axis(), pos - trafo_params->center()));
   }
 
+  // normalize by number of hits per bin
   auto hist_tmp = histogram.get_histogram();
   auto tot_count = histogram.get_tot_count();
-  for (size_t ind = 0; ind < hist_tmp.size(); ++ind) {
-    if (tot_count[ind] > 0) {
-      hist_tmp[ind] /= static_cast<double>(tot_count[ind]);
-    }
-  }
+  std::transform(hist_tmp.begin(), hist_tmp.end(), tot_count.begin(),
+                 hist_tmp.begin(), [](auto hi, auto ci) {
+                   return ci > 0 ? hi / static_cast<double>(ci) : 0.;
+                 });
   return hist_tmp;
 }