Loading [MathJax]/extensions/TeX/AMSmath.js
ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages Concepts
BoundaryHandling.hpp
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
20#pragma once
21
23
25
26#include <blockforest/StructuredBlockForest.h>
27#include <domain_decomposition/BlockDataID.h>
28#include <domain_decomposition/IBlock.h>
29#include <field/FlagField.h>
30#include <field/FlagUID.h>
31
32#include <utils/Vector.hpp>
33
34#include <cassert>
35#include <functional>
36#include <memory>
37#include <tuple>
38#include <unordered_map>
39
40namespace walberla {
41
42/**
43 * @brief Boundary class optimized for sparse data.
44 *
45 * Instead of storing the boundary data on a vector field,
46 * store individual vectors in a map.
47 * The global cell is used as key.
48 *
49 * Requires a custom communicator:
50 * @ref walberla::field::communication::BoundaryPackInfo.
51 */
52template <typename T, typename BoundaryClass> class BoundaryHandling {
53private:
54 /** Flag for domain cells, i.e. all cells. */
55 FlagUID const Domain_flag{"domain"};
56 /** Flag for boundary cells. */
57 FlagUID const Boundary_flag{"boundary"};
58
59 /** Container for the map between cells and values. */
60 class DynamicValueCallback {
61 public:
62 DynamicValueCallback() {
63 m_value_boundary = std::make_shared<std::unordered_map<Cell, T>>();
64 }
65
66 [[nodiscard]] T operator()(
67 Cell const &local,
68 std::shared_ptr<blockforest::StructuredBlockForest> const &blocks,
69 IBlock &block) const {
71 blocks->transformBlockLocalToGlobalCell(global, block, local);
72 return get_value(global);
73 }
74
75 void set_node_boundary_value(Utils::Vector3i const &node, T const &val) {
76 auto const global = Cell(node[0], node[1], node[2]);
77 (*m_value_boundary)[global] = val;
78 }
79
80 void unset_node_boundary_value(Utils::Vector3i const &node) {
81 auto const global = Cell(node[0], node[1], node[2]);
82 assert(m_value_boundary->count(global));
83 m_value_boundary->erase(global);
84 }
85
86 [[nodiscard]] auto &
87 get_node_boundary_value(Utils::Vector3i const &node) const {
88 auto const global = Cell(node[0], node[1], node[2]);
89 return get_value(global);
90 }
91
92 bool node_is_boundary(Utils::Vector3i const &node) const {
93 auto const global = Cell(node[0], node[1], node[2]);
94 return m_value_boundary->count(global) != 0;
95 }
96
97 private:
98 std::shared_ptr<std::unordered_map<Cell, T>> m_value_boundary;
99 static constexpr T default_value{};
100
101 [[nodiscard]] T const &get_value(Cell const &cell) const {
102 if (m_value_boundary->count(cell) == 0) {
103 return default_value;
104 }
105 return m_value_boundary->at(cell);
106 }
107 };
108
109 [[nodiscard]] inline auto get_flag_field_and_flag(IBlock *block) const {
110 auto const flag_field =
111 block->template uncheckedFastGetData<FlagField>(m_flag_field_id);
112 auto const boundary_flag = flag_field->getFlag(Boundary_flag);
113 return std::make_tuple(flag_field, boundary_flag);
114 }
115
116public:
117 using value_type = T;
118 using FlagField = field::FlagField<uint8_t>;
119
120 BoundaryHandling(std::shared_ptr<StructuredBlockForest> blocks,
122 : m_blocks(std::move(blocks)), m_flag_field_id(flag_field_id),
123 m_callback(DynamicValueCallback()), m_pending_changes(false) {
124 // reinitialize the flag field
125 for (auto block = m_blocks->begin(); block != m_blocks->end(); ++block) {
126 flag_reset_kernel(block->template getData<FlagField>(m_flag_field_id));
127 }
128 // instantiate the boundary sweep
129 std::function callback = m_callback;
130 m_boundary =
131 std::make_shared<BoundaryClass>(m_blocks, value_field_id, callback);
132 }
133
134 void operator()(IBlock *block) { (*m_boundary)(block); }
135
136 [[nodiscard]] bool node_is_boundary(Utils::Vector3i const &node) const {
137 return m_callback.node_is_boundary(node);
138 }
139
140 [[nodiscard]] auto &
142 return m_callback.get_node_boundary_value(node);
143 }
144
145 void set_node_value_at_boundary(Utils::Vector3i const &node, T const &v,
146 BlockAndCell const &bc) {
147 auto [flag_field, boundary_flag] = get_flag_field_and_flag(bc.block);
148 m_callback.set_node_boundary_value(node, v);
149 flag_field->addFlag(bc.cell, boundary_flag);
150 m_pending_changes = true;
151 }
152
153 void unpack_node(Utils::Vector3i const &node, T const &v) {
154 m_callback.set_node_boundary_value(node, v);
155 }
156
158 BlockAndCell const &bc) {
159 auto [flag_field, boundary_flag] = get_flag_field_and_flag(bc.block);
160 m_callback.unset_node_boundary_value(node);
161 flag_field->removeFlag(bc.cell, boundary_flag);
162 m_pending_changes = true;
163 }
164
165 /** Assign boundary conditions to boundary cells. */
167 if (m_pending_changes) {
168 m_boundary->template fillFromFlagField<FlagField>(
169 m_blocks, m_flag_field_id, Boundary_flag, Domain_flag);
170 m_pending_changes = false;
171 }
172 }
173
174 std::tuple<int, int, int> block_dims(IBlock const &block) const {
175 auto const field = block.template getData<FlagField>(m_flag_field_id);
176 return {static_cast<int>(field->xSize()), static_cast<int>(field->ySize()),
177 static_cast<int>(field->zSize())};
178 }
179
180private:
181 std::shared_ptr<StructuredBlockForest> m_blocks;
182 BlockDataID m_flag_field_id;
183 DynamicValueCallback m_callback;
184 std::shared_ptr<BoundaryClass> m_boundary;
185 bool m_pending_changes;
186
187 /** Register flags and reset all cells. */
188 void flag_reset_kernel(FlagField *flag_field) {
189 // register flags
190 if (!flag_field->flagExists(Domain_flag))
191 flag_field->registerFlag(Domain_flag);
192 if (!flag_field->flagExists(Boundary_flag))
193 flag_field->registerFlag(Boundary_flag);
194 // mark all cells as domain cells and fluid cells
195 auto domain_flag = flag_field->getFlag(Domain_flag);
196 auto boundary_flag = flag_field->getFlag(Boundary_flag);
197 for (auto it = flag_field->begin(); it != flag_field->end(); ++it) {
198 flag_field->addFlag(it.x(), it.y(), it.z(), domain_flag);
199 flag_field->removeFlag(it.x(), it.y(), it.z(), boundary_flag);
200 }
201 }
202};
203
204} // namespace walberla
Vector implementation and trait types for boost qvm interoperability.
Definition Cell.hpp:96
Boundary class optimized for sparse data.
bool node_is_boundary(Utils::Vector3i const &node) const
auto & get_node_value_at_boundary(Utils::Vector3i const &node) const
void boundary_update()
Assign boundary conditions to boundary cells.
void unpack_node(Utils::Vector3i const &node, T const &v)
field::FlagField< uint8_t > FlagField
void set_node_value_at_boundary(Utils::Vector3i const &node, T const &v, BlockAndCell const &bc)
void remove_node_from_boundary(Utils::Vector3i const &node, BlockAndCell const &bc)
BoundaryHandling(std::shared_ptr< StructuredBlockForest > blocks, BlockDataID value_field_id, BlockDataID flag_field_id)
std::tuple< int, int, int > block_dims(IBlock const &block) const
cudaStream_t stream[1]
CUDA streams for parallel computing on CPU and GPU.
static double * block(double *p, std::size_t index, std::size_t size)
Definition elc.cpp:172
\file PackInfoPdfDoublePrecision.cpp \author pystencils