dune/solvers/test/obstacletnnmgtest.cc

@@ -76,7 +76,7 @@ void solveObstacleProblemByTNNMG(const GridType& grid, const MatrixType& mat, Ve
     tnnmgStep.preprocess();
     tnnmgStep.nestedIteration();
 
-    // cretae loop solver
+    // create loop solver
     Solver solver(tnnmgStep, maxIterations, tolerance, norm, Solver::FULL);
 
     // solve problem

dune/solvers/test/proximalnewtonsolvertest.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+
+#include <config.h>
+
+#include <cmath>
+#include <iostream>
+#include <sstream>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/norms/twonorm.hh>
+#include <dune/solvers/solvers/cholmodsolver.hh>
+#include <dune/solvers/solvers/proximalnewtonsolver.hh>
+
+#include "common.hh"
+
+using namespace Dune;
+
+// grid dimension
+const int dim = 3;
+
+// f(x) = 1/4 * ||x||^4
+template<class Matrix, class Vector>
+struct F
+{
+  void assembleGradientAndHessian( const Vector& x, Vector& grad, Matrix& hess ) const
+  {
+    auto norm2 = x.two_norm2();
+
+    grad = x;
+    grad *= norm2;
+
+    for(size_t i=0; i<x.size(); i++)
+      for(size_t j=0; j<x.size(); j++)
+        hess[i][j] = (i==j)*norm2 + 2*x[i]*x[j];
+  }
+
+  double computeEnergy( const Vector& x ) const
+  {
+    return  0.25* x.two_norm2() * x.two_norm2();
+  }
+};
+
+
+// g(x) = ( -1, -1, -1, ... , -1)^T * x
+template<class Vector>
+struct G
+{
+  double operator()( const Vector& x ) const
+  {
+    auto realX = x;
+    realX += offset_;
+
+    double sum = 0.0;
+    for ( size_t i=0; i<realX.size(); i++ )
+      sum -= realX[i];
+
+    return sum;
+  }
+
+  void updateOffset( const Vector& offset )
+  {
+    offset_ = offset;
+  }
+
+  Vector offset_;
+};
+
+template<class Matrix, class Vector>
+struct SecondOrdersolver
+{
+  using MatrixType = Matrix;
+
+  template<class BitVector>
+  void minimize( const Matrix& hess, const Vector& grad, const G<Vector>& g, double reg, Vector& dx, const BitVector& /*ignore*/) const
+  {
+    // add reg
+    auto HH = hess;
+    for ( size_t i=0; i<dx.size(); i++ )
+      HH[i][i] += reg;
+
+    // add g to grad
+    auto gg = grad;
+    gg *= -1.0;
+    for ( size_t i=0; i<dx.size(); i++ )
+      gg[i] += 1.0;
+
+
+    Solvers::CholmodSolver cholmod( HH, dx, gg );
+    cholmod.solve();
+  }
+};
+
+
+
+int main (int argc, char *argv[])
+{
+  // initialize MPI, finalize is done automatically on exit
+  [[maybe_unused]] MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+  const int size = 10;
+
+  using Vector = FieldVector<double,size>;
+  using Matrix = FieldMatrix<double,size>;
+  using BitVector = std::bitset<size>;
+
+  // this is the analytical solution of the problem
+  Vector exactSol, zeroVector;
+  exactSol = std::pow( size, -1.0/3.0 );
+
+  zeroVector = 0.0;
+
+  F<Matrix,Vector> f;
+  G<Vector> g;
+  g.updateOffset(zeroVector);
+  SecondOrdersolver<Matrix,Vector> sos;
+
+
+  // choose some random initial vector
+  Vector sol;
+  for ( size_t i=0; i<sol.size(); i++ )
+    sol[i] = i;
+
+  TwoNorm<Vector> norm;
+
+  // create the solver
+  Solvers::ProximalNewton::SimpleRegUpdater regUpdater;
+  double reg = 42.0;
+  Solvers::ProximalNewtonSolver proximalNewtonSolver( f, g, sos, sol, norm, regUpdater, reg, 100, 1e-14, Solver::FULL);
+
+  // set some empty ignore field
+  BitVector ignore;
+  proximalNewtonSolver.setIgnore( ignore );
+
+  // go!
+  proximalNewtonSolver.solve();
+
+  // compute diff the exact solution
+  sol -= exactSol;
+
+  std::cout << "Difference to exact solution = " << sol.two_norm() << std::endl;
+
+  bool passed = sol.two_norm() < 1e-15;
+
+  return passed ? 0 : 1;
+}

dune/solvers/test/umfpacksolvertest.cc

@@ -9,6 +9,8 @@
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/tuplevector.hh>
+#include <dune/common/version.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
 
@@ -62,6 +64,99 @@ struct UMFPackSolverTestSuite
 
 
 
+/**  \brief test for the UMFPackSolver class with multitype matrix */
+template <size_t blocksize>
+struct UMFPackMultiTypeSolverTestSuite
+{
+  template <class GridType>
+  bool check(const GridType& grid)
+  {
+    double tol = 1e-11;
+
+    bool passed = true;
+
+    using Problem =
+    SymmetricSampleProblem<blocksize, typename GridType::LeafGridView>;
+    Problem p(grid.leafGridView());
+
+    // create a multitype problem
+    using BlockRow = MultiTypeBlockVector<typename Problem::Matrix,typename Problem::Matrix>;
+    using Matrix = MultiTypeBlockMatrix<BlockRow,BlockRow>;
+
+    using Vector = MultiTypeBlockVector<typename Problem::Vector,typename Problem::Vector>;
+
+    Matrix A;
+    Vector u, rhs;
+
+    using namespace Dune::Indices;
+    // initialize each block
+    A[_0][_0] = p.A;
+    A[_0][_1] = p.A;
+    A[_1][_0] = p.A;
+    A[_1][_1] = p.A;
+
+    u[_0] = p.u;
+    u[_1] = p.u;
+    rhs[_0] = p.rhs;
+    rhs[_1] = p.rhs;
+
+    using IgnoreBlock = decltype(p.ignore);
+    using Ignore = TupleVector<IgnoreBlock,IgnoreBlock>;
+
+    Ignore ignore;
+    ignore[_0] = p.ignore;
+    ignore[_1] = p.ignore;
+
+    typedef Solvers::UMFPackSolver<Matrix,Vector> Solver;
+    Solver solver(A,u,rhs);
+
+    solver.setIgnore(ignore);
+
+    // solve blockdiagonal
+    A[_0][_1] *= 0.0;
+    A[_1][_0] *= 0.0;
+    solver.solve();
+
+    // the solution is ( u_ex, u_ex )
+    if (p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) > tol)
+    {
+      std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) << std::endl;
+      passed = false;
+    }
+
+    // solve a problem more generic problem without ignore field
+
+    A[_0][_1] = p.A;
+    A[_1][_0] = p.A;
+
+    // make A regular
+    A[_0][_0] *= 2.0;
+
+    u[_0] = p.u;
+    u[_1] = p.u;
+
+    // solve problem
+    Solver solver2(A,u,rhs);
+    solver2.preprocess();
+    solver2.solve();
+
+    // check residual
+    auto resisual = rhs;
+    A.mmv(u,resisual);
+
+    if (p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) > tol)
+    {
+      std::cout << "The UMFPackSolver did not produce a satisfactory result. ||residual||=" << p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) << std::endl;
+      passed = false;
+    }
+
+
+    return passed;
+  }
+};
+
+
+
 int main(int argc, char** argv)
 {
     Dune::MPIHelper::instance(argc, argv);
@@ -71,9 +166,15 @@ int main(int argc, char** argv)
     UMFPackSolverTestSuite<0> testsuite0;
     UMFPackSolverTestSuite<1> testsuite1;
     UMFPackSolverTestSuite<2> testsuite2;
+
+    UMFPackMultiTypeSolverTestSuite<3> testsuite3;
+
     passed = checkWithStandardGrids(testsuite0);
     passed = checkWithStandardGrids(testsuite1);
     passed = checkWithStandardGrids(testsuite2);
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 10)
+    passed = checkWithStandardGrids(testsuite3);
+#endif
 
     return passed ? 0 : 1;
 }

dune/solvers/transferoperators/genericmultigridtransfer.hh

@@ -10,7 +10,7 @@
 
 #include <dune/geometry/referenceelements.hh>
 
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
 #include <dune/matrix-vector/axpy.hh>
 #include <dune/matrix-vector/transformmatrix.hh>
@@ -81,7 +81,7 @@ public:
         int cols = grid.size(cL, dim);
 
         // A factory for the shape functions
-        typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
+        typedef typename Dune::LagrangeLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
         typedef typename P1FECache::FiniteElementType FEType;
         P1FECache p1FECache;
 
@@ -228,7 +228,7 @@ public:
         int cols = gvCoarse.size(dim);
 
         // A factory for the shape functions
-        typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
+        typedef typename Dune::LagrangeLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
         typedef typename P1FECache::FiniteElementType FEType;
         P1FECache p1FECache;