lbSolver.h

Go to the documentation of this file.

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2012-2021 Mathias J. Krause, Benjamin Förster,
 *  Julius Jessberger, 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 LBSOLVER_H
#define LBSOLVER_H
 
#include "solverParameters.h"
#include "utilities/timer.h"
#include "utilities/typeIndexedContainers.h"
#include "utilities/typeMap.h"
 
 
namespace olb {
 
template <typename T, typename PARAMETERS>
class BaseSolver {
 
private:
  mutable OstreamManager                 clout {std::cout, "BaseSolver"};
 
public:
  using t = T;
  using Parameters_t = PARAMETERS;
  utilities::TypeIndexedSharedPtrTuple<PARAMETERS>           _parameters;
 
protected:
  bool                                   _isInitialized {false};
  std::size_t                            _iT {0};
  bool                                   _finishedTimeLoop {false};
 
public:
  BaseSolver(utilities::TypeIndexedSharedPtrTuple<PARAMETERS> params)
   : _parameters(params)
  { }
 
  virtual void initialize() = 0;
 
  virtual void preProcess() { }
 
  void solve()
  {
    if (! _isInitialized) {
      initialize();
    }
 
    prepareSimulation();
 
    do
    {
      timeStep(_iT);
      ++_iT;
    } while (! exitCondition(_iT));
 
    _finishedTimeLoop = true;
    postSimulation();
  }
 
  virtual void postProcess() { }
 
protected:
  virtual void prepareSimulation() = 0;
 
  virtual void timeStep(std::size_t iT) = 0;
 
  virtual bool exitCondition(std::size_t iT) const = 0;
 
  virtual void postSimulation() = 0;
 
public:
  template <typename KEY>
  auto& parameters(KEY = KEY()) {
    return *_parameters.template get<KEY>();
  }
  template <typename KEY>
  auto& parameters(KEY = KEY()) const {
    return *_parameters.template get<KEY>();
  }
  template <typename KEY>
  auto& parameters(meta::id<KEY>) {
    return *_parameters.template get<KEY>();
  }
  template <typename KEY>
  auto& parameters(meta::id<KEY>) const {
    return *_parameters.template get<KEY>();
  }
};
 
 
template<
  typename T,
  typename PARAMETERS,
  typename LATTICES
>
class LbSolver : public BaseSolver<T,PARAMETERS> {
 
private:
  mutable OstreamManager                             clout {std::cout, "LbSolver"};
 
protected:
  static constexpr unsigned dim = LATTICES::values_t::template get<0>::d;
 
  template<typename... DESCRIPTORS>
  using SuperLattices = std::tuple<std::shared_ptr<SuperLattice<T,DESCRIPTORS>>...>;
 
  std::shared_ptr<SuperGeometry<T,dim>>              _sGeometry;
  std::shared_ptr<CuboidGeometry<T,dim>>             _cGeometry;
  std::shared_ptr<LoadBalancer<T>>                   _loadBalancer;
 
  typename LATTICES::values_t::template decompose_into<SuperLattices>  _sLattices;
 
  std::unique_ptr<util::Timer<BaseType<T>>>          _timer;
 
  static constexpr bool isStationary = PARAMETERS::keys_t::template contains<names::Stationarity>();
  static constexpr unsigned getNumberStationaryLattices() {
    if constexpr (isStationary) {
      return PARAMETERS::template value<names::Stationarity>::numberOfStationaryLattices;
    } else {
      return 0;
    }
    __builtin_unreachable();
  }
  std::array<std::unique_ptr<util::ValueTracer<T>>,
    getNumberStationaryLattices()> _convergenceCheck;
 
  static constexpr bool outputGnuplot = PARAMETERS::keys_t::template contains<names::VisualizationGnuplot>();
  static constexpr bool outputImages  = PARAMETERS::keys_t::template contains<names::VisualizationImages>();
  static constexpr bool outputVTK     = PARAMETERS::keys_t::template contains<names::VisualizationVTK>();
 
  bool                                               _exitMaxU {false};
  BaseType<T>                                        _boundMaxU {1.0};
  std::size_t                                        _itCheckStability {1};
  std::size_t                                        _itBoundaryUpdate {1};
 
public:
  LbSolver(utilities::TypeIndexedSharedPtrTuple<PARAMETERS> params) : LbSolver::BaseSolver(params)
  { }
 
  // Allows to return the geometric/ lattice data without solving
  void buildAndReturn();
 
  void initialize() override;
 
protected:
  void prepareSimulation() override;
 
  void timeStep(std::size_t iT) override;
 
  void postSimulation() override;
 
 
  // Inheritant has to provide this method
  virtual void prepareGeometry() = 0;
 
  // Inheritant has to provide this method
  virtual void prepareLattices() = 0;
 
  // Inheritant has to provide this method
  virtual void setInitialValues() = 0;
 
  // Inheritant has to provide this method
  virtual void setBoundaryValues(std::size_t iT) = 0;
 
  // Inheritant may provide this method
  virtual void getResults(std::size_t iT) { };
 
  // Inheritant may provide these method
  // This function is called at every time step
  virtual void computeResults(std::size_t iT) { };
  // This function is called after the simulation
  virtual void computeResults() { };
 
  virtual bool exitCondition(std::size_t iT) const override;
 
  virtual bool checkStability(std::size_t iT);
 
private:
  void build();
  void renewLattices();
 
  // ------------ Output generation -------------------------------------------
protected:
  virtual void printLog(std::size_t iT) const;
 
  virtual void writeLogConverter() const;
 
  virtual void prepareVTK() const;
 
  // Inheritant may override this method for vtk output
  virtual void writeVTK(std::size_t iT) const { };
 
  // Inheritant may override this method for image output
  virtual void writeImages(std::size_t iT) const { };
 
  // Inheritant may override this method for gnuplot output
  virtual void writeGnuplot(std::size_t iT) const { };
 
  // ------------ Access to converters and lattices ---------------------------
  template<typename... ARGS>
  auto& converter(ARGS&&... args) {
    using SimulationParameters_t = typename PARAMETERS::template value<names::Simulation>;
    if constexpr(std::is_invocable_v<decltype(&SimulationParameters_t::converter),ARGS...>) {
      return this->parameters(names::Simulation()).converter(args...);
    }
    else {
      return *this->parameters(names::Simulation()).converter;
    }
    __builtin_unreachable();
  }
  template<typename... ARGS>
  auto& converter(ARGS&&... args) const {
    using SimulationParameters_t = typename PARAMETERS::template value<names::Simulation>;
    if constexpr(std::is_invocable_v<decltype(&SimulationParameters_t::converter),ARGS...>) {
      return this->parameters(names::Simulation()).converter(args...);
    }
    else {
      return *this->parameters(names::Simulation()).converter;
    }
    __builtin_unreachable();
  }
  template <typename KEY>
  auto& lattice(KEY = KEY()) {
    return *std::get<(LATTICES::keys_t::template index<KEY>())>(_sLattices);
  }
  template <typename KEY>
  auto& lattice(KEY = KEY()) const {
    return *std::get<(LATTICES::keys_t::template index<KEY>())>(_sLattices);
  }
  template <typename KEY>
  auto& lattice(meta::id<KEY>) {
    return *std::get<(LATTICES::keys_t::template index<KEY>())>(_sLattices);
  }
  template <typename KEY>
  auto& lattice(meta::id<KEY>) const {
    return *std::get<(LATTICES::keys_t::template index<KEY>())>(_sLattices);
  }
  auto& lattice() {
    static_assert((LATTICES::size == 1), "Lattice name must be provided");
    return *std::get<0>(_sLattices);
  }
  auto& lattice() const {
    static_assert((LATTICES::size == 1), "Lattice name must be provided");
    return *std::get<0>(_sLattices);
  }
  auto& geometry() {
    return *_sGeometry;
  }
  auto& geometry() const {
    return *_sGeometry;
  }
  /*auto& results() {
    return *(this->_simulationResults);
  }
  auto& results() const {
    return *(this->_simulationResults);
  }*/
};
 
 
template<typename T, typename SOLVER>
std::function<T (const std::vector<T>&, unsigned)> getCallable(std::shared_ptr<SOLVER> solver){
  return [=](const std::vector<T>& control, unsigned optiStep) -> T {
    solver->parameters(names::Opti()).applyControl(control);
    solver->parameters(names::OutputOpti()).counterOptiStep = optiStep;
    solver->solve();
    return solver->parameters(names::Results()).objective;
  };
}
 
template<typename SOLVER>
std::function<void ()> doPostProcess(std::shared_ptr<SOLVER> solver){
  return [=]() {
    solver->postProcess();
  };
}
 
 
// ------------------ Creator functions for XML interface ---------------------
 
template <class SOLVER>
std::shared_ptr<SOLVER> createLbSolver(XMLreader const& xml)
{
  using parameters_map = typename SOLVER::BaseSolver::Parameters_t;
  auto paramsTuple = createParametersTuple<parameters_map>(xml);
  return std::make_shared<SOLVER>(paramsTuple);
}
 
 
} // namespace olb
 
#endif