29#include <blockforest/Initialization.h>
30#include <blockforest/StructuredBlockForest.h>
31#include <blockforest/communication/UniformBufferedScheme.h>
32#include <domain_decomposition/BlockDataID.h>
33#include <domain_decomposition/IBlock.h>
34#include <field/AddToStorage.h>
35#include <field/GhostLayerField.h>
36#include <field/communication/PackInfo.h>
37#include <field/vtk/FlagFieldCellFilter.h>
38#include <field/vtk/VTKWriter.h>
39#include <stencil/D3Q19.h>
40#include <stencil/D3Q27.h>
41#if defined(__CUDACC__)
42#include <gpu/AddGPUFieldToStorage.h>
43#include <gpu/communication/MemcpyPackInfo.h>
44#include <gpu/communication/UniformGPUScheme.h>
47#include "../BoundaryHandling.hpp"
48#include "../BoundaryPackInfo.hpp"
49#include "../utils/boundary.hpp"
50#include "../utils/types_conversion.hpp"
54#if defined(__CUDACC__)
72#include <initializer_list>
87template <
typename FloatType, lbmpy::Arch Architecture>
91 typename detail::KernelTrait<FloatType,
94 typename detail::KernelTrait<FloatType,
102 typename detail::BoundaryHandlingTrait<
103 FloatType, Architecture>::Dynamic_UBB>;
105 std::variant<CollisionModelThermalized, CollisionModelLeesEdwards>;
116 template <
typename FT, lbmpy::Arch AT = lbmpy::Arch::CPU>
struct FieldTrait {
117 using PdfField = field::GhostLayerField<FT, Stencil::Size>;
119 template <
class Field>
120 using PackInfo = field::communication::PackInfo<Field>;
121 template <
class Stencil>
123 blockforest::communication::UniformBufferedScheme<Stencil>;
124 template <
class Stencil>
126 blockforest::communication::UniformBufferedScheme<Stencil>;
129#if defined(__CUDACC__)
133 template <
class Field>
134 using PackInfo = gpu::communication::MemcpyPackInfo<Field>;
135 template <
class Stencil>
137 template <
class Stencil>
139 blockforest::communication::UniformBufferedScheme<Stencil>;
151#if defined(__CUDACC__)
152 using GPUField = gpu::GPUField<FloatType>;
157 return numeric_cast<FloatType>(t);
161 return static_cast<std::size_t
>(Stencil::Size);
165 return std::is_same_v<FloatType, double>;
169 class CollideSweepVisitor {
174 cm.v_s_ =
static_cast<decltype(cm.v_s_)
>(
175 m_lees_edwards_callbacks->get_shear_velocity());
179 CollideSweepVisitor() =
default;
180 explicit CollideSweepVisitor(std::shared_ptr<LeesEdwardsPack> callbacks) {
181 m_lees_edwards_callbacks = std::move(callbacks);
185 std::shared_ptr<LeesEdwardsPack> m_lees_edwards_callbacks{};
187 CollideSweepVisitor m_run_collide_sweep{};
189 FloatType shear_mode_relaxation_rate()
const {
190 return FloatType{2} / (FloatType{6} *
m_viscosity + FloatType{1});
193 FloatType odd_mode_relaxation_rate(
194 FloatType shear_relaxation,
195 FloatType magic_number = FloatType{3} / FloatType{16})
const {
196 return (FloatType{4} - FloatType{2} * shear_relaxation) /
197 (FloatType{4} * magic_number * shear_relaxation + FloatType{2} -
201 void reset_boundary_handling() {
207 FloatType pressure_tensor_correction_factor()
const {
211 void pressure_tensor_correction(Matrix3<FloatType> &tensor)
const {
212 auto const revert_factor = pressure_tensor_correction_factor();
213 for (
auto const i : {1u, 2u, 3u, 5u, 6u, 7u}) {
214 tensor[i] *= revert_factor;
218 class interpolation_illegal_access :
public std::runtime_error {
220 explicit interpolation_illegal_access(std::string
const &field,
222 std::array<int, 3>
const &
node,
224 : std::runtime_error(
"Access to LB " + field +
" field failed") {
226 <<
"], weight " <<
weight <<
"\n";
233 std::string
const &reason)
234 : std::runtime_error(
"VTKOutput object '" + vtk_uid +
"' " + reason) {}
265 typename stencil::D3Q27>;
268 typename stencil::D3Q27>;
275 Architecture>::template RegularCommScheme<Stencil>;
276 template <
class Field>
293 std::shared_ptr<InterpolateAndShiftAtBoundary<_PdfField, FloatType>>
295 std::shared_ptr<InterpolateAndShiftAtBoundary<_VectorField, FloatType>>
297 std::shared_ptr<InterpolateAndShiftAtBoundary<_VectorField, FloatType>>
309 [[nodiscard]] std::optional<CellInterval>
313 auto const &cell_min = lower_corner;
317 if (not lower_bc or not upper_bc) {
320 assert(&(*(lower_bc->block)) == &(*(upper_bc->block)));
321 return {CellInterval(lower_bc->cell, upper_bc->cell)};
336 auto const &blocks =
m_lattice->get_blocks();
337 auto const n_ghost_layers =
m_lattice->get_ghost_layers();
339#ifdef ESPRESSO_BUILD_WITH_AVX_KERNELS
340#if defined(__AVX512F__)
341 constexpr uint_t alignment = 64;
342#elif defined(__AVX__)
343 constexpr uint_t alignment = 32;
344#elif defined(__SSE__)
345 constexpr uint_t alignment = 16;
347#error "Unsupported arch, check walberla src/field/allocation/FieldAllocator.h"
349 using value_type =
typename Field::value_type;
350 using Allocator = field::AllocateAligned<value_type, alignment>;
351 auto const allocator = std::make_shared<Allocator>();
352 auto const empty_set = Set<SUID>::emptySet();
353 return field::addToStorage<Field>(
354 blocks, tag, field::internal::defaultSize, FloatType{0}, field::fzyx,
355 n_ghost_layers,
false, {}, empty_set, empty_set, allocator);
357 return field::addToStorage<Field>(blocks, tag, FloatType{0}, field::fzyx,
361#if defined(__CUDACC__)
363 auto field_id = gpu::addGPUFieldToStorage<GPUField>(
364 blocks, tag, Field::F_SIZE, field::fzyx, n_ghost_layers);
365 if constexpr (std::is_same_v<Field, _VectorField>) {
367 auto field =
block->template getData<GPUField>(field_id);
370 }
else if constexpr (std::is_same_v<Field, _PdfField>) {
372 auto field =
block->template getData<GPUField>(field_id);
374 field, std::array<FloatType, Stencil::Size>{});
388 auto const &blocks =
m_lattice->get_blocks();
389 auto const n_ghost_layers =
m_lattice->get_ghost_layers();
390 if (n_ghost_layers == 0
u)
391 throw std::runtime_error(
"At least one ghost layer must be used");
405 for (
auto b = blocks->begin(); b != blocks->end(); ++b) {
411 blocks,
"flag field", n_ghost_layers);
413 reset_boundary_handling();
417 std::make_shared<PDFStreamingCommunicator>(blocks);
432 std::make_shared<BoundaryFullCommunicator>(blocks);
434 std::make_shared<field::communication::PackInfo<FlagField>>(
436 auto boundary_packinfo = std::make_shared<
444 m_reset_force = std::make_shared<ResetForce<PdfField, VectorField>>(
451 m_stream = std::make_shared<StreamSweep>(
456 void integrate_stream(std::shared_ptr<Lattice_T>
const &blocks) {
457 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
461 void integrate_collide(std::shared_ptr<Lattice_T>
const &blocks) {
463 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
464 std::visit(m_run_collide_sweep, cm_variant, std::variant<IBlock *>(&*b));
465 if (
auto *cm = std::get_if<CollisionModelThermalized>(&cm_variant)) {
470 auto has_lees_edwards_bc()
const {
471 return std::holds_alternative<CollisionModelLeesEdwards>(
475 void apply_lees_edwards_pdf_interpolation(
476 std::shared_ptr<Lattice_T>
const &blocks) {
477 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
481 void apply_lees_edwards_vel_interpolation_and_shift(
482 std::shared_ptr<Lattice_T>
const &blocks) {
483 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
487 void apply_lees_edwards_last_applied_force_interpolation(
488 std::shared_ptr<Lattice_T>
const &blocks) {
489 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
493 void integrate_reset_force(std::shared_ptr<Lattice_T>
const &blocks) {
494 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
498 void integrate_boundaries(std::shared_ptr<Lattice_T>
const &blocks) {
499 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
503 void integrate_push_scheme() {
506 integrate_reset_force(blocks);
508 integrate_collide(blocks);
511 integrate_boundaries(blocks);
513 integrate_stream(blocks);
518 void integrate_pull_scheme() {
520 integrate_reset_force(blocks);
522 integrate_boundaries(blocks);
524 integrate_stream(blocks);
526 integrate_collide(blocks);
534 auto &vtk_handle = it.second;
535 if (vtk_handle->enabled) {
536 vtk::writeFiles(vtk_handle->ptr)();
537 vtk_handle->execution_count++;
544 if (has_lees_edwards_bc()) {
545 integrate_pull_scheme();
547 integrate_push_scheme();
564 if (has_lees_edwards_bc()) {
566 apply_lees_edwards_pdf_interpolation(blocks);
567 apply_lees_edwards_vel_interpolation_and_shift(blocks);
568 apply_lees_edwards_last_applied_force_interpolation(blocks);
573 auto const omega = shear_mode_relaxation_rate();
574 auto const omega_odd = odd_mode_relaxation_rate(omega);
578 uint32_t{0
u}, uint32_t{0
u},
m_kT, omega, omega, omega_odd, omega, seed,
580 obj.block_offset_generator =
581 [
this](IBlock *
const block, uint32_t &block_offset_0,
582 uint32_t &block_offset_1, uint32_t &block_offset_2) {
584 auto const &ci = blocks->getBlockCellBB(*
block);
585 block_offset_0 =
static_cast<uint32_t
>(ci.xMin());
586 block_offset_1 =
static_cast<uint32_t
>(ci.yMin());
587 block_offset_2 =
static_cast<uint32_t
>(ci.zMin());
593 std::unique_ptr<LeesEdwardsPack> &&lees_edwards_pack)
override {
595#if defined(__CUDACC__)
597 throw std::runtime_error(
"Lees-Edwards LB doesn't support GPU yet");
600 auto const shear_direction = lees_edwards_pack->shear_direction;
601 auto const shear_plane_normal = lees_edwards_pack->shear_plane_normal;
602 auto const shear_vel =
FloatType_c(lees_edwards_pack->get_shear_velocity());
603 auto const omega = shear_mode_relaxation_rate();
604 if (shear_plane_normal != 1u) {
605 throw std::domain_error(
606 "Lees-Edwards LB only supports shear_plane_normal=\"y\"");
609 auto const n_ghost_layers = lattice.get_ghost_layers();
610 auto const blocks = lattice.get_blocks();
612 FloatType_c(lattice.get_grid_dimensions()[shear_plane_normal]);
619 std::make_shared<InterpolateAndShiftAtBoundary<_PdfField, FloatType>>(
621 shear_direction, shear_plane_normal,
626 shear_direction, shear_plane_normal,
632 n_ghost_layers, shear_direction, shear_plane_normal,
637 unsigned int shear_plane_normal)
const override {
641 throw std::runtime_error(
642 "MD and LB Lees-Edwards boundary conditions disagree");
660 std::optional<Utils::Vector3d>
662 bool consider_ghosts =
false)
const override {
671 auto field = bc->block->template uncheckedFastGetData<VectorField>(
684 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
689 auto const vel = to_vector3<FloatType>(v);
699 std::vector<double> out;
700 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
702 auto const &
block = *(lattice.get_blocks()->begin());
705 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
706 auto const lower_cell = ci->min();
707 auto const upper_cell = ci->max();
708 out.reserve(ci->numCells());
709 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
710 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
711 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
712 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
716 for (uint_t
f = 0
u;
f < 3u; ++
f) {
717 out.emplace_back(double_c(vec[
f]));
721 for (uint_t
f = 0
u;
f < 3u; ++
f) {
722 out.emplace_back(double_c(vec[
f]));
734 std::vector<double>
const &
velocity)
override {
735 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
737 auto &
block = *(lattice.get_blocks()->begin());
743 auto const lower_cell = ci->min();
744 auto const upper_cell = ci->max();
747 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
748 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
749 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
750 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
751 auto const cell =
Cell{x, y, z};
752 Vector3<FloatType> vec;
753 for (uint_t
f = 0
u;
f < 3u; ++
f) {
765 [[nodiscard]]
bool is_gpu() const noexcept
override {
770 std::vector<Utils::Vector3d>
const &forces)
override {
771 assert(
pos.size() == forces.size());
776 for (std::size_t i = 0ul; i <
pos.size(); ++i) {
780#if defined(__CUDACC__)
783 auto const &
block = *(lattice.get_blocks()->begin());
784 auto const origin =
block.getAABB().min();
785 std::vector<FloatType> host_pos;
786 std::vector<FloatType> host_force;
787 host_pos.reserve(3ul *
pos.size());
788 host_force.reserve(3ul * forces.size());
789 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
790 for (
auto const &vec :
pos) {
792 for (std::size_t i : {0ul, 1ul, 2ul}) {
793 host_pos.emplace_back(
static_cast<FloatType
>(vec[i] - origin[i]));
796 for (
auto const &vec : forces) {
798 for (std::size_t i : {0ul, 1ul, 2ul}) {
799 host_force.emplace_back(
static_cast<FloatType
>(vec[i]));
802 auto const gl = lattice.get_ghost_layers();
803 auto field =
block.template uncheckedFastGetData<VectorField>(
810 std::vector<Utils::Vector3d>
816 std::vector<Utils::Vector3d> vel{};
817 vel.reserve(
pos.size());
818 for (
auto const &vec :
pos) {
820 assert(
res.has_value());
821 vel.emplace_back(*
res);
825#if defined(__CUDACC__)
828 auto const &
block = *(lattice.get_blocks()->begin());
829 auto const origin =
block.getAABB().min();
830 std::vector<FloatType> host_pos;
831 host_pos.reserve(3ul *
pos.size());
832 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
833 for (
auto const &vec :
pos) {
835 for (std::size_t i : {0ul, 1ul, 2ul}) {
836 host_pos.emplace_back(
static_cast<FloatType
>(vec[i] - origin[i]));
839 auto const gl = lattice.get_ghost_layers();
843 std::vector<Utils::Vector3d> vel{};
844 vel.reserve(
res.size() / 3ul);
845 for (
auto it =
res.begin(); it !=
res.end(); it += 3) {
847 static_cast<double>(*(it + 1)),
848 static_cast<double>(*(it + 2))});
856 std::optional<Utils::Vector3d>
858 bool consider_points_in_halo =
false)
const override {
859 if (!consider_points_in_halo and !
m_lattice->pos_in_local_domain(
pos))
861 if (consider_points_in_halo and !
m_lattice->pos_in_local_halo(
pos))
871 throw interpolation_illegal_access(
"velocity",
pos,
node,
weight);
876 return {std::move(v)};
879 std::optional<double>
881 bool consider_points_in_halo =
false)
const override {
882 if (!consider_points_in_halo and !
m_lattice->pos_in_local_domain(
pos))
884 if (consider_points_in_halo and !
m_lattice->pos_in_local_halo(
pos))
894 throw interpolation_illegal_access(
"density",
pos,
node,
weight);
899 return {std::move(dens)};
907 auto const force_at_node = [
this, &
force](std::array<int, 3>
const node,
912 auto const weighted_force = to_vector3<FloatType>(
weight *
force);
914 bc->block->template uncheckedFastGetData<VectorField>(
923 std::optional<Utils::Vector3d>
935 std::optional<Utils::Vector3d>
937 bool consider_ghosts =
false)
const override {
956 auto const vec = to_vector3<FloatType>(
force);
965 std::vector<double> out;
966 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
968 auto const &
block = *(lattice.get_blocks()->begin());
971 auto const lower_cell = ci->min();
972 auto const upper_cell = ci->max();
973 auto const n_values = 3u * ci->numCells();
974 out.reserve(n_values);
975 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
976 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
977 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
978 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
980 for (uint_t
f = 0
u;
f < 3u; ++
f) {
981 out.emplace_back(double_c(vec[
f]));
986 assert(out.size() == n_values);
993 std::vector<double>
const &
force)
override {
994 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
996 auto &
block = *(lattice.get_blocks()->begin());
999 auto const lower_cell = ci->min();
1000 auto const upper_cell = ci->max();
1001 auto it =
force.begin();
1002 assert(
force.size() == 3u * ci->numCells());
1003 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1004 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1005 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1006 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1007 Vector3<FloatType> vec;
1008 for (uint_t
f = 0
u;
f < 3u; ++
f) {
1020 std::optional<std::vector<double>>
1022 bool consider_ghosts =
false)
const override {
1025 return std::nullopt;
1027 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1029 std::vector<double> population(Stencil::Size);
1030 for (uint_t
f = 0
u;
f < Stencil::Size; ++
f) {
1031 population[
f] = double_c(pop[
f]);
1034 return {std::move(population)};
1038 std::vector<double>
const &population)
override {
1043 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1044 std::array<FloatType, Stencil::Size> pop;
1045 for (uint_t
f = 0
u;
f < Stencil::Size; ++
f) {
1056 std::vector<double> out;
1057 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1059 auto const &
block = *(lattice.get_blocks()->begin());
1062 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1063 if constexpr (std::is_same_v<
typename decltype(values)::value_type,
1065 out = std::move(values);
1067 out = std::vector<double>(values.begin(), values.end());
1076 std::vector<double>
const &population)
override {
1077 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1079 auto &
block = *(lattice.get_blocks()->begin());
1081 assert(population.size() ==
stencil_size() * ci->numCells());
1082 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1083 std::vector<FloatType>
const values(population.begin(), population.end());
1089 std::optional<double>
1091 bool consider_ghosts =
false)
const override {
1094 return std::nullopt;
1097 bc->block->template uncheckedFastGetData<PdfField>(
m_pdf_field_id);
1107 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1116 std::vector<double> out;
1117 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1119 auto const &
block = *(lattice.get_blocks()->begin());
1122 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1123 if constexpr (std::is_same_v<
typename decltype(values)::value_type,
1125 out = std::move(values);
1127 out = std::vector<double>(values.begin(), values.end());
1129 assert(out.size() == ci->numCells());
1136 std::vector<double>
const &
density)
override {
1137 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1139 auto &
block = *(lattice.get_blocks()->begin());
1141 assert(
density.size() == ci->numCells());
1142 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1143 std::vector<FloatType>
const values(
density.begin(),
density.end());
1148 std::optional<Utils::Vector3d>
1150 bool consider_ghosts =
false)
const override {
1153 return std::nullopt;
1165 return bc.has_value();
1171 std::vector<std::optional<Utils::Vector3d>> out;
1172 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1174 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
1175 auto const lower_cell = ci->min();
1176 auto const upper_cell = ci->max();
1177 auto const n_values = ci->numCells();
1178 out.reserve(n_values);
1179 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1180 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1181 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1187 out.emplace_back(std::nullopt);
1192 assert(out.size() == n_values);
1199 std::vector<std::optional<Utils::Vector3d>>
const &
velocity)
override {
1200 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1202 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
1203 auto const lower_cell = ci->min();
1204 auto const upper_cell = ci->max();
1207 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1208 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1209 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1212 auto const &opt = *it;
1215 node, to_vector3<FloatType>(*opt), *bc);
1226 std::optional<Utils::Vector3d>
1230 return std::nullopt;
1240 return bc.has_value();
1245 bool consider_ghosts =
false)
const override {
1248 return std::nullopt;
1256 std::vector<bool> out;
1257 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1259 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
1260 auto const lower_cell = ci->min();
1261 auto const upper_cell = ci->max();
1262 auto const n_values = ci->numCells();
1263 out.reserve(n_values);
1264 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1265 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1266 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1272 assert(out.size() == n_values);
1280 reset_boundary_handling();
1286 std::vector<double>
const &data_flat)
override {
1295 std::optional<Utils::VectorXd<9>>
1299 return std::nullopt;
1301 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1303 pressure_tensor_correction(tensor);
1311 std::vector<double> out;
1312 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1314 auto const &
block = *(lattice.get_blocks()->begin());
1316 auto const lower_cell = ci->min();
1317 auto const upper_cell = ci->max();
1318 auto const n_values = 9u * ci->numCells();
1319 out.reserve(n_values);
1320 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1321 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1322 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1323 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1324 auto const cell =
Cell{x, y, z};
1326 pressure_tensor_correction(tensor);
1327 for (
auto i = 0
u; i < 9u; ++i) {
1328 out.emplace_back(tensor[i]);
1333 assert(out.size() == n_values);
1341 Matrix3<FloatType> tensor(FloatType{0});
1344 WALBERLA_FOR_ALL_CELLS_XYZ(pdf_field, {
1350 pressure_tensor_correction(tensor);
1351 return to_vector9d(tensor) * (1. /
static_cast<double>(number_of_nodes));
1357 Vector3<FloatType> mom(FloatType{0});
1376 [[nodiscard]]
double get_kT() const noexcept
override {
1377 return numeric_cast<double>(
m_kT);
1382 if (!cm or
m_kT == 0.) {
1383 return std::nullopt;
1385 return {
static_cast<uint64_t
>(cm->time_step_)};
1390 if (!cm or
m_kT == 0.) {
1391 throw std::runtime_error(
"This LB instance is unthermalized");
1394 static_cast<uint32_t
>(std::numeric_limits<uint_t>::max()));
1395 cm->time_step_ =
static_cast<uint32_t
>(counter);
1412 throw std::runtime_error(
"VTK output not supported for GPU");
1416 vtk_obj.addCellExclusionFilter(fluid_filter);
1421 template <
typename Field_T, u
int_t F_SIZE_ARG,
typename OutputType>
1422 class VTKWriter :
public vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG> {
1424 VTKWriter(ConstBlockDataID
const &block_id, std::string
const &
id,
1425 FloatType unit_conversion)
1426 : vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG>(id),
1432 WALBERLA_ASSERT_NOT_NULLPTR(this->block_);
1441 template <
typename OutputType =
float>
1447 OutputType
evaluate(cell_idx_t
const x, cell_idx_t
const y,
1448 cell_idx_t
const z, cell_idx_t
const)
override {
1449 WALBERLA_ASSERT_NOT_NULLPTR(this->
m_field);
1456 template <
typename OutputType =
float>
1462 OutputType
evaluate(cell_idx_t
const x, cell_idx_t
const y,
1463 cell_idx_t
const z, cell_idx_t
const f)
override {
1464 WALBERLA_ASSERT_NOT_NULLPTR(this->
m_field);
1471 template <
typename OutputType =
float>
1475 std::string
const &
id, FloatType unit_conversion,
1476 FloatType off_diag_factor)
1482 OutputType
evaluate(cell_idx_t
const x, cell_idx_t
const y,
1483 cell_idx_t
const z, cell_idx_t
const f)
override {
1484 WALBERLA_ASSERT_NOT_NULLPTR(this->
m_field);
1485 auto const pressure =
1487 auto const revert_factor =
1489 return numeric_cast<OutputType>(this->
m_conversion * revert_factor *
1490 pressure[uint_c(
f)]);
1498 int flag_observables)
override {
1500 throw std::runtime_error(
"VTK output not supported for GPU");
1503 auto const unit_conversion =
FloatType_c(units.at(
"density"));
1508 auto const unit_conversion =
FloatType_c(units.at(
"velocity"));
1513 auto const unit_conversion =
FloatType_c(units.at(
"pressure"));
1516 pressure_tensor_correction_factor()));
LBWalberlaBase provides the public interface of the LB waLBerla bridge.
Vector implementation and trait types for boost qvm interoperability.
__shared__ int pos[MAXDEPTH *THREADS5/WARPSIZE]
__shared__ int node[MAXDEPTH *THREADS5/WARPSIZE]
Interface of a lattice-based fluid 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.
walberla::blockforest::StructuredBlockForest Lattice_T
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.
field::FlagField< uint8_t > FlagField
OutputType evaluate(cell_idx_t const x, cell_idx_t const y, cell_idx_t const z, cell_idx_t const) override
OutputType evaluate(cell_idx_t const x, cell_idx_t const y, cell_idx_t const z, cell_idx_t const f) override
FloatType const m_off_diag_factor
PressureTensorVTKWriter(ConstBlockDataID const &block_id, std::string const &id, FloatType unit_conversion, FloatType off_diag_factor)
void configure() override
VTKWriter(ConstBlockDataID const &block_id, std::string const &id, FloatType unit_conversion)
FloatType const m_conversion
ConstBlockDataID const m_block_id
OutputType evaluate(cell_idx_t const x, cell_idx_t const y, cell_idx_t const z, cell_idx_t const f) override
Class that runs and controls the LB on waLBerla.
void add_forces_at_pos(std::vector< Utils::Vector3d > const &pos, std::vector< Utils::Vector3d > const &forces) override
typename FieldTrait< FloatType, Architecture >::template RegularCommScheme< typename stencil::D3Q27 > RegularFullCommunicator
Full communicator.
std::vector< double > get_slice_last_applied_force(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
LatticeWalberla::Lattice_T Lattice_T
Lattice model (e.g.
std::shared_ptr< RegularFullCommunicator > m_pdf_full_communicator
std::vector< std::optional< Utils::Vector3d > > get_slice_velocity_at_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
std::optional< Utils::Vector3d > get_node_last_applied_force(Utils::Vector3i const &node, bool consider_ghosts=false) const override
stencil::D3Q19 Stencil
Stencil for collision and streaming operations.
std::optional< Utils::Vector3d > get_node_velocity_at_boundary(Utils::Vector3i const &node, bool consider_ghosts=false) const override
void ghost_communication() override
std::optional< Utils::Vector3d > get_node_velocity(Utils::Vector3i const &node, bool consider_ghosts=false) const override
void reallocate_ubb_field() override
BlockDataID m_pdf_tmp_field_id
typename detail::KernelTrait< FloatType, Architecture >::CollisionModelLeesEdwards CollisionModelLeesEdwards
std::size_t get_force_field_id() const noexcept override
BlockDataID m_pdf_field_id
std::shared_ptr< CollisionModel > m_collision_model
typename FieldTrait< FloatType, Architecture >::VectorField VectorField
void set_slice_velocity(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< double > const &velocity) override
void integrate_vtk_writers() override
bool remove_node_from_boundary(Utils::Vector3i const &node) override
std::optional< Utils::Vector3d > get_node_boundary_force(Utils::Vector3i const &node) const override
std::optional< double > get_density_at_pos(Utils::Vector3d const &pos, bool consider_points_in_halo=false) const override
void set_slice_velocity_at_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< std::optional< Utils::Vector3d > > const &velocity) override
std::vector< double > get_slice_population(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
void set_rng_state(uint64_t counter) override
std::vector< double > get_slice_velocity(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
std::shared_ptr< LatticeWalberla > m_lattice
std::shared_ptr< LeesEdwardsPack > m_lees_edwards_callbacks
LBWalberlaImpl(std::shared_ptr< LatticeWalberla > lattice, double viscosity, double density)
bool set_node_last_applied_force(Utils::Vector3i const &node, Utils::Vector3d const &force) override
std::shared_ptr< InterpolateAndShiftAtBoundary< _VectorField, FloatType > > m_lees_edwards_vel_interpol_sweep
std::vector< bool > get_slice_is_boundary(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
std::vector< double > get_slice_density(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
void update_boundary_from_shape(std::vector< int > const &raster_flat, std::vector< double > const &data_flat) override
std::shared_ptr< BoundaryModel > m_boundary
std::vector< double > get_slice_pressure_tensor(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const override
bool set_node_density(Utils::Vector3i const &node, double density) override
void set_collision_model(double kT, unsigned int seed) override
std::variant< CollisionModelThermalized, CollisionModelLeesEdwards > CollisionModel
std::vector< Utils::Vector3d > get_velocities_at_pos(std::vector< Utils::Vector3d > const &pos) override
std::optional< bool > get_node_is_boundary(Utils::Vector3i const &node, bool consider_ghosts=false) const override
BlockDataID m_velocity_field_id
std::shared_ptr< ResetForce< PdfField, VectorField > > m_reset_force
virtual bool is_double_precision() const noexcept override
BlockDataID m_flag_field_id
std::shared_ptr< StreamSweep > m_stream
FlagUID const Boundary_flag
Flag for boundary cells.
~LBWalberlaImpl() override=default
std::optional< uint64_t > get_rng_state() const override
std::optional< std::vector< double > > get_node_population(Utils::Vector3i const &node, bool consider_ghosts=false) const override
FloatType m_density
kinematic viscosity
double get_viscosity() const noexcept override
double get_kT() const noexcept override
Utils::Vector3d get_momentum() const override
void set_viscosity(double viscosity) override
std::shared_ptr< BoundaryFullCommunicator > m_boundary_communicator
std::size_t stencil_size() const noexcept override
auto add_to_storage(std::string const tag)
Convenience function to add a field with a custom allocator.
typename detail::KernelTrait< FloatType, Architecture >::InitialPDFsSetter InitialPDFsSetter
bool set_node_velocity_at_boundary(Utils::Vector3i const &node, Utils::Vector3d const &velocity) override
stencil::D3Q27 StencilFull
Stencil for ghost communication (includes domain corners).
void register_vtk_field_writers(walberla::vtk::VTKOutput &vtk_obj, LatticeModel::units_map const &units, int flag_observables) override
typename BoundaryModel::FlagField FlagField
std::optional< double > get_node_density(Utils::Vector3i const &node, bool consider_ghosts=false) const override
void check_lebc(unsigned int shear_direction, unsigned int shear_plane_normal) const override
void set_slice_density(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< double > const &density) override
typename FieldTrait< FloatType >::PdfField _PdfField
Utils::Vector3d get_external_force() const noexcept override
BlockDataID m_force_to_be_applied_id
FloatType FloatType_c(T t) const
std::size_t get_velocity_field_id() const noexcept override
Utils::VectorXd< 9 > get_pressure_tensor() const override
std::shared_ptr< InterpolateAndShiftAtBoundary< _PdfField, FloatType > > m_lees_edwards_pdf_interpol_sweep
BlockDataID m_vec_tmp_field_id
void set_external_force(Utils::Vector3d const &ext_force) override
std::optional< Utils::Vector3d > get_velocity_at_pos(Utils::Vector3d const &pos, bool consider_points_in_halo=false) const override
void set_slice_last_applied_force(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< double > const &force) override
void set_collision_model(std::unique_ptr< LeesEdwardsPack > &&lees_edwards_pack) override
typename FieldTrait< FloatType, Architecture >::template BoundaryCommScheme< typename stencil::D3Q27 > BoundaryFullCommunicator
typename FieldTrait< FloatType, Architecture >::template RegularCommScheme< Stencil > PDFStreamingCommunicator
Regular communicator.
LatticeWalberla const & get_lattice() const noexcept override
void set_slice_population(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner, std::vector< double > const &population) override
typename FieldTrait< FloatType, Architecture >::template PackInfo< Field > PackInfo
double get_density() const noexcept override
typename detail::KernelTrait< FloatType, Architecture >::CollisionModelThermalized CollisionModelThermalized
typename FieldTrait< FloatType, Architecture >::PdfField PdfField
void register_vtk_field_filters(walberla::vtk::VTKOutput &vtk_obj) override
std::optional< Utils::Vector3d > get_node_force_to_be_applied(Utils::Vector3i const &node) const override
typename FieldTrait< FloatType >::VectorField _VectorField
bool set_node_velocity(Utils::Vector3i const &node, Utils::Vector3d const &v) override
void clear_boundaries() override
std::optional< CellInterval > get_interval(Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner) const
bool set_node_population(Utils::Vector3i const &node, std::vector< double > const &population) override
bool is_gpu() const noexcept override
std::shared_ptr< PDFStreamingCommunicator > m_pdf_streaming_communicator
BlockDataID m_last_applied_force_field_id
void integrate() override
std::optional< Utils::VectorXd< 9 > > get_node_pressure_tensor(Utils::Vector3i const &node) const override
std::shared_ptr< InterpolateAndShiftAtBoundary< _VectorField, FloatType > > m_lees_edwards_last_applied_force_interpol_sweep
typename detail::KernelTrait< FloatType, Architecture >::StreamSweep StreamSweep
void ghost_communication_boundary()
void ghost_communication_pdfs()
bool add_force_at_pos(Utils::Vector3d const &pos, Utils::Vector3d const &force) override
void setup_boundary_handle(std::shared_ptr< LatticeWalberla > lattice, std::shared_ptr< Boundary_T > boundary)
static double weight(int type, double r_cut, double k, double r)
static double * block(double *p, std::size_t index, std::size_t size)
T product(Vector< T, N > const &v)
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, double const rho_in, Cell const &cell)
double get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
std::vector< double > get(gpu::GPUField< double > const *vec_field, std::vector< double > const &pos, uint gl)
void set(gpu::GPUField< double > const *vec_field, std::vector< double > const &pos, std::vector< double > const &forces, uint gl)
Vector3< double > reduce(GhostLayerField< double, uint_t{19u}> const *pdf_field, GhostLayerField< double, uint_t{3u}> const *force_field)
std::array< double, 19u > get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, std::array< double, 19u > const &pop, Cell const &cell)
void initialize(GhostLayerField< double, uint_t{19u}> *pdf_field, std::array< double, 19u > const &pop)
Matrix3< double > get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
void initialize(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec)
void add(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec, Cell const &cell)
void set(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec, Cell const &cell)
Vector3< double > get(GhostLayerField< double, uint_t{3u}> const *vec_field, Cell const &cell)
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, GhostLayerField< double, uint_t{3u}> const *force_field, Vector3< double > const &u, Cell const &cell)
std::optional< BlockAndCell > get_block_and_cell(::LatticeWalberla const &lattice, Utils::Vector3i const &node, bool consider_ghost_layers)
void interpolate_bspline_at_pos(Utils::Vector3d const &pos, Function const &f)
Utils::VectorXd< 9 > to_vector9d(Matrix3< double > const &m)
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)
static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs)
Velocity of the virtual site.
DEVICE_QUALIFIER constexpr iterator begin() noexcept
DEVICE_QUALIFIER constexpr size_type size() const noexcept
field::GhostLayerField< FT, uint_t{3u}> VectorField
field::communication::PackInfo< Field > PackInfo
field::GhostLayerField< FT, Stencil::Size > PdfField
blockforest::communication::UniformBufferedScheme< Stencil > RegularCommScheme
blockforest::communication::UniformBufferedScheme< Stencil > BoundaryCommScheme