jacobian_iterative.hh

Go to the documentation of this file.

// -*- tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_PYTHON_MMESH_JACOBIAN_HH
#define DUNE_PYTHON_MMESH_JACOBIAN_HH
 
#include <dune/python/pybind11/functional.h>
#include <dune/python/pybind11/numpy.h>
#include <dune/python/pybind11/pybind11.h>
#include <dune/python/pybind11/stl.h>
 
#include <array>
#include <dune/fem/operator/linear/spoperator.hh>
#include <functional>
#include <list>
#include <map>
#include <memory>
#include <sstream>
#include <type_traits>
 
namespace Dune {
 
namespace Python {
 
namespace MMesh {
using Dune::Indices::_0;
using Dune::Indices::_1;
 
template <class GridPart, class Intersection, class Entity>
const typename GridPart::IntersectionType convert(
    const GridPart& gridPart, const Intersection& intersection,
    const Entity& inside, Dune::PriorityTag<0>) {
  const Entity outside = gridPart.convert(intersection.outside());
 
  for (auto is : intersections(gridPart, inside))
    if (is.outside() == outside) return is;
  DUNE_THROW(InvalidStateException, "Intersection not found!");
}
 
template <class GridPart, class Intersection, class Entity>
std::enable_if_t<
    std::is_convertible_v<Intersection, typename GridPart::IntersectionType>,
    const typename GridPart::IntersectionType>
convert(const GridPart& gridPart, const Intersection& intersection,
        const Entity& inside, Dune::PriorityTag<1>) {
  return intersection;
}
 
template <class GridPart, class Intersection, class Entity>
const typename GridPart::IntersectionType convert(
    const GridPart& gridPart, const Intersection& intersection,
    const Entity& inside) {
  return convert(gridPart, intersection, inside, Dune::PriorityTag<1>{});
}
 
template <class DomainSpace, class RangeSpace,
          class Partition = Dune::Partitions::All>
struct NeighborInterfaceStencil : public Fem::Stencil<DomainSpace, RangeSpace> {
  typedef Fem::Stencil<DomainSpace, RangeSpace> BaseType;
 
  typedef typename DomainSpace::GridPartType DomainGridPart;
  typedef typename RangeSpace::GridPartType RangeGridPart;
 
 public:
  NeighborInterfaceStencil(const DomainSpace& dSpace, const RangeSpace& rSpace)
      : BaseType(dSpace, rSpace),
        domainGridPart_(dSpace.gridPart()),
        rangeGridPart_(rSpace.gridPart()) {}
 
  void setupStencil() const {
    const auto& mmesh = domainGridPart_.grid().getMMesh();
    for (const auto& entity : elements(domainGridPart_, Partition{})) {
      const auto mmeshIntersection = mmesh.asIntersection(entity);
      const auto inside = rangeGridPart_.convert(mmeshIntersection.inside());
      const auto intersection =
          convert(rangeGridPart_, mmeshIntersection, inside);
 
      BaseType::fill(entity, intersection.inside());
      BaseType::fill(entity, intersection.outside());
    }
  }
 
 private:
  const DomainGridPart& domainGridPart_;
  const RangeGridPart& rangeGridPart_;
};
 
template <class DomainSpace, class RangeSpace,
          class Partition = Dune::Partitions::All>
struct InterfaceNeighborStencil : public Fem::Stencil<DomainSpace, RangeSpace> {
  typedef Fem::Stencil<DomainSpace, RangeSpace> BaseType;
 
  typedef typename DomainSpace::GridPartType DomainGridPart;
  typedef typename RangeSpace::GridPartType RangeGridPart;
 
 public:
  InterfaceNeighborStencil(const DomainSpace& dSpace, const RangeSpace& rSpace)
      : BaseType(dSpace, rSpace),
        domainGridPart_(dSpace.gridPart()),
        rangeGridPart_(rSpace.gridPart()) {}
 
  void setupStencil() const {
    for (const auto& entity : elements(domainGridPart_, Partition{})) {
      for (const auto& intersection : intersections(domainGridPart_, entity)) {
        if (intersection.neighbor() &&
            domainGridPart_.grid().isInterface(intersection)) {
          auto ientity = domainGridPart_.grid().asInterfaceEntity(intersection);
          BaseType::fill(entity, ientity);
        }
      }
    }
  }
 
 private:
  const DomainGridPart& domainGridPart_;
  const RangeGridPart& rangeGridPart_;
};
 
template <class X, class BulkDF, class InterfaceDF>
class ParallelizedScalarProduct : public ScalarProduct<X> {
 public:
  typedef X domain_type;
  typedef typename X::field_type field_type;
  typedef typename FieldTraits<field_type>::real_type real_type;
 
  ParallelizedScalarProduct(const BulkDF& uh, const InterfaceDF& th)
      : u(uh), v(uh), t(th), s(th) {}
 
  virtual field_type dot(const X& x, const X& y) const {
    u.blockVector() = x[_0];
    t.blockVector() = x[_1];
    v.blockVector() = y[_0];
    s.blockVector() = y[_1];
    return u.scalarProductDofs(v) + t.scalarProductDofs(s);  // already parallel
  }
 
  virtual real_type norm(const X& x) const {
    u.blockVector() = x[_0];
    t.blockVector() = x[_1];
    return std::sqrt(u.scalarProductDofs(u) +
                     t.scalarProductDofs(t));  // already parallel
  }
 
 private:
  mutable BulkDF u, v;
  mutable InterfaceDF t, s;
};
 
template <class M, class X, class Y, class BulkDF, class InterfaceDF>
class ParallelizedMatrixAdapter : public AssembledLinearOperator<M, X, Y> {
 public:
  typedef M matrix_type;
  typedef X domain_type;
  typedef Y range_type;
  typedef typename X::field_type field_type;
 
  explicit ParallelizedMatrixAdapter(const M& A, const BulkDF& uh,
                                     const InterfaceDF& th)
      : _A_(stackobject_to_shared_ptr(A)), u(uh), t(th) {}
 
  void apply(const X& x, Y& y) const override {
    _A_->mv(x, y);
    communicate(y);
  }
 
  void applyscaleadd(field_type alpha, const X& x, Y& y) const override {
    _A_->usmv(alpha, x, y);
    communicate(y);
  }
 
  const M& getmat() const override { return *_A_; }
 
  SolverCategory::Category category() const override {
    return SolverCategory::sequential;
  }
 
 private:
  void communicate(Y& y) const {
    u.blockVector() = y[_0];
    t.blockVector() = y[_1];
    u.communicate();
    t.communicate();
    y[_0] = u.blockVector();
    y[_1] = t.blockVector();
  }
 
  const std::shared_ptr<const M> _A_;
  mutable BulkDF u;
  mutable InterfaceDF t;
};
 
template <class Sch, class ISch, class Sol, class ISol>
class Jacobian {
  using ThisType = Jacobian<Sch, ISch, Sol, ISol>;
 
 public:
  using Scheme = Sch;
  using IScheme = ISch;
  using Solution = Sol;
  using ISolution = ISol;
 
  using AType = typename Scheme::JacobianOperatorType;
  using BType = Dune::Fem::ISTLLinearOperator<ISolution, Solution>;
  using CType = Dune::Fem::ISTLLinearOperator<Solution, ISolution>;
  using DType = typename IScheme::JacobianOperatorType;
 
  // Block Matrix
  using AMatrixBlock = Dune::BCRSMatrix<typename AType::MatrixBlockType>;
  using BMatrixBlock = Dune::BCRSMatrix<typename BType::MatrixBlockType>;
  using CMatrixBlock = Dune::BCRSMatrix<typename CType::MatrixBlockType>;
  using DMatrixBlock = Dune::BCRSMatrix<typename DType::MatrixBlockType>;
 
  using ABRowType = Dune::MultiTypeBlockVector<AMatrixBlock, BMatrixBlock>;
  using CDRowType = Dune::MultiTypeBlockVector<CMatrixBlock, DMatrixBlock>;
  using BlockMatrix = Dune::MultiTypeBlockMatrix<ABRowType, CDRowType>;
 
  // Block Vector
  using AVectorBlock = typename AType::ColumnBlockVectorType;
  using DVectorBlock = typename DType::ColumnBlockVectorType;
  using BlockVector = Dune::MultiTypeBlockVector<AVectorBlock, DVectorBlock>;
 
  typedef typename AType::DomainSpaceType BulkSpaceType;
  typedef typename DType::DomainSpaceType InterfaceSpaceType;
  using ParameterType = Fem::NewtonParameter<
      typename IScheme::LinearInverseOperatorType::SolverParameterType>;
 
  Jacobian(const Scheme& scheme, const IScheme& ischeme, const Solution& uh,
           const ISolution& th, const double eps,
           const std::function<void()>& callback)
      : scheme_(scheme),
        ischeme_(ischeme),
        A_("A", uh.space(), uh.space()),
        B_("B", th.space(), uh.space()),
        C_("C", uh.space(), th.space()),
        D_("D", th.space(), th.space()),
        eps_(eps),
        callback_(callback),
        n_(uh.size()),
        m_(th.size()) {
    x_[_0] = uh.blockVector();
    x_[_1] = th.blockVector();
  }
 
  const Scheme& scheme() const { return scheme_; };
  const BlockMatrix& M() const { return M_; };
 
  void init() {
    if (n_ == 0 || m_ == 0) return;
 
    NeighborInterfaceStencil<InterfaceSpaceType, BulkSpaceType> stencilB(
        B_.domainSpace(), B_.rangeSpace());
    stencilB.setupStencil();
    B_.reserve(stencilB);
 
    InterfaceNeighborStencil<BulkSpaceType, InterfaceSpaceType> stencilC(
        C_.domainSpace(), C_.rangeSpace());
    stencilC.setupStencil();
    C_.reserve(stencilC);
  }
 
  void update(const Solution& uh, const ISolution& th) {
    scheme_.jacobian(uh, A_);
    ischeme_.jacobian(th, D_);
 
    if (n_ > 0 && m_ > 0) {
      B_.clear();
      C_.clear();
      assembleB(scheme_, th, uh);
      assembleC(ischeme_, uh, th);
    }
  }
 
  void solve(const Solution& f, const ISolution& g, Solution& u, ISolution& t) {
    auto setBlock = [this](const auto& block, const auto blockrow,
                           const auto blockcol) {
      this->M_[blockrow][blockcol] = block.exportMatrix();
    };
 
    if (n_ > 0) setBlock(A_, _0, _0);
    if (n_ > 0 && m_ > 0) {
      setBlock(B_, _0, _1);
      setBlock(C_, _1, _0);
    }
    if (m_ > 0) setBlock(D_, _1, _1);
 
    b_[_0] = f.blockVector();
    b_[_1] = g.blockVector();
 
    auto params = std::make_shared<Fem::ISTLSolverParameter>(
        ParameterType(scheme_.parameter()).linear());
    const size_t numIterations = params->preconditionerIteration();
    const double relaxFactor = params->relaxation();
 
    using SolverAdapterType = Fem::ISTLSolverAdapter<-1, BlockVector>;
    using ReductionType = typename SolverAdapterType::ReductionType;
    SolverAdapterType solver(ReductionType(params), params);
    solver.setMaxIterations(params->maxIterations());
 
    ParallelizedMatrixAdapter<BlockMatrix, BlockVector, BlockVector, Solution,
                              ISolution>
        linop(M_, u, t);
    ParallelizedScalarProduct<BlockVector, Solution, ISolution> scp(u, t);
    Dune::Richardson<BlockVector, BlockVector> prec(1.0);
    InverseOperatorResult res;
 
    solver(linop, scp, prec, b_, x_, res);
 
    u.blockVector() = x_[_0];
    t.blockVector() = x_[_1];
    u.communicate();
    t.communicate();
  }
 
 private:
  template <class Scheme, class DomainGridFunction, class RangeGridFunction>
  void assembleB(const Scheme& scheme, const DomainGridFunction& t,
                 const RangeGridFunction& u) {
    typedef InterfaceSpaceType DomainSpaceType;
    typedef BulkSpaceType RangeSpaceType;
 
    typedef Fem::TemporaryLocalMatrix<DomainSpaceType, RangeSpaceType>
        TemporaryLocalMatrixType;
    typedef typename RangeGridFunction::DiscreteFunctionSpaceType
        DiscreteFunctionSpaceType;
 
    const auto& gridPart = t.gridPart();
    const auto& grid = gridPart.grid();
    const auto& mmesh = grid.getMMesh();
 
    auto dFIn = u;
    auto dFOut = u;
 
    Dune::Fem::TemporaryLocalFunction<DiscreteFunctionSpaceType> FmTmpIn(
        u.space());
    Dune::Fem::TemporaryLocalFunction<DiscreteFunctionSpaceType> FmTmpOut(
        u.space());
    Dune::Fem::TemporaryLocalFunction<DiscreteFunctionSpaceType> FpTmpIn(
        u.space());
    Dune::Fem::TemporaryLocalFunction<DiscreteFunctionSpaceType> FpTmpOut(
        u.space());
 
    Dune::Fem::MutableLocalFunction<DomainGridFunction> tLocal(t);
 
    Dune::Fem::ConstLocalFunction<RangeGridFunction> uInside(u);
    Dune::Fem::ConstLocalFunction<RangeGridFunction> uOutside(u);
 
    Dune::Fem::MutableLocalFunction<RangeGridFunction> dFLocalIn(dFIn);
    Dune::Fem::MutableLocalFunction<RangeGridFunction> dFLocalOut(dFOut);
 
    TemporaryLocalMatrixType localMatrixIn(B_.domainSpace(), B_.rangeSpace());
    TemporaryLocalMatrixType localMatrixOut(B_.domainSpace(), B_.rangeSpace());
 
    for (const auto& interface : elements(gridPart, Partitions::all)) {
      tLocal.bind(interface);
      auto& tDof = tLocal.localDofVector();
 
      const auto mmeshIntersection = mmesh.asIntersection(interface);
      const auto inside = u.gridPart().convert(mmeshIntersection.inside());
      const auto intersection =
          convert(u.gridPart(), mmeshIntersection, inside);
 
      const auto& outside = intersection.outside();
 
      FmTmpIn.bind(inside);
      FmTmpOut.bind(outside);
      FpTmpIn.bind(inside);
      FpTmpOut.bind(outside);
 
      uInside.bind(inside);
      uOutside.bind(outside);
 
      dFLocalIn.bind(inside);
      dFLocalOut.bind(outside);
 
      localMatrixIn.init(interface, inside);
      localMatrixOut.init(interface, outside);
 
      localMatrixIn.clear();
      localMatrixOut.clear();
 
      for (std::size_t i = 0; i < tDof.size(); ++i) {
        dFIn.clear();
        dFOut.clear();
 
        FmTmpIn.clear();
        FmTmpOut.clear();
        FpTmpIn.clear();
        FpTmpOut.clear();
 
        double h = std::max(tDof[i] * eps_, eps_);
        tDof[i] -= h;
        callback_();
        scheme.fullOperator().impl().addSkeletonIntegral(
            intersection, uInside, uOutside, FmTmpIn, FmTmpOut);
        tDof[i] += h;
 
        dFIn.addLocalDofs(inside, FmTmpIn.localDofVector());
        dFOut.addLocalDofs(outside, FmTmpOut.localDofVector());
 
        dFIn *= -1.;
        dFOut *= -1.;
 
        tDof[i] += h;
        callback_();
        scheme.fullOperator().impl().addSkeletonIntegral(
            intersection, uInside, uOutside, FpTmpIn, FpTmpOut);
        tDof[i] -= h;
 
        dFIn.addLocalDofs(inside, FpTmpIn.localDofVector());
        dFOut.addLocalDofs(outside, FpTmpOut.localDofVector());
 
        dFIn /= 2 * h;
        dFOut /= 2 * h;
 
        for (std::size_t j = 0; j < dFLocalIn.localDofVector().size(); ++j)
          localMatrixIn.set(j, i, dFLocalIn[j]);
 
        for (std::size_t j = 0; j < dFLocalOut.localDofVector().size(); ++j)
          localMatrixOut.set(j, i, dFLocalOut[j]);
      }
 
      B_.addLocalMatrix(interface, inside, localMatrixIn);
      B_.addLocalMatrix(interface, outside, localMatrixOut);
    }
    B_.compress();
  }
 
  template <class Scheme, class DomainGridFunction, class RangeGridFunction>
  void assembleC(const Scheme& ischeme, const DomainGridFunction& u,
                 const RangeGridFunction& t) {
    typedef BulkSpaceType DomainSpaceType;
    typedef InterfaceSpaceType RangeSpaceType;
 
    typedef Fem::TemporaryLocalMatrix<DomainSpaceType, RangeSpaceType>
        TemporaryLocalMatrixType;
    typedef typename RangeGridFunction::DiscreteFunctionSpaceType
        DiscreteFunctionSpaceType;
 
    const auto& gridPart = u.gridPart();
    const auto& grid = gridPart.grid();
 
    auto dG = t;
 
    Dune::Fem::TemporaryLocalFunction<DiscreteFunctionSpaceType> GmTmp(
        t.space());
    Dune::Fem::TemporaryLocalFunction<DiscreteFunctionSpaceType> GpTmp(
        t.space());
 
    Dune::Fem::MutableLocalFunction<DomainGridFunction> uLocal(u);
 
    Dune::Fem::ConstLocalFunction<RangeGridFunction> tInterface(t);
    Dune::Fem::MutableLocalFunction<RangeGridFunction> dGLocal(dG);
 
    TemporaryLocalMatrixType localMatrix(C_.domainSpace(), C_.rangeSpace());
 
    for (const auto& element : elements(gridPart, Partitions::all)) {
      for (const auto& intersection : intersections(gridPart, element)) {
        if (grid.isInterface(intersection)) {
          const auto& interface = grid.asInterfaceEntity(intersection);
 
          uLocal.bind(element);
          auto& uDof = uLocal.localDofVector();
 
          GmTmp.bind(interface);
          GpTmp.bind(interface);
          tInterface.bind(interface);
 
          dGLocal.bind(interface);
          localMatrix.init(element, interface);
          localMatrix.clear();
 
          for (std::size_t i = 0; i < uDof.size(); ++i) {
            dG.clear();
 
            GmTmp.clear();
            GpTmp.clear();
 
            double h = std::max(uDof[i] * eps_, eps_);
            uDof[i] -= h;
            callback_();
            ischeme.fullOperator().impl().addInteriorIntegral(tInterface,
                                                              GmTmp);
            uDof[i] += h;
 
            dG.addLocalDofs(interface, GmTmp.localDofVector());
            dG *= -1.;
 
            uDof[i] += h;
            callback_();
            ischeme.fullOperator().impl().addInteriorIntegral(tInterface,
                                                              GpTmp);
            uDof[i] -= h;
 
            dG.addLocalDofs(interface, GpTmp.localDofVector());
            dG /= 2 * h;
 
            for (std::size_t j = 0; j < dGLocal.localDofVector().size(); ++j)
              localMatrix.set(j, i, dGLocal[j]);
          }
 
          C_.addLocalMatrix(element, interface, localMatrix);
        }
      }
    }
    C_.compress();
  }
 
  const Scheme& scheme_;
  const IScheme& ischeme_;
  AType A_;
  BType B_;
  CType C_;
  DType D_;
  BlockMatrix M_;
  BlockVector x_, b_;
  std::size_t n_, m_;
  const double eps_;
  const std::function<void()> callback_;
};
 
template <class Jacobian, class... options>
inline static auto registerJacobian(
    pybind11::handle scope, pybind11::class_<Jacobian, options...> cls) {
  using Solution = typename Jacobian::Solution;
  using ISolution = typename Jacobian::ISolution;
 
  cls.def(
      "init", [](Jacobian& self) { self.init(); },
      R"doc(
            Initialize the mixed-dimensional jacobian.
          )doc");
 
  cls.def(
      "update",
      [](Jacobian& self, const Solution& uh, const ISolution& th) {
        self.update(uh, th);
      },
      R"doc(
            Update the mixed-dimensional jacobian.
          )doc");
 
  cls.def(
      "solve",
      [](Jacobian& self, const Solution& f, const ISolution& g, Solution& ux,
         ISolution& tx) { self.solve(f, g, ux, tx); },
      R"doc(
            Solve the mixed-dimensional jacobian.
          )doc");
}
 
}  // namespace MMesh
 
}  // namespace Python
 
}  // namespace Dune
 
#endif