22#include <blockforest/communication/UniformBufferedScheme.h>
23#include <field/AddToStorage.h>
24#include <field/FlagField.h>
25#include <field/FlagUID.h>
26#include <field/GhostLayerField.h>
27#include <field/communication/PackInfo.h>
28#include <field/vtk/FlagFieldCellFilter.h>
29#include <field/vtk/VTKWriter.h>
30#include <stencil/D3Q27.h>
32#include "../BoundaryHandling.hpp"
33#include "../utils/boundary.hpp"
34#include "../utils/types_conversion.hpp"
55template <std::
size_t FluxCount = 13,
typename FloatType =
double>
57 using ContinuityKernel =
59 using DiffusiveFluxKernel =
61 using AdvectiveFluxKernel =
63 using FrictionCouplingKernel =
65 using DiffusiveFluxKernelElectrostatic =
73 using FluxField = GhostLayerField<FloatType, FluxCount>;
74 using FlagField = walberla::FlagField<walberla::uint8_t>;
82 return numeric_cast<FloatType>(t);
90 return std::is_same<FloatType, double>::value;
94 FloatType m_diffusion;
99 bool m_friction_coupling;
129 [[nodiscard]] std::optional<CellInterval>
133 auto const &cell_min = lower_corner;
137 if (not lower_bc or not upper_bc) {
140 assert(&(*(lower_bc->block)) == &(*(upper_bc->block)));
141 return {CellInterval(lower_bc->cell, upper_bc->cell)};
157 typename stencil::D3Q27>;
163 double density,
bool advection,
bool friction_coupling)
165 m_valency(
FloatType_c(valency)), m_ext_efield(ext_efield),
166 m_advection(advection), m_friction_coupling(friction_coupling),
170 field::fzyx,
m_lattice->get_ghost_layers());
172 field::addFlattenedShallowCopyToStorage<DensityField>(
174 "flattened density field");
176 m_lattice->get_blocks(),
"flux field", FloatType{0}, field::fzyx,
179 field::addFlattenedShallowCopyToStorage<FluxField>(
187 m_lattice->get_blocks(),
"flag field density",
192 m_lattice->get_blocks(),
"flag field flux",
197 std::make_shared<FullCommunicator>(
m_lattice->get_blocks());
199 std::make_shared<field::communication::PackInfo<DensityField>>(
210 [[nodiscard]]
double get_kT() const noexcept
override {
return m_kT; }
218 return m_friction_coupling;
233 m_friction_coupling = friction_coupling;
236 m_ext_efield = field;
242 void kernel_boundary_density() {
244 (*m_boundary_density)(&
block);
248 void kernel_boundary_flux() {
250 (*m_boundary_flux)(&
block);
254 void kernel_continuity() {
256 (*m_continuity).run(&
block);
260 void kernel_diffusion() {
270 void kernel_advection(
const std::size_t &velocity_id) {
273 BlockDataID(velocity_id));
279 void kernel_friction_coupling(
const std::size_t &force_id) {
280 auto kernel = FrictionCouplingKernel(
288 void kernel_diffusion_electrostatic(
const std::size_t &potential_id) {
290 auto kernel = DiffusiveFluxKernelElectrostatic(
301 void kernel_migration() {}
303 void updated_boundary_fields() {
311 auto &vtk_handle = it.second;
312 if (vtk_handle->enabled) {
313 vtk::writeFiles(vtk_handle->ptr)();
314 vtk_handle->execution_count++;
320 void integrate(std::size_t potential_id, std::size_t velocity_id,
321 std::size_t force_id)
override {
323 updated_boundary_fields();
329 if (potential_id == walberla::BlockDataID{}) {
330 throw std::runtime_error(
"Walberla EK: electrostatic potential enabled "
331 "but no field accessible. potential id is " +
332 std::to_string(potential_id));
334 kernel_diffusion_electrostatic(potential_id);
340 kernel_boundary_flux();
343 if (force_id == walberla::BlockDataID{}) {
344 throw std::runtime_error(
"Walberla EK: friction coupling enabled but "
345 "no force field accessible. force_id is " +
346 std::to_string(force_id) +
347 ". Hint: LB may be inactive.");
349 kernel_friction_coupling(force_id);
353 if (velocity_id == walberla::BlockDataID{}) {
354 throw std::runtime_error(
"Walberla EK: advection enabled but no "
355 "velocity field accessible. velocity_id is " +
356 std::to_string(velocity_id) +
357 ". Hint: LB may be inactive.");
359 kernel_advection(velocity_id);
364 kernel_boundary_density();
371 static_assert(std::is_same_v<std::size_t, walberla::uint_t>);
387 [[nodiscard]] std::optional<double>
389 bool consider_ghosts =
false)
const override {
395 auto const density_field =
398 return {double_c(density_field->get(bc->cell))};
401 [[nodiscard]] std::vector<double>
404 std::vector<double> out;
405 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
407 auto const &
block = *(lattice.get_blocks()->begin());
408 auto const density_field =
410 auto const lower_cell = ci->min();
411 auto const upper_cell = ci->max();
412 auto const n_values = ci->numCells();
413 out.reserve(n_values);
414 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
415 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
416 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
417 out.emplace_back(density_field->get(
Cell{x, y, z}));
421 assert(out.size() == n_values);
428 std::vector<double>
const &
density)
override {
429 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
431 auto &
block = *(lattice.get_blocks()->begin());
435 auto const lower_cell = ci->min();
436 auto const upper_cell = ci->max();
437 assert(
density.size() == ci->numCells());
438 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
439 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
440 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
462 node, to_vector3<FloatType>(flux), *bc);
467 [[nodiscard]] std::optional<Utils::Vector3d>
469 bool consider_ghosts =
false)
const override {
499 [[nodiscard]] std::optional<double>
501 bool consider_ghosts =
false)
const override {
511 std::vector<std::optional<double>>
const &
density)
override {
512 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
514 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
515 auto const lower_cell = ci->min();
516 auto const upper_cell = ci->max();
518 assert(
density.size() == ci->numCells());
519 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
520 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
521 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
524 auto const &opt = *it;
538 [[nodiscard]] std::vector<std::optional<double>>
542 std::vector<std::optional<double>> out;
543 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
545 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
546 auto const lower_cell = ci->min();
547 auto const upper_cell = ci->max();
548 auto const n_values = ci->numCells();
549 out.reserve(n_values);
550 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
551 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
552 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
555 out.emplace_back(double_c(
558 out.emplace_back(std::nullopt);
563 assert(out.size() == n_values);
570 std::vector<std::optional<Utils::Vector3d>>
const &flux)
override {
571 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
573 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
574 auto const lower_cell = ci->min();
575 auto const upper_cell = ci->max();
576 auto it = flux.begin();
577 assert(flux.size() == ci->numCells());
578 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
579 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
580 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
583 auto const &opt = *it;
586 node, to_vector3<FloatType>(*opt), *bc);
597 [[nodiscard]] std::vector<std::optional<Utils::Vector3d>>
601 std::vector<std::optional<Utils::Vector3d>> out;
602 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
604 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
605 auto const lower_cell = ci->min();
606 auto const upper_cell = ci->max();
607 auto const n_values = ci->numCells();
608 out.reserve(n_values);
609 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
610 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
611 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
617 out.emplace_back(std::nullopt);
622 assert(out.size() == n_values);
627 [[nodiscard]] std::vector<bool>
630 std::vector<bool> out;
631 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
633 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
634 auto const lower_cell = ci->min();
635 auto const upper_cell = ci->max();
636 auto const n_values = ci->numCells();
637 out.reserve(n_values);
638 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
639 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
640 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
647 assert(out.size() == n_values);
662 [[nodiscard]] std::optional<bool>
664 bool consider_ghosts)
const override {
672 [[nodiscard]] std::optional<bool>
674 bool consider_ghosts)
const override {
682 [[nodiscard]] std::optional<bool>
684 bool consider_ghosts =
false)
const override {
694 const std::vector<int> &raster_flat,
695 const std::vector<double> &data_flat)
override {
703 const std::vector<int> &raster_flat,
704 const std::vector<double> &data_flat)
override {
725 vtk_obj.addCellExclusionFilter(fluid_filter);
729 template <
typename Field_T, u
int_t F_SIZE_ARG,
typename OutputType>
730 class VTKWriter :
public vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG> {
732 VTKWriter(ConstBlockDataID
const &block_id, std::string
const &
id,
733 FloatType unit_conversion)
734 : vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG>(id),
740 WALBERLA_ASSERT_NOT_NULLPTR(this->block_);
749 template <
typename OutputType = float,
756 OutputType
evaluate(cell_idx_t
const x, cell_idx_t
const y,
757 cell_idx_t
const z, cell_idx_t
const)
override {
758 WALBERLA_ASSERT_NOT_NULLPTR(this->
m_field);
759 auto const density = VectorTrait<typename DensityField::value_type>::get(
760 this->
m_field->get(x, y, z, 0), uint_c(0));
768 int flag_observables)
override {
770 auto const unit_conversion =
FloatType_c(units.at(
"density"));
Vector implementation and trait types for boost qvm interoperability.
__shared__ int node[MAXDEPTH *THREADS5/WARPSIZE]
Interface of a lattice-based electrokinetic model.
std::map< std::string, std::shared_ptr< VTKHandle > > m_vtk_auto
VTK writers that are executed automatically.
std::unordered_map< std::string, double > units_map
Class that runs and controls the BlockForest in waLBerla.
auto get_grid_dimensions() const
static Vector< T, N > broadcast(T const &s)
Create a vector that has all entries set to one value.
Boundary class optimized for sparse data.
OutputType evaluate(cell_idx_t const x, cell_idx_t const y, cell_idx_t const z, cell_idx_t const) override
FloatType const m_conversion
void configure() override
VTKWriter(ConstBlockDataID const &block_id, std::string const &id, FloatType unit_conversion)
ConstBlockDataID const m_block_id
Class that runs and controls the EK on waLBerla.
void update_density_boundary_from_shape(const std::vector< int > &raster_flat, const std::vector< double > &data_flat) override
double get_kT() const noexcept override
std::optional< CellInterval > get_interval(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const
blockforest::communication::UniformBufferedScheme< typename stencil::D3Q27 > FullCommunicator
std::shared_ptr< LatticeWalberla > m_lattice
Block forest.
bool get_friction_coupling() const noexcept override
void set_friction_coupling(bool friction_coupling) override
void register_vtk_field_writers(walberla::vtk::VTKOutput &vtk_obj, LatticeModel::units_map const &units, int flag_observables) override
std::vector< std::optional< double > > get_slice_density_at_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
void set_slice_density_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< std::optional< double > > const &density) override
~EKinWalberlaImpl() override=default
double get_valency() const noexcept override
void register_vtk_field_filters(walberla::vtk::VTKOutput &vtk_obj) override
walberla::FlagField< walberla::uint8_t > FlagField
void set_slice_flux_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< std::optional< Utils::Vector3d > > const &flux) override
bool remove_node_from_flux_boundary(Utils::Vector3i const &node) override
void clear_flux_boundaries() override
bool remove_node_from_density_boundary(Utils::Vector3i const &node) override
std::optional< Utils::Vector3d > get_node_flux_at_boundary(Utils::Vector3i const &node, bool consider_ghosts=false) const override
void set_diffusion(double diffusion) override
EKinWalberlaImpl(std::shared_ptr< LatticeWalberla > lattice, double diffusion, double kT, double valency, Utils::Vector3d const &ext_efield, double density, bool advection, bool friction_coupling)
std::size_t stencil_size() const noexcept override
BlockDataID m_flag_field_density_id
void update_flux_boundary_from_shape(const std::vector< int > &raster_flat, const std::vector< double > &data_flat) override
bool set_node_density_boundary(Utils::Vector3i const &node, double density) override
void set_ext_efield(Utils::Vector3d const &field) override
std::vector< double > get_slice_density(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
std::unique_ptr< BoundaryModelDensity > m_boundary_density
std::optional< bool > get_node_is_flux_boundary(Utils::Vector3i const &node, bool consider_ghosts) const override
BlockDataID m_flux_field_id
BlockDataID m_density_field_flattened_id
bool get_advection() const noexcept override
std::size_t get_density_id() const noexcept override
GhostLayerField< FloatType, FluxCount > FluxField
FloatType FloatType_c(T t)
void integrate_vtk_writers() override
double get_diffusion() const noexcept override
FlagUID const Boundary_flag
Flag for boundary cells.
void reset_density_boundary_handling()
BlockDataID m_density_field_id
bool set_node_flux_boundary(Utils::Vector3i const &node, Utils::Vector3d const &flux) override
std::optional< double > get_node_density(Utils::Vector3i const &node, bool consider_ghosts=false) const override
void reset_flux_boundary_handling()
std::shared_ptr< FullCommunicator > m_full_communication
std::unique_ptr< BoundaryModelFlux > m_boundary_flux
void set_slice_density(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< double > const &density) override
FlagUID const Domain_flag
Flag for domain cells, i.e.
bool is_double_precision() const noexcept override
void ghost_communication() override
void reallocate_density_boundary_field()
void set_kT(double kT) override
void set_valency(double valency) override
std::optional< double > get_node_density_at_boundary(Utils::Vector3i const &node, bool consider_ghosts=false) const override
void integrate(std::size_t potential_id, std::size_t velocity_id, std::size_t force_id) override
std::optional< bool > get_node_is_density_boundary(Utils::Vector3i const &node, bool consider_ghosts) const override
void reallocate_flux_boundary_field()
GhostLayerField< FloatType, 1 > DensityField
Utils::Vector3d get_ext_efield() const noexcept override
LatticeWalberla const & get_lattice() const noexcept override
bool set_node_density(Utils::Vector3i const &node, double density) override
std::unique_ptr< ContinuityKernel > m_continuity
BlockDataID m_flux_field_flattened_id
void clear_density_boundaries() override
std::optional< bool > get_node_is_boundary(Utils::Vector3i const &node, bool consider_ghosts=false) const override
void set_advection(bool advection) override
BlockDataID m_flag_field_flux_id
std::vector< std::optional< Utils::Vector3d > > get_slice_flux_at_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
std::vector< bool > get_slice_is_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
static double * block(double *p, std::size_t index, std::size_t size)
std::optional< BlockAndCell > get_block_and_cell(::LatticeWalberla const &lattice, Utils::Vector3i const &node, bool consider_ghost_layers)
std::vector< double > fill_3D_scalar_array(std::vector< double > const &vec_flat, Utils::Vector3i const &grid_size)
void set_boundary_from_grid(BoundaryModel &boundary, LatticeWalberla const &lattice, std::vector< int > const &raster_flat, std::vector< DataType > const &data_flat)
std::vector< Utils::Vector3d > fill_3D_vector_array(std::vector< double > const &vec_flat, Utils::Vector3i const &grid_size)
Utils::Vector3d to_vector3d(Vector3< float > const &v)