Es walberla

src/core/grid_based_algorithms/lb_interface.cpp

@@ -158,8 +158,38 @@ double lb_lbfluid_get_viscosity() {
 }
 
 double lb_lbfluid_get_bulk_viscosity() {
-  if (lattice_switch == ActiveLB::WALBERLA)
-    throw std::runtime_error("Getting bulk viscosity not implemented.");
+#ifdef LB_WALBERLA
+  if (lattice_switch == ActiveLB::WALBERLA) {
+    return lb_walberla()->get_bulk_viscosity();
+  }
+#endif
+  throw NoLBActive();
+}
+
+double lb_lbfluid_get_shear_viscosity() {
+#ifdef LB_WALBERLA
+  if (lattice_switch == ActiveLB::WALBERLA) {
+    return lb_walberla()->get_shear_viscosity();
+  }
+#endif
+  throw NoLBActive();
+}
+
+double lb_lbfluid_get_gamma_odd() {
+#ifdef LB_WALBERLA
+  if (lattice_switch == ActiveLB::WALBERLA) {
+    return lb_walberla()->get_gamma_odd();
+  }
+#endif
+  throw NoLBActive();
+}
+
+double lb_lbfluid_get_gamma_even() {
+#ifdef LB_WALBERLA
+  if (lattice_switch == ActiveLB::WALBERLA) {
+    return lb_walberla()->get_gamma_even();
+  }
+#endif
   throw NoLBActive();
 }
 

src/core/grid_based_algorithms/lb_interface.hpp

@@ -125,11 +125,21 @@ double lb_lbfluid_get_agrid();
  */
 double lb_lbfluid_get_gamma_even();
 
+/**
+ * @brief Get the global LB relaxation parameter for odd modes.
+ */
+double lb_lbfluid_get_gamma_odd();
+
 /**
  * @brief Get the global LB bulk viscosity.
  */
 double lb_lbfluid_get_bulk_viscosity();
 
+/**
+ * @brief Get the global LB shear viscosity.
+ */
+double lb_lbfluid_get_shear_viscosity();
+
 /**
  * @brief Get the global LB viscosity.
  */

src/core/grid_based_algorithms/lb_walberla_instance.cpp

@@ -55,23 +55,45 @@ LBWalberlaParams *lb_walberla_params() {
   return lb_walberla_params_instance;
 }
 
-void init_lb_walberla(double viscosity, double density, double agrid,
-                      double tau, const Utils::Vector3d &box_dimensions,
-                      const Utils::Vector3i &node_grid, double kT,
-                      unsigned int seed) {
+void init_lb_walberla_magic_number(double viscosity, double magic_number,
+                                   double density, double agrid, double tau,
+                                   const Utils::Vector3d &box_dimensions,
+                                   const Utils::Vector3i &node_grid, double kT,
+                                   unsigned int seed) {
   // Exceptions need to be converted to runtime errors so they can be
   // handled from Python in a parallel simulation
   try {
-    lb_walberla_instance = new_lb_walberla(viscosity, density, agrid, tau,
-                                           box_dimensions, node_grid, kT, seed);
+    lb_walberla_instance =
+        new_lb_walberla(viscosity, magic_number, density, agrid, tau,
+                        box_dimensions, node_grid, kT, seed);
     lb_walberla_params_instance = new LBWalberlaParams{agrid, tau};
   } catch (const std::exception &e) {
     runtimeErrorMsg() << "Error during Walberla initialization: " << e.what();
     lb_walberla_instance = nullptr;
     lb_walberla_params_instance = nullptr;
   }
 }
-REGISTER_CALLBACK(init_lb_walberla)
+REGISTER_CALLBACK(init_lb_walberla_magic_number)
+
+void init_lb_walberla_relaxation_rates(LBRelaxationRates relaxation_rates,
+                                       double density, double agrid, double tau,
+                                       const Utils::Vector3d &box_dimensions,
+                                       const Utils::Vector3i &node_grid,
+                                       double kT, unsigned int seed) {
+  // Exceptions need to be converted to runtime errors so they can be
+  // handled from Python in a parallel simulation
+  try {
+    lb_walberla_instance =
+        new_lb_walberla(relaxation_rates, density, agrid, tau, box_dimensions,
+                        node_grid, kT, seed);
+    lb_walberla_params_instance = new LBWalberlaParams{agrid, tau};
+  } catch (const std::exception &e) {
+    runtimeErrorMsg() << "Error during Walberla initialization: " << e.what();
+    lb_walberla_instance = nullptr;
+    lb_walberla_params_instance = nullptr;
+  }
+}
+REGISTER_CALLBACK(init_lb_walberla_relaxation_rates)
 
 void destruct_lb_walberla() {
   delete lb_walberla_instance;
@@ -81,11 +103,24 @@ void destruct_lb_walberla() {
 }
 REGISTER_CALLBACK(destruct_lb_walberla)
 
-void mpi_init_lb_walberla(double viscosity, double density, double agrid,
-                          double tau, double kT, unsigned int seed) {
-  Communication::mpiCallbacks().call_all(init_lb_walberla, viscosity, density,
-                                         agrid, tau, box_geo.length(),
-                                         node_grid, kT, seed);
+void mpi_init_lb_walberla(double viscosity, double magic_number, double density,
+                          double agrid, double tau, double kT,
+                          unsigned int seed) {
+  Communication::mpiCallbacks().call_all(
+      init_lb_walberla_magic_number, viscosity, magic_number, density, agrid,
+      tau, box_geo.length(), node_grid, kT, seed);
+  if (lb_walberla_instance) {
+    lb_lbfluid_set_lattice_switch(ActiveLB::WALBERLA);
+    lb_lbfluid_sanity_checks();
+  }
+}
+
+void mpi_init_lb_walberla(LBRelaxationRates relaxation_rates, double density,
+                          double agrid, double tau, double kT,
+                          unsigned int seed) {
+  Communication::mpiCallbacks().call_all(init_lb_walberla_relaxation_rates,
+                                         relaxation_rates, density, agrid, tau,
+                                         box_geo.length(), node_grid, kT, seed);
   if (lb_walberla_instance) {
     lb_lbfluid_set_lattice_switch(ActiveLB::WALBERLA);
     lb_lbfluid_sanity_checks();

src/core/grid_based_algorithms/lb_walberla_instance.hpp

@@ -47,12 +47,26 @@ LBWalberlaParams *lb_walberla_params();
  *  WALBERLA
  *
  *  @param viscosity Fluid viscosity
+ *  @param magic_number Fluid magic number
  *  @param density Fluid density
  *  @param agrid  Size of one lb cell
  *  @param tau    LB time step
  */
-void mpi_init_lb_walberla(double viscosity, double density, double agrid,
-                          double tau, double kT, unsigned int seed);
+void mpi_init_lb_walberla(double viscosity, double magic_number, double density,
+                          double agrid, double tau, double kT,
+                          unsigned int seed);
+
+/** @brief Create the LBWalberla instance and sets the lattice switch to
+ *  WALBERLA
+ *
+ *  @param relaxation_rates Fluid relaxation rates
+ *  @param density Fluid density
+ *  @param agrid  Size of one lb cell
+ *  @param tau    LB time step
+ */
+void mpi_init_lb_walberla(LBRelaxationRates relaxation_rates, double density,
+                          double agrid, double tau, double kT,
+                          unsigned int seed);
 
 /** @brief Destruct the LBWalberla instance and set lattice switch to NONE */
 void mpi_destruct_lb_walberla();

src/python/espressomd/lb.pxd

@@ -38,6 +38,14 @@ cdef class HydrodynamicInteraction(Actor):
     #
 
 IF LB_WALBERLA:
+    cdef extern from "LBWalberlaBase.hpp":
+        cdef cppclass LBRelaxationRates:
+            double omega_bulk;
+            double omega_shear;
+            double omega_odd;
+            double omega_even;
+
+
     cdef extern from "grid_based_algorithms/lb_interface.hpp":
 
         cdef enum OutputVTK:
@@ -64,6 +72,9 @@ cdef extern from "grid_based_algorithms/lb_interface.hpp":
     void lb_lbfluid_set_ext_force_density(const Vector3d forcedensity) except +
     const Vector3d lb_lbfluid_get_ext_force_density() except +
     double lb_lbfluid_get_bulk_viscosity() except +
+    double lb_lbfluid_get_shear_viscosity() except +
+    double lb_lbfluid_get_gamma_odd() except +
+    double lb_lbfluid_get_gamma_even() except +
     void lb_lbfluid_sanity_checks() except +
     void lb_lbfluid_save_checkpoint(string filename, bool binary) except +
     void lb_lbfluid_load_checkpoint(string filename, bool binary) except +
@@ -111,7 +122,8 @@ cdef extern from "grid_based_algorithms/lb_interpolation.hpp" namespace "Interpo
 
 IF LB_WALBERLA:
     cdef extern from "grid_based_algorithms/lb_walberla_instance.hpp":
-        void mpi_init_lb_walberla(double viscosity, double density, double agrid, double tau, double kT, unsigned int seed) except +
+        void mpi_init_lb_walberla(double viscosity, double magic_number, double density, double agrid, double tau, double kT, unsigned int seed) except +
+        void mpi_init_lb_walberla(LBRelaxationRates relaxation_rates, double density, double agrid, double tau, double kT, unsigned int seed) except +
         void mpi_destruct_lb_walberla() except +
 
 
@@ -142,6 +154,15 @@ cdef inline double python_lbfluid_get_viscosity(p_agrid, p_tau) except +:
 cdef inline double python_lbfluid_get_bulk_viscosity(p_agrid, p_tau) except +:
     return lb_lbfluid_get_bulk_viscosity() / p_tau * (p_agrid * p_agrid)
 
+cdef inline double python_lbfluid_get_shear_viscosity(p_agrid, p_tau) except +:
+    return lb_lbfluid_get_shear_viscosity() / p_tau * (p_agrid * p_agrid)
+
+cdef inline double python_lbfluid_get_gamma_odd(p_agrid, p_tau) except +:
+    return lb_lbfluid_get_gamma_odd() / p_tau * (p_agrid * p_agrid)
+
+cdef inline double python_lbfluid_get_gamma_even(p_agrid, p_tau) except +:
+    return lb_lbfluid_get_gamma_even() / p_tau * (p_agrid * p_agrid)
+
 cdef inline double python_lbfluid_get_gamma() except +:
     return lb_lbcoupling_get_gamma()
 

src/python/espressomd/lb.pyx

@@ -183,8 +183,12 @@ cdef class HydrodynamicInteraction(Actor):
         Fluid density.
     visc : :obj:`float`
         Fluid kinematic viscosity.
+    magic_number : :obj:`float`
+        Magic number to define the relaxation time from viscosity.
     bulk_visc : :obj:`float`, optional
         Fluid bulk viscosity.
+    shear_visc : :obj:`float`, optional
+        Fluid shear viscosity.
     gamma_odd : :obj:`int`, optional
         Relaxation parameter :math:`\\gamma_{\\textrm{odd}}` for kinetic modes.
     gamma_even : :obj:`int`, optional
@@ -231,10 +235,10 @@ cdef class HydrodynamicInteraction(Actor):
             raise ValueError("tau has to be a positive double")
 
     def valid_keys(self):
-        return "agrid", "dens", "ext_force_density", "visc", "tau", "bulk_visc", "gamma_odd", "gamma_even", "kT", "seed"
+        return "agrid", "dens", "ext_force_density", "visc", "magic_number", "tau", "bulk_visc", "shear_visc", "gamma_odd", "gamma_even", "kT", "seed"
 
     def required_keys(self):
-        return ["dens", "agrid", "visc", "tau"]
+        return ["dens", "agrid", "tau"]
 
     def default_params(self):
         return {"agrid": -1.0,
@@ -365,6 +369,22 @@ cdef class HydrodynamicInteraction(Actor):
         def __get__(self):
             return python_lbfluid_get_viscosity(self.agrid, self.tau)
 
+    property bulk_viscosity:
+        def __get__(self):
+            return python_lbfluid_get_bulk_viscosity(self.agrid, self.tau)
+
+    property shear_viscosity:
+        def __get__(self):
+            return python_lbfluid_get_shear_viscosity(self.agrid, self.tau)
+
+    property gamma_odd:
+        def __get__(self):
+            return python_lbfluid_get_gamma_odd(self.agrid, self.tau)
+
+    property gamma_even:
+        def __get__(self):
+            return python_lbfluid_get_gamma_even(self.agrid, self.tau)
+
     property tau:
         def __get__(self):
             return lb_lbfluid_get_tau()
@@ -399,6 +419,10 @@ IF LB_WALBERLA:
                     "ext_force_density": [0.0, 0.0, 0.0],
                     "visc": -1.0,
                     "bulk_visc": -1.0,
+                    "shear_visc": -1.0,
+                    "gamma_odd": -1.0,
+                    "gamma_even": -1.0,
+                    "magic_number": 3.0/16.0,
                     "tau": -1.0,
                     "kT": 0.0,
                     "seed": 0}
@@ -407,17 +431,47 @@ IF LB_WALBERLA:
             raise Exception("This may not be called")
 
         def _activate_method(self):
+            cdef LBRelaxationRates relaxation_rates
+
             self.validate_params()
 
             # unit conversion
             lb_visc = self._params["visc"] * \
                 self._params['tau'] / self._params['agrid']**2
+            lb_magic_number = self._params["magic_number"]
+
+            lb_bulk_visc = self._params["bulk_visc"] * \
+                self._params['tau'] / self._params['agrid']**2
+            lb_shear_visc = self._params["shear_visc"] * \
+                self._params['tau'] / self._params['agrid']**2
+            lb_gamma_odd = self._params["gamma_odd"] * \
+                self._params['tau'] / self._params['agrid']**2
+            lb_gamma_even = self._params["gamma_even"] * \
+                self._params['tau'] / self._params['agrid']**2
+
             lb_dens = self._params['dens'] * self._params['agrid']**3
             lb_kT = self._params['kT'] * \
-                self._params['tau']**2 / self._params['agrid']**2 
-            mpi_init_lb_walberla(
-                lb_visc, lb_dens, self._params["agrid"], self._params["tau"],
-                lb_kT, self._params['seed'])
+                self._params['tau']**2 / self._params['agrid']**2
+
+            default_params = self.default_params()
+
+            if ((not lb_bulk_visc == default_params["bulk_visc"]) and
+               (not lb_shear_visc == default_params["shear_visc"]) and
+               (not lb_gamma_odd == default_params["gamma_odd"]) and
+               (not lb_gamma_even == default_params["gamma_even"])):
+
+                relaxation_rates.omega_bulk = lb_bulk_visc
+                relaxation_rates.omega_shear = lb_shear_visc
+                relaxation_rates.omega_odd = lb_gamma_odd
+                relaxation_rates.omega_even = lb_gamma_even
+                mpi_init_lb_walberla(
+                    relaxation_rates, lb_dens, self._params["agrid"],
+                    self._params["tau"], lb_kT, self._params['seed'])
+            else:
+                mpi_init_lb_walberla(
+                    lb_visc, lb_magic_number, lb_dens,
+                    self._params["agrid"], self._params["tau"],
+                    lb_kT, self._params['seed'])
             utils.handle_errors("LB fluid activation")
             self.ext_force_density = self._params["ext_force_density"]
 

src/walberla_bridge/LBWalberlaBase.hpp

@@ -105,8 +105,12 @@ class LBWalberlaBase {
   virtual Utils::Vector3d get_external_force() const = 0;
 
   // Global parameters
-  virtual void set_viscosity(double viscosity) = 0;
+  virtual void set_viscosity(double viscosity, double magic_number) = 0;
   virtual double get_viscosity() const = 0;
+  virtual double get_bulk_viscosity() const = 0;
+  virtual double get_shear_viscosity() const = 0;
+  virtual double get_gamma_odd() const = 0;
+  virtual double get_gamma_even() const = 0;
   virtual double get_kT() const = 0;
 
   // Grid, domain, halo
@@ -159,4 +163,18 @@ enum class OutputVTK : unsigned {
   pressure_tensor = 1u << 2u,
 };
 
+struct LBRelaxationRates {
+  double omega_bulk;
+  double omega_shear;
+  double omega_odd;
+  double omega_even;
+  template <class Archive>
+  void serialize(Archive &ar, const unsigned int version) {
+    ar &omega_bulk;
+    ar &omega_shear;
+    ar &omega_odd;
+    ar &omega_even;
+  }
+};
+
 #endif // LB_WALBERLA_H