Dirichlet_double_precision.h

Go to the documentation of this file.

//======================================================================================================================
//
//  This file is part of waLBerla. waLBerla is free software: you can
//  redistribute it and/or modify it under the terms of the GNU General Public
//  License as published by the Free Software Foundation, either version 3 of
//  the License, or (at your option) any later version.
//
//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
//  for more details.
//
//  You should have received a copy of the GNU General Public License along
//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \\file Dirichlet_double_precision.h
//! \\author pystencils
//======================================================================================================================
 
// kernel generated with pystencils v1.4+1.ge851f4e, lbmpy v1.4+1.ge9efe34,
// sympy v1.12.1, lbmpy_walberla/pystencils_walberla from waLBerla commit
// 007e77e077ad9d22b5eed6f3d3118240993e553c
 
#pragma once
#include "core/DataTypes.h"
#include "core/logging/Logging.h"
 
#include "blockforest/StructuredBlockForest.h"
#include "core/debug/Debug.h"
#include "domain_decomposition/BlockDataID.h"
#include "domain_decomposition/IBlock.h"
#include "field/FlagField.h"
#include "field/GhostLayerField.h"
 
#include <functional>
#include <memory>
#include <vector>
 
#ifdef __GNUC__
#define RESTRICT __restrict__
#else
#define RESTRICT
#endif
 
#ifdef WALBERLA_BUILD_WITH_HALF_PRECISION_SUPPORT
using walberla::half;
#endif
 
namespace walberla {
namespace pystencils {
 
class Dirichlet_double_precision {
public:
  struct IndexInfo {
    int32_t x;
    int32_t y;
    int32_t z;
    int32_t dir;
    double value;
    IndexInfo(int32_t x_, int32_t y_, int32_t z_, int32_t dir_)
        : x(x_), y(y_), z(z_), dir(dir_), value() {}
    bool operator==(const IndexInfo &o) const {
      return x == o.x && y == o.y && z == o.z && dir == o.dir &&
             floatIsEqual(value, o.value);
    }
  };
 
  class IndexVectors {
  public:
    using CpuIndexVector = std::vector<IndexInfo>;
 
    enum Type { ALL = 0, INNER = 1, OUTER = 2, NUM_TYPES = 3 };
 
    IndexVectors() = default;
    bool operator==(IndexVectors const &other) const {
      return other.cpuVectors_ == cpuVectors_;
    }
 
    auto &indexVector(Type t) { return cpuVectors_[t]; }
    auto const &indexVector(Type t) const { return cpuVectors_[t]; }
    IndexInfo *pointerCpu(Type t) {
      return cpuVectors_[t].empty() ? nullptr : cpuVectors_[t].data();
    }
 
    void syncGPU() {}
 
  private:
    std::vector<CpuIndexVector> cpuVectors_{NUM_TYPES};
  };
 
  Dirichlet_double_precision(
      const std::shared_ptr<StructuredBlockForest> &blocks,
      BlockDataID fieldID_,
      std::function<double(const Cell &,
                           const shared_ptr<StructuredBlockForest> &, IBlock &)>
          &dirichletCallback)
      : elementInitaliser(dirichletCallback), fieldID(fieldID_) {
    auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const) {
      return new IndexVectors();
    };
    indexVectorID = blocks->addStructuredBlockData<IndexVectors>(
        createIdxVector, "IndexField_Dirichlet_double_precision");
  }
 
  void run(IBlock *block);
 
  void operator()(IBlock *block) { run(block); }
 
  void inner(IBlock *block);
 
  void outer(IBlock *block);
 
  Vector3<double> getForce(IBlock * /*block*/) {
 
    WALBERLA_ABORT(
        "Boundary condition was not generated including force calculation.")
    return Vector3<double>(double_c(0.0));
  }
 
  std::function<void(IBlock *)> getSweep() {
    return [this](IBlock *b) { this->run(b); };
  }
 
  std::function<void(IBlock *)> getInnerSweep() {
    return [this](IBlock *b) { this->inner(b); };
  }
 
  std::function<void(IBlock *)> getOuterSweep() {
    return [this](IBlock *b) { this->outer(b); };
  }
 
  template <typename FlagField_T>
  void fillFromFlagField(const std::shared_ptr<StructuredBlockForest> &blocks,
                         ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID,
                         FlagUID domainFlagUID) {
    for (auto &block : *blocks)
      fillFromFlagField<FlagField_T>(blocks, &block, flagFieldID,
                                     boundaryFlagUID, domainFlagUID);
  }
 
  template <typename FlagField_T>
  void fillFromFlagField(const shared_ptr<StructuredBlockForest> &blocks,
                         IBlock *block, ConstBlockDataID flagFieldID,
                         FlagUID boundaryFlagUID, FlagUID domainFlagUID) {
    auto *indexVectors = block->getData<IndexVectors>(indexVectorID);
    auto &indexVectorAll = indexVectors->indexVector(IndexVectors::ALL);
    auto &indexVectorInner = indexVectors->indexVector(IndexVectors::INNER);
    auto &indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER);
 
    auto *flagField = block->getData<FlagField_T>(flagFieldID);
 
    if (!(flagField->flagExists(boundaryFlagUID) &&
          flagField->flagExists(domainFlagUID)))
      return;
 
    auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
    auto domainFlag = flagField->getFlag(domainFlagUID);
 
    auto inner = flagField->xyzSize();
    inner.expand(cell_idx_t(-1));
 
    indexVectorAll.clear();
    indexVectorInner.clear();
    indexVectorOuter.clear();
 
    auto flagWithGLayers = flagField->xyzSizeWithGhostLayer();
    for (auto it = flagField->beginWithGhostLayerXYZ(); it != flagField->end();
         ++it) {
 
      if (!isFlagSet(it, boundaryFlag))
        continue;
      if (flagWithGLayers.contains(it.x() + cell_idx_c(0),
                                   it.y() + cell_idx_c(0),
                                   it.z() + cell_idx_c(0)) &&
          isFlagSet(it.neighbor(0, 0, 0, 0), domainFlag)) {
 
        auto element = IndexInfo(it.x(), it.y(), it.z(), 0);
        double InitialisatonAdditionalData =
            elementInitaliser(Cell(it.x(), it.y(), it.z()), blocks, *block);
        element.value = InitialisatonAdditionalData;
        indexVectorAll.emplace_back(element);
        if (inner.contains(it.x(), it.y(), it.z()))
          indexVectorInner.emplace_back(element);
        else
          indexVectorOuter.emplace_back(element);
      }
    }
 
    indexVectors->syncGPU();
  }
 
private:
  void run_impl(IBlock *block, IndexVectors::Type type);
 
  BlockDataID indexVectorID;
 
  std::function<double(const Cell &, const shared_ptr<StructuredBlockForest> &,
                       IBlock &)>
      elementInitaliser;
 
public:
  BlockDataID fieldID;
};
 
} // namespace pystencils
} // namespace walberla