ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
EKSpecies.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 "EKSpecies.hpp"
26
28
29#include <boost/mpi.hpp>
30#include <boost/mpi/collectives/all_reduce.hpp>
31#include <boost/mpi/collectives/broadcast.hpp>
32
33#include <algorithm>
34#include <cassert>
35#include <functional>
36#include <memory>
37#include <optional>
38#include <sstream>
39#include <stdexcept>
40#include <string>
41#include <vector>
42
44
45std::unordered_map<std::string, int> const EKVTKHandle::obs_map = {
46 {"density", static_cast<int>(EKOutputVTK::density)},
47};
48
50 VariantMap const &parameters) {
51 if (method == "update_flux_boundary_from_shape") {
53 std::transform(values.begin(), values.end(), values.begin(),
54 [this](double v) { return v * m_conv_flux; });
55
56 m_instance->update_flux_boundary_from_shape(
57 get_value<std::vector<int>>(parameters, "raster"), values);
58 return {};
59 }
60 if (method == "update_density_boundary_from_shape") {
62 std::transform(values.begin(), values.end(), values.begin(),
63 [this](double v) { return v * m_conv_density; });
64 m_instance->update_density_boundary_from_shape(
65 get_value<std::vector<int>>(parameters, "raster"), values);
66 return {};
67 }
68 if (method == "clear_flux_boundaries") {
69 m_instance->clear_flux_boundaries();
70 return {};
71 }
72 if (method == "clear_density_boundaries") {
73 m_instance->clear_density_boundaries();
74 return {};
75 }
76 if (method == "save_checkpoint") {
77 auto const path = get_value<std::string>(parameters, "path");
78 auto const mode = get_value<int>(parameters, "mode");
79 save_checkpoint(path, mode);
80 return {};
81 }
82 if (method == "load_checkpoint") {
83 auto const path = get_value<std::string>(parameters, "path");
84 auto const mode = get_value<int>(parameters, "mode");
85 load_checkpoint(path, mode);
86 return {};
87 }
89}
90
92 auto const diffusion = get_value<double>(params, "diffusion");
93 auto const ext_efield = get_value<Utils::Vector3d>(params, "ext_efield");
94 auto const density = get_value<double>(params, "density");
95 auto const kT = get_value<double>(params, "kT");
98 auto const ek_density = density * m_conv_density;
99 auto const ek_kT = kT * m_conv_energy;
101 m_lattice->lattice(), ek_diffusion, ek_kT,
103 get_value<bool>(params, "advection"),
104 get_value<bool>(params, "friction_coupling"),
105 get_value<bool>(params, "single_precision"),
106 get_value_or<bool>(params, "thermalized", false),
107 static_cast<uint>(get_value_or<int>(params, "seed", 0)));
108}
109
114 auto const agrid = get_value<double>(m_lattice->get_parameter("agrid"));
115 auto const density = get_value<double>(params, "density");
116 auto const kT = get_value<double>(params, "kT");
117 auto const tau = m_tau = get_value<double>(params, "tau");
118 auto const seed = get_value_or<int>(params, "seed", 0);
119 context()->parallel_try_catch([&]() {
120 if (get_value<bool>(params, "thermalized") and
121 not params.contains("seed")) {
122 throw std::invalid_argument(
123 "Parameter 'seed' is required for thermalized EKSpecies");
124 }
125 if (seed < 0) {
126 throw std::domain_error("Parameter 'seed' must be >= 0");
127 }
128 if (tau <= 0.) {
129 throw std::domain_error("Parameter 'tau' must be > 0");
130 }
131 if (kT < 0.) {
132 throw std::domain_error("Parameter 'kT' must be >= 0");
133 }
134 if (density < 0.) {
135 throw std::domain_error("Parameter 'density' must be >= 0");
136 }
137 m_conv_energy = Utils::int_pow<2>(tau) / Utils::int_pow<2>(agrid);
138 m_conv_diffusion = tau / Utils::int_pow<2>(agrid);
139 m_conv_ext_efield = Utils::int_pow<2>(tau) / agrid;
140 m_conv_density = Utils::int_pow<3>(agrid);
141 m_conv_flux = tau * Utils::int_pow<2>(agrid);
144 for (auto &vtk : m_vtk_writers) {
145 vtk->attach_to_lattice(m_instance, get_latice_to_md_units_conversion());
146 }
147 });
148}
149
150void EKSpecies::load_checkpoint(std::string const &filename, int mode) {
151 auto &ek_obj = *m_instance;
152
153 auto const read_metadata = [&ek_obj](CheckpointFile &cpfile) {
154 auto const expected_grid_size = ek_obj.get_lattice().get_grid_dimensions();
158 std::stringstream message;
159 message << "grid dimensions mismatch, read [" << read_grid_size << "], "
160 << "expected [" << expected_grid_size << "].";
161 throw std::runtime_error(message.str());
162 }
163 };
164
165 auto const read_data = [&ek_obj](CheckpointFile &cpfile) {
166 auto const grid_size = ek_obj.get_lattice().get_grid_dimensions();
167 auto const i_max = grid_size[0];
168 auto const j_max = grid_size[1];
169 auto const k_max = grid_size[2];
171 for (int i = 0; i < i_max; i++) {
172 for (int j = 0; j < j_max; j++) {
173 for (int k = 0; k < k_max; k++) {
174 auto const ind = Utils::Vector3i{{i, j, k}};
175 cpfile.read(cpnode.density);
176 cpfile.read(cpnode.is_boundary_density);
177 if (cpnode.is_boundary_density) {
178 cpfile.read(cpnode.density_boundary);
179 }
180 cpfile.read(cpnode.is_boundary_flux);
181 if (cpnode.is_boundary_flux) {
182 cpfile.read(cpnode.flux_boundary);
183 }
184 ek_obj.set_node_density(ind, cpnode.density);
185 if (cpnode.is_boundary_density) {
186 ek_obj.set_node_density_boundary(ind, cpnode.density_boundary);
187 }
188 if (cpnode.is_boundary_flux) {
189 ek_obj.set_node_flux_boundary(ind, cpnode.flux_boundary);
190 }
191 }
192 }
193 }
194 };
195
196 auto const on_success = [&ek_obj]() { ek_obj.ghost_communication(); };
197
200}
201
202void EKSpecies::save_checkpoint(std::string const &filename, int mode) {
203 auto &ek_obj = *m_instance;
204
205 auto const write_metadata = [&ek_obj,
206 mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
207 Context const &context) {
208 auto const grid_size = ek_obj.get_lattice().get_grid_dimensions();
209 if (context.is_head_node()) {
210 cpfile_ptr->write(grid_size);
212 }
213 };
214
215 auto const on_failure = [](std::shared_ptr<CheckpointFile> const &,
216 Context const &context) {
217 if (context.is_head_node()) {
218 auto failure = true;
219 boost::mpi::broadcast(context.get_comm(), failure, 0);
220 }
221 };
222
223 auto const write_data = [&ek_obj,
224 mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
225 Context const &context) {
226 auto const get_node_checkpoint =
227 [&](Utils::Vector3i const &ind) -> std::optional<EKWalberlaNodeState> {
228 auto const density = ek_obj.get_node_density(ind);
229 auto const is_b_d = ek_obj.get_node_is_density_boundary(ind);
230 auto const dens_b = ek_obj.get_node_density_at_boundary(ind);
231 auto const is_b_f = ek_obj.get_node_is_flux_boundary(ind);
232 auto const flux_b = ek_obj.get_node_flux_at_boundary(ind);
234 ((*is_b_d) ? dens_b.has_value() : true) and
235 ((*is_b_f) ? flux_b.has_value() : true)) {
238 cpnode.is_boundary_density = *is_b_d;
239 if (*is_b_d) {
240 cpnode.density_boundary = *dens_b;
241 }
242 cpnode.is_boundary_flux = *is_b_f;
243 if (*is_b_f) {
244 cpnode.flux_boundary = *flux_b;
245 }
246 return cpnode;
247 }
248 return std::nullopt;
249 };
250
251 auto failure = false;
252 auto const &comm = context.get_comm();
253 auto const is_head_node = context.is_head_node();
254 auto const unit_test_mode = (mode != static_cast<int>(CptMode::ascii)) and
255 (mode != static_cast<int>(CptMode::binary));
256 auto const grid_size = ek_obj.get_lattice().get_grid_dimensions();
257 auto const i_max = grid_size[0];
258 auto const j_max = grid_size[1];
259 auto const k_max = grid_size[2];
261 for (int i = 0; i < i_max; i++) {
262 for (int j = 0; j < j_max; j++) {
263 for (int k = 0; k < k_max; k++) {
264 auto const ind = Utils::Vector3i{{i, j, k}};
265 auto const result = get_node_checkpoint(ind);
266 if (!unit_test_mode) {
267 assert(1 == boost::mpi::all_reduce(comm, static_cast<int>(!!result),
268 std::plus<>()) &&
269 "Incorrect number of return values");
270 }
271 if (is_head_node) {
272 if (result) {
273 cpnode = *result;
274 } else {
275 comm.recv(boost::mpi::any_source, 42, cpnode);
276 }
277 auto &cpfile = *cpfile_ptr;
278 cpfile.write(cpnode.density);
279 cpfile.write(cpnode.is_boundary_density);
280 if (cpnode.is_boundary_density) {
281 cpfile.write(cpnode.density_boundary);
282 }
283 cpfile.write(cpnode.is_boundary_flux);
284 if (cpnode.is_boundary_flux) {
285 cpfile.write(cpnode.flux_boundary);
286 }
287 boost::mpi::broadcast(comm, failure, 0);
288 } else {
289 if (result) {
290 comm.send(0, 42, *result);
291 }
292 boost::mpi::broadcast(comm, failure, 0);
293 if (failure) {
294 return;
295 }
296 }
297 }
298 }
299 }
300 };
301
304}
305
306} // namespace ScriptInterface::walberla
307
308#endif // WALBERLA
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
Context * context() const
Responsible context.
void make_instance(VariantMap const &params) override
Definition EKSpecies.cpp:91
Variant do_call_method(std::string const &method, VariantMap const &parameters) override
Definition EKSpecies.cpp:49
::LatticeModel::units_map get_latice_to_md_units_conversion() const override
void do_construct(VariantMap const &params) override
Variant do_call_method(std::string const &method_name, VariantMap const &params) override
This file contains the defaults for ESPResSo.
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
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
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
std::shared_ptr< EKinWalberlaBase > new_ek_walberla(std::shared_ptr< LatticeWalberla > const &lattice, double diffusion, double kT, double valency, Utils::Vector3d ext_efield, double density, bool advection, bool friction_coupling, bool single_precision, bool thermalized, unsigned int seed)
static SteepestDescentParameters params
Currently active steepest descent instance.
Checkpoint data for a EK node.