ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
LBWalberlaImpl.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2019-2024 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
22/**
23 * @file
24 * @ref walberla::LBWalberlaImpl implements the interface of the LB
25 * waLBerla bridge using sweeps generated by lbmpy
26 * (see <tt>maintainer/walberla_kernels</tt>).
27 */
28
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>
46#endif
47
48#include "../BoundaryHandling.hpp"
49#include "../BoundaryPackInfo.hpp"
50#include "../utils/boundary.hpp"
51#include "../utils/types_conversion.hpp"
53#include "ResetForce.hpp"
54#include "lb_kernels.hpp"
55#if defined(__CUDACC__)
56#include "lb_kernels.cuh"
57#endif
58
64
65#include <utils/Vector.hpp>
68
69#include <array>
70#include <bitset>
71#include <cmath>
72#include <cstddef>
73#include <functional>
74#include <initializer_list>
75#include <limits>
76#include <memory>
77#include <optional>
78#include <stdexcept>
79#include <string>
80#include <tuple>
81#include <type_traits>
82#include <utility>
83#include <variant>
84#include <vector>
85
86namespace walberla {
87
88/** @brief Class that runs and controls the LB on waLBerla. */
89template <typename FloatType, lbmpy::Arch Architecture>
91protected:
93 typename detail::KernelTrait<FloatType,
94 Architecture>::CollisionModelLeesEdwards;
96 typename detail::KernelTrait<FloatType,
97 Architecture>::CollisionModelThermalized;
104 typename detail::BoundaryHandlingTrait<
105 FloatType, Architecture>::Dynamic_UBB>;
107 std::variant<CollisionModelThermalized, CollisionModelLeesEdwards>;
108
109public:
110 /** @brief Stencil for collision and streaming operations. */
111 using Stencil = stencil::D3Q19;
112 /** @brief Stencil for ghost communication (includes domain corners). */
113 using StencilFull = stencil::D3Q27;
114 /** @brief Lattice model (e.g. blockforest). */
116
117protected:
118 template <typename FT, lbmpy::Arch AT = lbmpy::Arch::CPU> struct FieldTrait {
119 using PdfField = field::GhostLayerField<FT, Stencil::Size>;
120 using VectorField = field::GhostLayerField<FT, uint_t{3u}>;
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>;
133 };
134
135#if defined(__CUDACC__)
136 template <typename FT> struct FieldTrait<FT, lbmpy::Arch::GPU> {
137 private:
138 static auto constexpr AT = lbmpy::Arch::GPU;
139 template <class Field>
140 using MemcpyPackInfo = gpu::communication::MemcpyPackInfo<Field>;
141
142 public:
143 template <typename Stencil>
144 class UniformGPUScheme
145 : public gpu::communication::UniformGPUScheme<Stencil> {
146 public:
147 explicit UniformGPUScheme(auto const &bf)
148 : gpu::communication::UniformGPUScheme<Stencil>(
149 bf, /* sendDirectlyFromGPU */ false,
150 /* useLocalCommunication */ false) {}
151 };
152 using PdfField = gpu::GPUField<FT>;
153 using VectorField = gpu::GPUField<FT>;
154 template <class Field> using PackInfo = MemcpyPackInfo<Field>;
159 template <class Stencil>
160 using RegularCommScheme = UniformGPUScheme<Stencil>;
161 template <class Stencil>
162 using BoundaryCommScheme =
163 blockforest::communication::UniformBufferedScheme<Stencil>;
164 };
165#endif
166
167 // "underlying" field types (`GPUField` has no f-size info at compile time)
170
171public:
175#if defined(__CUDACC__)
176 using GPUField = gpu::GPUField<FloatType>;
177 using PdfFieldCpu =
179 using VectorFieldCpu =
181#endif
182
183 struct GhostComm {
184 /** @brief Ghost communication operations. */
185 enum GhostCommFlags : unsigned {
186 PDF, ///< PDFs communication
187 VEL, ///< velocities communication
188 LAF, ///< last applied forces communication
189 UBB, ///< boundaries communication
190 SIZE
191 };
192 };
193
194public:
195 template <typename T> FloatType FloatType_c(T t) const {
196 return numeric_cast<FloatType>(t);
197 }
198
199 [[nodiscard]] std::size_t stencil_size() const noexcept override {
200 return static_cast<std::size_t>(Stencil::Size);
201 }
202
203 [[nodiscard]] virtual bool is_double_precision() const noexcept override {
204 return std::is_same_v<FloatType, double>;
205 }
206
207private:
208 class CollideSweepVisitor {
209 public:
210 using StructuredBlockStorage = LatticeWalberla::Lattice_T;
211
212 void operator()(CollisionModelThermalized &cm, IBlock *b) {
213 cm.configure(m_storage, b);
214 cm(b);
215 }
216
217 void operator()(CollisionModelLeesEdwards &cm, IBlock *b) {
218 cm.v_s_ = static_cast<decltype(cm.v_s_)>(
219 m_lees_edwards_callbacks->get_shear_velocity());
220 cm(b);
221 }
222
223 CollideSweepVisitor() = default;
224 CollideSweepVisitor(std::shared_ptr<StructuredBlockStorage> storage) {
225 m_storage = std::move(storage);
226 }
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);
231 }
232
233 private:
234 std::shared_ptr<StructuredBlockStorage> m_storage{};
235 std::shared_ptr<LeesEdwardsPack> m_lees_edwards_callbacks{};
236 };
237 CollideSweepVisitor m_run_collide_sweep{};
238
239 FloatType shear_mode_relaxation_rate() const {
240 return FloatType{2} / (FloatType{6} * m_viscosity + FloatType{1});
241 }
242
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} -
248 shear_relaxation);
249 }
250
251 void reset_boundary_handling() {
252 auto const &blocks = get_lattice().get_blocks();
253 m_boundary = std::make_shared<BoundaryModel>(blocks, m_pdf_field_id,
255 }
256
257 FloatType pressure_tensor_correction_factor() const {
258 return m_viscosity / (m_viscosity + FloatType{1} / FloatType{6});
259 }
260
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;
265 }
266 }
267
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;
272 }
273 }
274
275 class interpolation_illegal_access : public std::runtime_error {
276 public:
277 explicit interpolation_illegal_access(std::string const &field,
278 Utils::Vector3d const &pos,
279 std::array<int, 3> const &node,
280 double weight)
281 : std::runtime_error("Access to LB " + field + " field failed") {
282 std::cerr << "pos [" << pos << "], node [" << Utils::Vector3i(node)
283 << "], weight " << weight << "\n";
284 }
285 };
286
287 class vtk_runtime_error : public std::runtime_error {
288 public:
289 explicit vtk_runtime_error(std::string const &vtk_uid,
290 std::string const &reason)
291 : std::runtime_error("VTKOutput object '" + vtk_uid + "' " + reason) {}
292 };
293
294protected:
295 // Member variables
296 FloatType m_viscosity; /// kinematic viscosity
297 FloatType m_density;
298 FloatType m_kT;
299 unsigned int m_seed;
300
301 // Block data access handles
302 BlockDataID m_pdf_field_id;
304 BlockDataID m_flag_field_id;
305
308
311
312#if defined(__CUDACC__)
313 std::optional<BlockDataID> m_pdf_cpu_field_id;
314 std::optional<BlockDataID> m_vel_cpu_field_id;
315#endif
316
317 /** Flag for boundary cells. */
318 FlagUID const Boundary_flag{"boundary"};
319 bool m_has_boundaries{false};
320
321 /**
322 * @brief Full communicator.
323 * We use the D3Q27 directions to update cells along the diagonals during
324 * a full ghost communication. This is needed to properly update the corners
325 * of the ghost layer when setting cell velocities or populations.
326 */
328 typename FieldTrait<FloatType, Architecture>::template RegularCommScheme<
329 typename stencil::D3Q27>;
331 typename FieldTrait<FloatType, Architecture>::template BoundaryCommScheme<
332 typename stencil::D3Q27>;
333 /**
334 * @brief Regular communicator.
335 * We use the same directions as the stencil during integration.
336 */
338 typename FieldTrait<FloatType,
339 Architecture>::template RegularCommScheme<Stencil>;
340 template <class Field>
341 using PackInfo =
343
344 // communicators
345 std::shared_ptr<BoundaryFullCommunicator> m_boundary_communicator;
346 std::shared_ptr<RegularFullCommunicator> m_full_communicator;
347 std::shared_ptr<RegularFullCommunicator> m_pdf_communicator;
348 std::shared_ptr<RegularFullCommunicator> m_vel_communicator;
349 std::shared_ptr<RegularFullCommunicator> m_laf_communicator;
350 std::shared_ptr<PDFStreamingCommunicator> m_pdf_streaming_communicator;
351 std::bitset<GhostComm::SIZE> m_pending_ghost_comm;
352
353 // ResetForce sweep + external force handling
354 std::shared_ptr<ResetForce<PdfField, VectorField>> m_reset_force;
355
356 // Stream sweep
357 std::shared_ptr<StreamSweep> m_stream;
358
359 // Lees Edwards boundary interpolation
360 std::shared_ptr<LeesEdwardsPack> m_lees_edwards_callbacks;
361 std::shared_ptr<InterpolateAndShiftAtBoundary<_PdfField, FloatType>>
363 std::shared_ptr<InterpolateAndShiftAtBoundary<_VectorField, FloatType>>
365 std::shared_ptr<InterpolateAndShiftAtBoundary<_VectorField, FloatType>>
367
368 // Collision sweep
369 std::shared_ptr<CollisionModel> m_collision_model;
370
371 // boundaries
372 std::shared_ptr<BoundaryModel> m_boundary;
373
374 // lattice
375 std::shared_ptr<LatticeWalberla> m_lattice;
376
377#if defined(__CUDACC__)
378 std::shared_ptr<gpu::HostFieldAllocator<FloatType>> m_host_field_allocator;
379#endif
380
381 [[nodiscard]] std::optional<CellInterval>
382 get_interval(Utils::Vector3i const &lower_corner,
383 Utils::Vector3i const &upper_corner) const {
384 auto const &lattice = get_lattice();
385 auto const &cell_min = lower_corner;
386 auto const cell_max = upper_corner - Utils::Vector3i::broadcast(1);
387 auto const lower_bc = get_block_and_cell(lattice, cell_min, true);
388 auto const upper_bc = get_block_and_cell(lattice, cell_max, true);
389 if (not lower_bc or not upper_bc) {
390 return std::nullopt;
391 }
392 assert(&(*(lower_bc->block)) == &(*(upper_bc->block)));
393 return {CellInterval(lower_bc->cell, upper_bc->cell)};
394 }
395
396 /**
397 * @brief Convenience function to add a field with a custom allocator.
398 *
399 * When vectorization is off, let waLBerla decide which memory allocator
400 * to use. When vectorization is on, the aligned memory allocator is
401 * required, otherwise <tt>cpu_vectorize_info["assume_aligned"]</tt> will
402 * trigger assertions. That is because for single-precision kernels the
403 * waLBerla heuristic in <tt>src/field/allocation/FieldAllocator.h</tt>
404 * will fall back to @c StdFieldAlloc, yet @c AllocateAligned is needed
405 * for intrinsics to work.
406 */
407 template <typename Field> auto add_to_storage(std::string const tag) {
408 auto const &blocks = m_lattice->get_blocks();
409 auto const n_ghost_layers = m_lattice->get_ghost_layers();
410 if constexpr (Architecture == lbmpy::Arch::CPU) {
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;
418#else
419#error "Unsupported arch, check walberla src/field/allocation/FieldAllocator.h"
420#endif
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);
428#else // ESPRESSO_BUILD_WITH_AVX_KERNELS
429 return field::addToStorage<Field>(blocks, tag, FloatType{0}, field::fzyx,
430 n_ghost_layers);
431#endif // ESPRESSO_BUILD_WITH_AVX_KERNELS
432 }
433#if defined(__CUDACC__)
434 else {
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>) {
438 for (auto block = blocks->begin(); block != blocks->end(); ++block) {
439 auto field = block->template getData<GPUField>(field_id);
440 lbm::accessor::Vector::initialize(field, Vector3<FloatType>{0});
441 }
442 } else if constexpr (std::is_same_v<Field, _PdfField>) {
443 for (auto block = blocks->begin(); block != blocks->end(); ++block) {
444 auto field = block->template getData<GPUField>(field_id);
446 field, std::array<FloatType, Stencil::Size>{});
447 }
448 }
449 return field_id;
450 }
451#endif
452 }
453
455 auto const setup = [this]<typename PackInfoPdf, typename PackInfoVec>() {
456 auto const &blocks = m_lattice->get_blocks();
458 std::make_shared<PDFStreamingCommunicator>(blocks);
459 m_pdf_streaming_communicator->addPackInfo(
460 std::make_shared<PackInfoPdf>(m_pdf_field_id));
461 m_pdf_streaming_communicator->addPackInfo(
462 std::make_shared<PackInfoVec>(m_last_applied_force_field_id));
463 };
465 using PackInfoPdf = typename FieldTrait::PackInfoStreamingPdf;
466 using PackInfoVec = typename FieldTrait::PackInfoStreamingVec;
467 if (m_has_boundaries or (m_collision_model and has_lees_edwards_bc())) {
468 setup.template operator()<PackInfo<PdfField>, PackInfoVec>();
469 } else {
470 setup.template operator()<PackInfoPdf, PackInfoVec>();
471 }
472 }
473
474public:
475 LBWalberlaImpl(std::shared_ptr<LatticeWalberla> lattice, double viscosity,
476 double density)
478 m_kT(FloatType{0}), m_seed(0u), m_lattice(std::move(lattice)) {
479
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");
484
485 // Initialize and register fields (must use the "underlying" types)
486 m_pdf_field_id = add_to_storage<_PdfField>("pdfs");
487 m_pdf_tmp_field_id = add_to_storage<_PdfField>("pdfs_tmp");
488 m_last_applied_force_field_id = add_to_storage<_VectorField>("force last");
489 m_force_to_be_applied_id = add_to_storage<_VectorField>("force next");
490 m_velocity_field_id = add_to_storage<_VectorField>("velocity");
491 m_vel_tmp_field_id = add_to_storage<_VectorField>("velocity_tmp");
492#if defined(__CUDACC__)
493 m_host_field_allocator =
494 std::make_shared<gpu::HostFieldAllocator<FloatType>>();
495#endif
496
497 // Initialize and register pdf field
498 auto pdf_setter =
501 for (auto b = blocks->begin(); b != blocks->end(); ++b) {
502 pdf_setter(&*b);
503 }
504
505 // Initialize and register flag field (fluid/boundary)
506 m_flag_field_id = field::addFlagFieldToStorage<FlagField>(
507 blocks, "flag field", n_ghost_layers);
508 // Initialize boundary sweep
509 reset_boundary_handling();
510
511 // Set up the communication and register fields
513
514 m_full_communicator = std::make_shared<RegularFullCommunicator>(blocks);
515 m_full_communicator->addPackInfo(
516 std::make_shared<PackInfo<PdfField>>(m_pdf_field_id));
517 m_full_communicator->addPackInfo(
519 m_full_communicator->addPackInfo(
520 std::make_shared<PackInfo<VectorField>>(m_velocity_field_id));
521
522 m_pdf_communicator = std::make_shared<RegularFullCommunicator>(blocks);
523 m_vel_communicator = std::make_shared<RegularFullCommunicator>(blocks);
524 m_laf_communicator = std::make_shared<RegularFullCommunicator>(blocks);
525 m_pdf_communicator->addPackInfo(
526 std::make_shared<PackInfo<PdfField>>(m_pdf_field_id));
527 m_vel_communicator->addPackInfo(
528 std::make_shared<PackInfo<VectorField>>(m_velocity_field_id));
529 m_laf_communicator->addPackInfo(
531
533 std::make_shared<BoundaryFullCommunicator>(blocks);
534 m_boundary_communicator->addPackInfo(
535 std::make_shared<field::communication::PackInfo<FlagField>>(
537 auto boundary_packinfo = std::make_shared<
540 boundary_packinfo->setup_boundary_handle(m_lattice, m_boundary);
541 m_boundary_communicator->addPackInfo(boundary_packinfo);
542
544
545 // Instantiate the sweep responsible for force double buffering and
546 // external forces
547 m_reset_force = std::make_shared<ResetForce<PdfField, VectorField>>(
549
550 // Prepare LB sweeps
551 // Note: For now, combined collide-stream sweeps cannot be used,
552 // because the collide-push variant is not supported by lbmpy.
553 // The following functors are individual in-place collide and stream steps
554 m_stream = std::make_shared<StreamSweep>(
556 }
557
558private:
559 void integrate_stream(std::shared_ptr<Lattice_T> const &blocks) {
560 for (auto b = blocks->begin(); b != blocks->end(); ++b)
561 (*m_stream)(&*b);
562 }
563
564 void integrate_collide(std::shared_ptr<Lattice_T> const &blocks) {
565 auto &cm_variant = *m_collision_model;
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)) {
569 cm->time_step_++;
570 }
571 }
572
573 auto has_lees_edwards_bc() const {
574 return std::holds_alternative<CollisionModelLeesEdwards>(
576 }
577
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)
582 }
583
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)
588 }
589
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)
594 }
595
596 void integrate_reset_force(std::shared_ptr<Lattice_T> const &blocks) {
597 for (auto b = blocks->begin(); b != blocks->end(); ++b)
598 (*m_reset_force)(&*b);
599 }
600
601 void integrate_boundaries(std::shared_ptr<Lattice_T> const &blocks) {
602 for (auto b = blocks->begin(); b != blocks->end(); ++b)
603 (*m_boundary)(&*b);
604 }
605
606 void integrate_push_scheme() {
607 auto const &blocks = get_lattice().get_blocks();
608 // Reset force fields
609 integrate_reset_force(blocks);
610 // LB collide
611 integrate_collide(blocks);
612 m_pdf_streaming_communicator->communicate();
613 // Handle boundaries
614 if (m_has_boundaries) {
615 integrate_boundaries(blocks);
616 }
617 // LB stream
618 integrate_stream(blocks);
619 // Mark pending ghost layer updates
623 // Refresh ghost layers
625 }
626
627 void integrate_pull_scheme() {
628 auto const &blocks = get_lattice().get_blocks();
629 // Handle boundaries
630 if (m_has_boundaries) {
631 integrate_boundaries(blocks);
632 }
633 // LB stream
634 integrate_stream(blocks);
635 // LB collide
636 integrate_collide(blocks);
637 // Reset force fields
638 integrate_reset_force(blocks);
639 // Mark pending ghost layer updates
643 // Refresh ghost layers
645 }
646
647protected:
648 void integrate_vtk_writers() override {
649 for (auto const &it : m_vtk_auto) {
650 auto &vtk_handle = it.second;
651 if (vtk_handle->enabled) {
652 vtk::writeFiles(vtk_handle->ptr)();
653 vtk_handle->execution_count++;
654 }
655 }
656 }
657
658public:
659 void integrate() override {
660 if (has_lees_edwards_bc()) {
661 integrate_pull_scheme();
662 } else {
663 integrate_push_scheme();
664 }
665 // Handle VTK writers
667 }
668
669 void ghost_communication() override {
670 if (m_pending_ghost_comm.any()) {
673 }
674 }
675
676 void ghost_communication_pdf() override {
678 m_pdf_communicator->communicate();
679 if (has_lees_edwards_bc()) {
680 auto const &blocks = get_lattice().get_blocks();
681 apply_lees_edwards_pdf_interpolation(blocks);
682 }
684 }
685 }
686
687 void ghost_communication_vel() override {
689 m_vel_communicator->communicate();
690 if (has_lees_edwards_bc()) {
691 auto const &blocks = get_lattice().get_blocks();
692 apply_lees_edwards_vel_interpolation_and_shift(blocks);
693 }
695 }
696 }
697
700 m_laf_communicator->communicate();
701 if (has_lees_edwards_bc()) {
702 auto const &blocks = get_lattice().get_blocks();
703 apply_lees_edwards_last_applied_force_interpolation(blocks);
704 }
706 }
707 }
708
715
717 m_full_communicator->communicate();
718 if (has_lees_edwards_bc()) {
719 auto const &blocks = get_lattice().get_blocks();
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);
723 }
727 }
728
730 if (has_lees_edwards_bc()) {
731 m_full_communicator->communicate();
732 auto const &blocks = get_lattice().get_blocks();
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);
739 }
740 }
741
742 void set_collision_model(double kT, unsigned int seed) override {
743 auto const omega = shear_mode_relaxation_rate();
744 auto const omega_odd = odd_mode_relaxation_rate(omega);
745 auto const blocks = get_lattice().get_blocks();
746 m_kT = FloatType_c(kT);
747 m_seed = seed;
749 m_pdf_field_id, m_kT, omega, omega,
750 omega_odd, omega, seed, uint32_t{0u});
751 m_collision_model = std::make_shared<CollisionModel>(std::move(obj));
752 m_run_collide_sweep = CollideSweepVisitor(blocks);
754 }
755
757 std::unique_ptr<LeesEdwardsPack> &&lees_edwards_pack) override {
758 assert(m_kT == 0.);
759#if defined(__CUDACC__)
760 if constexpr (Architecture == lbmpy::Arch::GPU) {
761 throw std::runtime_error("Lees-Edwards LB doesn't support GPU yet");
762 }
763#endif
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\"");
771 }
772 auto const &lattice = get_lattice();
773 auto const n_ghost_layers = lattice.get_ghost_layers();
774 auto const blocks = lattice.get_blocks();
775 auto const agrid =
776 FloatType_c(lattice.get_grid_dimensions()[shear_plane_normal]);
777 auto obj = CollisionModelLeesEdwards(
778 m_last_applied_force_field_id, m_pdf_field_id, agrid, omega, shear_vel);
779 m_collision_model = std::make_shared<CollisionModel>(std::move(obj));
780 m_lees_edwards_callbacks = std::move(lees_edwards_pack);
781 m_run_collide_sweep = CollideSweepVisitor(blocks, m_lees_edwards_callbacks);
783 std::make_shared<InterpolateAndShiftAtBoundary<_PdfField, FloatType>>(
784 blocks, m_pdf_field_id, m_pdf_tmp_field_id, n_ghost_layers,
785 shear_direction, shear_plane_normal,
786 m_lees_edwards_callbacks->get_pos_offset);
787 m_lees_edwards_vel_interpol_sweep = std::make_shared<
789 blocks, m_velocity_field_id, m_vel_tmp_field_id, n_ghost_layers,
790 shear_direction, shear_plane_normal,
791 m_lees_edwards_callbacks->get_pos_offset,
792 m_lees_edwards_callbacks->get_shear_velocity);
796 n_ghost_layers, shear_direction, shear_plane_normal,
797 m_lees_edwards_callbacks->get_pos_offset);
799 }
800
801 void check_lebc(unsigned int shear_direction,
802 unsigned int shear_plane_normal) const override {
804 if (m_lees_edwards_callbacks->shear_direction != shear_direction or
805 m_lees_edwards_callbacks->shear_plane_normal != shear_plane_normal) {
806 throw std::runtime_error(
807 "MD and LB Lees-Edwards boundary conditions disagree");
808 }
809 }
810 }
811
812 void set_viscosity(double viscosity) override {
813 m_viscosity = FloatType_c(viscosity);
814 }
815
816 [[nodiscard]] double get_viscosity() const noexcept override {
817 return numeric_cast<double>(m_viscosity);
818 }
819
820 [[nodiscard]] double get_density() const noexcept override {
821 return numeric_cast<double>(m_density);
822 }
823
824 // Velocity
825 std::optional<Utils::Vector3d>
827 bool consider_ghosts = false) const override {
828 assert(not(consider_ghosts and m_pending_ghost_comm.test(GhostComm::VEL)));
829 assert(not(consider_ghosts and m_pending_ghost_comm.test(GhostComm::UBB)));
830 if (m_has_boundaries) {
831 auto const is_boundary = get_node_is_boundary(node, consider_ghosts);
832 if (is_boundary and *is_boundary) {
833 return get_node_velocity_at_boundary(node, consider_ghosts);
834 }
835 }
836 auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
837 if (!bc)
838 return std::nullopt;
839
840 auto field = bc->block->template uncheckedFastGetData<VectorField>(
842 auto const vec = lbm::accessor::Vector::get(field, bc->cell);
843 return to_vector3d(vec);
844 }
845
847 Utils::Vector3d const &v) override {
850 auto bc = get_block_and_cell(get_lattice(), node, false);
851 if (!bc)
852 return false;
853
854 // We have to set both, the pdf and the stored velocity field
855 auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
856 auto vel_field =
857 bc->block->template getData<VectorField>(m_velocity_field_id);
858 auto force_field =
859 bc->block->template getData<VectorField>(m_last_applied_force_field_id);
860 auto const vel = to_vector3<FloatType>(v);
861 lbm::accessor::Velocity::set(pdf_field, vel_field, force_field, vel,
862 bc->cell);
863
864 return true;
865 }
866
867 std::vector<double>
869 Utils::Vector3i const &upper_corner) const override {
870 std::vector<double> out;
871 if (auto const ci = get_interval(lower_corner, upper_corner)) {
872 auto const &lattice = get_lattice();
873 auto const &block = *(lattice.get_blocks()->begin());
874 auto const field =
875 block.template getData<VectorField>(m_velocity_field_id);
876 auto const values = lbm::accessor::Vector::get(field, *ci);
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,
880 double>) {
881 out = std::move(values);
882 } else {
883 out = std::vector<double>(values.begin(), values.end());
884 }
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) {
892 auto const node = local_offset + Utils::Vector3i{{x, y, z}};
893 if (m_boundary->node_is_boundary(node)) {
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);
898 }
899 } else {
900 std::advance(it, 3l);
901 }
902 }
903 }
904 }
905 }
906 return out;
907 }
908
909 void set_slice_velocity(Utils::Vector3i const &lower_corner,
910 Utils::Vector3i const &upper_corner,
911 std::vector<double> const &velocity) override {
914 if (auto const ci = get_interval(lower_corner, upper_corner)) {
915 auto const &lattice = get_lattice();
916 auto &block = *(lattice.get_blocks()->begin());
917 auto pdf_field = block.template getData<PdfField>(m_pdf_field_id);
918 auto force_field =
919 block.template getData<VectorField>(m_last_applied_force_field_id);
920 auto vel_field = block.template getData<VectorField>(m_velocity_field_id);
921 assert(++(lattice.get_blocks()->begin()) == lattice.get_blocks()->end());
922 assert(velocity.size() == 3u * ci->numCells());
923 std::vector<FloatType> const values(velocity.begin(), velocity.end());
924 lbm::accessor::Velocity::set(pdf_field, vel_field, force_field, values,
925 *ci);
926 }
927 }
928
929 [[nodiscard]] bool is_gpu() const noexcept override {
930 return Architecture == lbmpy::Arch::GPU;
931 }
932
933 void add_forces_at_pos(std::vector<Utils::Vector3d> const &pos,
934 std::vector<Utils::Vector3d> const &forces) override {
935 assert(pos.size() == forces.size());
936 if (pos.empty()) {
937 return;
938 }
939 if constexpr (Architecture == lbmpy::Arch::CPU) {
940 for (std::size_t i = 0ul; i < pos.size(); ++i) {
941 add_force_at_pos(pos[i], forces[i]);
942 }
943 }
944#if defined(__CUDACC__)
945 if constexpr (Architecture == lbmpy::Arch::GPU) {
946 auto const &lattice = get_lattice();
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) {
955#pragma unroll
956 for (std::size_t i : {0ul, 1ul, 2ul}) {
957 host_pos.emplace_back(static_cast<FloatType>(vec[i] - origin[i]));
958 }
959 }
960 for (auto const &vec : forces) {
961#pragma unroll
962 for (std::size_t i : {0ul, 1ul, 2ul}) {
963 host_force.emplace_back(static_cast<FloatType>(vec[i]));
964 }
965 }
966 auto const gl = lattice.get_ghost_layers();
967 auto field = block.template uncheckedFastGetData<VectorField>(
969 lbm::accessor::Interpolation::set(field, host_pos, host_force, gl);
970 }
971#endif
972 }
973
974 std::vector<Utils::Vector3d>
975 get_velocities_at_pos(std::vector<Utils::Vector3d> const &pos) override {
976 if (pos.empty()) {
977 return {};
978 }
979 if constexpr (Architecture == lbmpy::Arch::CPU) {
980 std::vector<Utils::Vector3d> vel{};
981 vel.reserve(pos.size());
982 for (auto const &vec : pos) {
983 auto res = get_velocity_at_pos(vec, true);
984 assert(res.has_value());
985 vel.emplace_back(*res);
986 }
987 return vel;
988 }
989#if defined(__CUDACC__)
990 if constexpr (Architecture == lbmpy::Arch::GPU) {
991 auto const &lattice = get_lattice();
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) {
998#pragma unroll
999 for (std::size_t i : {0ul, 1ul, 2ul}) {
1000 host_pos.emplace_back(static_cast<FloatType>(vec[i] - origin[i]));
1001 }
1002 }
1003 auto const gl = lattice.get_ghost_layers();
1004 auto field =
1005 block.template uncheckedFastGetData<VectorField>(m_velocity_field_id);
1006 auto const res = lbm::accessor::Interpolation::get(field, host_pos, gl);
1007 std::vector<Utils::Vector3d> vel{};
1008 vel.reserve(res.size() / 3ul);
1009 for (auto it = res.begin(); it != res.end(); it += 3) {
1010 vel.emplace_back(Utils::Vector3d{static_cast<double>(*(it + 0)),
1011 static_cast<double>(*(it + 1)),
1012 static_cast<double>(*(it + 2))});
1013 }
1014 return vel;
1015 }
1016#endif
1017 return {};
1018 }
1019
1020 std::optional<Utils::Vector3d>
1022 bool consider_points_in_halo = false) const override {
1023 assert(not m_pending_ghost_comm.test(GhostComm::VEL));
1024 assert(not m_pending_ghost_comm.test(GhostComm::UBB));
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;
1029 Utils::Vector3d v{0., 0., 0.};
1031 pos, [this, &v, &pos](std::array<int, 3> const node, double weight) {
1032 // Nodes with zero weight might not be accessible, because they can be
1033 // outside ghost layers
1034 if (weight != 0.) {
1035 auto const res = get_node_velocity(Utils::Vector3i(node), true);
1036 if (!res) {
1037 throw interpolation_illegal_access("velocity", pos, node, weight);
1038 }
1039 v += *res * weight;
1040 }
1041 });
1042 return {std::move(v)};
1043 }
1044
1045 std::optional<double>
1047 bool consider_points_in_halo = false) const override {
1048 assert(not m_pending_ghost_comm.test(GhostComm::PDF));
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;
1053 double dens = 0.;
1055 pos, [this, &dens, &pos](std::array<int, 3> const node, double weight) {
1056 // Nodes with zero weight might not be accessible, because they can be
1057 // outside ghost layers
1058 if (weight != 0.) {
1059 auto const res = get_node_density(Utils::Vector3i(node), true);
1060 if (!res) {
1061 throw interpolation_illegal_access("density", pos, node, weight);
1062 }
1063 dens += *res * weight;
1064 }
1065 });
1066 return {std::move(dens)};
1067 }
1068
1069 // Local force
1071 Utils::Vector3d const &force) override {
1072 if (!m_lattice->pos_in_local_halo(pos))
1073 return false;
1074 auto const force_at_node = [this, &force](std::array<int, 3> const node,
1075 double weight) {
1076 auto const bc =
1078 if (bc) {
1079 auto const weighted_force = to_vector3<FloatType>(weight * force);
1080 auto force_field =
1081 bc->block->template uncheckedFastGetData<VectorField>(
1083 lbm::accessor::Vector::add(force_field, weighted_force, bc->cell);
1084 }
1085 };
1086 interpolate_bspline_at_pos(pos, force_at_node);
1087 return true;
1088 }
1089
1090 std::optional<Utils::Vector3d>
1091 get_node_force_to_be_applied(Utils::Vector3i const &node) const override {
1092 auto const bc = get_block_and_cell(get_lattice(), node, true);
1093 if (!bc)
1094 return std::nullopt;
1095
1096 auto field =
1097 bc->block->template getData<VectorField>(m_force_to_be_applied_id);
1098 auto const vec = lbm::accessor::Vector::get(field, bc->cell);
1099 return to_vector3d(vec);
1100 }
1101
1102 std::optional<Utils::Vector3d>
1104 bool consider_ghosts = false) const override {
1105 assert(not(consider_ghosts and m_pending_ghost_comm.test(GhostComm::LAF)));
1106 auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
1107 if (!bc)
1108 return std::nullopt;
1109
1110 auto const field =
1111 bc->block->template getData<VectorField>(m_last_applied_force_field_id);
1112 auto const vec = lbm::accessor::Vector::get(field, bc->cell);
1113 return to_vector3d(vec);
1114 }
1115
1117 Utils::Vector3d const &force) override {
1120 auto bc = get_block_and_cell(get_lattice(), node, false);
1121 if (!bc)
1122 return false;
1123
1124 auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
1125 auto force_field =
1126 bc->block->template getData<VectorField>(m_last_applied_force_field_id);
1127 auto vel_field =
1128 bc->block->template getData<VectorField>(m_velocity_field_id);
1129 auto const vec = to_vector3<FloatType>(force);
1130 lbm::accessor::Force::set(pdf_field, vel_field, force_field, vec, bc->cell);
1131
1132 return true;
1133 }
1134
1135 std::vector<double> get_slice_last_applied_force(
1136 Utils::Vector3i const &lower_corner,
1137 Utils::Vector3i const &upper_corner) const override {
1138 std::vector<double> out;
1139 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1140 auto const &lattice = get_lattice();
1141 auto const &block = *(lattice.get_blocks()->begin());
1142 auto const field =
1143 block.template getData<VectorField>(m_last_applied_force_field_id);
1144 auto const values = lbm::accessor::Vector::get(field, *ci);
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,
1148 double>) {
1149 out = std::move(values);
1150 } else {
1151 out = std::vector<double>(values.begin(), values.end());
1152 }
1153 }
1154 return out;
1155 }
1156
1158 Utils::Vector3i const &upper_corner,
1159 std::vector<double> const &force) override {
1162 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1163 auto const &lattice = get_lattice();
1164 auto &block = *(lattice.get_blocks()->begin());
1165 auto pdf_field = block.template getData<PdfField>(m_pdf_field_id);
1166 auto force_field =
1167 block.template getData<VectorField>(m_last_applied_force_field_id);
1168 auto vel_field = block.template getData<VectorField>(m_velocity_field_id);
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());
1172 lbm::accessor::Force::set(pdf_field, vel_field, force_field, values, *ci);
1173 }
1174 }
1175
1176 // Population
1177 std::optional<std::vector<double>>
1179 bool consider_ghosts = false) const override {
1180 assert(not(consider_ghosts and m_pending_ghost_comm.test(GhostComm::PDF)));
1181 auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
1182 if (!bc)
1183 return std::nullopt;
1184
1185 auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
1186 auto const pop = lbm::accessor::Population::get(pdf_field, bc->cell);
1187 std::vector<double> population(Stencil::Size);
1188 for (uint_t f = 0u; f < Stencil::Size; ++f) {
1189 population[f] = double_c(pop[f]);
1190 }
1191
1192 return {std::move(population)};
1193 }
1194
1196 std::vector<double> const &population) override {
1199 auto bc = get_block_and_cell(get_lattice(), node, false);
1200 if (!bc)
1201 return false;
1202
1203 auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
1204 auto force_field =
1205 bc->block->template getData<VectorField>(m_last_applied_force_field_id);
1206 auto vel_field =
1207 bc->block->template getData<VectorField>(m_velocity_field_id);
1208 std::array<FloatType, Stencil::Size> pop;
1209 for (uint_t f = 0u; f < Stencil::Size; ++f) {
1210 pop[f] = FloatType_c(population[f]);
1211 }
1212 lbm::accessor::Population::set(pdf_field, vel_field, force_field, pop,
1213 bc->cell);
1214
1215 return true;
1216 }
1217
1218 std::vector<double>
1220 Utils::Vector3i const &upper_corner) const override {
1221 std::vector<double> out;
1222 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1223 auto const &lattice = get_lattice();
1224 auto const &block = *(lattice.get_blocks()->begin());
1225 auto const pdf_field = block.template getData<PdfField>(m_pdf_field_id);
1226 auto const values = lbm::accessor::Population::get(pdf_field, *ci);
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,
1230 double>) {
1231 out = std::move(values);
1232 } else {
1233 out = std::vector<double>(values.begin(), values.end());
1234 }
1235 }
1236 return out;
1237 }
1238
1239 void set_slice_population(Utils::Vector3i const &lower_corner,
1240 Utils::Vector3i const &upper_corner,
1241 std::vector<double> const &population) override {
1242 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1243 auto const &lattice = get_lattice();
1244 auto &block = *(lattice.get_blocks()->begin());
1245 auto pdf_field = block.template getData<PdfField>(m_pdf_field_id);
1246 auto force_field =
1247 block.template getData<VectorField>(m_last_applied_force_field_id);
1248 auto vel_field = block.template getData<VectorField>(m_velocity_field_id);
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());
1252 lbm::accessor::Population::set(pdf_field, vel_field, force_field, values,
1253 *ci);
1254 }
1255 }
1256
1257 // Density
1258 std::optional<double>
1260 bool consider_ghosts = false) const override {
1261 assert(not(consider_ghosts and m_pending_ghost_comm.test(GhostComm::PDF)));
1262 auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
1263 if (!bc)
1264 return std::nullopt;
1265
1266 auto pdf_field =
1267 bc->block->template uncheckedFastGetData<PdfField>(m_pdf_field_id);
1268 auto const density = lbm::accessor::Density::get(pdf_field, bc->cell);
1269 return {double_c(density)};
1270 }
1271
1272 bool set_node_density(Utils::Vector3i const &node, double density) override {
1274 auto bc = get_block_and_cell(get_lattice(), node, false);
1275 if (!bc)
1276 return false;
1277
1278 auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
1279 lbm::accessor::Density::set(pdf_field, FloatType_c(density), bc->cell);
1280
1281 return true;
1282 }
1283
1284 std::vector<double>
1286 Utils::Vector3i const &upper_corner) const override {
1287 std::vector<double> out;
1288 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1289 auto const &lattice = get_lattice();
1290 auto const &block = *(lattice.get_blocks()->begin());
1291 auto const pdf_field = block.template getData<PdfField>(m_pdf_field_id);
1292 auto const values = lbm::accessor::Density::get(pdf_field, *ci);
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,
1296 double>) {
1297 out = std::move(values);
1298 } else {
1299 out = std::vector<double>(values.begin(), values.end());
1300 }
1301 }
1302 return out;
1303 }
1304
1305 void set_slice_density(Utils::Vector3i const &lower_corner,
1306 Utils::Vector3i const &upper_corner,
1307 std::vector<double> const &density) override {
1309 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1310 auto const &lattice = get_lattice();
1311 auto &block = *(lattice.get_blocks()->begin());
1312 auto pdf_field = block.template getData<PdfField>(m_pdf_field_id);
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());
1316 lbm::accessor::Density::set(pdf_field, values, *ci);
1317 }
1318 }
1319
1320 std::optional<Utils::Vector3d>
1322 bool consider_ghosts = false) const override {
1323 assert(not(consider_ghosts and m_pending_ghost_comm.test(GhostComm::UBB)));
1324 auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
1325 if (!bc or !m_boundary->node_is_boundary(node))
1326 return std::nullopt;
1327
1328 return {to_vector3d(m_boundary->get_node_value_at_boundary(node))};
1329 }
1330
1332 Utils::Vector3d const &velocity) override {
1335 auto bc = get_block_and_cell(get_lattice(), node, true);
1336 if (bc) {
1337 m_boundary->set_node_value_at_boundary(
1338 node, to_vector3<FloatType>(velocity), *bc);
1339 }
1340 return bc.has_value();
1341 }
1342
1343 std::vector<std::optional<Utils::Vector3d>> get_slice_velocity_at_boundary(
1344 Utils::Vector3i const &lower_corner,
1345 Utils::Vector3i const &upper_corner) const override {
1346 std::vector<std::optional<Utils::Vector3d>> out;
1347 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1348 auto const &lattice = get_lattice();
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) {
1357 auto const node = local_offset + Utils::Vector3i{{x, y, z}};
1358 if (m_boundary->node_is_boundary(node)) {
1359 out.emplace_back(
1360 to_vector3d(m_boundary->get_node_value_at_boundary(node)));
1361 } else {
1362 out.emplace_back(std::nullopt);
1363 }
1364 }
1365 }
1366 }
1367 assert(out.size() == n_values);
1368 }
1369 return out;
1370 }
1371
1373 Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner,
1374 std::vector<std::optional<Utils::Vector3d>> const &velocity) override {
1377 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1378 auto const &lattice = get_lattice();
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();
1382 auto it = velocity.begin();
1383 assert(velocity.size() == ci->numCells());
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) {
1387 auto const node = local_offset + Utils::Vector3i{{x, y, z}};
1388 auto const bc = get_block_and_cell(lattice, node, false);
1389 auto const &opt = *it;
1390 if (opt) {
1391 m_boundary->set_node_value_at_boundary(
1392 node, to_vector3<FloatType>(*opt), *bc);
1393 } else {
1394 m_boundary->remove_node_from_boundary(node, *bc);
1395 }
1396 ++it;
1397 }
1398 }
1399 }
1400 }
1401 }
1402
1403 std::optional<Utils::Vector3d>
1404 get_node_boundary_force(Utils::Vector3i const &node) const override {
1405 auto const bc = get_block_and_cell(get_lattice(), node, true);
1406 if (!bc or !m_boundary->node_is_boundary(node))
1407 return std::nullopt;
1408
1409 return get_node_last_applied_force(node, true);
1410 }
1411
1412 bool remove_node_from_boundary(Utils::Vector3i const &node) override {
1413 auto bc = get_block_and_cell(get_lattice(), node, true);
1414 if (bc) {
1415 m_boundary->remove_node_from_boundary(node, *bc);
1416 }
1417 return bc.has_value();
1418 }
1419
1420 std::optional<bool>
1422 bool consider_ghosts = false) const override {
1423 assert(not(consider_ghosts and m_pending_ghost_comm.test(GhostComm::UBB)));
1424 auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
1425 if (!bc)
1426 return std::nullopt;
1427
1428 return {m_boundary->node_is_boundary(node)};
1429 }
1430
1431 std::vector<bool>
1433 Utils::Vector3i const &upper_corner) const override {
1434 std::vector<bool> out;
1435 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1436 auto const &lattice = get_lattice();
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) {
1445 auto const node = local_offset + Utils::Vector3i{x, y, z};
1446 out.emplace_back(m_boundary->node_is_boundary(node));
1447 }
1448 }
1449 }
1450 assert(out.size() == n_values);
1451 }
1452 return out;
1453 }
1454
1455 void reallocate_ubb_field() override { m_boundary->boundary_update(); }
1456
1458 if (not m_has_boundaries) {
1459 m_has_boundaries = true;
1461 }
1462 m_has_boundaries = true;
1463 }
1464
1465 void clear_boundaries() override {
1466 reset_boundary_handling();
1469 m_has_boundaries = false;
1471 }
1472
1473 void
1474 update_boundary_from_shape(std::vector<int> const &raster_flat,
1475 std::vector<double> const &data_flat) override {
1478 auto const grid_size = get_lattice().get_grid_dimensions();
1479 auto const data = fill_3D_vector_array(data_flat, grid_size);
1480 set_boundary_from_grid(*m_boundary, get_lattice(), raster_flat, data);
1483 }
1484
1485 // Pressure tensor
1486 std::optional<Utils::VectorXd<9>>
1487 get_node_pressure_tensor(Utils::Vector3i const &node) const override {
1488 auto bc = get_block_and_cell(get_lattice(), node, false);
1489 if (!bc)
1490 return std::nullopt;
1491
1492 auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
1493 auto tensor = lbm::accessor::PressureTensor::get(pdf_field, bc->cell);
1494 pressure_tensor_correction(tensor);
1495
1496 return to_vector9d(tensor);
1497 }
1498
1499 std::vector<double> get_slice_pressure_tensor(
1500 Utils::Vector3i const &lower_corner,
1501 Utils::Vector3i const &upper_corner) const override {
1502 std::vector<double> out;
1503 if (auto const ci = get_interval(lower_corner, upper_corner)) {
1504 auto const &lattice = get_lattice();
1505 auto const &block = *(lattice.get_blocks()->begin());
1506 auto const pdf_field = block.template getData<PdfField>(m_pdf_field_id);
1507 auto values = lbm::accessor::PressureTensor::get(pdf_field, *ci);
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));
1512 }
1513 if constexpr (std::is_same_v<typename decltype(values)::value_type,
1514 double>) {
1515 out = std::move(values);
1516 } else {
1517 out = std::vector<double>(values.begin(), values.end());
1518 }
1519 }
1520 return out;
1521 }
1522
1523 // Global pressure tensor
1524 [[nodiscard]] Utils::VectorXd<9> get_pressure_tensor() const override {
1525 auto const &blocks = get_lattice().get_blocks();
1526 Matrix3<FloatType> tensor(FloatType{0});
1527 for (auto block = blocks->begin(); block != blocks->end(); ++block) {
1528 auto pdf_field = block->template getData<PdfField>(m_pdf_field_id);
1529 tensor += lbm::accessor::PressureTensor::reduce(pdf_field);
1530 }
1531 auto const grid_size = get_lattice().get_grid_dimensions();
1532 auto const number_of_nodes = Utils::product(grid_size);
1533 pressure_tensor_correction(tensor);
1534 return to_vector9d(tensor) * (1. / static_cast<double>(number_of_nodes));
1535 }
1536
1537 // Global momentum
1538 [[nodiscard]] Utils::Vector3d get_momentum() const override {
1539 auto const &blocks = get_lattice().get_blocks();
1540 Vector3<FloatType> mom(FloatType{0});
1541 for (auto block = blocks->begin(); block != blocks->end(); ++block) {
1542 auto pdf_field = block->template getData<PdfField>(m_pdf_field_id);
1543 auto force_field =
1544 block->template getData<VectorField>(m_last_applied_force_field_id);
1545 mom += lbm::accessor::MomentumDensity::reduce(pdf_field, force_field);
1546 }
1547 return to_vector3d(mom);
1548 }
1549
1550 // Global external force
1551 void set_external_force(Utils::Vector3d const &ext_force) override {
1552 m_reset_force->set_ext_force(ext_force);
1553 }
1554
1555 [[nodiscard]] Utils::Vector3d get_external_force() const noexcept override {
1556 return m_reset_force->get_ext_force();
1557 }
1558
1559 [[nodiscard]] double get_kT() const noexcept override {
1560 return numeric_cast<double>(m_kT);
1561 }
1562
1563 [[nodiscard]] unsigned int get_seed() const noexcept override {
1564 return m_seed;
1565 }
1566
1567 [[nodiscard]] std::optional<uint64_t> get_rng_state() const override {
1568 auto const cm = std::get_if<CollisionModelThermalized>(&*m_collision_model);
1569 if (!cm or m_kT == 0.) {
1570 return std::nullopt;
1571 }
1572 return {static_cast<uint64_t>(cm->time_step_)};
1573 }
1574
1575 void set_rng_state(uint64_t counter) override {
1576 auto const cm = std::get_if<CollisionModelThermalized>(&*m_collision_model);
1577 if (!cm or m_kT == 0.) {
1578 throw std::runtime_error("This LB instance is unthermalized");
1579 }
1580 assert(counter <=
1581 static_cast<uint32_t>(std::numeric_limits<uint_t>::max()));
1582 cm->time_step_ = static_cast<uint32_t>(counter);
1583 }
1584
1585 [[nodiscard]] LatticeWalberla const &get_lattice() const noexcept override {
1586 return *m_lattice;
1587 }
1588
1589 [[nodiscard]] std::size_t get_velocity_field_id() const noexcept override {
1590 return m_velocity_field_id;
1591 }
1592
1593 [[nodiscard]] std::size_t get_force_field_id() const noexcept override {
1595 }
1596
1597 void register_vtk_field_filters(walberla::vtk::VTKOutput &vtk_obj) override {
1598 field::FlagFieldCellFilter<FlagField> fluid_filter(m_flag_field_id);
1599 fluid_filter.addFlag(Boundary_flag);
1600 vtk_obj.addCellExclusionFilter(fluid_filter);
1601 }
1602
1603protected:
1604 template <typename Field_T, uint_t F_SIZE_ARG, typename OutputType>
1605 class VTKWriter : public vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG> {
1606 public:
1607 VTKWriter(ConstBlockDataID const &block_id, std::string const &id,
1608 FloatType unit_conversion)
1609 : vtk::BlockCellDataWriter<OutputType, F_SIZE_ARG>(id),
1610 m_block_id(block_id), m_field(nullptr),
1611 m_conversion(unit_conversion) {}
1612
1613 protected:
1614 void configure() override {
1615 WALBERLA_ASSERT_NOT_NULLPTR(this->block_);
1616 m_field = this->block_->template getData<Field_T>(m_block_id);
1617 }
1618
1619 ConstBlockDataID const m_block_id;
1620 Field_T const *m_field;
1621 FloatType const m_conversion;
1622 };
1623
1624 template <typename OutputType = float>
1625 class DensityVTKWriter : public VTKWriter<PdfField, 1u, OutputType> {
1626 public:
1628 using Base::Base;
1629 using Base::evaluate;
1630
1631 protected:
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);
1635 auto const density =
1636 lbm::accessor::Density::get(this->m_field, {x, y, z});
1637 return numeric_cast<OutputType>(this->m_conversion * density);
1638 }
1639 };
1640
1641 template <typename OutputType = float>
1642 class VelocityVTKWriter : public VTKWriter<VectorField, 3u, OutputType> {
1643 public:
1645 using Base::Base;
1646 using Base::evaluate;
1647
1648 protected:
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);
1652 auto const velocity =
1653 lbm::accessor::Vector::get(this->m_field, {x, y, z});
1654 return numeric_cast<OutputType>(this->m_conversion * velocity[uint_c(f)]);
1655 }
1656 };
1657
1658 template <typename OutputType = float>
1659 class PressureTensorVTKWriter : public VTKWriter<PdfField, 9u, OutputType> {
1660 public:
1662 using Base::Base;
1663 using Base::evaluate;
1664
1665 PressureTensorVTKWriter(ConstBlockDataID const &block_id,
1666 std::string const &id, FloatType unit_conversion,
1667 FloatType off_diag_factor)
1668 : Base(block_id, id, unit_conversion),
1669 m_off_diag_factor(off_diag_factor) {}
1670
1671 protected:
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 =
1678 (f == 0 or f == 4 or f == 8) ? FloatType{1} : m_off_diag_factor;
1679 return numeric_cast<OutputType>(this->m_conversion * revert_factor *
1680 pressure[uint_c(f)]);
1681 }
1682 FloatType const m_off_diag_factor;
1683 };
1684
1685public:
1686 void register_vtk_field_writers(walberla::vtk::VTKOutput &vtk_obj,
1687 LatticeModel::units_map const &units,
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);
1697 }
1698 };
1699#endif
1700 if (flag_observables & static_cast<int>(OutputVTK::density)) {
1701 auto const unit_conversion = FloatType_c(units.at("density"));
1702#if defined(__CUDACC__)
1703 if constexpr (Architecture == lbmpy::Arch::GPU) {
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>(
1708 blocks, *m_pdf_cpu_field_id, m_pdf_field_id));
1709 }
1710#endif
1711 vtk_obj.addCellDataWriter(std::make_shared<DensityVTKWriter<float>>(
1712 m_pdf_field_id, "density", unit_conversion));
1713 }
1714 if (flag_observables & static_cast<int>(OutputVTK::velocity_vector)) {
1715 auto const unit_conversion = FloatType_c(units.at("velocity"));
1716#if defined(__CUDACC__)
1717 if constexpr (Architecture == lbmpy::Arch::GPU) {
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>(
1723 blocks, *m_vel_cpu_field_id, m_velocity_field_id));
1724 }
1725#endif
1726 vtk_obj.addCellDataWriter(std::make_shared<VelocityVTKWriter<float>>(
1727 m_velocity_field_id, "velocity_vector", unit_conversion));
1728 }
1729 if (flag_observables & static_cast<int>(OutputVTK::pressure_tensor)) {
1730 auto const unit_conversion = FloatType_c(units.at("pressure"));
1731#if defined(__CUDACC__)
1732 if constexpr (Architecture == lbmpy::Arch::GPU) {
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>(
1737 blocks, *m_pdf_cpu_field_id, m_pdf_field_id));
1738 }
1739#endif
1740 vtk_obj.addCellDataWriter(
1741 std::make_shared<PressureTensorVTKWriter<float>>(
1742 m_pdf_field_id, "pressure_tensor", unit_conversion,
1743 pressure_tensor_correction_factor()));
1744 }
1745 }
1746
1747 ~LBWalberlaImpl() override = default;
1748};
1749
1750} // namespace walberla
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
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.
Definition Vector.hpp:110
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
PressureTensorVTKWriter(ConstBlockDataID const &block_id, std::string const &id, FloatType unit_conversion, FloatType off_diag_factor)
VTKWriter(ConstBlockDataID const &block_id, std::string const &id, FloatType unit_conversion)
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
typename detail::KernelTrait< FloatType, Architecture >::CollisionModelLeesEdwards CollisionModelLeesEdwards
std::size_t get_force_field_id() const noexcept override
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
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
std::shared_ptr< ResetForce< PdfField, VectorField > > m_reset_force
virtual bool is_double_precision() const noexcept override
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_vel() override
typename FieldTrait< FloatType >::PdfField _PdfField
Utils::Vector3d get_external_force() const noexcept override
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
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
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
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)
Definition dpd.cpp:78
static double * block(double *p, std::size_t index, std::size_t size)
Definition elc.cpp:172
T product(Vector< T, N > const &v)
Definition Vector.hpp:374
VectorXi< 3 > Vector3i
Definition Vector.hpp:173
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)
Definition boundary.hpp:81
std::vector< Utils::Vector3d > fill_3D_vector_array(std::vector< double > const &vec_flat, Utils::Vector3i const &grid_size)
Definition boundary.hpp:36
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.
Definition relative.cpp:64
DEVICE_QUALIFIER constexpr iterator begin() noexcept
Definition Array.hpp:128
DEVICE_QUALIFIER constexpr size_type size() const noexcept
Definition Array.hpp:154
DEVICE_QUALIFIER constexpr iterator end() noexcept
Definition Array.hpp:140
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.
@ LAF
last applied forces communication