ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
LBFluid.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2021-2023 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19#include "config/config.hpp"
20
21#ifdef WALBERLA
22
23#include "LBFluid.hpp"
26
27#include "core/BoxGeometry.hpp"
32
34
38
39#include <utils/Vector.hpp>
40#include <utils/matrix.hpp>
42
43#include <boost/mpi.hpp>
44#include <boost/mpi/collectives/all_reduce.hpp>
45#include <boost/mpi/collectives/broadcast.hpp>
46
47#include <algorithm>
48#include <cassert>
49#include <cstddef>
50#include <functional>
51#include <memory>
52#include <optional>
53#include <sstream>
54#include <stdexcept>
55#include <string>
56#include <vector>
57
59
60std::unordered_map<std::string, int> const LBVTKHandle::obs_map = {
61 {"density", static_cast<int>(OutputVTK::density)},
62 {"velocity_vector", static_cast<int>(OutputVTK::velocity_vector)},
63 {"pressure_tensor", static_cast<int>(OutputVTK::pressure_tensor)},
64};
65
66Variant LBFluid::do_call_method(std::string const &name,
67 VariantMap const &params) {
68 if (name == "activate") {
69 context()->parallel_try_catch([this]() {
71 });
72 m_is_active = true;
73 return {};
74 }
75 if (name == "deactivate") {
76 if (m_is_active) {
78 m_is_active = false;
79 }
80 return {};
81 }
82 if (name == "add_force_at_pos") {
83 auto const &box_geo = *::System::get_system().box_geo;
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);
87 m_instance->add_force_at_pos(folded_pos * m_conv_dist, f * m_conv_force);
88 return {};
89 }
90 if (name == "get_interpolated_velocity") {
91 auto const pos = get_value<Utils::Vector3d>(params, "pos");
92 return get_interpolated_velocity(pos);
93 }
94 if (name == "get_pressure_tensor") {
95 return get_average_pressure_tensor();
96 }
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);
101 return {};
102 }
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);
107 return {};
108 }
109 if (name == "clear_boundaries") {
110 m_instance->clear_boundaries();
112 return {};
113 }
114 if (name == "add_boundary_from_shape") {
115 m_instance->update_boundary_from_shape(
116 get_value<std::vector<int>>(params, "raster"),
117 get_value<std::vector<double>>(params, "values"));
118 return {};
119 }
120 if (name == "get_lattice_speed") {
121 return 1. / m_conv_speed;
122 }
123
125}
126
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();
132 auto const lb_visc = m_conv_visc * visc;
133 auto const lb_dens = m_conv_dens * dens;
135}
136
137#ifdef CUDA
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();
143 auto const lb_visc = m_conv_visc * visc;
144 auto const lb_dens = m_conv_dens * dens;
146}
147#endif // CUDA
148
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);
160 m_is_active = false;
161 auto const seed = get_value<int>(params, "seed");
162 context()->parallel_try_catch([&]() {
163 if (tau <= 0.) {
164 throw std::domain_error("Parameter 'tau' must be > 0");
165 }
166 m_conv_dist = 1. / agrid;
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);
169 m_conv_dens = Utils::int_pow<3>(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);
173 m_conv_force_dens = Utils::int_pow<2>(tau) * Utils::int_pow<2>(agrid);
174 auto const lb_visc = m_conv_visc * visc;
175 auto const lb_dens = m_conv_dens * dens;
176 auto const lb_kT = m_conv_energy * kT;
177 auto const lb_ext_f = m_conv_force_dens * ext_f;
178 if (seed < 0) {
179 throw std::domain_error("Parameter 'seed' must be >= 0");
180 }
181 if (lb_kT < 0.) {
182 throw std::domain_error("Parameter 'kT' must be >= 0");
183 }
184 if (lb_dens <= 0.) {
185 throw std::domain_error("Parameter 'density' must be > 0");
186 }
187 if (lb_visc < 0.) {
188 throw std::domain_error("Parameter 'kinematic_viscosity' must be >= 0");
189 }
192 static_cast<unsigned int>(seed));
193 m_instance->set_external_force(lb_ext_f);
194 m_instance->ghost_communication();
195 for (auto &vtk : m_vtk_writers) {
196 vtk->attach_to_lattice(m_instance, get_latice_to_md_units_conversion());
197 }
198 });
199}
200
201std::vector<Variant> LBFluid::get_average_pressure_tensor() const {
202 auto const local = m_instance->get_pressure_tensor() / m_conv_press;
203 auto const tensor_flat = mpi_reduce_sum(context()->get_comm(), local);
205 std::copy(tensor_flat.begin(), tensor_flat.end(), tensor.m_data.begin());
206 return std::vector<Variant>{tensor.row<0>().as_vector(),
207 tensor.row<1>().as_vector(),
208 tensor.row<2>().as_vector()};
209}
210
211Variant LBFluid::get_interpolated_velocity(Utils::Vector3d const &pos) const {
212 auto const &box_geo = *::System::get_system().box_geo;
213 auto const lb_pos = box_geo.folded_position(pos) * m_conv_dist;
214 auto const result = m_instance->get_velocity_at_pos(lb_pos);
215 return Utils::Mpi::reduce_optional(context()->get_comm(), result) /
217}
218
219void LBFluid::load_checkpoint(std::string const &filename, int mode) {
220 auto &lb_obj = *m_instance;
221
222 auto const read_metadata = [&lb_obj](CheckpointFile &cpfile) {
223 auto const expected_grid_size = lb_obj.get_lattice().get_grid_dimensions();
224 auto const expected_pop_size = lb_obj.stencil_size();
226 std::size_t read_pop_size;
228 cpfile.read(read_pop_size);
230 std::stringstream message;
231 message << "grid dimensions mismatch, read [" << read_grid_size << "], "
232 << "expected [" << expected_grid_size << "].";
233 throw std::runtime_error(message.str());
234 }
236 throw std::runtime_error("population size mismatch, read " +
237 std::to_string(read_pop_size) + ", expected " +
238 std::to_string(expected_pop_size) + ".");
239 }
240 };
241
242 auto const read_data = [&lb_obj](CheckpointFile &cpfile) {
243 auto const grid_size = lb_obj.get_lattice().get_grid_dimensions();
244 auto const i_max = grid_size[0];
245 auto const j_max = grid_size[1];
246 auto const k_max = grid_size[2];
248 cpnode.populations.resize(lb_obj.stencil_size());
249 for (int i = 0; i < i_max; i++) {
250 for (int j = 0; j < j_max; j++) {
251 for (int k = 0; k < k_max; k++) {
252 auto const ind = Utils::Vector3i{{i, j, k}};
253 cpfile.read(cpnode.populations);
254 cpfile.read(cpnode.last_applied_force);
255 cpfile.read(cpnode.is_boundary);
256 if (cpnode.is_boundary) {
257 cpfile.read(cpnode.slip_velocity);
258 }
259 lb_obj.set_node_population(ind, cpnode.populations);
260 lb_obj.set_node_last_applied_force(ind, cpnode.last_applied_force);
261 if (cpnode.is_boundary) {
262 lb_obj.set_node_velocity_at_boundary(ind, cpnode.slip_velocity);
263 }
264 }
265 }
266 }
267 };
268
269 auto const on_success = [&lb_obj]() {
270 lb_obj.ghost_communication();
271 lb_obj.reallocate_ubb_field();
272 };
273
276}
277
278void LBFluid::save_checkpoint(std::string const &filename, int mode) {
279 auto &lb_obj = *m_instance;
280
281 auto const write_metadata = [&lb_obj,
282 mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
283 Context const &context) {
284 auto const grid_size = lb_obj.get_lattice().get_grid_dimensions();
285 auto const pop_size = lb_obj.stencil_size();
286 if (context.is_head_node()) {
287 cpfile_ptr->write(grid_size);
288 cpfile_ptr->write(pop_size);
290 }
291 };
292
293 auto const on_failure = [](std::shared_ptr<CheckpointFile> const &,
294 Context const &context) {
295 if (context.is_head_node()) {
296 auto failure = true;
297 boost::mpi::broadcast(context.get_comm(), failure, 0);
298 }
299 };
300
301 auto const write_data = [&lb_obj,
302 mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
303 Context const &context) {
304 auto const get_node_checkpoint =
305 [&](Utils::Vector3i const &ind) -> std::optional<LBWalberlaNodeState> {
306 auto const pop = lb_obj.get_node_population(ind);
307 auto const laf = lb_obj.get_node_last_applied_force(ind);
308 auto const lbb = lb_obj.get_node_is_boundary(ind);
309 auto const vbb = lb_obj.get_node_velocity_at_boundary(ind);
310 if (pop and laf and lbb and ((*lbb) ? vbb.has_value() : true)) {
313 cpnode.last_applied_force = *laf;
314 cpnode.is_boundary = *lbb;
315 if (*lbb) {
316 cpnode.slip_velocity = *vbb;
317 }
318 return {cpnode};
319 }
320 return std::nullopt;
321 };
322
323 auto failure = false;
324 auto const &comm = context.get_comm();
325 auto const is_head_node = context.is_head_node();
326 auto const unit_test_mode = (mode != static_cast<int>(CptMode::ascii)) and
327 (mode != static_cast<int>(CptMode::binary));
328 auto const grid_size = lb_obj.get_lattice().get_grid_dimensions();
329 auto const i_max = grid_size[0];
330 auto const j_max = grid_size[1];
331 auto const k_max = grid_size[2];
333 for (int i = 0; i < i_max; i++) {
334 for (int j = 0; j < j_max; j++) {
335 for (int k = 0; k < k_max; k++) {
336 auto const ind = Utils::Vector3i{{i, j, k}};
337 auto const result = get_node_checkpoint(ind);
338 if (!unit_test_mode) {
339 assert(1 == boost::mpi::all_reduce(comm, static_cast<int>(!!result),
340 std::plus<>()) &&
341 "Incorrect number of return values");
342 }
343 if (is_head_node) {
344 if (result) {
345 cpnode = *result;
346 } else {
347 comm.recv(boost::mpi::any_source, 42, cpnode);
348 }
349 auto &cpfile = *cpfile_ptr;
350 cpfile.write(cpnode.populations);
351 cpfile.write(cpnode.last_applied_force);
352 cpfile.write(cpnode.is_boundary);
353 if (cpnode.is_boundary) {
354 cpfile.write(cpnode.slip_velocity);
355 }
356 boost::mpi::broadcast(comm, failure, 0);
357 } else {
358 if (result) {
359 comm.send(0, 42, *result);
360 }
361 boost::mpi::broadcast(comm, failure, 0);
362 if (failure) {
363 return;
364 }
365 }
366 }
367 }
368 }
369 };
370
373}
374
375} // namespace ScriptInterface::walberla
376
377#endif // WALBERLA
Vector implementation and trait types for boost qvm interoperability.
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 &params) override
Definition LBFluid.cpp:127
void make_instance(VariantMap const &params) override
Definition LBFluid.cpp:138
void do_construct(VariantMap const &params) override
Definition LBFluid.cpp:149
Variant do_call_method(std::string const &name, VariantMap const &params) override
Definition LBFluid.cpp:66
std::shared_ptr<::LB::LBWalberlaParams > m_lb_params
Definition LBFluid.hpp:52
::LatticeModel::units_map get_latice_to_md_units_conversion() const override
Definition LBFluid.hpp:124
Variant do_call_method(std::string const &method_name, VariantMap const &params) override
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
Definition Variant.hpp:69
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.
Definition Variant.hpp:67
System & get_system()
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.
std::vector< double > populations
void update_collision_model()
void reset()
Remove the LB solver.
Definition lb/Solver.cpp:66
void set(Args... args)
Set the LB solver.
Matrix representation with static size.
Definition matrix.hpp:65