Running Allen-Cahn example

src/CMakeLists.txt

 add_executable("binary-allen-cahn" binary-allen-cahn.cc)
 target_link_dune_default_libraries("binary-allen-cahn")
+
+dune_add_copy_dependency(binary-allen-cahn "binary-allen-cahn.parset")

src/binary-allen-cahn.cc

@@ -5,8 +5,293 @@
 # include "config.h"
 #endif
 #include <iostream>
+
+
 #include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
 #include <dune/common/exceptions.hh> // We use exceptions
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/stringutility.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
+#include <dune/fufem/assemblers/dunefunctionsfunctionalassembler.hh>
+#include <dune/fufem/assemblers/istlbackend.hh>
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+
+#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/dunefunctionsl2functionalassembler.hh>
+
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/transferoperators/truncateddensemgtransfer.hh>
+
+#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
+#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
+#include <dune/tnnmg/projections/obstacledefectprojection.hh>
+
+
+
+
+
+template<class Grid, class GridView, class Basis>
+class BinaryAllenCahnProblem
+{
+public:
+  static const unsigned int dim = Grid::dimension;
+
+  using Field = double;
+  using Vector = Dune::BlockVector<Dune::FieldVector<Field,1>>;
+  using BitVector = Dune::Solvers::DefaultBitVector_t<Vector>;
+
+  using LocalMatrix =  Dune::FieldMatrix<Field,1,1>;
+  using Matrix = Dune::BCRSMatrix<LocalMatrix>;
+
+  using GlobalCoordinate = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+  using Function = std::function<Field(GlobalCoordinate)>;
+  using GridFunction = Dune::Functions::GridViewFunction<Field(GlobalCoordinate), GridView>;
+
+  BinaryAllenCahnProblem(Grid& grid, const GridView& gridView, const Basis& basis, const Dune::ParameterTree& parameterSet) :
+    grid_(grid),
+    gridView_(gridView),
+    basis_(basis),
+    parameterSet_(parameterSet)
+  {
+    epsilon_ = parameterSet_.get<Field>("epsilon");
+    timeStepSize_ = parameterSet_.get<Field>("timestepsize");
+    timeSteps_ = parameterSet_.get<std::size_t>("timesteps");
+
+    solverMaxIterations_ = parameterSet_.get<std::size_t>("solver.maxIterations");
+    solverTolerance_ = parameterSet_.get<Field>("solver.tolerance");
+
+    outputFilenameFormat_ = parameterSet_.get<std::string>("outputFilenameFormat");
+
+    u_.resize(basis_.size());
+    u_ = 0.0;
+
+    Dune::Solvers::resizeInitialize(dirichletNodes_, u_, false);
+
+    u0Function_ = [](auto x) {
+      x -= .5;
+      return (x.infinity_norm() < .2)*2.0 - 1.0;
+    };
+  }
+
+
+  void timeLoop()
+  {
+    using VectorBackend = typename Dune::Functions::HierarchicVectorWrapper<Vector, Field>;
+    Dune::Functions::interpolate(basis_, VectorBackend(u_), u0Function_);
+    writeSolution(0);
+    for(std::size_t t=1; t<=timeSteps_; ++t)
+    {
+      std::cout << "Timestep:" << t << std::endl;
+      if (t==1)
+        assembleMatrix();
+
+      assembleRHS(t==1);
+
+      solveObstacleProblem();
+
+      writeSolution(t);
+
+      uOld_ = u_;
+    }
+  }
+
+  void assembleRHS(bool firstTimeStep)
+  {
+    std::cout << "Assembling right hand side" << std::endl;
+    if (firstTimeStep)
+    {
+      using Assembler = Dune::Fufem::DuneFunctionsFunctionalAssembler<Basis>;
+      using FiniteElement = std::decay_t<decltype(basis_.localView().tree().finiteElement())>;
+      using VectorBackend = typename Dune::Functions::HierarchicVectorWrapper<Vector, Field>;
+      
+      auto assembler = Assembler{basis_};
+
+      auto uOldGridFunction = Dune::Functions::makeAnalyticGridViewFunction(u0Function_, basis_.gridView());
+      auto localAssembler = Dune::Fufem::DuneFunctionsL2FunctionalAssembler<Grid, FiniteElement, GridFunction>(uOldGridFunction, QuadratureRuleKey(dim,1));
+
+      assembler.assembleBulk(VectorBackend(rhs_), localAssembler);
+      rhs_ *= 1.0/timeStepSize_;
+    }
+    else
+    {
+      massMatrix_.mv(uOld_, rhs_);
+      rhs_ *= 1.0/timeStepSize_;
+    }
+  }
+
+  void assembleMatrix()
+  {
+    std::cout << "Assembling matrix" << std::endl;
+    using MatrixBuilder = Dune::Fufem::MatrixBuilder<Matrix>;
+    using Assembler = Dune::Fufem::DuneFunctionsOperatorAssembler<Basis, Basis>;
+    using FiniteElement = std::decay_t<decltype(basis_.localView().tree().finiteElement())>;
+
+
+    auto assembler = Assembler{basis_, basis_};
+    {
+      auto matrixBackend = Dune::Fufem::istlMatrixBackend(stiffnessMatrix_);
+      auto localAssembler = LaplaceAssembler<Grid, FiniteElement, FiniteElement>();
+      assembler.assembleBulk(matrixBackend, localAssembler);
+    }
+    {
+      auto matrixBackend = Dune::Fufem::istlMatrixBackend(massMatrix_);
+      auto localAssembler = MassAssembler<Grid, FiniteElement, FiniteElement>();
+      assembler.assembleBulk(matrixBackend, localAssembler);
+    }
+
+    systemMatrix_ = stiffnessMatrix_;
+    systemMatrix_ *= epsilon_;
+    systemMatrix_.axpy(epsilon_/timeStepSize_ - 1.0/epsilon_, massMatrix_);
+  }
+
+  void solveObstacleProblem()
+  {
+    auto linearBaseSolverStep = TruncatedBlockGSStep<Matrix,Vector>();
+    auto baseEnergyNorm = EnergyNorm<Matrix,Vector>(linearBaseSolverStep);
+    auto linearBaseSolver = Dune::Solvers::LoopSolver<Vector>(&linearBaseSolverStep,
+                                                                               100,
+                                                                               1e-8,
+                                                                               &baseEnergyNorm,
+                                                                               Solver::QUIET);
+    auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+
+    // setup transfer operators
+    auto transferOperators = std::vector<TruncatedDenseMGTransfer<Vector>>{grid_.maxLevel()};
+    auto transferOperatorAssembler =  TransferOperatorAssembler<Grid>(grid_);
+    transferOperatorAssembler.assembleOperatorHierarchy(transferOperators);
+
+    // setup multigrid step
+    auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+    linearMultigridStep->setMGType(1, 3, 3);
+    linearMultigridStep->setBaseSolver(linearBaseSolver);
+    linearMultigridStep->setSmoother(&smoother);
+    linearMultigridStep->setTransferOperators(transferOperators);
+    
+    // setup upper and lower obstacle
+    Vector lower(rhs_.size());
+    Vector upper(rhs_.size());
+    lower = -1.0;
+    upper = 1.0;
+
+    using Functional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix, Vector, Vector, Vector, double>;
+    auto J = Functional(systemMatrix_, rhs_, lower, upper);
+
+    auto localSolver = gaussSeidelLocalSolver(Dune::TNNMG::ScalarObstacleSolver());
+
+    using Linearization = Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<Functional, BitVector>;
+    using DefectProjection = Dune::TNNMG::ObstacleDefectProjection;
+    using LineSearchSolver = Dune::TNNMG::ScalarObstacleSolver;
+
+    using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, decltype(localSolver), BitVector>;
+    auto nonlinearSmoother = std::make_shared<NonlinearSmoother>(J, u_, localSolver);
+
+    using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+
+    auto step = Step(J, u_, nonlinearSmoother, linearMultigridStep, 1, DefectProjection(), LineSearchSolver());
+
+    auto norm = EnergyNorm<Matrix, Vector>(systemMatrix_);
+    auto solver = Dune::Solvers::LoopSolver<Vector>(&step, solverMaxIterations_, solverTolerance_, &norm, Solver::FULL);
+
+
+    step.setIgnore(dirichletNodes_);
+    step.setPreSmoothingSteps(3);
+
+    solver.addCriterion(
+              [&](){
+                  return Dune::formatString("   % 12d", step.linearization().truncated().count());
+              },
+              "   truncated   ");
+
+//    double initialEnergy = J(u_);
+//    solver.addCriterion(
+//              [&](){
+//                  static double oldEnergy=initialEnergy;
+//                  double currentEnergy = J(u_);
+//                  double decrease = currentEnergy - oldEnergy;
+//                  oldEnergy = currentEnergy;
+//                  return Dune::formatString("   % 12.5e", decrease);
+//              },
+//              "   decrease    ");
+
+    solver.addCriterion(
+              [&](){
+                  return Dune::formatString("   % 12.5e", step.lastDampingFactor());
+              },
+              "   step size     ");
+
+    solver.addCriterion(
+              [&](){
+                  return Dune::formatString("   % 12.5e", J(u_));
+              },
+              "   energy      ");
+
+    solver.addCriterion(Dune::Solvers::correctionNormCriterion(step, norm, 1e-10));
+
+    solver.preprocess();
+    solver.solve();
+  }
+
+  template<class... Arg>
+  void writeSolution(Arg&&... args) const
+  {
+    using VectorBackend = typename Dune::Functions::HierarchicVectorWrapper<const Vector, Field>;
+    
+    //  Make a discrete function from the FE basis and the coefficient vector
+    auto uFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(basis_, VectorBackend(u_));
+
+    auto vtkWriter = Dune::VTKWriter<GridView>(gridView_);
+    vtkWriter.addVertexData(uFunction, Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
+    vtkWriter.write(Dune::formatString(outputFilenameFormat_, args...));
+  }
+
+
+private:
+  Grid& grid_;
+  GridView gridView_;
+  const Basis& basis_;
+  const Dune::ParameterTree& parameterSet_;
+
+  Field epsilon_;
+  Field timeStepSize_;
+  std::size_t timeSteps_;
+  std::string outputFilenameFormat_;
+
+  std::size_t solverMaxIterations_;
+  Field solverTolerance_;
+
+  Matrix stiffnessMatrix_;
+  Matrix massMatrix_;
+  Matrix systemMatrix_;
+  Vector u_;
+  Vector uOld_;
+  Vector rhs_;
+
+  BitVector dirichletNodes_;
+
+  Function u0Function_;
+};
+
+
 
 int main(int argc, char** argv)
 {
@@ -19,6 +304,47 @@ int main(int argc, char** argv)
     else
       std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
         <<" processes!"<<std::endl;
+
+    // Read parameter file
+    Dune::ParameterTree parameterSet;
+    Dune::ParameterTreeParser::readINITree("binary-allen-cahn.parset", parameterSet);
+
+    const int dim = 2;
+    const int refine = parameterSet.get<std::size_t>("grid.global_refinements");
+
+    // Build a test grid
+    using Grid = Dune::YaspGrid<dim>;
+
+    Dune::FieldVector<double,dim> h(1);
+    std::array<int,dim> n;
+    n.fill(1);
+    n[0] = 2;
+
+    Grid grid(h,n);
+    grid.globalRefine(refine);
+
+    auto gridView = grid.leafGridView();
+
+    using GridView = decltype(gridView);
+
+    using namespace Dune::Functions::BasisBuilder;
+
+    auto basis = makeBasis(gridView, lagrange<1>());
+    std::cout << refine << " , " << basis.size() << std::endl;
+
+    using Basis = decltype(basis);
+
+
+
+    auto problem = BinaryAllenCahnProblem<Grid, GridView, Basis>(grid, gridView, basis, parameterSet.sub("allen-cahn"));
+
+    problem.timeLoop();
+
+
+
+
+
+
     return 0;
   }
   catch (Dune::Exception &e){

src/binary-allen-cahn.parset

+[grid]
+global_refinements = 7
+
+[allen-cahn]
+epsilon = 1e-1
+timestepsize = 5e-3
+timesteps = 100
+outputFilenameFormat = allen-cahn-%05d
+
+[allen-cahn.solver]
+maxIterations = 100
+tolerance = 1e-10
+