24#ifndef SIMD_OPERATOR_H
25#define SIMD_OPERATOR_H
45 operator bool()
const {
56template <
typename T,
unsigned D>
64 std::array<Pack<T>,D> _packs;
80 meta::call_n_times<D>([&](
unsigned iD) {
81 _packs[iD] = &_data[iD][_index];
88 _packs(rhs._packs) { }
90 template <
typename U,
typename IMPL>
93 for (
unsigned iDim=0; iDim < D; ++iDim) {
108template <
typename T,
typename DESCRIPTOR,
typename V,
typename... RW_FIELDS>
121 std::tuple<FieldD<V,DESCRIPTOR,RW_FIELDS>...> _fields;
134 using FIELD =
typename decltype(field)::type;
135 auto& array = _lattice.template getField<FIELD>();
136 auto& pack = std::get<(rw_fields::template index<FIELD>())>(_fields);
138 for (
unsigned iD=0; iD < DESCRIPTOR::template size<FIELD>(); ++iD) {
148 using FIELD =
typename decltype(field)::type;
149 auto& array = _lattice.template getField<FIELD>();
150 auto& pack = std::get<(rw_fields::template index<FIELD>())>(_fields);
152 for (
unsigned iD=0; iD < DESCRIPTOR::template size<FIELD>(); ++iD) {
160 return std::get<(rw_fields::template index<descriptors::POPULATION>())>(_fields)[iPop];
164 template <
typename FIELD>
166 if constexpr (rw_fields::template contains<FIELD>()) {
167 if constexpr (DESCRIPTOR::template size<FIELD>() == 1) {
168 return std::get<(rw_fields::template index<FIELD>())>(_fields)[0];
170 return std::get<(rw_fields::template index<FIELD>())>(_fields);
173 auto& fieldArray = _lattice.template getField<FIELD>();
174 if constexpr (DESCRIPTOR::template size<FIELD>() == 1) {
178 return &fieldArray[iD][_iCell];
182 __builtin_unreachable();
186 template <
typename FIELD>
188 if constexpr (rw_fields::template contains<FIELD>()) {
189 std::get<(rw_fields::template index<FIELD>())>(_fields) = value;
191 auto& array = _lattice.template getField<FIELD>();
192 for (
unsigned iD=0; iD < DESCRIPTOR::template size<FIELD>(); ++iD) {
199 template <
typename FIELD>
202 return std::get<(rw_fields::template index<FIELD>())>(_fields);
206 template <
typename FIELD>
209 return getField<FIELD>();
213 template <
typename FIELD>
214 std::enable_if_t<rw_fields::template contains<FIELD>(),V&>
216 return std::get<(rw_fields::template index<FIELD>())>(_fields)[iD];
220 template <
typename FIELD>
221 std::enable_if_t<!rw_fields::template contains<FIELD>(),V>
223 return &_lattice.template getField<FIELD>()[iD][_iCell];
230template <
typename T,
typename DESCRIPTOR,
typename DYNAMICS>
237 _parameters{parameters} {
241 return DYNAMICS().apply(cell, *_parameters);
245 return typename DYNAMICS::MomentaF().computeRho(cell);
248 typename DYNAMICS::MomentaF().computeU(cell, u);
251 typename DYNAMICS::MomentaF().computeJ(cell, j);
254 typename DYNAMICS::MomentaF().computeRhoU(cell, rho, u);
257 typename DYNAMICS::MomentaF().computeStress(cell, rho, u, pi);
260 typename DYNAMICS::MomentaF().computeAllMomenta(cell, rho, u, pi);
264 if constexpr (DYNAMICS::parameters::template contains<descriptors::OMEGA>()) {
265 return _parameters->template get<descriptors::OMEGA>();
269 __builtin_unreachable();
273 return DYNAMICS().computeEquilibrium(iPop, rho, u);
277 typename DYNAMICS::MomentaF().defineRho(cell, rho);
281 typename DYNAMICS::MomentaF().defineU(cell, u);
285 typename DYNAMICS::MomentaF().defineRhoU(cell, rho, u);
289 typename DYNAMICS::MomentaF().defineAllMomenta(cell, rho, u, pi);
293 typename DYNAMICS::MomentaF().inverseShiftRhoU(cell, rho, u);
309template <
typename T,
typename DESCRIPTOR,
typename DYNAMICS>
313 std::unique_ptr<DYNAMICS> _dynamics;
314 std::unique_ptr<cpu::Dynamics<T,DESCRIPTOR,Platform::CPU_SIMD>> _concreteDynamics;
327 if (
auto cellStatistic = _dynamicsOfCells[iCell]->collide(cell)) {
328 statistics.
increment(cellStatistic.rho, cellStatistic.uSqr);
340 if constexpr (dynamics::is_vectorizable_v<DYNAMICS>) {
343 auto cellStatistic = DYNAMICS().apply(cell, parameters);
344 for (
unsigned i=0; i < cpu::simd::Pack<T>::size; ++i) {
346 if (cellStatistic.rho[i] != T{-1}) {
347 statistics.
increment(cellStatistic.rho[i], cellStatistic.uSqr[i]);
349 }
else if (subdomain[iCell+i]) {
350 applyOther(block, statistics, iCell+i);
354 for (std::size_t i=iCell; i < iCell+cpu::simd::Pack<T>::size; ++i) {
356 applyOther(block, statistics, i);
365 _dynamics(new DYNAMICS()),
366 _parameters(nullptr),
370 std::type_index
id()
const override
372 return typeid(DYNAMICS);
377 return _mask->weight();
380 void set(
CellID iCell,
bool state,
bool overlap)
override
383 if constexpr (!std::is_same_v<DYNAMICS,NoDynamics<T,DESCRIPTOR>>) {
385 _mask->set(iCell, state);
389 _dynamicsOfCells[iCell] = _concreteDynamics.get();
395 return _dynamics.get();
400 _parameters = &block.template getData<OperatorParameters<DYNAMICS>>().parameters;
401 _mask = &block.template getData<DynamicsMask<DYNAMICS>>();
402 if constexpr (dynamics::has_parametrized_momenta_v<DYNAMICS>) {
403 _dynamics->setMomentaParameters(_parameters);
408 _dynamicsOfCells = block.template getField<cpu::DYNAMICS<T,DESCRIPTOR,Platform::CPU_SIMD>>()[0].data();
417 throw std::runtime_error(
"Platform::CPU_SIMD currently only support CollisionDispatchStrategy::Dominant");
422 #ifdef PARALLEL_MODE_OMP
423 #pragma omp declare reduction(+ : typename LatticeStatistics<T>::Aggregatable : omp_out += omp_in) initializer (omp_priv={})
426 if constexpr (dynamics::is_vectorizable_v<DYNAMICS>) {
430 #ifdef PARALLEL_MODE_OMP
431 #pragma omp parallel for schedule(static) reduction(+ : statistics)
434 apply(block, subdomain, mask, *_parameters, statistics, iCell);
437 #ifdef PARALLEL_MODE_OMP
438 #pragma omp parallel for schedule(static) reduction(+ : statistics)
440 for (std::size_t iCell=0; iCell < block.
getNcells(); ++iCell) {
443 if (
auto cellStatistic = DYNAMICS().apply(cell, *_parameters)) {
444 statistics.
increment(cellStatistic.rho, cellStatistic.uSqr);
446 }
else if (subdomain[iCell]) {
447 applyOther(block, statistics, iCell);
459template <
typename T,
typename DESCRIPTOR, CONCEPT(CellOperator) OPERATOR>
461 :
public BlockO<T,DESCRIPTOR,Platform::CPU_SIMD> {
463 std::vector<CellID> _cells;
469 std::type_index
id()
const override
471 return typeid(OPERATOR);
477 _cells.emplace_back(iCell);
488 std::sort(_cells.begin(), _cells.end());
489 _cells.erase(std::unique(_cells.begin(), _cells.end()), _cells.end());
492 if (_cells.size() > 0) {
494 #ifdef PARALLEL_MODE_OMP
495 #pragma omp parallel for schedule(static) firstprivate(cell)
497 for (
CellID iCell : _cells) {
499 OPERATOR().apply(cell);
507template <
typename T,
typename DESCRIPTOR, CONCEPT(CellOperator) OPERATOR>
509 :
public BlockO<T,DESCRIPTOR,Platform::CPU_SIMD> {
511 std::vector<CellID> _cells;
519 std::type_index
id()
const override
521 return typeid(OPERATOR);
527 _cells.emplace_back(iCell);
534 _parameters = &block.template getData<OperatorParameters<OPERATOR>>().parameters;
540 std::sort(_cells.begin(), _cells.end());
541 _cells.erase(std::unique(_cells.begin(), _cells.end()), _cells.end());
544 if (_cells.size() > 0) {
546 #ifdef PARALLEL_MODE_OMP
547 #pragma omp parallel for schedule(static) firstprivate(cell)
549 for (
CellID iCell : _cells) {
551 OPERATOR().apply(cell, *_parameters);
563template <
typename T,
typename DESCRIPTOR, CONCEPT(BlockOperator) OPERATOR>
565 :
public BlockO<T,DESCRIPTOR,Platform::CPU_SIMD> {
567 std::type_index
id()
const override
569 return typeid(OPERATOR);
574 throw std::logic_error(
"BlockO::set not supported for OperatorScope::PerBlock");
579 OPERATOR().setup(block);
584 OPERATOR().apply(block);
LatticeStatistics< T > & getStatistics()
Return a handle to the LatticeStatistics object.
std::size_t getNcells() const
Get number of cells.
Collision operation of concrete DYNAMICS on concrete block lattices of PLATFORM.
Implementation of BlockLattice on a concrete PLATFORM.
Block application of concrete OPERATOR called using SCOPE on PLATFORM.
Cell concept for concrete block lattices on CPU platforms.
void setCellId(std::size_t iCell)
Implementation of the Cell concept for vectorized collision operators.
auto getField() const
Return pack-valued copy of FIELD.
V & operator[](unsigned iPop)
Return reference to iPop population pack.
std::enable_if_t< rw_fields::template contains< FIELD >(), FieldD< V, DESCRIPTOR, FIELD > & > getFieldPointer()
Return reference to pack-valued interim storage vector of r/w field.
std::enable_if_t<!rw_fields::template contains< FIELD >(), FieldD< V, DESCRIPTOR, FIELD > > getFieldPointer()
Return pack-valued copy of non r/w field.
void setField(FieldD< V, DESCRIPTOR, FIELD > &&value)
Set compoents of FIELD from pack-valued vector.
std::enable_if_t<!rw_fields::template contains< FIELD >(), V > getFieldComponent(unsigned iD)
Return pack-valued copy of non r/w field component.
std::enable_if_t< rw_fields::template contains< FIELD >(), V & > getFieldComponent(unsigned iD)
Return reference to pack-valued interim storage component of r/w field.
Cell(ConcreteBlockLattice< T, DESCRIPTOR, Platform::CPU_SIMD > &lattice, std::size_t iCell, Mask< T > &mask)
Load r/w fields into SIMD packs.
~Cell()
Store modified r/w fields back into lattice taking into account the mask.
Plain column for SIMD CPU targets.
Implementation of cpu::Dynamics for concrete DYNAMICS on SIMD blocks.
void computeJ(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T *j) override
CellStatistic< T > collide(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell) override
void inverseShiftRhoU(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T &rho, T *u) override
T computeRho(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell) override
void computeRhoU(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T &rho, T *u) override
void computeStress(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T &rho, T *u, T *pi) override
T computeEquilibrium(int iPop, T rho, T *u) override
void defineAllMomenta(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T &rho, T *u, T *pi) override
void defineU(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T *u) override
void computeAllMomenta(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T &rho, T *u, T *pi) override
ConcreteDynamics(ParametersOfOperatorD< T, DESCRIPTOR, DYNAMICS > *parameters)
T getOmegaOrFallback(T fallback) override
void defineRhoU(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T &rho, T *u) override
void computeU(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T *u) override
void defineRho(cpu::Cell< T, DESCRIPTOR, Platform::CPU_SIMD > &cell, T &rho) override
SIMD-specific pointer to a pack of rows of a D-dimensional field.
Pack< T > * getComponentPointer(unsigned iDim)
FieldPtr< T, D > & operator=(const GenericVector< U, D, IMPL > &rhs)
FieldPtr(ColumnVector< Column< T >, D > &columns, std::size_t index)
FieldPtr(FieldPtr< T, D > &&rhs)
const Pack< T > * getComponentPointer(unsigned iDim) const
Interface for post-processing steps – header file.
void maskstore(T *target, Mask< T > mask, Pack< T > value)
Top level namespace for all of OpenLB.
std::uint32_t CellID
Type for sequential block-local cell indices.
@ Simulation
Data available on host for e.g. functor evaluation.
Platform
OpenLB execution targets.
@ CPU_SIMD
Basic scalar CPU.
CollisionDispatchStrategy
Collision dispatch strategy.
@ Dominant
Apply dominant dynamics using mask and fallback to virtual dispatch for others.
typename ParametersD< T, DESCRIPTOR >::template include< typename OPERATOR::parameters > ParametersOfOperatorD
Deduce ParametersD of OPERATOR w.r.t. T and DESCRIPTOR.
@ PerBlock
Per-block application, i.e. OPERATOR::apply is passed a ConcreteBlockLattice.
@ PerCell
Per-cell application, i.e. OPERATOR::apply is passed a CELL concept implementation.
@ PerCellWithParameters
Per-cell application with parameters, i.e. OPERATOR::apply is passed a CELL concept implementation an...
Collision operation on concrete blocks of PLATFORM.
Base of block-wide operators such as post processors.
cpu::simd::Pack< T > uSqr
Return value of any collision.
Interface for per-cell dynamics.
Generic vector of values supporting basic arithmetic.
constexpr const T & operator[](unsigned iDim) const any_platform
void increment(T rho, T uSqr)
GenericVector< T, D, FieldPtr< T, D > > type
CPU specific field mirroring BlockDynamicsMap.
Virtual interface for dynamically-dispatched dynamics access on CPU targets.