interface.h

Go to the documentation of this file.

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2021 Adrian Kummerlaender
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program 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 2
 *  of the License, or (at your option) any later version.
 *
 *  This program 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 this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/
 
#ifndef DYNAMICS_INTERFACE_H
#define DYNAMICS_INTERFACE_H
 
#include "lbm.h"
#include "core/concepts.h"
#include "momenta/interface.h"
 
#include <type_traits>
 
namespace olb {
 
template <typename T, typename DESCRIPTOR> class Cell;
template <typename T, typename DESCRIPTOR> class ConstCell;
 
 
template <typename T>
struct CellStatistic {
  T rho;
  T uSqr;
 
  operator bool() const any_platform {
    return rho != T{-1} && uSqr != T{-1};
  }
};
 
template <typename T, typename DESCRIPTOR>
struct Dynamics {
  using value_t = T;
  using descriptor_t = DESCRIPTOR;
 
  virtual ~Dynamics() any_platform { }
 
  virtual std::type_index id() = 0;
  virtual std::string getName() const {
    return "Dynamics";
  };
 
  virtual void initialize(Cell<T,DESCRIPTOR>& cell) { };
 
  virtual AbstractParameters<T,DESCRIPTOR>& getParameters(BlockLattice<T,DESCRIPTOR>& block) = 0;
 
  virtual CellStatistic<T> collide(Cell<T,DESCRIPTOR>& cell) {
    throw std::bad_function_call();
  };
 
  virtual T computeRho(ConstCell<T,DESCRIPTOR>& cell) const = 0;
  virtual void computeU(ConstCell<T,DESCRIPTOR>& cell, T u[DESCRIPTOR::d]) const = 0;
  virtual void computeJ(ConstCell<T,DESCRIPTOR>& cell, T j[DESCRIPTOR::d]) const = 0;
  virtual void computeStress(ConstCell<T,DESCRIPTOR>& cell,
                             T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n]) const = 0;
  virtual void computeRhoU(ConstCell<T,DESCRIPTOR>& cell, T& rho, T u[DESCRIPTOR::d]) const = 0;
  virtual void computeAllMomenta(ConstCell<T,DESCRIPTOR>& cell,
                                 T& rho, T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n]) const = 0;
 
  virtual void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) = 0;
  virtual void defineU(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d]) = 0;
  virtual void defineRhoU(Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d]) = 0;
  virtual void defineAllMomenta(Cell<T,DESCRIPTOR>& cell,
                                T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal<DESCRIPTOR >::n]) = 0;
 
  virtual T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d]) const any_platform = 0;
  void iniEquilibrium(Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d]) {
    T rhoShift(rho);
    T uShift[DESCRIPTOR::d];
    util::copyN(uShift, u, DESCRIPTOR::d);
    inverseShiftRhoU(cell, rhoShift, uShift);
    for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
      cell[iPop] = computeEquilibrium(iPop, rhoShift, uShift);
    }
  };
  void iniRegularized(Cell<T,DESCRIPTOR>& cell,
                      T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal<DESCRIPTOR>::n])  {
    iniEquilibrium(cell, rho, u);
    for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
      cell[iPop] += equilibrium<DESCRIPTOR>::template fromPiToFneq<T>(iPop, pi);
    }
  };
  // This is relevant for correct initialization for some forcing/ sourcing schemes
  virtual void inverseShiftRhoU(ConstCell<T,DESCRIPTOR>& cell, T& rho, T u[DESCRIPTOR::d]) const { };
 
};
 
namespace dynamics {
 
 
struct DefaultCombination {
  static std::string getName() {
    return "Default";
  }
 
  template <typename DESCRIPTOR, typename MOMENTA>
  using combined_momenta = typename MOMENTA::template type<DESCRIPTOR>;
 
  template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM>
  using combined_equilibrium = typename EQUILIBRIUM::template type<DESCRIPTOR,MOMENTA>;
 
  template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM, typename COLLISION>
  using combined_collision = typename COLLISION::template type<DESCRIPTOR,MOMENTA,EQUILIBRIUM>;
 
  template <typename DESCRIPTOR, typename MOMENTA, typename EQUILIBRIUM, typename COLLISION>
  using combined_parameters = typename COLLISION::parameters;
};
 
 
template <typename DYNAMICS, typename = void>
struct is_vectorizable : std::true_type { };
 
 
template <typename DYNAMICS>
struct is_vectorizable<
  DYNAMICS,
  std::enable_if_t<!DYNAMICS::is_vectorizable>
> : std::false_type { };
 
template <typename DYNAMICS>
static constexpr bool is_vectorizable_v = is_vectorizable<DYNAMICS>::value;
 
 
 
template <
  typename T, typename DESCRIPTOR,
  typename MOMENTA, typename EQUILIBRIUM, typename COLLISION,
  typename COMBINATION_RULE = DefaultCombination
>
struct Tuple final : public Dynamics<T,DESCRIPTOR> {
  using MomentaF     = typename COMBINATION_RULE::template combined_momenta<DESCRIPTOR,MOMENTA>;
  using EquilibriumF = typename COMBINATION_RULE::template combined_equilibrium<DESCRIPTOR,MOMENTA,EQUILIBRIUM>;
  using CollisionO   = typename COMBINATION_RULE::template combined_collision<DESCRIPTOR,MOMENTA,EQUILIBRIUM,COLLISION>;
 
  using parameters = typename COMBINATION_RULE::template combined_parameters<DESCRIPTOR,MOMENTA,EQUILIBRIUM,COLLISION>;
 
  constexpr static bool is_vectorizable = is_vectorizable_v<CollisionO>;
 
  template <typename NEW_MOMENTA>
  using exchange_momenta = Tuple<
    T, DESCRIPTOR,
    NEW_MOMENTA, EQUILIBRIUM, COLLISION,
    COMBINATION_RULE
  >;
 
  template <typename NEW_RULE>
  using exchange_combination_rule = Tuple<
    T, DESCRIPTOR,
    MOMENTA, EQUILIBRIUM, COLLISION,
    NEW_RULE
  >;
 
  template <template<typename> typename WRAPPER>
  using wrap_collision = Tuple<
    T, DESCRIPTOR,
    MOMENTA, EQUILIBRIUM, WRAPPER<COLLISION>,
    COMBINATION_RULE
  >;
 
  std::type_index id() override {
    return typeid(Tuple);
  }
 
  std::string getName() const override {
    return "dynamics::Tuple<"
      + MomentaF().getName() + ","
      + EQUILIBRIUM::getName() + ","
      + COLLISION::getName() + ","
      + COMBINATION_RULE::getName() +
    ">";
  };
 
  template <typename FIELD>
  constexpr bool hasParameter() const {
    return parameters::template contains<FIELD>();
  };
 
  AbstractParameters<T,DESCRIPTOR>& getParameters(BlockLattice<T,DESCRIPTOR>& block) override {
    return block.template getData<OperatorParameters<Tuple>>();
  }
 
  void initialize(Cell<T,DESCRIPTOR>& cell) override {
    MomentaF().initialize(cell);
  };
 
  template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
  CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform {
    // Copy parameters to enable changes in composed dynamics
    auto params = static_cast<
      ParametersOfOperatorD<T,DESCRIPTOR,Tuple>&
    >(parameters).template copyAs<V>();
    return CollisionO().apply(cell, params);
  };
 
  T computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d]) const override any_platform {
    return EquilibriumF().compute(iPop, rho, u);
  };
 
  T computeRho(ConstCell<T,DESCRIPTOR>& cell) const override {
    return MomentaF().computeRho(cell);
  };
  void computeU(ConstCell<T,DESCRIPTOR>& cell, T u[DESCRIPTOR::d]) const override {
    MomentaF().computeU(cell, u);
  };
  void computeJ(ConstCell<T,DESCRIPTOR>& cell, T j[DESCRIPTOR::d]) const override {
    MomentaF().computeJ(cell, j);
  };
  void computeStress(ConstCell<T,DESCRIPTOR>& cell,
                     T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n]) const override {
    MomentaF().computeStress(cell, rho, u, pi);
  };
  void computeRhoU(ConstCell<T,DESCRIPTOR>& cell, T& rho, T u[DESCRIPTOR::d]) const override {
    MomentaF().computeRhoU(cell, rho, u);
  };
  void computeAllMomenta(ConstCell<T,DESCRIPTOR>& cell,
                         T& rho, T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n]) const override {
    MomentaF().computeAllMomenta(cell, rho, u, pi);
  };
 
  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override {
    MomentaF().defineRho(cell, rho);
  };
  void defineU(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d]) override {
    MomentaF().defineU(cell, u);
  };
  void defineRhoU(Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d]) override {
    MomentaF().defineRhoU(cell, rho, u);
  };
  void defineAllMomenta(Cell<T,DESCRIPTOR>& cell,
                        T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal<DESCRIPTOR >::n]) override {
    MomentaF().defineAllMomenta(cell, rho, u, pi);
  };
  void inverseShiftRhoU(ConstCell<T,DESCRIPTOR>& cell, T& rho, T u[DESCRIPTOR::d]) const override {
    MomentaF().inverseShiftRhoU(cell, rho, u);
  }
};
 
template <typename T, typename DESCRIPTOR, typename MOMENTA>
struct CustomCollision : public Dynamics<T,DESCRIPTOR> {
  using value_t = T;
  using descriptor_t = DESCRIPTOR;
 
  using MomentaF = typename MOMENTA::template type<DESCRIPTOR>;
 
  void initialize(Cell<T,DESCRIPTOR>& cell) override {
    MomentaF().initialize(cell);
  };
 
  T computeRho(ConstCell<T,DESCRIPTOR>& cell) const override {
    return MomentaF().computeRho(cell);
  };
  void computeU(ConstCell<T,DESCRIPTOR>& cell, T u[DESCRIPTOR::d]) const override {
    MomentaF().computeU(cell, u);
  };
  void computeJ(ConstCell<T,DESCRIPTOR>& cell, T j[DESCRIPTOR::d]) const override {
    MomentaF().computeJ(cell, j);
  };
  void computeStress(ConstCell<T,DESCRIPTOR>& cell,
                     T rho, const T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n]) const override {
    MomentaF().computeStress(cell, rho, u, pi);
  };
  void computeRhoU(ConstCell<T,DESCRIPTOR>& cell, T& rho, T u[DESCRIPTOR::d]) const override {
    MomentaF().computeRhoU(cell, rho, u);
  };
  void computeAllMomenta(ConstCell<T,DESCRIPTOR>& cell,
                         T& rho, T u[DESCRIPTOR::d], T pi[util::TensorVal<DESCRIPTOR >::n]) const override {
    MomentaF().computeAllMomenta(cell, rho, u, pi);
  };
 
  void defineRho(Cell<T,DESCRIPTOR>& cell, T rho) override {
    MomentaF().defineRho(cell, rho);
  };
  void defineU(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d]) override {
    MomentaF().defineU(cell, u);
  };
  void defineRhoU(Cell<T,DESCRIPTOR>& cell, T rho, const T u[DESCRIPTOR::d]) override {
    MomentaF().defineRhoU(cell, rho, u);
  };
  void defineAllMomenta(Cell<T,DESCRIPTOR>& cell,
                        T rho, const T u[DESCRIPTOR::d], const T pi[util::TensorVal<DESCRIPTOR >::n]) override {
    MomentaF().defineAllMomenta(cell, rho, u, pi);
  };
  void inverseShiftRhoU(ConstCell<T,DESCRIPTOR>& cell, T& rho, T u[DESCRIPTOR::d]) const override {
    MomentaF().inverseShiftRhoU(cell, rho, u);
  }
};
 
 
template <typename PARAMETER, typename DYNAMICS>
struct ParameterFromCell final : public CustomCollision<
  typename DYNAMICS::value_t,
  typename DYNAMICS::descriptor_t,
  typename DYNAMICS::MomentaF::abstract
> {
  using value_t = typename DYNAMICS::value_t;
  using descriptor_t = typename DYNAMICS::descriptor_t;
 
  using MomentaF = typename DYNAMICS::MomentaF;
 
  using parameters = typename DYNAMICS::parameters;
 
  template<typename M>
  using exchange_momenta = ParameterFromCell<PARAMETER,typename DYNAMICS::template exchange_momenta<M>>;
 
  std::type_index id() override {
    return typeid(ParameterFromCell);
  };
 
  AbstractParameters<value_t,descriptor_t>&
  getParameters(BlockLattice<value_t,descriptor_t>& block) override {
    return block.template getData<OperatorParameters<ParameterFromCell>>();
  }
 
  template <CONCEPT(MinimalCell) CELL, typename PARAMETERS, typename V=typename CELL::value_t>
  CellStatistic<V> apply(CELL& cell, PARAMETERS& parameters) any_platform {
    parameters.template set<PARAMETER>(
      cell.template getField<PARAMETER>());
    return DYNAMICS().apply(cell, parameters);
  };
 
  value_t computeEquilibrium(int iPop, value_t rho, const value_t u[descriptor_t::d]) const override {
    return DYNAMICS().computeEquilibrium(iPop, rho, u);
  };
 
  std::string getName() const override {
    return "ParameterFromCell<" + DYNAMICS().getName() + ">";
  };
 
};
 
template <typename DYNAMICS, typename CELL, typename PARAMETERS, typename = void>
struct is_generic : std::false_type { };
 
template <typename DYNAMICS, typename CELL, typename PARAMETERS>
struct is_generic<
  DYNAMICS, CELL, PARAMETERS,
  std::enable_if_t<std::is_member_function_pointer_v<decltype(&DYNAMICS::template apply<CELL,PARAMETERS>)>>
> : std::true_type { };
 
 
template <typename T, typename DESCRIPTOR, typename DYNAMICS>
static constexpr bool is_generic_v = is_generic<DYNAMICS, Cell<T,DESCRIPTOR>, AbstractParameters<T,DESCRIPTOR>>::value;
 
template <typename DYNAMICS, typename = void>
struct has_parametrized_momenta : std::false_type { };
 
 
template <typename DYNAMICS>
struct has_parametrized_momenta<
  DYNAMICS,
  std::enable_if_t<DYNAMICS::has_parametrized_momenta>
> : std::true_type { };
 
template <typename DYNAMICS>
static constexpr bool has_parametrized_momenta_v = has_parametrized_momenta<DYNAMICS>::value;
 
}
 
}
 
#endif