30#include "system/System.hpp"
31#include "thermostat.hpp"
54static auto is_solver_set(std::unique_ptr<Solver::Implementation>
const &ptr) {
55 return ptr !=
nullptr and ptr->solver.has_value();
58static void check_solver(std::unique_ptr<Solver::Implementation>
const &ptr) {
73 std::visit([](
auto &ptr) { ptr->propagate(); }, *impl->solver);
79 std::visit([&](
auto &ptr) { ptr->sanity_checks(system); }, *impl->solver);
85 std::visit([=](
auto &ptr) { ptr->veto_time_step(time_step); },
92 std::visit([=](
auto &ptr) { ptr->veto_kT(kT); }, *impl->solver);
97 unsigned int shear_plane_normal)
const {
99 auto const callback = [=](
auto &ptr) {
100 ptr->lebc_sanity_checks(shear_direction, shear_plane_normal);
102 std::visit(callback, *impl->solver);
108 auto &solver = *impl->solver;
109 std::visit([](
auto &ptr) { ptr->on_cell_structure_change(); }, solver);
115 std::visit([](
auto const &ptr) { ptr->veto_boxl_change(); }, *impl->solver);
121 std::visit([](
auto &ptr) { ptr->on_boxl_change(); }, *impl->solver);
127 std::visit([](
auto &ptr) { ptr->on_node_grid_change(); }, *impl->solver);
133 std::visit([](
auto &ptr) { ptr->on_timestep_change(); }, *impl->solver);
139 std::visit([](
auto &ptr) { ptr->on_temperature_change(); }, *impl->solver);
145 return std::visit([](
auto &ptr) {
return ptr->is_gpu(); }, *impl->solver);
150 return std::visit([](
auto &ptr) {
return ptr->get_agrid(); }, *impl->solver);
155 return std::visit([](
auto &ptr) {
return ptr->get_tau(); }, *impl->solver);
160 return std::visit([](
auto &ptr) {
return ptr->get_kT(); }, *impl->solver);
165 return std::visit([](
auto &ptr) {
return ptr->get_pressure_tensor(); },
169std::optional<Utils::Vector3d>
177 auto const lb_pos = box_geo.folded_position(
pos) * m_conv.
pos_to_lb;
178 return ptr->get_velocity_at_pos(lb_pos,
false);
188 auto const lb_pos = box_geo.folded_position(
pos) * m_conv.
pos_to_lb;
189 return ptr->get_density_at_pos(lb_pos,
false);
198 auto const res = ptr->get_velocity_at_pos(
pos * m_conv.
pos_to_lb,
true);
206 std::vector<Utils::Vector3d>
const &
pos)
const {
209 std::vector<Utils::Vector3d> pos_lb;
210 pos_lb.reserve(
pos.size());
211 for (
auto const &pos_md :
pos) {
212 pos_lb.emplace_back(pos_md * m_conv.
pos_to_lb);
214 auto res = ptr->get_velocities_at_pos(pos_lb);
215 for (
auto &v :
res) {
224 std::vector<Utils::Vector3d>
const &forces) {
227 std::vector<Utils::Vector3d> pos_lb;
228 std::vector<Utils::Vector3d> force_lb;
229 pos_lb.reserve(
pos.size());
230 force_lb.reserve(
pos.size());
231 for (
auto const &pos_md :
pos) {
232 pos_lb.emplace_back(pos_md * m_conv.
pos_to_lb);
234 for (
auto const &force_md : forces) {
235 force_lb.emplace_back(force_md * m_conv.
force_to_lb);
237 ptr->add_forces_at_pos(pos_lb, force_lb);
248 throw std::runtime_error(
"Cannot apply force to LB");
256 return std::visit([](
auto const &ptr) {
return ptr->get_momentum(); },
260template <>
void Solver::set<LBNone>(std::shared_ptr<LBNone> lb_instance) {
262 assert(not impl->solver.has_value());
263 impl->solver = lb_instance;
268void Solver::set<LBWalberla>(std::shared_ptr<LBWalberlaBase> lb_fluid,
269 std::shared_ptr<LBWalberlaParams> lb_params) {
271 assert(not impl->solver.has_value());
273 auto lb_instance = std::make_shared<LBWalberla>(lb_fluid, lb_params);
274 lb_instance->sanity_checks(system);
275 auto const &lebc = system.box_geo->lees_edwards_bc();
276 lb_fluid->check_lebc(lebc.shear_direction, lebc.shear_plane_normal);
277 impl->solver = lb_instance;
278 auto const agrid = lb_instance->get_agrid();
279 auto const tau = lb_instance->get_tau();
280 m_conv =
Conversions{1. / agrid, agrid / tau, tau * tau / agrid};
LBWalberlaBase provides the public interface of the LB waLBerla bridge.
Vector implementation and trait types for boost qvm interoperability.
__shared__ int pos[MAXDEPTH *THREADS5/WARPSIZE]
std::shared_ptr< BoxGeometry > box_geo
This file contains the defaults for ESPResSo.
static auto is_solver_set(std::unique_ptr< Solver::Implementation > const &ptr)
static void check_solver(std::unique_ptr< Solver::Implementation > const &ptr)
void veto_boxl_change() const
void lebc_sanity_checks(unsigned int shear_direction, unsigned int shear_plane_normal) const
Perform LB LEbc parameter checks.
void add_force_density(Utils::Vector3d const &pos, Utils::Vector3d const &force_density)
Add a force density to the fluid at the given position.
Utils::Vector3d get_coupling_interpolated_velocity(Utils::Vector3d const &pos) const
Calculate the interpolated fluid velocity in MD units.
void sanity_checks() const
Perform LB parameter and boundary velocity checks.
void on_timestep_change()
void propagate()
Propagate the LB fluid.
void reset()
Remove the LB solver.
std::optional< Utils::Vector3d > get_interpolated_velocity(Utils::Vector3d const &pos) const
Calculate the interpolated fluid velocity in LB units.
bool is_solver_set() const
Return true if a LB solver is active.
double get_tau() const
Get the LB time step.
void on_temperature_change()
std::optional< double > get_interpolated_density(Utils::Vector3d const &pos) const
Calculate the interpolated fluid density in LB units.
void add_forces_at_pos(std::vector< Utils::Vector3d > const &pos, std::vector< Utils::Vector3d > const &forces)
void veto_kT(double kT) const
Check if a thermostat is compatible with the LB temperature.
double get_agrid() const
Get the LB grid spacing.
Utils::VectorXd< 9 > get_pressure_tensor() const
void veto_time_step(double time_step) const
Check if a MD time step is compatible with the LB tau.
void on_node_grid_change()
Utils::Vector3d get_momentum() const
void on_cell_structure_change()
std::vector< Utils::Vector3d > get_coupling_interpolated_velocities(std::vector< Utils::Vector3d > const &pos) const
double get_kT() const
Get the thermal energy parameter of the LB fluid.