diff --git a/maintainer/benchmarks/lb.py b/maintainer/benchmarks/lb.py index 6472e14442a..46efb3abc35 100644 --- a/maintainer/benchmarks/lb.py +++ b/maintainer/benchmarks/lb.py @@ -62,7 +62,7 @@ args = parser.parse_args() # process and check arguments -n_iterations = 30 +n_iterations = 200 assert args.volume_fraction > 0, "--volume_fraction must be a positive number" assert args.volume_fraction < np.pi / (3 * np.sqrt(2)), \ "--volume_fraction exceeds the physical limit of sphere packing (~0.74)" diff --git a/maintainer/benchmarks/lj.py b/maintainer/benchmarks/lj.py index 6aed78b64eb..6f72a7bc105 100644 --- a/maintainer/benchmarks/lj.py +++ b/maintainer/benchmarks/lj.py @@ -44,7 +44,7 @@ # process and check arguments measurement_steps = int(np.round(5e6 / args.particles_per_core, -2)) -n_iterations = 30 +n_iterations = 2000 assert args.volume_fraction > 0, "volume_fraction must be a positive number" assert args.volume_fraction < np.pi / (3 * np.sqrt(2)), \ "volume_fraction exceeds the physical limit of sphere packing (~0.74)" @@ -119,10 +119,10 @@ system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=42) # tuning and equilibration -min_skin = 0.2 +min_skin = 0.1 max_skin = 1.0 print("Tune skin: {:.3f}".format(system.cell_system.tune_skin( - min_skin=min_skin, max_skin=max_skin, tol=0.05, int_steps=100))) + min_skin=min_skin, max_skin=max_skin, tol=0.005, int_steps=300))) print("Equilibration") system.integrator.run(min(5 * measurement_steps, 60000)) print("Tune skin: {:.3f}".format(system.cell_system.tune_skin( diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 957b37b92b3..298ad41c7a7 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -51,8 +51,10 @@ #include "short_range_loop.hpp" #include "system/System.hpp" #include "thermostat.hpp" + #include "thermostats/langevin_inline.hpp" #include "virtual_sites/relative.hpp" +#include #include #include @@ -184,6 +186,9 @@ void System::System::calculate_forces() { #endif // Use combined function instead of two separate calls + // Prepare LB coupling state (initiates velocity interpolation) + std::optional lb_coupling_state; + auto const elc_kernel = coulomb.pair_force_elc_kernel(); auto const coulomb_kernel = coulomb.pair_force_kernel(); auto const dipoles_kernel = dipoles.pair_force_kernel(); @@ -347,11 +352,6 @@ void System::System::calculate_forces() { // Must be done here. Forces need to be ghost-communicated immersed_boundaries->volume_conservation(*cell_structure); - if (thermostat->lb and (propagation->used_propagations & - PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE)) { - lb_couple_particles(); - } - #ifdef CUDA #ifdef CALIPER CALI_MARK_BEGIN("copy_forces_from_GPU"); @@ -362,6 +362,13 @@ void System::System::calculate_forces() { #endif #endif // CUDA + // Apply LB forces after GPU copy (hides interpolation latency) + if (thermostat->lb and (propagation->used_propagations & + PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE)) { + lb_coupling_state = lb_prepare_particle_coupling(); + lb_apply_particle_forces(*lb_coupling_state); + } + #ifdef VIRTUAL_SITES_RELATIVE if (propagation->used_propagations & (PropagationMode::TRANS_VS_RELATIVE | PropagationMode::ROT_VS_RELATIVE)) { diff --git a/src/core/lb/Solver.cpp b/src/core/lb/Solver.cpp index 300bbb0db0d..30de8d52157 100644 --- a/src/core/lb/Solver.cpp +++ b/src/core/lb/Solver.cpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include #include @@ -277,6 +278,14 @@ std::vector Solver::get_coupling_interpolated_velocities( *impl->solver); } +std::future> +Solver::get_coupling_interpolated_velocities_async( + std::vector const &pos) const { + return std::async(std::launch::async, [this, pos]() { + return get_coupling_interpolated_velocities(pos); + }); +} + void Solver::add_forces_at_pos(std::vector const &pos, std::vector const &forces) { std::visit( diff --git a/src/core/lb/Solver.hpp b/src/core/lb/Solver.hpp index 101b69ee23c..05233adef75 100644 --- a/src/core/lb/Solver.hpp +++ b/src/core/lb/Solver.hpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -190,6 +191,10 @@ struct Solver : public System::Leaf { std::vector get_coupling_interpolated_velocities( std::vector const &pos) const; + std::future> + get_coupling_interpolated_velocities_async( + std::vector const &pos) const; + void add_forces_at_pos(std::vector const &pos, std::vector const &forces); diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index 57e8e42c3fb..f35884430eb 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -207,125 +208,167 @@ Utils::Vector3d ParticleCoupling::get_noise_term(Particle const &p) const { return m_noise_pref_wo_gamma * Utils::hadamard_product(sqrt(gamma), noise); } -void ParticleCoupling::kernel(std::vector const &particles) { +/** @brief Collect particles to couple to LB and initiate obtaining Lb + * interpolated velocities */ +ParticleCouplingState +ParticleCoupling::prepare_coupling(std::vector const &particles) { +#ifdef CALIPER + CALI_CXX_MARK_FUNCTION; +#endif + ParticleCouplingState state; + state.coupled_particle_data.reserve(particles.size()); + if (particles.empty()) { - return; + return state; } - enum coupling_modes { none, particle_force, swimmer_force_on_fluid }; + auto const halo = 0.5 * m_lb.get_agrid(); auto const halo_vec = Utils::Vector3d::broadcast(halo); auto const fully_inside_lower = m_local_box.my_left() + 2. * halo_vec; auto const fully_inside_upper = m_local_box.my_right() - 2. * halo_vec; auto const halo_lower_corner = m_local_box.my_left() - halo_vec; auto const halo_upper_corner = m_local_box.my_right() + halo_vec; - std::vector positions_velocity_coupling; - std::vector positions_force_coupling; - std::vector force_coupling_forces; - std::vector positions_force_coupling_counter; - std::vector coupled_particles; + + // First pass: determine positions and coupling modes for (auto ptr : particles) { auto &p = *ptr; - auto span_size = uint8_t{1u}; auto const folded_pos = m_box_geo.folded_position(p.pos()); + + ParticleCouplingState::CoupledParticleData data; + data.particle = ptr; + data.mode = ParticleCouplingState::none; + data.force_positions_start = state.all_force_positions.size(); + data.force_positions_count = 0; + + // Temporary vector to collect positions for this particle + std::vector temp_positions; + temp_positions.reserve(8); + if (in_box(folded_pos, fully_inside_lower, fully_inside_upper)) { - positions_force_coupling.emplace_back(folded_pos); + // Determine force coupling positions + // If the folded position of the particle is further than 1/2 agrid away + // from any wall, ghosts of the particle (shifted by +- box_l in any + // coordinate) are not within the domain of this MPi rank including halos + // (-agrid/2 +local_box_min..local_box_max+agrid/2) + temp_positions.emplace_back(folded_pos); } else { - auto const old_size = positions_force_coupling.size(); + // the particle is close to a boundary. We might have to also couple + // halos of the particle shifted by +- box_l in any coordinate positions_in_halo_impl(folded_pos, halo_lower_corner, halo_upper_corner, - m_box_geo, positions_force_coupling); - auto const new_size = positions_force_coupling.size(); - span_size = static_cast(new_size - old_size); + m_box_geo, temp_positions); } - auto coupling_mode = none; + + // Check if particle should be coupled #ifdef ENGINE if (p.swimming().is_engine_force_on_fluid) { - coupling_mode = swimmer_force_on_fluid; + // swimmers excert a force on the fluid, but there is no + // coupling force based on the velocity difference between particle and + // fluid + data.mode = ParticleCouplingState::swimmer_force_on_fluid; } #endif - if (coupling_mode == none) { - for (auto end = positions_force_coupling.end(); - auto const &pos : std::views::counted(end - span_size, span_size)) { + if (data.mode == ParticleCouplingState::none) { + // Check if any position is within velocity coupling region + for (auto const &pos : temp_positions) { if (pos >= halo_lower_corner and pos < halo_upper_corner) { - positions_velocity_coupling.emplace_back(pos); - coupling_mode = particle_force; + // Coupling force on particle and fluid based on velocity difference. + // I.e., flwo velocity has to be obtained at this position. + // Note: Only one interpolation per particle is needed + // since velocities in the halo layers are identical + data.velocity_coupling_index = + state.positions_velocity_coupling.size(); + state.positions_velocity_coupling.emplace_back(pos); + data.mode = ParticleCouplingState::particle_force; break; } } } - if (coupling_mode == none) { - positions_force_coupling.erase(positions_force_coupling.end() - span_size, - positions_force_coupling.end()); - } else { - coupled_particles.emplace_back(ptr); - positions_force_coupling_counter.emplace_back(span_size); + + // Only keep particles that will be coupled + if (data.mode != ParticleCouplingState::none) { + // Add positions directly to global vector + state.all_force_positions.insert(state.all_force_positions.end(), + temp_positions.begin(), + temp_positions.end()); + data.force_positions_count = temp_positions.size(); + state.coupled_particle_data.emplace_back(std::move(data)); } } - if (coupled_particles.empty()) { + if (!state.coupled_particle_data.empty()) { + state.interpolated_velocities = m_lb.get_coupling_interpolated_velocities( + state.positions_velocity_coupling); + } + + return state; +} + +void ParticleCoupling::apply_forces(ParticleCouplingState &state) { + if (state.coupled_particle_data.empty()) { return; } - auto interpolated_velocities = - m_lb.get_coupling_interpolated_velocities(positions_velocity_coupling); auto const &domain_lower_corner = m_local_box.my_left(); auto const &domain_upper_corner = m_local_box.my_right(); - auto it_interpolated_velocities = interpolated_velocities.begin(); - auto it_positions_force_coupling = positions_force_coupling.begin(); - auto it_positions_velocity_coupling = positions_velocity_coupling.begin(); - auto it_positions_force_coupling_counter = - positions_force_coupling_counter.begin(); - for (auto ptr : coupled_particles) { - auto &p = *ptr; - auto coupling_mode = particle_force; -#ifdef ENGINE - if (p.swimming().is_engine_force_on_fluid) { - coupling_mode = swimmer_force_on_fluid; - } -#endif + + // Pre-allocate forces vector with same size as positions + std::vector all_forces; + all_forces.reserve(state.all_force_positions.size()); + + // Calculate and apply forces + for (auto &data : state.coupled_particle_data) { + auto &p = *data.particle; Utils::Vector3d force_on_particle = {}; - if (coupling_mode == particle_force) { + + if (data.mode == ParticleCouplingState::particle_force) { #ifndef THERMOSTAT_PER_PARTICLE if (m_thermostat.gamma > 0.) #endif { - auto &v_fluid = *it_interpolated_velocities; + auto v_fluid = + state.interpolated_velocities[*data.velocity_coupling_index]; + auto const &vel_pos = + state.positions_velocity_coupling[*data.velocity_coupling_index]; + if (m_box_geo.type() == BoxType::LEES_EDWARDS) { // Account for the case where the interpolated velocity has been read // from a ghost of the particle across the LE boundary (or vice versa) // Then the particle velocity is shifted by +,- the LE shear velocity - auto const vel_correction = lees_edwards_vel_shift( - *it_positions_velocity_coupling, p.pos(), m_box_geo); + auto const vel_correction = + lees_edwards_vel_shift(vel_pos, p.pos(), m_box_geo); v_fluid += vel_correction; } auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid); auto const random_force = get_noise_term(p); force_on_particle = drag_force + random_force; } - ++it_interpolated_velocities; - ++it_positions_velocity_coupling; } auto force_on_fluid = -force_on_particle; #ifdef ENGINE - if (coupling_mode == swimmer_force_on_fluid) { + if (data.mode == ParticleCouplingState::swimmer_force_on_fluid) { force_on_fluid = p.calc_director() * p.swimming().f_swim; } #endif - auto const span_size = *it_positions_force_coupling_counter; - ++it_positions_force_coupling_counter; - for (uint8_t i{0u}; i < span_size; ++i) { - auto &pos = *it_positions_force_coupling; + // Apply forces to particle and collect forces for LB + bool particle_force_applied = false; + for (size_t i = 0; i < data.force_positions_count; ++i) { + auto const &pos = + state.all_force_positions[data.force_positions_start + i]; if (pos >= domain_lower_corner and pos < domain_upper_corner) { /* Particle is in our LB volume, so this node * is responsible to adding its force */ - p.force() += force_on_particle; + if (!particle_force_applied) { + p.force() += force_on_particle; + particle_force_applied = true; + } } - force_coupling_forces.emplace_back(force_on_fluid); - ++it_positions_force_coupling; + all_forces.emplace_back(force_on_fluid); } } - m_lb.add_forces_at_pos(positions_force_coupling, force_coupling_forces); + + m_lb.add_forces_at_pos(state.all_force_positions, all_forces); } #if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY) @@ -344,15 +387,16 @@ static void lb_coupling_sanity_checks(Particle const &p) { } // namespace LB -void System::System::lb_couple_particles() { +LB::ParticleCouplingState System::System::lb_prepare_particle_coupling() { #ifdef CALIPER CALI_CXX_MARK_FUNCTION; #endif assert(thermostat->lb != nullptr); + if (thermostat->lb->couple_to_md) { if (not lb.is_solver_set()) { runtimeErrorMsg() << "The LB thermostat requires a LB fluid"; - return; + return {}; } auto const real_particles = cell_structure->local_particles(); auto const ghost_particles = cell_structure->ghost_particles(); @@ -370,6 +414,32 @@ void System::System::lb_couple_particles() { } } } - coupling.kernel(particles); + + return coupling.prepare_coupling(particles); } + return {}; +} + +void System::System::lb_apply_particle_forces( + LB::ParticleCouplingState &state) { +#ifdef CALIPER + CALI_CXX_MARK_FUNCTION; +#endif + assert(thermostat->lb != nullptr); + if (thermostat->lb->couple_to_md && !state.coupled_particle_data.empty()) { + if (not lb.is_solver_set()) { + runtimeErrorMsg() << "The LB thermostat requires a LB fluid"; + return; + } + LB::ParticleCoupling coupling{*thermostat->lb, lb, *box_geo, *local_geo}; + coupling.apply_forces(state); + } +} + +void System::System::lb_couple_particles() { +#ifdef CALIPER + CALI_CXX_MARK_FUNCTION; +#endif + auto coupling_state = lb_prepare_particle_coupling(); + lb_apply_particle_forces(coupling_state); } diff --git a/src/core/lb/particle_coupling.hpp b/src/core/lb/particle_coupling.hpp index ff6f29178fd..6e30ca8caa1 100644 --- a/src/core/lb/particle_coupling.hpp +++ b/src/core/lb/particle_coupling.hpp @@ -66,6 +66,34 @@ Utils::Vector3d lb_drag_force(LB::Solver const &lb, double lb_gamma, namespace LB { +// State structure to hold particle coupling data between phases +struct ParticleCouplingState { + enum coupling_modes { none, particle_force, swimmer_force_on_fluid }; + + // Structure to hold all data related to a coupled particle + struct CoupledParticleData { + Particle *particle; + size_t + force_positions_start; // Start index in global force positions vector + size_t force_positions_count; // Number of force positions for this particle + std::optional velocity_coupling_index; + coupling_modes mode; + }; + + std::vector coupled_particle_data; + std::vector positions_velocity_coupling; + std::vector + all_force_positions; // Direct storage for all force positions + std::vector interpolated_velocities; + + void clear() { + coupled_particle_data.clear(); + positions_velocity_coupling.clear(); + all_force_positions.clear(); + interpolated_velocities.clear(); + } +}; + class ParticleCoupling { LBThermostat const &m_thermostat; LB::Solver &m_lb; @@ -91,7 +119,11 @@ class ParticleCoupling { } Utils::Vector3d get_noise_term(Particle const &p) const; - void kernel(std::vector const &particles); + + // Split kernel into two phases + ParticleCouplingState + prepare_coupling(std::vector const &particles); + void apply_forces(ParticleCouplingState &state); }; /** diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index 4797fde87eb..dcb1d103a78 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -65,6 +65,9 @@ class AutoUpdateAccumulators; namespace Constraints { class Constraints; } +namespace LB { +struct ParticleCouplingState; +} struct NptIsoParameters; struct InstantaneousPressure; @@ -245,6 +248,12 @@ class System : public std::enable_shared_from_this { /** @brief Calculate particle-lattice interactions. */ void lb_couple_particles(); + /** @brief Prepare particle-lattice coupling (phase 1). */ + LB::ParticleCouplingState lb_prepare_particle_coupling(); + + /** @brief Apply particle-lattice forces (phase 2). */ + void lb_apply_particle_forces(LB::ParticleCouplingState &state); + /** \name Hook procedures * These procedures are called if several significant changes to * the system happen which may make a reinitialization of subsystems diff --git a/src/core/unit_tests/lb_particle_coupling_test.cpp b/src/core/unit_tests/lb_particle_coupling_test.cpp index 4041e885e4c..d2bbaf80f04 100644 --- a/src/core/unit_tests/lb_particle_coupling_test.cpp +++ b/src/core/unit_tests/lb_particle_coupling_test.cpp @@ -275,7 +275,8 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, swimmer_force, bdata::make(kTs), kT) { { if (in_local_halo(local_box, p.pos(), params.agrid)) { LB::ParticleCoupling coupling{thermostat, lb, box_geo, local_box}; - coupling.kernel({&p}); + auto state = coupling.prepare_coupling({&p}); + coupling.apply_forces(state); auto const interpolated = LB::get_force_to_be_applied(p.pos()); auto const expected = params.force_md_to_lb(Utils::Vector3d{0., 0., p.swimming().f_swim}); @@ -339,7 +340,8 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, particle_coupling, bdata::make(kTs), // coupling { if (in_local_halo(local_box, p.pos(), params.agrid)) { - coupling.kernel({&p}); + auto state = coupling.prepare_coupling({&p}); + coupling.apply_forces(state); BOOST_CHECK_SMALL((p.force() - expected).norm(), eps); auto const interpolated = -LB::get_force_to_be_applied(p.pos());