external force density in first integration step

samples/billiard.py

@@ -186,7 +186,7 @@ def main():
     mass = np.array([0.17])
 
     num_types = len(wca_sig)
-    
+
     def mix_eps(eps1, eps2, rule='LB'):
         return math.sqrt(eps1 * eps2)
 

src/core/grid_based_algorithms/lb.cpp

@@ -172,6 +172,51 @@ static double fluidstep = 0.0;
 #include "grid.hpp"
 
 /********************** The Main LB Part *************************************/
+
+/**
+ * @brief Initialize fluid nodes.
+ * @param[out] fields         Vector containing the fluid nodes
+ * @param[in]  lb_parameters  Parameters for the LB
+ * @param[in]  lb_lattice     Lattice instance
+ */
+void lb_initialize_fields(std::vector<LB_FluidNode> &fields,
+                          LB_Parameters const &lb_parameters,
+                          Lattice const &lb_lattice) {
+  fields.resize(lb_lattice.halo_grid_volume);
+  for (auto &field : fields) {
+    field.force_density = lb_parameters.ext_force_density;
+#ifdef LB_BOUNDARIES
+    field.boundary = false;
+#endif // LB_BOUNDARIES
+  }
+}
+
+/** (Re-)allocate memory for the fluid and initialize pointers. */
+void lb_realloc_fluid(LB_FluidData &lb_fluid_a, LB_FluidData &lb_fluid_b,
+                      const Lattice::index_t halo_grid_volume,
+                      LB_Fluid &lb_fluid, LB_Fluid &lb_fluid_post) {
+  const std::array<int, 2> size = {{D3Q19::n_vel, halo_grid_volume}};
+
+  lb_fluid_a.resize(size);
+  lb_fluid_b.resize(size);
+
+  using Utils::Span;
+  for (int i = 0; i < size[0]; i++) {
+    lb_fluid[i] = Span<double>(lb_fluid_a[i].origin(), size[1]);
+    lb_fluid_post[i] = Span<double>(lb_fluid_b[i].origin(), size[1]);
+  }
+}
+
+void lb_set_equilibrium_populations(const Lattice &lb_lattice,
+                                    const LB_Parameters &lb_parameters) {
+  for (Lattice::index_t index = 0; index < lb_lattice.halo_grid_volume;
+       ++index) {
+    lb_set_population_from_density_momentum_density_stress(
+        index, lb_parameters.density, Utils::Vector3d{} /*momentum density*/,
+        Utils::Vector6d{} /*stress*/);
+  }
+}
+
 void lb_init(const LB_Parameters &lb_parameters) {
   if (lb_parameters.agrid <= 0.0) {
     runtimeErrorMsg()
@@ -181,7 +226,6 @@ void lb_init(const LB_Parameters &lb_parameters) {
     return;
 
   /* initialize the local lattice domain */
-
   try {
     lblattice = Lattice(lb_parameters.agrid, 0.5 /*offset*/, 1 /*halo size*/,
                         local_geo.length(), local_geo.my_right(),
@@ -193,40 +237,28 @@ void lb_init(const LB_Parameters &lb_parameters) {
 
   /* allocate memory for data structures */
   lb_realloc_fluid(lbfluid_a, lbfluid_b, lblattice.halo_grid_volume, lbfluid,
-                   lbfluid_post, lbfields);
+                   lbfluid_post);
+
+  lb_initialize_fields(lbfields, lbpar, lblattice);
 
   /* prepare the halo communication */
   lb_prepare_communication(update_halo_comm, lblattice);
 
   /* initialize derived parameters */
   lb_reinit_parameters(lbpar);
 
-  /* setup the initial populations */
-  lb_reinit_fluid(lbfields, lblattice, lbpar);
-}
-
-void lb_reinit_fluid(std::vector<LB_FluidNode> &lb_fields,
-                     const Lattice &lb_lattice,
-                     const LB_Parameters &lb_parameters) {
-  std::fill(lb_fields.begin(), lb_fields.end(), LB_FluidNode());
-  /* default values for fields in lattice units */
-  Utils::Vector3d momentum_density{};
-  Utils::Vector6d stress{};
-
-  for (Lattice::index_t index = 0; index < lb_lattice.halo_grid_volume;
-       ++index) {
-    // sets equilibrium distribution
-    lb_set_population_from_density_momentum_density_stress(
-        index, lb_parameters.density, momentum_density, stress);
-
-#ifdef LB_BOUNDARIES
-    lb_fields[index].boundary = 0;
-#endif // LB_BOUNDARIES
-  }
+  lb_set_equilibrium_populations(lblattice, lbpar);
 
 #ifdef LB_BOUNDARIES
   LBBoundaries::lb_init_boundaries();
-#endif // LB_BOUNDARIES
+#endif
+}
+
+void lb_reinit_fluid(std::vector<LB_FluidNode> &lb_fields,
+                     Lattice const &lb_lattice,
+                     LB_Parameters const &lb_parameters) {
+  lb_set_equilibrium_populations(lb_lattice, lb_parameters);
+  lb_initialize_fields(lbfields, lb_parameters, lb_lattice);
 }
 
 void lb_reinit_parameters(LB_Parameters &lb_parameters) {
@@ -587,25 +619,6 @@ void lb_fluid_set_rng_state(uint64_t counter) {
 
 /***********************************************************************/
 
-/** (Re-)allocate memory for the fluid and initialize pointers. */
-void lb_realloc_fluid(LB_FluidData &lb_fluid_a, LB_FluidData &lb_fluid_b,
-                      const Lattice::index_t halo_grid_volume,
-                      LB_Fluid &lb_fluid, LB_Fluid &lb_fluid_post,
-                      std::vector<LB_FluidNode> &lb_fields) {
-  const std::array<int, 2> size = {{D3Q19::n_vel, halo_grid_volume}};
-
-  lb_fluid_a.resize(size);
-  lb_fluid_b.resize(size);
-
-  using Utils::Span;
-  for (int i = 0; i < size[0]; i++) {
-    lb_fluid[i] = Span<double>(lb_fluid_a[i].origin(), size[1]);
-    lb_fluid_post[i] = Span<double>(lb_fluid_b[i].origin(), size[1]);
-  }
-
-  lb_fields.resize(halo_grid_volume);
-}
-
 /** Set up the structures for exchange of the halo regions.
  *  See also \ref halo.cpp
  */
@@ -938,18 +951,6 @@ inline void lb_collide_stream() {
   }
 #endif // LB_BOUNDARIES
 
-#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
-  // Safeguard the node forces so that we can later use them for the IBM
-  // particle update
-  // In the following loop the lbfields[XX].force are reset to zero
-  // Safeguard the node forces so that we can later use them for the IBM
-  // particle update In the following loop the lbfields[XX].force are reset to
-  // zero
-  for (int i = 0; i < lblattice.halo_grid_volume; ++i) {
-    lbfields[i].force_density_buf = lbfields[i].force_density;
-  }
-#endif
-
   Lattice::index_t index = lblattice.halo_offset;
   for (int z = 1; z <= lblattice.grid[2]; z++) {
     for (int y = 1; y <= lblattice.grid[1]; y++) {
@@ -974,6 +975,12 @@ inline void lb_collide_stream() {
           auto const modes_with_forces =
               lb_apply_forces(index, thermalized_modes, lbpar, lbfields);
 
+#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
+          // Safeguard the node forces so that we can later use them for the IBM
+          // particle update
+          lbfields[index].force_density_buf = lbfields[index].force_density;
+#endif
+
           /* reset the force density */
           lbfields[index].force_density = lbpar.ext_force_density;
 

src/core/grid_based_algorithms/lb.hpp

@@ -143,13 +143,6 @@ extern Lattice lblattice;
 
 extern HaloCommunicator update_halo_comm;
 
-void lb_realloc_fluid(boost::multi_array<double, 2> &lb_fluid_a,
-                      boost::multi_array<double, 2> &lb_fluid_b,
-                      Lattice::index_t halo_grid_volume,
-                      std::array<Utils::Span<double>, 19> &lb_fluid,
-                      std::array<Utils::Span<double>, 19> &lb_fluid_post,
-                      std::vector<LB_FluidNode> &lb_fields);
-
 void lb_init(const LB_Parameters &lb_parameters);
 
 void lb_reinit_fluid(std::vector<LB_FluidNode> &lb_fields,
@@ -287,6 +280,9 @@ void lb_calc_fluid_momentum(double *result, const LB_Parameters &lb_parameters,
                             const std::vector<LB_FluidNode> &lb_fields,
                             const Lattice &lb_lattice);
 void lb_collect_boundary_forces(double *result);
+void lb_initialize_fields(std::vector<LB_FluidNode> &fields,
+                          LB_Parameters const &lb_parameters,
+                          Lattice const &lb_lattice);
 
 /*@}*/
 

src/core/grid_based_algorithms/lb_interface.cpp

@@ -89,6 +89,17 @@ void mpi_lb_set_population(Utils::Vector3i const &index,
 
 REGISTER_CALLBACK(mpi_lb_set_population)
 
+void mpi_lb_set_force_density(Utils::Vector3i const &index,
+                              Utils::Vector3d const &force_density) {
+  lb_set(index, [&](auto index) {
+    auto const linear_index =
+        get_linear_index(lblattice.local_index(index), lblattice.halo_grid);
+    lbfields[linear_index].force_density = force_density;
+  });
+}
+
+REGISTER_CALLBACK(mpi_lb_set_force_density)
+
 auto mpi_lb_get_momentum_density(Utils::Vector3i const &index) {
   return lb_calc_fluid_kernel(index, [&](auto modes, auto force_density) {
     return lb_calc_momentum_density(modes, force_density);
@@ -1265,6 +1276,7 @@ void lb_lbnode_set_velocity(const Utils::Vector3i &ind,
         lb_get_population_from_density_momentum_density_stress(
             density, momentum_density, stress);
     mpi_call_all(mpi_lb_set_population, ind, population);
+    mpi_call_all(mpi_lb_set_force_density, ind, Utils::Vector3d{});
   }
 }
 
@@ -1308,6 +1320,7 @@ void lb_lbfluid_on_lb_params_change(LBParam field) {
     break;
   case LBParam::VISCOSITY:
   case LBParam::EXT_FORCE_DENSITY:
+    lb_initialize_fields(lbfields, lbpar, lblattice);
   case LBParam::BULKVISC:
   case LBParam::KT:
   case LBParam::GAMMA_ODD:

src/core/grid_based_algorithms/lbgpu.hpp

@@ -213,10 +213,6 @@ void lb_init_boundaries_GPU(int n_lb_boundaries, int number_of_boundnodes,
                             int *host_boundary_index_list,
                             float *lb_bounday_velocity);
 #endif
-void lb_init_extern_nodeforcedensities_GPU(
-    int n_extern_node_force_densities,
-    LB_extern_nodeforcedensity_gpu *host_extern_node_force_densities,
-    LB_parameters_gpu *lbpar_gpu);
 
 void lb_set_agrid_gpu(double agrid);