Pressure tensor with generated kernel

src/core/grid_based_algorithms/lb_interface.cpp

@@ -531,11 +531,12 @@ lb_lbnode_get_pressure_tensor(const Utils::Vector3i &ind) {
         ::Communication::Result::one_rank, Walberla::get_node_pressure_tensor,
         ind);
 
-    // reverts the correction done by walberla
-    //    auto const revert_factor =
-    //        lb_lbfluid_get_viscosity() / (lb_lbfluid_get_viscosity() + 1.0
-    //        / 6.0); stress[1] /= revert_factor; stress[3] /= revert_factor;
-    //        stress[4] /= revert_factor;
+    // correction factor
+    auto const adjust_factor =
+            lb_lbfluid_get_viscosity() / (lb_lbfluid_get_viscosity() + 1.0 / 6.0);
+    stress[1] *= adjust_factor;
+    stress[3] *= adjust_factor;
+    stress[4] *= adjust_factor;
 
     return stress;
   }

src/core/grid_based_algorithms/lb_interface.hpp

@@ -247,3 +247,4 @@ void lb_lbfluid_add_force_at_pos(const Utils::Vector3d &pos,
 double lb_lbfluid_get_interpolated_density(const Utils::Vector3d &pos);
 
 #endif
+

src/core/grid_based_algorithms/lb_walberla_instance.cpp

@@ -110,3 +110,4 @@ void mpi_destruct_lb_walberla() {
   Communication::mpiCallbacks().call_all(destruct_lb_walberla);
 }
 #endif
+

src/core/grid_based_algorithms/lb_walberla_instance.hpp

@@ -64,3 +64,4 @@ void mpi_destruct_lb_walberla();
 #endif // LB_WALBERLA
 
 #endif
+

src/python/espressomd/lb.pxd

@@ -210,3 +210,4 @@ cdef inline python_lbfluid_set_gamma(p_gamma) except +:
         c_gamma = p_gamma
         # call c-function
         lb_lbcoupling_set_gamma(c_gamma)
+

src/python/espressomd/lb.pyx

@@ -396,7 +396,6 @@ IF LB_WALBERLA:
         """
         Initialize the lattice-Boltzmann method for hydrodynamic flow using waLBerla.
         See :class:`HydrodynamicInteraction` for the list of parameters.
-
         """
 
         def _set_params_in_es_core(self):
@@ -615,7 +614,6 @@ def _add_lb_slice_properties():
         """
         Setter function that sets attribute on every member of lb_slice.
         If values contains only one element, all members are set to it.
-
         """
 
         indices = [lb_slice.x_indices, lb_slice.y_indices, lb_slice.z_indices]
@@ -669,3 +667,4 @@ def _add_lb_slice_properties():
 
 
 _add_lb_slice_properties()
+

src/walberla_bridge/LBWalberlaBase.hpp

@@ -179,3 +179,4 @@ enum class OutputVTK : unsigned {
 };
 
 #endif // LB_WALBERLA_H
+

src/walberla_bridge/LBWalberlaD3Q19FluctuatingMRT.hpp

@@ -75,3 +75,4 @@ class LBWalberlaD3Q19FluctuatingMRT : public LBWalberlaImpl<LatticeModelName> {
 
 } // namespace walberla
 #undef LatticeModelName
+

src/walberla_bridge/LBWalberlaD3Q19MRT.hpp

@@ -49,3 +49,4 @@ class LBWalberlaD3Q19MRT : public LBWalberlaImpl<LatticeModelName> {
 } // namespace walberla
 
 #undef LatticeModelName
+

src/walberla_bridge/generated_kernels/Dynamic_UBB.cpp

@@ -19,9 +19,11 @@
 
 #include <cmath>
 
-#include "Dynamic_UBB.h"
 #include "core/DataTypes.h"
 #include "core/Macros.h"
+#include "Dynamic_UBB.h"
+
+
 
 #define FUNC_PREFIX
 
@@ -43,102 +45,38 @@ namespace lbm {
 #endif
 
 namespace internal_dynamic_ubb_boundary_Dynamic_UBB {
-static FUNC_PREFIX void dynamic_ubb_boundary_Dynamic_UBB(
-    uint8_t *RESTRICT const _data_indexVector, double *RESTRICT _data_pdfs,
-    int64_t const _stride_pdfs_0, int64_t const _stride_pdfs_1,
-    int64_t const _stride_pdfs_2, int64_t const _stride_pdfs_3,
-    int64_t indexVectorSize) {
-
-  const int32_t f_in_inv_dir_idx[] = {0, 2,  1,  4,  3,  6,  5,  10, 9, 8,
-                                      7, 16, 15, 18, 17, 12, 11, 14, 13};
-  const int32_t f_out_offsets_x[] = {0, 0, 0, -1, 1, 0, 0, -1, 1, -1,
-                                     1, 0, 0, -1, 1, 0, 0, -1, 1};
-  const int32_t f_out_offsets_y[] = {0,  1, -1, 0, 0, 0, 0,  1, 1, -1,
-                                     -1, 1, -1, 0, 0, 1, -1, 0, 0};
-  const int32_t f_out_offsets_z[] = {0, 0, 0, 0, 0, 1,  -1, 0,  0, 0,
-                                     0, 1, 1, 1, 1, -1, -1, -1, -1};
-
-  const double weights[] = {
-      0.333333333333333,  0.0555555555555556, 0.0555555555555556,
-      0.0555555555555556, 0.0555555555555556, 0.0555555555555556,
-      0.0555555555555556, 0.0277777777777778, 0.0277777777777778,
-      0.0277777777777778, 0.0277777777777778, 0.0277777777777778,
-      0.0277777777777778, 0.0277777777777778, 0.0277777777777778,
-      0.0277777777777778, 0.0277777777777778, 0.0277777777777778,
-      0.0277777777777778};
-
-  const int32_t neighbour_offset_x[] = {0, 0, 0, -1, 1, 0, 0, -1, 1, -1,
-                                        1, 0, 0, -1, 1, 0, 0, -1, 1};
-  const int32_t neighbour_offset_y[] = {0,  1, -1, 0, 0, 0, 0,  1, 1, -1,
-                                        -1, 1, -1, 0, 0, 1, -1, 0, 0};
-  const int32_t neighbour_offset_z[] = {0, 0, 0, 0, 0, 1,  -1, 0,  0, 0,
-                                        0, 1, 1, 1, 1, -1, -1, -1, -1};
-
-  for (int64_t ctr_0 = 0; ctr_0 < indexVectorSize; ctr_0 += 1) {
-    const int32_t x = *((int32_t *)(&_data_indexVector[40 * ctr_0]));
-    const int32_t y = *((int32_t *)(&_data_indexVector[40 * ctr_0 + 4]));
-    const int32_t z = *((int32_t *)(&_data_indexVector[40 * ctr_0 + 8]));
-    const int32_t dir = *((int32_t *)(&_data_indexVector[40 * ctr_0 + 12]));
-    const double vel0Term =
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_1 + _stride_pdfs_2 * z + 8 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z + 4 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z + _stride_pdfs_2 + 14 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z - _stride_pdfs_2 + 18 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_0 + _stride_pdfs_1 * y -
-                   _stride_pdfs_1 + _stride_pdfs_2 * z + 10 * _stride_pdfs_3];
-    const double vel1Term =
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y + _stride_pdfs_1 +
-                   _stride_pdfs_2 * z + _stride_pdfs_2 + 11 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y + _stride_pdfs_1 +
-                   _stride_pdfs_2 * z + _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y + _stride_pdfs_1 +
-                   _stride_pdfs_2 * z - _stride_pdfs_2 + 15 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x - _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_1 + _stride_pdfs_2 * z + 7 * _stride_pdfs_3];
-    const double vel2Term =
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z + _stride_pdfs_2 + 5 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y - _stride_pdfs_1 +
-                   _stride_pdfs_2 * z + _stride_pdfs_2 + 12 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x - _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z + _stride_pdfs_2 + 13 * _stride_pdfs_3];
-    const double rho =
-        vel0Term + vel1Term + vel2Term +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z - _stride_pdfs_2 + 6 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y - _stride_pdfs_1 +
-                   _stride_pdfs_2 * z + 2 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y - _stride_pdfs_1 +
-                   _stride_pdfs_2 * z - _stride_pdfs_2 + 16 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x - _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z + 3 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x - _stride_pdfs_0 + _stride_pdfs_1 * y +
-                   _stride_pdfs_2 * z - _stride_pdfs_2 + 17 * _stride_pdfs_3] +
-        _data_pdfs[_stride_pdfs_0 * x - _stride_pdfs_0 + _stride_pdfs_1 * y -
-                   _stride_pdfs_1 + _stride_pdfs_2 * z + 9 * _stride_pdfs_3];
-    _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_1 * y + _stride_pdfs_2 * z +
-               _stride_pdfs_3 * f_in_inv_dir_idx[dir]] =
-        -rho *
-            (6.0 * *((double *)(&_data_indexVector[40 * ctr_0 + 16])) *
-                 neighbour_offset_x[dir] +
-             6.0 * *((double *)(&_data_indexVector[40 * ctr_0 + 24])) *
-                 neighbour_offset_y[dir] +
-             6.0 * *((double *)(&_data_indexVector[40 * ctr_0 + 32])) *
-                 neighbour_offset_z[dir]) *
-            weights[dir] +
-        _data_pdfs[_stride_pdfs_0 * x + _stride_pdfs_0 * f_out_offsets_x[dir] +
-                   _stride_pdfs_1 * y + _stride_pdfs_1 * f_out_offsets_y[dir] +
-                   _stride_pdfs_2 * z + _stride_pdfs_2 * f_out_offsets_z[dir] +
-                   _stride_pdfs_3 * dir];
-  }
+static FUNC_PREFIX void dynamic_ubb_boundary_Dynamic_UBB(uint8_t * RESTRICT const _data_indexVector, double * RESTRICT _data_pdfs, int64_t const _stride_pdfs_0, int64_t const _stride_pdfs_1, int64_t const _stride_pdfs_2, int64_t const _stride_pdfs_3, int64_t indexVectorSize)
+{
+
+   const int32_t f_in_inv_dir_idx [] = { 0,2,1,4,3,6,5,10,9,8,7,16,15,18,17,12,11,14,13 }; 
+   const int32_t f_out_offsets_x [] = { 0,0,0,-1,1,0,0,-1,1,-1,1,0,0,-1,1,0,0,-1,1 }; 
+   const int32_t f_out_offsets_y [] = { 0,1,-1,0,0,0,0,1,1,-1,-1,1,-1,0,0,1,-1,0,0 }; 
+   const int32_t f_out_offsets_z [] = { 0,0,0,0,0,1,-1,0,0,0,0,1,1,1,1,-1,-1,-1,-1 }; 
+
+
+   const double weights [] = { 0.333333333333333,0.0555555555555556,0.0555555555555556,0.0555555555555556,0.0555555555555556,0.0555555555555556,0.0555555555555556,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778,0.0277777777777778 };
+
+
+
+   const int32_t neighbour_offset_x [] = { 0,0,0,-1,1,0,0,-1,1,-1,1,0,0,-1,1,0,0,-1,1 }; 
+   const int32_t neighbour_offset_y [] = { 0,1,-1,0,0,0,0,1,1,-1,-1,1,-1,0,0,1,-1,0,0 }; 
+   const int32_t neighbour_offset_z [] = { 0,0,0,0,0,1,-1,0,0,0,0,1,1,1,1,-1,-1,-1,-1 }; 
+
+   for (int64_t ctr_0 = 0; ctr_0 < indexVectorSize; ctr_0 += 1)
+   {
+      const int32_t x = *((int32_t *)(& _data_indexVector[40*ctr_0]));
+      const int32_t y = *((int32_t *)(& _data_indexVector[40*ctr_0 + 4]));
+      const int32_t z = *((int32_t *)(& _data_indexVector[40*ctr_0 + 8]));
+      const int32_t dir = *((int32_t *)(& _data_indexVector[40*ctr_0 + 12]));
+      const double vel0Term = _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_1 + _stride_pdfs_2*z + 8*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_2*z + 4*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_2*z + _stride_pdfs_2 + 14*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_2*z - _stride_pdfs_2 + 18*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_0 + _stride_pdfs_1*y - _stride_pdfs_1 + _stride_pdfs_2*z + 10*_stride_pdfs_3];
+      const double vel1Term = _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y + _stride_pdfs_1 + _stride_pdfs_2*z + _stride_pdfs_2 + 11*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y + _stride_pdfs_1 + _stride_pdfs_2*z + _stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y + _stride_pdfs_1 + _stride_pdfs_2*z - _stride_pdfs_2 + 15*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x - _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_1 + _stride_pdfs_2*z + 7*_stride_pdfs_3];
+      const double vel2Term = _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y + _stride_pdfs_2*z + _stride_pdfs_2 + 5*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y - _stride_pdfs_1 + _stride_pdfs_2*z + _stride_pdfs_2 + 12*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x - _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_2*z + _stride_pdfs_2 + 13*_stride_pdfs_3];
+      const double rho = vel0Term + vel1Term + vel2Term + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y + _stride_pdfs_2*z - _stride_pdfs_2 + 6*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y + _stride_pdfs_2*z] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y - _stride_pdfs_1 + _stride_pdfs_2*z + 2*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y - _stride_pdfs_1 + _stride_pdfs_2*z - _stride_pdfs_2 + 16*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x - _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_2*z + 3*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x - _stride_pdfs_0 + _stride_pdfs_1*y + _stride_pdfs_2*z - _stride_pdfs_2 + 17*_stride_pdfs_3] + _data_pdfs[_stride_pdfs_0*x - _stride_pdfs_0 + _stride_pdfs_1*y - _stride_pdfs_1 + _stride_pdfs_2*z + 9*_stride_pdfs_3];
+      _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_1*y + _stride_pdfs_2*z + _stride_pdfs_3*f_in_inv_dir_idx[dir]] = -rho*(6.0**((double *)(& _data_indexVector[40*ctr_0 + 16]))*neighbour_offset_x[dir] + 6.0**((double *)(& _data_indexVector[40*ctr_0 + 24]))*neighbour_offset_y[dir] + 6.0**((double *)(& _data_indexVector[40*ctr_0 + 32]))*neighbour_offset_z[dir])*weights[dir] + _data_pdfs[_stride_pdfs_0*x + _stride_pdfs_0*f_out_offsets_x[dir] + _stride_pdfs_1*y + _stride_pdfs_1*f_out_offsets_y[dir] + _stride_pdfs_2*z + _stride_pdfs_2*f_out_offsets_z[dir] + _stride_pdfs_3*dir];
+   }
+}
 }
-} // namespace internal_dynamic_ubb_boundary_Dynamic_UBB
+
 
 #ifdef __GNUC__
 #pragma GCC diagnostic pop
@@ -148,34 +86,47 @@ static FUNC_PREFIX void dynamic_ubb_boundary_Dynamic_UBB(
 #pragma pop
 #endif
 
-void Dynamic_UBB::run_impl(IBlock *block, IndexVectors::Type type) {
-  auto *indexVectors = block->getData<IndexVectors>(indexVectorID);
-  int64_t indexVectorSize = int64_c(indexVectors->indexVector(type).size());
-  if (indexVectorSize == 0)
-    return;
 
-  auto pointer = indexVectors->pointerCpu(type);
+void Dynamic_UBB::run_impl(IBlock * block, IndexVectors::Type type)
+{
+   auto * indexVectors = block->getData<IndexVectors>(indexVectorID);
+   int64_t indexVectorSize = int64_c( indexVectors->indexVector(type).size() );
+   if( indexVectorSize == 0)
+      return;
 
-  uint8_t *_data_indexVector = reinterpret_cast<uint8_t *>(pointer);
+
+   auto pointer = indexVectors->pointerCpu(type);
+
 
-  auto pdfs = block->getData<field::GhostLayerField<double, 19>>(pdfsID);
+   uint8_t * _data_indexVector = reinterpret_cast<uint8_t*>(pointer);
 
-  WALBERLA_ASSERT_GREATER_EQUAL(0, -int_c(pdfs->nrOfGhostLayers()));
-  double *RESTRICT _data_pdfs = pdfs->dataAt(0, 0, 0, 0);
-  const int64_t _stride_pdfs_0 = int64_t(pdfs->xStride());
-  const int64_t _stride_pdfs_1 = int64_t(pdfs->yStride());
-  const int64_t _stride_pdfs_2 = int64_t(pdfs->zStride());
-  const int64_t _stride_pdfs_3 = int64_t(1 * int64_t(pdfs->fStride()));
-  internal_dynamic_ubb_boundary_Dynamic_UBB::dynamic_ubb_boundary_Dynamic_UBB(
-      _data_indexVector, _data_pdfs, _stride_pdfs_0, _stride_pdfs_1,
-      _stride_pdfs_2, _stride_pdfs_3, indexVectorSize);
+   auto pdfs = block->getData< field::GhostLayerField<double, 19> >(pdfsID);
+
+
+   WALBERLA_ASSERT_GREATER_EQUAL(0, -int_c(pdfs->nrOfGhostLayers()));
+    double * RESTRICT _data_pdfs = pdfs->dataAt(0, 0, 0, 0);
+    const int64_t _stride_pdfs_0 = int64_t(pdfs->xStride());
+    const int64_t _stride_pdfs_1 = int64_t(pdfs->yStride());
+    const int64_t _stride_pdfs_2 = int64_t(pdfs->zStride());
+    const int64_t _stride_pdfs_3 = int64_t(1 * int64_t(pdfs->fStride()));
+    internal_dynamic_ubb_boundary_Dynamic_UBB::dynamic_ubb_boundary_Dynamic_UBB(_data_indexVector, _data_pdfs, _stride_pdfs_0, _stride_pdfs_1, _stride_pdfs_2, _stride_pdfs_3, indexVectorSize);
 }
 
-void Dynamic_UBB::run(IBlock *block) { run_impl(block, IndexVectors::ALL); }
+void Dynamic_UBB::run(IBlock * block)
+{
+   run_impl(block, IndexVectors::ALL);
+}
 
-void Dynamic_UBB::inner(IBlock *block) { run_impl(block, IndexVectors::INNER); }
+void Dynamic_UBB::inner(IBlock * block)
+{
+   run_impl(block, IndexVectors::INNER);
+}
 
-void Dynamic_UBB::outer(IBlock *block) { run_impl(block, IndexVectors::OUTER); }
+void Dynamic_UBB::outer(IBlock * block)
+{
+   run_impl(block, IndexVectors::OUTER);
+}
 
 } // namespace lbm
 } // namespace walberla
+