ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
EKReactionImplIndexed.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2022-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
28
29#include <blockforest/StructuredBlockForest.h>
30#include <domain_decomposition/BlockDataID.h>
31#include <domain_decomposition/IBlock.h>
32#include <field/AddToStorage.h>
33#include <field/FlagField.h>
34#include <field/FlagUID.h>
35
36#include <utils/Vector.hpp>
37
38#include <cassert>
39#include <cstddef>
40#include <memory>
41#include <optional>
42#include <tuple>
43#include <vector>
44
45namespace walberla {
46
48private:
49 BlockDataID m_flagfield_id;
50 BlockDataID m_indexvector_id;
51 bool m_pending_changes;
52
53public:
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 using FlagField = field::FlagField<uint8_t>;
60 using IndexVectors = detail::ReactionKernelIndexedSelector::KernelTrait<>::
61 ReactionKernelIndexed::IndexVectors;
62 using IndexInfo = detail::ReactionKernelIndexedSelector::KernelTrait<>::
63 ReactionKernelIndexed::IndexInfo;
64
65private:
66 auto get_flag_field_and_flag(IBlock *block, BlockDataID const &flagfield_id) {
67 auto const flag_field =
68 block->template uncheckedFastGetData<FlagField>(flagfield_id);
69 auto const boundary_flag = flag_field->getFlag(Boundary_flag);
70 return std::make_tuple(flag_field, boundary_flag);
71 }
72
73public:
74 EKReactionImplIndexed(std::shared_ptr<LatticeWalberla> const &lattice,
75 reactants_type const &reactants, double coefficient)
76 : EKReactionBaseIndexed(lattice, reactants, coefficient),
77 m_pending_changes(false) {
78 m_flagfield_id = field::addFlagFieldToStorage<FlagField>(
79 get_lattice()->get_blocks(), "flag field reaction",
80 get_lattice()->get_ghost_layers());
81
82 auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const) {
83 return new IndexVectors();
84 };
85 m_indexvector_id = get_lattice()
86 ->get_blocks()
87 ->template addStructuredBlockData<IndexVectors>(
88 createIdxVector, "IndexField");
89
90 for (auto &block : *get_lattice()->get_blocks()) {
91 auto flag_field = block.template getData<FlagField>(m_flagfield_id);
92 // register flags
93 flag_field->registerFlag(Domain_flag);
94 flag_field->registerFlag(Boundary_flag);
95 // mark all cells as domain cells and fluid cells
96 auto domain_flag = flag_field->getFlag(Domain_flag);
97 auto boundary_flag = flag_field->getFlag(Boundary_flag);
98 for (auto it = flag_field->begin(); it != flag_field->end(); ++it) {
99 flag_field->addFlag(it.x(), it.y(), it.z(), domain_flag);
100 flag_field->removeFlag(it.x(), it.y(), it.z(), boundary_flag);
101 }
102 }
103 }
104 ~EKReactionImplIndexed() override = default;
105
109
110 void perform_reaction() override {
112 auto kernel = detail::ReactionKernelIndexedSelector::get_kernel(
113 get_reactants(), get_coefficient(), m_indexvector_id);
114 for (auto &block : *get_lattice()->get_blocks()) {
115 kernel(&block);
116 }
117 }
118
120 bool is_boundary) override {
121 if (auto bc = get_block_and_cell(*get_lattice(), node, true)) {
122 auto const [flag_field, boundary_flag] =
123 get_flag_field_and_flag(bc->block, m_flagfield_id);
124 if (is_boundary) {
125 flag_field->addFlag(bc->cell, boundary_flag);
126 } else {
127 flag_field->removeFlag(bc->cell, boundary_flag);
128 }
129 m_pending_changes = true;
130 }
131 }
132
133 [[nodiscard]] std::optional<bool>
134 get_node_is_boundary(Utils::Vector3i const &node) override {
135 if (auto bc = get_block_and_cell(*get_lattice(), node, true)) {
136 auto const [flag_field, boundary_flag] =
137 get_flag_field_and_flag(bc->block, m_flagfield_id);
138 return {flag_field->isFlagSet(bc->cell, boundary_flag)};
139 }
140 return std::nullopt;
141 }
142
144 if (m_pending_changes) {
145 for (auto &block : *get_lattice()->get_blocks()) {
146 fillFromFlagField(block);
147 }
148 m_pending_changes = false;
149 }
150 }
151
152private:
153 void fillFromFlagField(IBlock &block) {
154 auto *indexVectors =
155 block.uncheckedFastGetData<IndexVectors>(m_indexvector_id);
156 auto &indexVectorAll = indexVectors->indexVector(IndexVectors::ALL);
157 auto &indexVectorInner = indexVectors->indexVector(IndexVectors::INNER);
158 auto &indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER);
159
160 auto *flagField = block.getData<FlagField>(m_flagfield_id);
161
162 assert(flagField->flagExists(Boundary_flag) and
163 flagField->flagExists(Domain_flag));
164
165 auto boundaryFlag = flagField->getFlag(Boundary_flag);
166 auto domainFlag = flagField->getFlag(Domain_flag);
167
168 auto inner = flagField->xyzSize();
169 inner.expand(cell_idx_t(-1));
170
171 indexVectorAll.clear();
172 indexVectorInner.clear();
173 indexVectorOuter.clear();
174
175 auto flagWithGLayers = flagField->xyzSizeWithGhostLayer();
176 for (auto it = flagField->beginWithGhostLayerXYZ(); it != flagField->end();
177 ++it) {
178 if (!isFlagSet(it, boundaryFlag))
179 continue;
180
181 if (flagWithGLayers.contains(it.x(), it.y(), it.z()) &&
182 isFlagSet(it.neighbor(0, 0, 0, 0), domainFlag)) {
183 auto element = IndexInfo(it.x(), it.y(), it.z());
184 indexVectorAll.push_back(element);
185 if (inner.contains(it.x(), it.y(), it.z())) {
186 indexVectorInner.push_back(element);
187 } else {
188 indexVectorOuter.push_back(element);
189 }
190 }
191 }
192
193 indexVectors->syncGPU();
194 }
195};
196
197} // namespace walberla
Vector implementation and trait types for boost qvm interoperability.
std::vector< std::shared_ptr< EKReactant > > reactants_type
auto get_lattice() const noexcept
double get_coefficient() const noexcept
auto const & get_reactants() const noexcept
std::optional< bool > get_node_is_boundary(Utils::Vector3i const &node) override
detail::ReactionKernelIndexedSelector::KernelTrait<>::ReactionKernelIndexed::IndexVectors IndexVectors
EKReactionImplIndexed(std::shared_ptr< LatticeWalberla > const &lattice, reactants_type const &reactants, double coefficient)
FlagUID const Boundary_flag
Flag for boundary cells.
detail::ReactionKernelIndexedSelector::KernelTrait<>::ReactionKernelIndexed::IndexInfo IndexInfo
~EKReactionImplIndexed() override=default
field::FlagField< uint8_t > FlagField
FlagUID const Domain_flag
Flag for domain cells, i.e.
void set_node_is_boundary(Utils::Vector3i const &node, bool is_boundary) override
static double * block(double *p, std::size_t index, std::size_t size)
Definition elc.cpp:172
\file PackInfoPdfDoublePrecision.cpp \author pystencils
std::optional< BlockAndCell > get_block_and_cell(::LatticeWalberla const &lattice, Utils::Vector3i const &node, bool consider_ghost_layers)