43#include <boost/mpi.hpp>
44#include <boost/mpi/collectives/all_reduce.hpp>
45#include <boost/mpi/collectives/broadcast.hpp>
60std::unordered_map<std::string, int>
const LBVTKHandle::obs_map = {
68 if (
name ==
"activate") {
75 if (
name ==
"deactivate") {
82 if (
name ==
"add_force_at_pos") {
84 auto const pos = get_value<Utils::Vector3d>(
params,
"pos");
85 auto const f = get_value<Utils::Vector3d>(
params,
"force");
86 auto const folded_pos = box_geo.folded_position(
pos);
90 if (
name ==
"get_interpolated_velocity") {
91 auto const pos = get_value<Utils::Vector3d>(
params,
"pos");
92 return get_interpolated_velocity(
pos);
94 if (
name ==
"get_pressure_tensor") {
95 return get_average_pressure_tensor();
97 if (
name ==
"load_checkpoint") {
98 auto const path = get_value<std::string>(
params,
"path");
99 auto const mode = get_value<int>(
params,
"mode");
100 load_checkpoint(path, mode);
103 if (
name ==
"save_checkpoint") {
104 auto const path = get_value<std::string>(
params,
"path");
105 auto const mode = get_value<int>(
params,
"mode");
106 save_checkpoint(path, mode);
109 if (
name ==
"clear_boundaries") {
114 if (
name ==
"add_boundary_from_shape") {
120 if (
name ==
"get_lattice_speed") {
128 auto const visc = get_value<double>(
params,
"kinematic_viscosity");
129 auto const dens = get_value<double>(
params,
"density");
130 auto const precision = get_value<bool>(
params,
"single_precision");
131 auto const lb_lattice =
m_lattice->lattice();
139 auto const visc = get_value<double>(
params,
"kinematic_viscosity");
140 auto const dens = get_value<double>(
params,
"density");
141 auto const precision = get_value<bool>(
params,
"single_precision");
142 auto const lb_lattice =
m_lattice->lattice();
150 m_lattice = get_value<std::shared_ptr<LatticeWalberla>>(
params,
"lattice");
152 get_value_or<decltype(m_vtk_writers)>(
params,
"vtk_writers", {});
153 auto const tau = get_value<double>(
params,
"tau");
154 auto const agrid = get_value<double>(
m_lattice->get_parameter(
"agrid"));
155 auto const visc = get_value<double>(
params,
"kinematic_viscosity");
156 auto const dens = get_value<double>(
params,
"density");
157 auto const kT = get_value<double>(
params,
"kT");
158 auto const ext_f = get_value<Utils::Vector3d>(
params,
"ext_force_density");
159 m_lb_params = std::make_shared<::LB::LBWalberlaParams>(agrid, tau);
164 throw std::domain_error(
"Parameter 'tau' must be > 0");
167 m_conv_visc = Utils::int_pow<1>(tau) / Utils::int_pow<2>(agrid);
168 m_conv_energy = Utils::int_pow<2>(tau) / Utils::int_pow<2>(agrid);
170 m_conv_speed = Utils::int_pow<1>(tau) / Utils::int_pow<1>(agrid);
171 m_conv_press = Utils::int_pow<2>(tau) * Utils::int_pow<1>(agrid);
172 m_conv_force = Utils::int_pow<2>(tau) / Utils::int_pow<1>(agrid);
179 throw std::domain_error(
"Parameter 'seed' must be >= 0");
182 throw std::domain_error(
"Parameter 'kT' must be >= 0");
185 throw std::domain_error(
"Parameter 'density' must be > 0");
188 throw std::domain_error(
"Parameter 'kinematic_viscosity' must be >= 0");
192 if (
auto le_protocol = system.lees_edwards->get_protocol()) {
194 throw std::runtime_error(
195 "Lees-Edwards LB doesn't support thermalization");
197 auto const &le_bc = system.box_geo->lees_edwards_bc();
198 auto lees_edwards_object = std::make_unique<LeesEdwardsPack>(
199 le_bc.shear_direction, le_bc.shear_plane_normal,
200 [
this, le_protocol, &system]() {
201 return get_pos_offset(system.get_sim_time(), *le_protocol) /
202 m_lb_params->get_agrid();
204 [
this, le_protocol, &system]() {
205 return get_shear_velocity(system.get_sim_time(), *le_protocol) *
206 (m_lb_params->get_tau() / m_lb_params->get_agrid());
208 m_instance->set_collision_model(std::move(lees_edwards_object));
220std::vector<Variant> LBFluid::get_average_pressure_tensor()
const {
224 std::copy(tensor_flat.begin(), tensor_flat.end(), tensor.m_data.begin());
225 return std::vector<Variant>{tensor.row<0>().as_vector(),
226 tensor.row<1>().as_vector(),
227 tensor.row<2>().as_vector()};
233 auto const result =
m_instance->get_velocity_at_pos(lb_pos);
238void LBFluid::load_checkpoint(std::string
const &filename,
int mode) {
241 auto const read_metadata = [&lb_obj](CheckpointFile &cpfile) {
242 auto const expected_grid_size = lb_obj.get_lattice().get_grid_dimensions();
243 auto const expected_pop_size = lb_obj.stencil_size();
245 std::size_t read_pop_size;
246 cpfile.read(read_grid_size);
247 cpfile.read(read_pop_size);
248 if (read_grid_size != expected_grid_size) {
249 std::stringstream message;
250 message <<
"grid dimensions mismatch, read [" << read_grid_size <<
"], "
251 <<
"expected [" << expected_grid_size <<
"].";
252 throw std::runtime_error(message.str());
254 if (read_pop_size != expected_pop_size) {
255 throw std::runtime_error(
"population size mismatch, read " +
256 std::to_string(read_pop_size) +
", expected " +
257 std::to_string(expected_pop_size) +
".");
261 auto const read_data = [&lb_obj](CheckpointFile &cpfile) {
262 auto const grid_size = lb_obj.get_lattice().get_grid_dimensions();
268 for (
int i = 0; i < i_max; i++) {
269 for (
int j = 0; j < j_max; j++) {
270 for (
int k = 0; k < k_max; k++) {
278 lb_obj.set_node_population(ind, cpnode.
populations);
281 lb_obj.set_node_velocity_at_boundary(ind, cpnode.
slip_velocity);
288 auto const on_success = [&lb_obj]() {
289 lb_obj.ghost_communication();
290 lb_obj.reallocate_ubb_field();
294 read_data, on_success);
297void LBFluid::save_checkpoint(std::string
const &filename,
int mode) {
300 auto const write_metadata = [&lb_obj,
301 mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
303 auto const grid_size = lb_obj.get_lattice().get_grid_dimensions();
304 auto const pop_size = lb_obj.stencil_size();
306 cpfile_ptr->write(grid_size);
307 cpfile_ptr->write(pop_size);
312 auto const on_failure = [](std::shared_ptr<CheckpointFile>
const &,
320 auto const write_data = [&lb_obj,
321 mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
323 auto const get_node_checkpoint =
324 [&](
Utils::Vector3i const &ind) -> std::optional<LBWalberlaNodeState> {
325 auto const pop = lb_obj.get_node_population(ind);
326 auto const laf = lb_obj.get_node_last_applied_force(ind);
327 auto const lbb = lb_obj.get_node_is_boundary(ind);
328 auto const vbb = lb_obj.get_node_velocity_at_boundary(ind);
329 if (pop and laf and lbb and ((*lbb) ? vbb.has_value() :
true)) {
342 auto failure =
false;
345 auto const unit_test_mode = (mode !=
static_cast<int>(
CptMode::ascii)) and
347 auto const grid_size = lb_obj.get_lattice().get_grid_dimensions();
352 for (
int i = 0; i < i_max; i++) {
353 for (
int j = 0; j < j_max; j++) {
354 for (
int k = 0; k < k_max; k++) {
356 auto const result = get_node_checkpoint(ind);
357 if (!unit_test_mode) {
358 assert(1 == boost::mpi::all_reduce(comm,
static_cast<int>(!!result),
360 "Incorrect number of return values");
366 comm.recv(boost::mpi::any_source, 42, cpnode);
368 auto &cpfile = *cpfile_ptr;
375 boost::mpi::broadcast(comm, failure, 0);
378 comm.send(0, 42, *result);
380 boost::mpi::broadcast(comm, failure, 0);
391 write_data, on_failure);
Vector implementation and trait types for boost qvm interoperability.
__shared__ int pos[MAXDEPTH *THREADS5/WARPSIZE]
virtual void parallel_try_catch(std::function< void()> const &cb) const =0
virtual bool is_head_node() const =0
virtual boost::mpi::communicator const & get_comm() const =0
boost::string_ref name() const
Context * context() const
Responsible context.
void make_instance(VariantMap const ¶ms) override
void make_instance(VariantMap const ¶ms) override
void do_construct(VariantMap const ¶ms) override
Variant do_call_method(std::string const &name, VariantMap const ¶ms) override
std::shared_ptr<::LB::LBWalberlaParams > m_lb_params
::LatticeModel::units_map get_latice_to_md_units_conversion() const override
std::shared_ptr< ::LBWalberlaBase > m_instance
Variant do_call_method(std::string const &method_name, VariantMap const ¶ms) override
std::vector< std::shared_ptr< LBVTKHandle > > m_vtk_writers
std::shared_ptr< LatticeWalberla > m_lattice
virtual void make_instance(VariantMap const ¶ms)=0
void on_lb_boundary_conditions_change()
Called when the LB boundary conditions change (geometry, slip velocity, or both).
std::shared_ptr< BoxGeometry > box_geo
This file contains the defaults for ESPResSo.
std::shared_ptr< LBWalberlaBase > new_lb_walberla_cpu(std::shared_ptr< LatticeWalberla > const &lattice, double viscosity, double density, bool single_precision)
std::shared_ptr< LBWalberlaBase > new_lb_walberla_gpu(std::shared_ptr< LatticeWalberla > const &lattice, double viscosity, double density, bool single_precision)
Matrix implementation and trait types for boost qvm interoperability.
void save_checkpoint_common(Context const &context, std::string const classname, std::string const &filename, int mode, F1 const write_metadata, F2 const write_data, F3 const on_failure)
void unit_test_handle(int mode)
Inject code for unit tests.
void load_checkpoint_common(Context const &context, std::string const classname, std::string const &filename, int mode, F1 const read_metadata, F2 const read_data, F3 const on_success)
T get_value(Variant const &v)
Extract value of specific type T from a Variant.
std::unordered_map< std::string, Variant > VariantMap
T mpi_reduce_sum(boost::mpi::communicator const &comm, T const &result)
Reduce object by sum on the head node.
boost::make_recursive_variant< None, bool, int, std::size_t, double, std::string, ObjectRef, Utils::Vector3b, Utils::Vector3i, Utils::Vector2d, Utils::Vector3d, Utils::Vector4d, std::vector< int >, std::vector< double >, std::vector< boost::recursive_variant_ >, std::unordered_map< int, boost::recursive_variant_ >, std::unordered_map< std::string, boost::recursive_variant_ > >::type Variant
Possible types for parameters.
T reduce_optional(boost::mpi::communicator const &comm, std::optional< T > const &result)
Reduce an optional on the head node.
static FUNC_PREFIX double *RESTRICT int64_t const int64_t const int64_t const int64_t const int64_t const int64_t const int64_t const int64_t const int64_t const int64_t const int64_t const double grid_size
static SteepestDescentParameters params
Currently active steepest descent instance.
Checkpoint data for a LB node.
Utils::Vector3d last_applied_force
std::vector< double > populations
Utils::Vector3d slip_velocity
void reset()
Remove the LB solver.
void set(Args... args)
Set the LB solver.
Matrix representation with static size.