diff --git a/src/core/event.cpp b/src/core/event.cpp index e423fa5e3fe..b809cffa65e 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -110,12 +110,14 @@ void on_integration_start(double time_step) { integrator_npt_sanity_checks(); #endif interactions_sanity_checks(); - lb_lbfluid_on_integration_start(); + lb_lbfluid_sanity_checks(time_step); /********************************************/ /* end sanity checks */ /********************************************/ + lb_lbfluid_on_integration_start(); + #ifdef CUDA MPI_Bcast(gpu_get_global_particle_vars_pointer_host(), sizeof(CUDA_global_part_vars), MPI_BYTE, 0, comm_cart); diff --git a/src/core/grid_based_algorithms/lb.cpp b/src/core/grid_based_algorithms/lb.cpp index 465c7524dcd..4b2a5d9c559 100644 --- a/src/core/grid_based_algorithms/lb.cpp +++ b/src/core/grid_based_algorithms/lb.cpp @@ -35,7 +35,6 @@ #include "grid.hpp" #include "grid_based_algorithms/lb_boundaries.hpp" #include "halo.hpp" -#include "integrate.hpp" #include "lb-d3q19.hpp" #include "random.hpp" @@ -983,13 +982,10 @@ void lb_collide_stream() { * for the lattice dynamics can be coarser than the MD time step, we * monitor the time since the last lattice update. */ -void lattice_boltzmann_update(double time_step) { - auto const factor = static_cast(round(lbpar.tau / time_step)); - +void lattice_boltzmann_update(int lb_steps_per_md_step) { fluidstep += 1; - if (fluidstep >= factor) { + if (fluidstep >= lb_steps_per_md_step) { fluidstep = 0; - lb_collide_stream(); } } diff --git a/src/core/grid_based_algorithms/lb.hpp b/src/core/grid_based_algorithms/lb.hpp index 0169498ea74..37de29d5c30 100644 --- a/src/core/grid_based_algorithms/lb.hpp +++ b/src/core/grid_based_algorithms/lb.hpp @@ -183,7 +183,7 @@ extern std::vector lbfields; * collisions. If boundaries are present, it also applies the boundary * conditions. */ -void lattice_boltzmann_update(double time_step); +void lattice_boltzmann_update(int lb_steps_per_md_step); void lb_sanity_checks(const LB_Parameters &lb_parameters); diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp index 4cc7211e543..b785ecd1f7e 100644 --- a/src/core/grid_based_algorithms/lb_interface.cpp +++ b/src/core/grid_based_algorithms/lb_interface.cpp @@ -52,9 +52,9 @@ struct NoLBActive : public std::exception { const char *what() const noexcept override { return "LB not activated"; } }; -void lb_lbfluid_update(double time_step) { +void lb_lbfluid_update(int lb_steps_per_md_step) { if (lattice_switch == ActiveLB::CPU) { - lattice_boltzmann_update(time_step); + lattice_boltzmann_update(lb_steps_per_md_step); } else if (lattice_switch == ActiveLB::GPU and this_node == 0) { #ifdef CUDA #ifdef ELECTROKINETICS @@ -62,7 +62,7 @@ void lb_lbfluid_update(double time_step) { ek_integrate(); } else { #endif - lattice_boltzmann_update_gpu(time_step); + lattice_boltzmann_update_gpu(lb_steps_per_md_step); #ifdef ELECTROKINETICS } #endif @@ -70,9 +70,9 @@ void lb_lbfluid_update(double time_step) { } } -void lb_lbfluid_propagate(double time_step) { +void lb_lbfluid_propagate(int lb_steps_per_md_step) { if (lattice_switch != ActiveLB::NONE) { - lb_lbfluid_update(time_step); + lb_lbfluid_update(lb_steps_per_md_step); if (lb_lbfluid_get_kT() > 0.0) { if (lattice_switch == ActiveLB::GPU) { #ifdef CUDA @@ -104,8 +104,7 @@ void lb_boundary_mach_check() { } } -void lb_lbfluid_sanity_checks() { - extern double time_step; +void lb_lbfluid_sanity_checks(double time_step) { if (lattice_switch == ActiveLB::GPU && this_node == 0) { #ifdef CUDA lb_GPU_sanity_checks(); @@ -123,7 +122,6 @@ void lb_lbfluid_sanity_checks() { } void lb_lbfluid_on_integration_start() { - lb_lbfluid_sanity_checks(); if (lattice_switch == ActiveLB::CPU) { halo_communication(update_halo_comm, reinterpret_cast(lbfluid[0].data())); @@ -418,18 +416,18 @@ void lb_lbfluid_set_tau(double tau) { } } -void check_tau_time_step_consistency(double tau, double time_s) { +void check_tau_time_step_consistency(double tau, double time_step) { auto const eps = std::numeric_limits::epsilon(); - if ((tau - time_s) / (tau + time_s) < -eps) + if ((tau - time_step) / (tau + time_step) < -eps) throw std::invalid_argument("LB tau (" + std::to_string(tau) + ") must be >= MD time_step (" + - std::to_string(time_s) + ")"); - auto const factor = tau / time_s; + std::to_string(time_step) + ")"); + auto const factor = tau / time_step; if (fabs(round(factor) - factor) / factor > eps) throw std::invalid_argument("LB tau (" + std::to_string(tau) + ") must be integer multiple of " "MD time_step (" + - std::to_string(time_s) + "). Factor is " + + std::to_string(time_step) + "). Factor is " + std::to_string(factor)); } diff --git a/src/core/grid_based_algorithms/lb_interface.hpp b/src/core/grid_based_algorithms/lb_interface.hpp index 009729315a5..af33d1d456e 100644 --- a/src/core/grid_based_algorithms/lb_interface.hpp +++ b/src/core/grid_based_algorithms/lb_interface.hpp @@ -37,7 +37,7 @@ extern ActiveLB lattice_switch; /** * @brief Propagate the LB fluid. */ -void lb_lbfluid_propagate(double time_step); +void lb_lbfluid_propagate(int lb_steps_per_md_step); /** * @brief Event handler for integration start. @@ -133,7 +133,7 @@ void lb_lbfluid_set_kT(double kT); /** * @brief Perform LB parameter and boundary velocity checks. */ -void lb_lbfluid_sanity_checks(); +void lb_lbfluid_sanity_checks(double time_step); /** * @brief Set the LB density for a single node. diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index 58adeb69992..80eaaec19e7 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -122,6 +122,7 @@ namespace { * @brief Add a force to the lattice force density. * @param pos Position of the force * @param force Force in MD units. + * @param time_step MD time step. */ void add_md_force(Utils::Vector3d const &pos, Utils::Vector3d const &force, double time_step) { @@ -137,9 +138,10 @@ void add_md_force(Utils::Vector3d const &pos, Utils::Vector3d const &force, * Section II.C. @cite ahlrichs99a * * @param[in] p The coupled particle. - * @param[in] f_random Additional force to be included. + * @param[in] f_random Additional force to be included. + * @param[in] time_step MD time step. * - * @return The viscous coupling force plus f_random. + * @return The viscous coupling force plus @p f_random. */ Utils::Vector3d lb_viscous_coupling(Particle const &p, Utils::Vector3d const &f_random, diff --git a/src/core/grid_based_algorithms/lbgpu.cpp b/src/core/grid_based_algorithms/lbgpu.cpp index 97addd4fdd0..e698785b412 100644 --- a/src/core/grid_based_algorithms/lbgpu.cpp +++ b/src/core/grid_based_algorithms/lbgpu.cpp @@ -100,14 +100,9 @@ bool ek_initialized = false; /*-----------------------------------------------------------*/ /** %Lattice Boltzmann update gpu called from integrate.cpp */ -void lattice_boltzmann_update_gpu(double time_step) { - - auto factor = static_cast(round(lbpar_gpu.tau / time_step)); - +void lattice_boltzmann_update_gpu(int lb_steps_per_md_step) { fluidstep += 1; - - if (fluidstep >= factor) { - + if (fluidstep >= lb_steps_per_md_step) { fluidstep = 0; lb_integrate_GPU(); } diff --git a/src/core/grid_based_algorithms/lbgpu.hpp b/src/core/grid_based_algorithms/lbgpu.hpp index c6c2c0cfe6c..ffcbd10768d 100644 --- a/src/core/grid_based_algorithms/lbgpu.hpp +++ b/src/core/grid_based_algorithms/lbgpu.hpp @@ -150,7 +150,7 @@ void lb_GPU_sanity_checks(); void lb_get_device_values_pointer(LB_rho_v_gpu **pointer_address); void lb_get_boundary_force_pointer(float **pointer_address); void lb_get_para_pointer(LB_parameters_gpu **pointer_address); -void lattice_boltzmann_update_gpu(double time_step); +void lattice_boltzmann_update_gpu(int lb_steps_per_md_step); /** Perform a full initialization of the lattice Boltzmann system. * All derived parameters and the fluid are reset to their default values. diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index dccfc906879..675da603bff 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -57,6 +57,10 @@ #include +#include +#include +#include +#include #include #ifdef VALGRIND_INSTRUMENTATION @@ -284,8 +288,12 @@ int integrate(int n_steps, int reuse_forces) { // propagate one-step functionalities if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) { - lb_lbfluid_propagate(time_step); - lb_lbcoupling_propagate(); + if (lb_lbfluid_get_lattice_switch() != ActiveLB::NONE) { + auto const lb_steps_per_md_step = + static_cast(std::round(lb_lbfluid_get_tau() / time_step)); + lb_lbfluid_propagate(lb_steps_per_md_step); + lb_lbcoupling_propagate(); + } #ifdef VIRTUAL_SITES virtual_sites()->after_lb_propagation(time_step); diff --git a/src/python/espressomd/lb.pxd b/src/python/espressomd/lb.pxd index c9be6fddbfd..19a466e95fa 100644 --- a/src/python/espressomd/lb.pxd +++ b/src/python/espressomd/lb.pxd @@ -66,7 +66,6 @@ cdef extern from "grid_based_algorithms/lb_interface.hpp": const Vector3d lb_lbfluid_get_ext_force_density() except + void lb_lbfluid_set_bulk_viscosity(double c_bulk_visc) except + double lb_lbfluid_get_bulk_viscosity() except + - void lb_lbfluid_sanity_checks() except + void lb_lbfluid_print_vtk_velocity(string filename) except + void lb_lbfluid_print_vtk_velocity(string filename, vector[int] bb1, vector[int] bb2) except + void lb_lbfluid_print_vtk_boundary(string filename) except + diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx index f950da4d2fc..5e4df1178df 100644 --- a/src/python/espressomd/lb.pyx +++ b/src/python/espressomd/lb.pyx @@ -54,7 +54,7 @@ cdef class HydrodynamicInteraction(Actor): Lattice constant. The box size in every direction must be an integer multiple of ``agrid``. tau : :obj:`float` - LB time step. The MD time step must be an integer multiple of ``tau``. + LB time step, must be an integer multiple of the MD time step. dens : :obj:`float` Fluid density. visc : :obj:`float` @@ -157,7 +157,6 @@ cdef class HydrodynamicInteraction(Actor): if "gamma_even" in self._params: python_lbfluid_set_gamma_even(self._params["gamma_even"]) - lb_lbfluid_sanity_checks() utils.handle_errors("LB fluid activation") def _get_params_from_es_core(self):