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/HostFieldAllocator.h>
44#include <gpu/communication/MemcpyPackInfo.h>
45#include <gpu/communication/UniformGPUScheme.h>
48#include "../BoundaryHandling.hpp"
49#include "../BoundaryPackInfo.hpp"
50#include "../utils/boundary.hpp"
51#include "../utils/types_conversion.hpp"
55#if defined(__CUDACC__)
74#include <initializer_list>
89template <
typename FloatType, lbmpy::Arch Architecture>
93 typename detail::KernelTrait<FloatType,
96 typename detail::KernelTrait<FloatType,
104 typename detail::BoundaryHandlingTrait<
105 FloatType, Architecture>::Dynamic_UBB>;
107 std::variant<CollisionModelThermalized, CollisionModelLeesEdwards>;
118 template <
typename FT, lbmpy::Arch AT = lbmpy::Arch::CPU>
struct FieldTrait {
119 using PdfField = field::GhostLayerField<FT, Stencil::Size>;
121 template <
class Field>
122 using PackInfo = field::communication::PackInfo<Field>;
127 template <
class Stencil>
129 blockforest::communication::UniformBufferedScheme<Stencil>;
130 template <
class Stencil>
132 blockforest::communication::UniformBufferedScheme<Stencil>;
135#if defined(__CUDACC__)
139 template <
class Field>
140 using MemcpyPackInfo = gpu::communication::MemcpyPackInfo<Field>;
143 template <
typename Stencil>
144 class UniformGPUScheme
145 :
public gpu::communication::UniformGPUScheme<Stencil> {
147 explicit UniformGPUScheme(
auto const &bf)
148 : gpu::communication::UniformGPUScheme<Stencil>(
154 template <
class Field>
using PackInfo = MemcpyPackInfo<Field>;
159 template <
class Stencil>
161 template <
class Stencil>
163 blockforest::communication::UniformBufferedScheme<Stencil>;
175#if defined(__CUDACC__)
176 using GPUField = gpu::GPUField<FloatType>;
179 using VectorFieldCpu =
196 return numeric_cast<FloatType>(t);
200 return static_cast<std::size_t
>(Stencil::Size);
204 return std::is_same_v<FloatType, double>;
208 class CollideSweepVisitor {
213 cm.configure(m_storage, b);
218 cm.v_s_ =
static_cast<decltype(cm.v_s_)
>(
219 m_lees_edwards_callbacks->get_shear_velocity());
223 CollideSweepVisitor() =
default;
224 CollideSweepVisitor(std::shared_ptr<StructuredBlockStorage> storage) {
225 m_storage = std::move(storage);
227 CollideSweepVisitor(std::shared_ptr<StructuredBlockStorage> storage,
228 std::shared_ptr<LeesEdwardsPack> callbacks) {
229 m_storage = std::move(storage);
230 m_lees_edwards_callbacks = std::move(callbacks);
234 std::shared_ptr<StructuredBlockStorage> m_storage{};
235 std::shared_ptr<LeesEdwardsPack> m_lees_edwards_callbacks{};
237 CollideSweepVisitor m_run_collide_sweep{};
239 FloatType shear_mode_relaxation_rate()
const {
240 return FloatType{2} / (FloatType{6} *
m_viscosity + FloatType{1});
243 FloatType odd_mode_relaxation_rate(
244 FloatType shear_relaxation,
245 FloatType magic_number = FloatType{3} / FloatType{16})
const {
246 return (FloatType{4} - FloatType{2} * shear_relaxation) /
247 (FloatType{4} * magic_number * shear_relaxation + FloatType{2} -
251 void reset_boundary_handling() {
257 FloatType pressure_tensor_correction_factor()
const {
261 void pressure_tensor_correction(Matrix3<FloatType> &tensor)
const {
262 auto const revert_factor = pressure_tensor_correction_factor();
263 for (
auto const i : {1u, 2u, 3u, 5u, 6u, 7u}) {
264 tensor[i] *= revert_factor;
268 void pressure_tensor_correction(std::span<FloatType, 9ul> tensor)
const {
269 auto const revert_factor = pressure_tensor_correction_factor();
270 for (
auto const i : {1u, 2u, 3u, 5u, 6u, 7u}) {
271 tensor[i] *= revert_factor;
275 class interpolation_illegal_access :
public std::runtime_error {
277 explicit interpolation_illegal_access(std::string
const &field,
279 std::array<int, 3>
const &node,
281 : std::runtime_error(
"Access to LB " + field +
" field failed") {
283 <<
"], weight " <<
weight <<
"\n";
290 std::string
const &reason)
291 : std::runtime_error(
"VTKOutput object '" + vtk_uid +
"' " + reason) {}
312#if defined(__CUDACC__)
313 std::optional<BlockDataID> m_pdf_cpu_field_id;
314 std::optional<BlockDataID> m_vel_cpu_field_id;
329 typename stencil::D3Q27>;
332 typename stencil::D3Q27>;
339 Architecture>::template RegularCommScheme<Stencil>;
340 template <
class Field>
361 std::shared_ptr<InterpolateAndShiftAtBoundary<_PdfField, FloatType>>
363 std::shared_ptr<InterpolateAndShiftAtBoundary<_VectorField, FloatType>>
365 std::shared_ptr<InterpolateAndShiftAtBoundary<_VectorField, FloatType>>
377#if defined(__CUDACC__)
378 std::shared_ptr<gpu::HostFieldAllocator<FloatType>> m_host_field_allocator;
381 [[nodiscard]] std::optional<CellInterval>
385 auto const &cell_min = lower_corner;
389 if (not lower_bc or not upper_bc) {
392 assert(&(*(lower_bc->block)) == &(*(upper_bc->block)));
393 return {CellInterval(lower_bc->cell, upper_bc->cell)};
408 auto const &blocks =
m_lattice->get_blocks();
409 auto const n_ghost_layers =
m_lattice->get_ghost_layers();
411#ifdef ESPRESSO_BUILD_WITH_AVX_KERNELS
412#if defined(__AVX512F__)
413 constexpr uint_t alignment = 64u;
414#elif defined(__AVX__)
415 constexpr uint_t alignment = 32u;
416#elif defined(__SSE__)
417 constexpr uint_t alignment = 16u;
419#error "Unsupported arch, check walberla src/field/allocation/FieldAllocator.h"
421 using value_type =
typename Field::value_type;
422 using Allocator = field::AllocateAligned<value_type, alignment>;
423 auto const allocator = std::make_shared<Allocator>();
424 auto const empty_set = Set<SUID>::emptySet();
425 return field::addToStorage<Field>(
426 blocks, tag, field::internal::defaultSize, FloatType{0}, field::fzyx,
427 n_ghost_layers,
false, {}, empty_set, empty_set, allocator);
429 return field::addToStorage<Field>(blocks, tag, FloatType{0}, field::fzyx,
433#if defined(__CUDACC__)
435 auto field_id = gpu::addGPUFieldToStorage<GPUField>(
436 blocks, tag, Field::F_SIZE, field::fzyx, n_ghost_layers);
437 if constexpr (std::is_same_v<Field, _VectorField>) {
439 auto field =
block->template getData<GPUField>(field_id);
442 }
else if constexpr (std::is_same_v<Field, _PdfField>) {
444 auto field =
block->template getData<GPUField>(field_id);
446 field, std::array<FloatType, Stencil::Size>{});
455 auto const setup = [
this]<
typename PackInfoPdf,
typename PackInfoVec>() {
456 auto const &blocks =
m_lattice->get_blocks();
458 std::make_shared<PDFStreamingCommunicator>(blocks);
470 setup.template operator()<PackInfoPdf, PackInfoVec>();
480 auto const &blocks =
m_lattice->get_blocks();
481 auto const n_ghost_layers =
m_lattice->get_ghost_layers();
482 if (n_ghost_layers == 0u)
483 throw std::runtime_error(
"At least one ghost layer must be used");
492#if defined(__CUDACC__)
493 m_host_field_allocator =
494 std::make_shared<gpu::HostFieldAllocator<FloatType>>();
501 for (
auto b = blocks->begin(); b != blocks->end(); ++b) {
507 blocks,
"flag field", n_ghost_layers);
509 reset_boundary_handling();
533 std::make_shared<BoundaryFullCommunicator>(blocks);
535 std::make_shared<field::communication::PackInfo<FlagField>>(
537 auto boundary_packinfo = std::make_shared<
547 m_reset_force = std::make_shared<ResetForce<PdfField, VectorField>>(
554 m_stream = std::make_shared<StreamSweep>(
559 void integrate_stream(std::shared_ptr<Lattice_T>
const &blocks) {
560 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
564 void integrate_collide(std::shared_ptr<Lattice_T>
const &blocks) {
566 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
567 std::visit(m_run_collide_sweep, cm_variant, std::variant<IBlock *>(&*b));
568 if (
auto *cm = std::get_if<CollisionModelThermalized>(&cm_variant)) {
573 auto has_lees_edwards_bc()
const {
574 return std::holds_alternative<CollisionModelLeesEdwards>(
578 void apply_lees_edwards_pdf_interpolation(
579 std::shared_ptr<Lattice_T>
const &blocks) {
580 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
584 void apply_lees_edwards_vel_interpolation_and_shift(
585 std::shared_ptr<Lattice_T>
const &blocks) {
586 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
590 void apply_lees_edwards_last_applied_force_interpolation(
591 std::shared_ptr<Lattice_T>
const &blocks) {
592 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
596 void integrate_reset_force(std::shared_ptr<Lattice_T>
const &blocks) {
597 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
601 void integrate_boundaries(std::shared_ptr<Lattice_T>
const &blocks) {
602 for (
auto b = blocks->begin(); b != blocks->end(); ++b)
606 void integrate_push_scheme() {
609 integrate_reset_force(blocks);
611 integrate_collide(blocks);
615 integrate_boundaries(blocks);
618 integrate_stream(blocks);
627 void integrate_pull_scheme() {
631 integrate_boundaries(blocks);
634 integrate_stream(blocks);
636 integrate_collide(blocks);
638 integrate_reset_force(blocks);
650 auto &vtk_handle = it.second;
651 if (vtk_handle->enabled) {
652 vtk::writeFiles(vtk_handle->ptr)();
653 vtk_handle->execution_count++;
660 if (has_lees_edwards_bc()) {
661 integrate_pull_scheme();
663 integrate_push_scheme();
679 if (has_lees_edwards_bc()) {
681 apply_lees_edwards_pdf_interpolation(blocks);
690 if (has_lees_edwards_bc()) {
692 apply_lees_edwards_vel_interpolation_and_shift(blocks);
701 if (has_lees_edwards_bc()) {
703 apply_lees_edwards_last_applied_force_interpolation(blocks);
718 if (has_lees_edwards_bc()) {
720 apply_lees_edwards_pdf_interpolation(blocks);
721 apply_lees_edwards_vel_interpolation_and_shift(blocks);
722 apply_lees_edwards_last_applied_force_interpolation(blocks);
730 if (has_lees_edwards_bc()) {
733 apply_lees_edwards_pdf_interpolation(blocks);
734 apply_lees_edwards_vel_interpolation_and_shift(blocks);
735 apply_lees_edwards_last_applied_force_interpolation(blocks);
743 auto const omega = shear_mode_relaxation_rate();
744 auto const omega_odd = odd_mode_relaxation_rate(omega);
750 omega_odd, omega, seed, uint32_t{0u});
752 m_run_collide_sweep = CollideSweepVisitor(blocks);
757 std::unique_ptr<LeesEdwardsPack> &&lees_edwards_pack)
override {
759#if defined(__CUDACC__)
761 throw std::runtime_error(
"Lees-Edwards LB doesn't support GPU yet");
764 auto const shear_direction = lees_edwards_pack->shear_direction;
765 auto const shear_plane_normal = lees_edwards_pack->shear_plane_normal;
766 auto const shear_vel =
FloatType_c(lees_edwards_pack->get_shear_velocity());
767 auto const omega = shear_mode_relaxation_rate();
768 if (shear_plane_normal != 1u) {
769 throw std::domain_error(
770 "Lees-Edwards LB only supports shear_plane_normal=\"y\"");
773 auto const n_ghost_layers = lattice.get_ghost_layers();
774 auto const blocks = lattice.get_blocks();
776 FloatType_c(lattice.get_grid_dimensions()[shear_plane_normal]);
783 std::make_shared<InterpolateAndShiftAtBoundary<_PdfField, FloatType>>(
785 shear_direction, shear_plane_normal,
790 shear_direction, shear_plane_normal,
796 n_ghost_layers, shear_direction, shear_plane_normal,
802 unsigned int shear_plane_normal)
const override {
806 throw std::runtime_error(
807 "MD and LB Lees-Edwards boundary conditions disagree");
825 std::optional<Utils::Vector3d>
827 bool consider_ghosts =
false)
const override {
832 if (is_boundary and *is_boundary) {
840 auto field = bc->block->template uncheckedFastGetData<VectorField>(
855 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
860 auto const vel = to_vector3<FloatType>(v);
870 std::vector<double> out;
871 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
873 auto const &
block = *(lattice.get_blocks()->begin());
877 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
878 assert(values.size() == 3u * ci->numCells());
879 if constexpr (std::is_same_v<
typename decltype(values)::value_type,
881 out = std::move(values);
883 out = std::vector<double>(values.begin(), values.end());
885 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
886 auto const lower_cell = ci->min();
887 auto const upper_cell = ci->max();
888 auto it = out.begin();
889 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
890 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
891 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
894 auto const &vec =
m_boundary->get_node_value_at_boundary(node);
895 for (uint_t f = 0u; f < 3u; ++f) {
896 (*it) = double_c(vec[f]);
897 std::advance(it, 1l);
900 std::advance(it, 3l);
911 std::vector<double>
const &
velocity)
override {
914 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
916 auto &
block = *(lattice.get_blocks()->begin());
921 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
929 [[nodiscard]]
bool is_gpu() const noexcept
override {
934 std::vector<Utils::Vector3d>
const &forces)
override {
935 assert(pos.size() == forces.size());
940 for (std::size_t i = 0ul; i < pos.size(); ++i) {
944#if defined(__CUDACC__)
947 auto const &
block = *(lattice.get_blocks()->begin());
948 auto const origin =
block.getAABB().min();
949 std::vector<FloatType> host_pos;
950 std::vector<FloatType> host_force;
951 host_pos.reserve(3ul * pos.size());
952 host_force.reserve(3ul * forces.size());
953 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
954 for (
auto const &vec : pos) {
956 for (std::size_t i : {0ul, 1ul, 2ul}) {
957 host_pos.emplace_back(
static_cast<FloatType
>(vec[i] - origin[i]));
960 for (
auto const &vec : forces) {
962 for (std::size_t i : {0ul, 1ul, 2ul}) {
963 host_force.emplace_back(
static_cast<FloatType
>(vec[i]));
966 auto const gl = lattice.get_ghost_layers();
967 auto field =
block.template uncheckedFastGetData<VectorField>(
974 std::vector<Utils::Vector3d>
980 std::vector<Utils::Vector3d> vel{};
981 vel.reserve(pos.size());
982 for (
auto const &vec : pos) {
984 assert(res.has_value());
985 vel.emplace_back(*res);
989#if defined(__CUDACC__)
992 auto const &
block = *(lattice.get_blocks()->begin());
993 auto const origin =
block.getAABB().min();
994 std::vector<FloatType> host_pos;
995 host_pos.reserve(3ul * pos.size());
996 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
997 for (
auto const &vec : pos) {
999 for (std::size_t i : {0ul, 1ul, 2ul}) {
1000 host_pos.emplace_back(
static_cast<FloatType
>(vec[i] - origin[i]));
1003 auto const gl = lattice.get_ghost_layers();
1007 std::vector<Utils::Vector3d> vel{};
1008 vel.reserve(res.size() / 3ul);
1009 for (
auto it = res.begin(); it != res.end(); it += 3) {
1011 static_cast<double>(*(it + 1)),
1012 static_cast<double>(*(it + 2))});
1020 std::optional<Utils::Vector3d>
1022 bool consider_points_in_halo =
false)
const override {
1025 if (!consider_points_in_halo and !
m_lattice->pos_in_local_domain(pos))
1026 return std::nullopt;
1027 if (consider_points_in_halo and !
m_lattice->pos_in_local_halo(pos))
1028 return std::nullopt;
1031 pos, [
this, &v, &pos](std::array<int, 3>
const node,
double weight) {
1037 throw interpolation_illegal_access(
"velocity", pos, node,
weight);
1042 return {std::move(v)};
1045 std::optional<double>
1047 bool consider_points_in_halo =
false)
const override {
1049 if (!consider_points_in_halo and !
m_lattice->pos_in_local_domain(pos))
1050 return std::nullopt;
1051 if (consider_points_in_halo and !
m_lattice->pos_in_local_halo(pos))
1052 return std::nullopt;
1055 pos, [
this, &dens, &pos](std::array<int, 3>
const node,
double weight) {
1061 throw interpolation_illegal_access(
"density", pos, node,
weight);
1066 return {std::move(dens)};
1074 auto const force_at_node = [
this, &force](std::array<int, 3>
const node,
1079 auto const weighted_force = to_vector3<FloatType>(
weight * force);
1081 bc->block->template uncheckedFastGetData<VectorField>(
1090 std::optional<Utils::Vector3d>
1094 return std::nullopt;
1102 std::optional<Utils::Vector3d>
1104 bool consider_ghosts =
false)
const override {
1108 return std::nullopt;
1124 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1129 auto const vec = to_vector3<FloatType>(force);
1138 std::vector<double> out;
1139 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1141 auto const &
block = *(lattice.get_blocks()->begin());
1145 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1146 assert(values.size() == 3u * ci->numCells());
1147 if constexpr (std::is_same_v<
typename decltype(values)::value_type,
1149 out = std::move(values);
1151 out = std::vector<double>(values.begin(), values.end());
1159 std::vector<double>
const &force)
override {
1162 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1164 auto &
block = *(lattice.get_blocks()->begin());
1169 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1170 assert(force.size() == 3u * ci->numCells());
1171 std::vector<FloatType>
const values(force.begin(), force.end());
1177 std::optional<std::vector<double>>
1179 bool consider_ghosts =
false)
const override {
1183 return std::nullopt;
1185 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1187 std::vector<double> population(Stencil::Size);
1188 for (uint_t f = 0u; f < Stencil::Size; ++f) {
1189 population[f] = double_c(pop[f]);
1192 return {std::move(population)};
1196 std::vector<double>
const &population)
override {
1203 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1208 std::array<FloatType, Stencil::Size> pop;
1209 for (uint_t f = 0u; f < Stencil::Size; ++f) {
1221 std::vector<double> out;
1222 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1224 auto const &
block = *(lattice.get_blocks()->begin());
1227 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1228 assert(values.size() ==
stencil_size() * ci->numCells());
1229 if constexpr (std::is_same_v<
typename decltype(values)::value_type,
1231 out = std::move(values);
1233 out = std::vector<double>(values.begin(), values.end());
1241 std::vector<double>
const &population)
override {
1242 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1244 auto &
block = *(lattice.get_blocks()->begin());
1249 assert(population.size() ==
stencil_size() * ci->numCells());
1250 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1251 std::vector<FloatType>
const values(population.begin(), population.end());
1258 std::optional<double>
1260 bool consider_ghosts =
false)
const override {
1264 return std::nullopt;
1267 bc->block->template uncheckedFastGetData<PdfField>(
m_pdf_field_id);
1278 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1287 std::vector<double> out;
1288 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1290 auto const &
block = *(lattice.get_blocks()->begin());
1293 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1294 assert(values.size() == ci->numCells());
1295 if constexpr (std::is_same_v<
typename decltype(values)::value_type,
1297 out = std::move(values);
1299 out = std::vector<double>(values.begin(), values.end());
1307 std::vector<double>
const &
density)
override {
1309 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1311 auto &
block = *(lattice.get_blocks()->begin());
1313 assert(
density.size() == ci->numCells());
1314 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1315 std::vector<FloatType>
const values(
density.begin(),
density.end());
1320 std::optional<Utils::Vector3d>
1322 bool consider_ghosts =
false)
const override {
1325 if (!bc or !
m_boundary->node_is_boundary(node))
1326 return std::nullopt;
1338 node, to_vector3<FloatType>(
velocity), *bc);
1340 return bc.has_value();
1346 std::vector<std::optional<Utils::Vector3d>> out;
1347 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1349 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
1350 auto const lower_cell = ci->min();
1351 auto const upper_cell = ci->max();
1352 auto const n_values = ci->numCells();
1353 out.reserve(n_values);
1354 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1355 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1356 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1362 out.emplace_back(std::nullopt);
1367 assert(out.size() == n_values);
1374 std::vector<std::optional<Utils::Vector3d>>
const &
velocity)
override {
1377 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1379 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
1380 auto const lower_cell = ci->min();
1381 auto const upper_cell = ci->max();
1384 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1385 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1386 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1389 auto const &opt = *it;
1392 node, to_vector3<FloatType>(*opt), *bc);
1394 m_boundary->remove_node_from_boundary(node, *bc);
1403 std::optional<Utils::Vector3d>
1406 if (!bc or !
m_boundary->node_is_boundary(node))
1407 return std::nullopt;
1415 m_boundary->remove_node_from_boundary(node, *bc);
1417 return bc.has_value();
1422 bool consider_ghosts =
false)
const override {
1426 return std::nullopt;
1434 std::vector<bool> out;
1435 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1437 auto const local_offset = std::get<0>(lattice.get_local_grid_range());
1438 auto const lower_cell = ci->min();
1439 auto const upper_cell = ci->max();
1440 auto const n_values = ci->numCells();
1441 out.reserve(n_values);
1442 for (
auto x = lower_cell.x(); x <= upper_cell.x(); ++x) {
1443 for (
auto y = lower_cell.y(); y <= upper_cell.y(); ++y) {
1444 for (
auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
1446 out.emplace_back(
m_boundary->node_is_boundary(node));
1450 assert(out.size() == n_values);
1466 reset_boundary_handling();
1475 std::vector<double>
const &data_flat)
override {
1486 std::optional<Utils::VectorXd<9>>
1490 return std::nullopt;
1492 auto pdf_field = bc->block->template getData<PdfField>(
m_pdf_field_id);
1494 pressure_tensor_correction(tensor);
1502 std::vector<double> out;
1503 if (
auto const ci =
get_interval(lower_corner, upper_corner)) {
1505 auto const &
block = *(lattice.get_blocks()->begin());
1508 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
1509 assert(values.size() == 9u * ci->numCells());
1510 for (
auto it = values.begin(); it != values.end(); std::advance(it, 9l)) {
1511 pressure_tensor_correction(std::span<FloatType, 9ul>(it, 9ul));
1513 if constexpr (std::is_same_v<
typename decltype(values)::value_type,
1515 out = std::move(values);
1517 out = std::vector<double>(values.begin(), values.end());
1526 Matrix3<FloatType> tensor(FloatType{0});
1533 pressure_tensor_correction(tensor);
1534 return to_vector9d(tensor) * (1. /
static_cast<double>(number_of_nodes));
1540 Vector3<FloatType> mom(FloatType{0});
1559 [[nodiscard]]
double get_kT() const noexcept
override {
1560 return numeric_cast<double>(
m_kT);
1563 [[nodiscard]]
unsigned int get_seed() const noexcept
override {
1569 if (!cm or
m_kT == 0.) {
1570 return std::nullopt;
1572 return {
static_cast<uint64_t
>(cm->time_step_)};
1577 if (!cm or
m_kT == 0.) {
1578 throw std::runtime_error(
"This LB instance is unthermalized");
1581 static_cast<uint32_t
>(std::numeric_limits<uint_t>::max()));
1582 cm->time_step_ =
static_cast<uint32_t
>(counter);
1600 vtk_obj.addCellExclusionFilter(fluid_filter);
1604 template <
typename Field_T, u
int_t F_SIZE_ARG,
typename OutputType>
1605 class VTKWriter :
public vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG> {
1607 VTKWriter(ConstBlockDataID
const &block_id, std::string
const &
id,
1608 FloatType unit_conversion)
1609 : vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG>(id),
1615 WALBERLA_ASSERT_NOT_NULLPTR(this->block_);
1624 template <
typename OutputType =
float>
1629 using Base::evaluate;
1632 OutputType
evaluate(cell_idx_t
const x, cell_idx_t
const y,
1633 cell_idx_t
const z, cell_idx_t
const)
override {
1634 WALBERLA_ASSERT_NOT_NULLPTR(this->
m_field);
1641 template <
typename OutputType =
float>
1646 using Base::evaluate;
1649 OutputType
evaluate(cell_idx_t
const x, cell_idx_t
const y,
1650 cell_idx_t
const z, cell_idx_t
const f)
override {
1651 WALBERLA_ASSERT_NOT_NULLPTR(this->
m_field);
1658 template <
typename OutputType =
float>
1663 using Base::evaluate;
1666 std::string
const &
id, FloatType unit_conversion,
1667 FloatType off_diag_factor)
1668 :
Base(block_id, id, unit_conversion),
1672 OutputType
evaluate(cell_idx_t
const x, cell_idx_t
const y,
1673 cell_idx_t
const z, cell_idx_t
const f)
override {
1674 WALBERLA_ASSERT_NOT_NULLPTR(this->
m_field);
1675 auto const pressure =
1677 auto const revert_factor =
1679 return numeric_cast<OutputType>(this->
m_conversion * revert_factor *
1680 pressure[uint_c(f)]);
1688 int flag_observables)
override {
1689#if defined(__CUDACC__)
1690 auto const allocate_cpu_field_if_empty =
1691 [&]<
typename Field>(
auto const &blocks, std::string name,
1692 std::optional<BlockDataID> &cpu_field) {
1693 if (not cpu_field) {
1694 cpu_field = field::addToStorage<Field>(
1695 blocks, name, FloatType{0}, field::fzyx,
1696 m_lattice->get_ghost_layers(), m_host_field_allocator);
1701 auto const unit_conversion =
FloatType_c(units.at(
"density"));
1702#if defined(__CUDACC__)
1704 auto const &blocks =
m_lattice->get_blocks();
1705 allocate_cpu_field_if_empty.template operator()<PdfFieldCpu>(
1706 blocks,
"pdfs_cpu", m_pdf_cpu_field_id);
1707 vtk_obj.addBeforeFunction(gpu::fieldCpyFunctor<PdfFieldCpu, PdfField>(
1715 auto const unit_conversion =
FloatType_c(units.at(
"velocity"));
1716#if defined(__CUDACC__)
1718 auto const &blocks =
m_lattice->get_blocks();
1719 allocate_cpu_field_if_empty.template operator()<VectorFieldCpu>(
1720 blocks,
"vel_cpu", m_vel_cpu_field_id);
1721 vtk_obj.addBeforeFunction(
1722 gpu::fieldCpyFunctor<VectorFieldCpu, VectorField>(
1730 auto const unit_conversion =
FloatType_c(units.at(
"pressure"));
1731#if defined(__CUDACC__)
1733 auto const &blocks =
m_lattice->get_blocks();
1734 allocate_cpu_field_if_empty.template operator()<PdfFieldCpu>(
1735 blocks,
"pdfs_cpu", m_pdf_cpu_field_id);
1736 vtk_obj.addBeforeFunction(gpu::fieldCpyFunctor<PdfFieldCpu, PdfField>(
1740 vtk_obj.addCellDataWriter(
1743 pressure_tensor_correction_factor()));
LBWalberlaBase provides the public interface of the LB waLBerla bridge.
Vector implementation and trait types for boost qvm interoperability.
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 DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same 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_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
std::shared_ptr< RegularFullCommunicator > m_full_communicator
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
unsigned int get_seed() const noexcept 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
BlockDataID m_vel_tmp_field_id
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
void setup_streaming_communicator()
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
void ghost_communication_laf()
void ghost_communication_vel() 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::shared_ptr< RegularFullCommunicator > m_laf_communicator
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
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
std::shared_ptr< RegularFullCommunicator > m_vel_communicator
typename FieldTrait< FloatType, Architecture >::template BoundaryCommScheme< typename stencil::D3Q27 > BoundaryFullCommunicator
std::bitset< GhostComm::SIZE > m_pending_ghost_comm
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 ghost_communication_pdf() override
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
void ghost_communication_push_scheme()
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)
void set(GhostLayerField< double, uint_t{19u}> const *pdf_field, GhostLayerField< double, uint_t{3u}> *velocity_field, GhostLayerField< double, uint_t{3u}> *force_field, Vector3< double > const &force, 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)
auto reduce(GhostLayerField< double, uint_t{19u}> const *pdf_field, GhostLayerField< double, uint_t{3u}> const *force_field)
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, std::array< double, 19u > const &pop, Cell const &cell)
auto get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
void initialize(GhostLayerField< double, uint_t{19u}> *pdf_field, std::array< double, 19u > const &pop)
auto get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
auto reduce(GhostLayerField< double, uint_t{19u}> const *pdf_field)
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)
auto 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}> *velocity_field, GhostLayerField< double, uint_t{3u}> const *force_field, Vector3< double > const &u, Cell const &cell)
\file PackInfoPdfDoublePrecision.cpp \author pystencils
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
DEVICE_QUALIFIER constexpr iterator end() noexcept
typename detail::KernelTrait< FT, AT >::PackInfoVec PackInfoStreamingVec
typename detail::KernelTrait< FT, AT >::PackInfoPdf PackInfoStreamingPdf
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
GhostCommFlags
Ghost communication operations.
@ VEL
velocities communication
@ LAF
last applied forces communication
@ UBB
boundaries communication