Implementation of a simple Newton iteration scheme

amdis/CMakeLists.txt

@@ -80,6 +80,7 @@ add_subdirectory("interpolators")
 add_subdirectory("io")
 add_subdirectory("linearalgebra")
 add_subdirectory("localoperators")
+add_subdirectory("nonlin")
 add_subdirectory("operations")
 add_subdirectory("typetree")
 add_subdirectory("utility")

amdis/DirichletBC.inc.hpp

@@ -44,8 +44,10 @@ apply(Mat& matrix, Sol& solution, Rhs& rhs)
 {
   assert(to_string(row_) == to_string(col_));
   auto interpolator = InterpolatorFactory<tag::assign>::create(*basis_, col_);
-  interpolator(solution, valueGridFct_, dirichletNodes_);
-  dirichletBC(matrix, solution, rhs, dirichletNodes_);
+  Sol boundarySolution{solution};
+  interpolator(boundarySolution, valueGridFct_, dirichletNodes_);
+  dirichletBC(matrix, boundarySolution, rhs, dirichletNodes_);
+  solution = std::move(boundarySolution);
 }
 
 } // end namespace AMDiS

amdis/nonlin/CMakeLists.txt

+install(FILES
+    NewtonIteration.hpp
+    NewtonIteration.inc.hpp
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/nonlin)

amdis/nonlin/NewtonIteration.hpp

+#pragma once
+
+#include <memory>
+#include <string>
+
+#include <amdis/AdaptInfo.hpp>
+#include <amdis/Flag.hpp>
+#include <amdis/Initfile.hpp>
+#include <amdis/ProblemIterationInterface.hpp>
+#include <amdis/ProblemStatBase.hpp>
+
+namespace AMDiS
+{
+  template <class ProblemType>
+  class NewtonIteration
+      : public ProblemIterationInterface
+  {
+    using Self    = NewtonIteration;
+    using Problem = ProblemType;
+
+    using SolutionVector = typename Problem::SolutionVector;
+
+  public:
+    /// Constructs a nonlinear iteration scheme
+    NewtonIteration (std::string name, Problem& prob)
+      : name_(std::move(name))
+      , prob_(&prob)
+    {
+      stepSolution_.reset(new SolutionVector(*prob_->solutionVector()));
+
+      Parameters::get(name_ + "->build cycle", buildCycle_);
+      Parameters::get(name_ + "->norm", norm_);
+    }
+
+    /// Implementation of \ref ProblemIterationInterface::beginIteration
+    void beginIteration(AdaptInfo& /*adaptInfo*/) override { /* do nothing */ }
+
+    /// Implementation of \ref ProblemIterationInterface::oneIteration
+    Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo) override;
+
+    /// Implementation of \ref ProblemIterationInterface::endIteration
+    void endIteration(AdaptInfo& adaptInfo) override
+    {
+      msg("{:>5} | {:>12} | {:>8} | {:>8} ", "iter.", "error", "red.", "tol");
+      if (errOld_ <= 0)
+        msg("{:5d} | {:12.5e} | {:->8} | {:8.2} ",
+          adaptInfo.spaceIteration()+1, err_, "-",adaptInfo.spaceTolerance(""));
+      else
+        msg("{:5d} | {:12.5e} | {:8.2e} | {:8.2} ",
+           adaptInfo.spaceIteration()+1, err_, err_/errOld_,adaptInfo.spaceTolerance(""));
+
+      errOld_ = err_;
+    }
+
+    /// Returns \ref problemStat.
+    ProblemStatBase& problem([[maybe_unused]] int n = 0) override
+    {
+      assert(n == 0); return *prob_;
+    }
+
+    /// Implementation of \ref ProblemIterationInterface::numProblems
+    int numProblems() const override { return 1; }
+
+    /// Returns the problem with the given name.
+    ProblemStatBase& problem([[maybe_unused]] std::string const& name) override
+    {
+      assert(prob_->name() == name);
+      return *prob_;
+    }
+
+    /// Returns the name of the problem.
+    std::string const& name() const override { return name_; }
+
+
+    /// Returns const-ref of \ref stepSolution_.
+    std::shared_ptr<SolutionVector const> stepSolutionVector() const
+    {
+      return stepSolution_;
+    }
+
+    /// Return a const view to a stepSolution component
+    /**
+     * \tparam Range  The range type return by evaluating the view in coordinates. If not specified,
+     *                it is automatically selected using \ref RangeType_t template.
+     **/
+    template <class Range = void, class... Indices>
+    auto stepSolution(Indices... ii) const
+    {
+      return valueOf<Range>(*stepSolution_, ii...);
+    }
+
+  protected:
+    /// Name of the newton problem, used to access initfile parameters
+    std::string name_;
+
+    /// Problem containing the Jacobian operators and objective function
+    Problem* prob_;
+
+    /// Solution of the update step
+    std::shared_ptr<SolutionVector> stepSolution_;
+
+    /// Build matrix every i'th iteration
+    /**
+     * 0 ........ build matrix only once
+     * i >= 1 ... rebuild matrix in i'th solver iteration
+     * [default: 1]
+     **/
+    int buildCycle_ = 1;
+
+    /// Type of norm to use for estimating the error
+    /**
+     * 1... L2-norm
+     * 2... H1-norm
+     * [default: 1]
+     **/
+    int norm_ = 1;
+
+  private:
+    double err_ = 0.0;
+    double errOld_ = -1.0;
+  };
+
+} // end namespace AMDiS
+
+#include <amdis/nonlin/NewtonIteration.inc.hpp>

amdis/nonlin/NewtonIteration.inc.hpp

+#pragma once
+
+#include <amdis/Integrate.hpp>
+#include <amdis/operations/Arithmetic.hpp>
+
+namespace AMDiS {
+
+template <class Problem>
+Flag NewtonIteration<Problem>::oneIteration(AdaptInfo& adaptInfo, Flag toDo)
+{
+  if (!toDo.isSet(SOLVE))
+    return Flag(1);
+
+  int iter = adaptInfo.spaceIteration()+1;
+  if (iter <= 0 || (buildCycle_ > 0 && (iter-1) % buildCycle_ == 0)) {
+    prob_->buildAfterAdapt(adaptInfo, 0, true, true);  // assemble DF(u0) and F(u)
+    prob_->solver()->init(prob_->systemMatrix()->impl());
+  }
+  else {
+    prob_->buildAfterAdapt(adaptInfo, 0, false, true); // assemble F(u)
+  }
+
+  // Initial guess is zero
+  stepSolution_->resizeZero();
+
+  // solve the linear system
+  Dune::InverseOperatorResult stat;
+  prob_->solver()->apply(stepSolution_->impl(), prob_->rhsVector()->impl(), stat);
+
+  // u = u + d
+  Recursive::transform(prob_->solutionVector()->impl(),
+    Operation::Plus{}, prob_->solutionVector()->impl(), stepSolution_->impl());
+
+  // L2-norm
+  err_ = integrate(sqr(this->stepSolution()), prob_->gridView());
+  if (norm_ > 1) // H1-norm
+    err_ += integrate(unary_dot(gradientOf(this->stepSolution())), prob_->gridView());
+  err_ = std::sqrt(err_);
+
+  adaptInfo.setEstSum(err_, "");
+  return Flag(1);
+}
+
+} // end namespace AMDiS

examples/CMakeLists.txt

@@ -18,6 +18,7 @@ add_amdis_executable(NAME stokes3.2d              SOURCES stokes3.cc
 add_amdis_executable(NAME traversal               SOURCES traversal.cc)
 add_amdis_executable(NAME treecontainer           SOURCES treecontainer.cc)
 add_amdis_executable(NAME vecellipt.2d            SOURCES vecellipt.cc            DIM 2 DOW 2)
+add_amdis_executable(NAME nonlin-poisson.2d       SOURCES nonlin-poisson.cc       DIM 2 DOW 2)
 
 add_dependencies(examples
   boundary.2d
@@ -35,6 +36,7 @@ add_dependencies(examples
   traversal
   treecontainer
   vecellipt.2d
+  nonlin-poisson.2d
 )
 
 if (ALBERTA_FOUND)

examples/init/nonlin-poisson.dat.2d

+mesh->global refinements: 4
+mesh->num cells: 8 8
+mesh->overlap: 0
+
+poisson->mesh: mesh
+poisson->symmetry: spd
+
+% ISTL backend parameters
+% =======================
+poisson->solver: direct
+poisson->solver->info: 0
+poisson->solver->max iteration: 10000
+poisson->solver->relative tolerance: 1.e-10
+poisson->solver->precon: ilu
+
+poisson->output->format: vtk
+poisson->output->filename: poisson.2d
+poisson->output->name: u
+poisson->output->output directory: ./output
+
+adapt->max iteration: 20
+adapt[]->tolerance: 1.e-10
+
+newton->build cycle: 2
+newton->norm: 1
\ No newline at end of file

examples/nonlin-poisson.cc

+#include <amdis/AMDiS.hpp>
+#include <amdis/AdaptStationary.hpp>
+#include <amdis/LocalOperators.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/common/Literals.hpp>
+#include <amdis/nonlin/NewtonIteration.hpp>
+
+#include <dune/grid/yaspgrid.hh>
+
+using namespace AMDiS;
+using namespace Dune::Functions::BasisFactory;
+
+int main(int argc, char** argv)
+{
+  Environment env(argc, argv);
+
+  auto f = [](auto const& x) {
+    double r2 = dot(x,x);
+    double ux = std::exp(-10.0 * r2);
+    return -(400.0 * r2 - 20.0 * x.size()) * ux;
+  };
+
+  using Grid = Dune::YaspGrid<GRIDDIM>;
+  auto gridPtr = MeshCreator<Grid>{"mesh"}.create();
+
+  ProblemStat prob("poisson", *gridPtr, lagrange<2>());
+  prob.initialize(INIT_ALL);
+
+  NewtonIteration newton{"newton", prob};
+
+  // <grad(du),grad(theta)> + <4*u_k^3 du,theta> = -<grad(u_k),grad(theta)> + <-u_k^4 + f(x),theta>
+  // du = -u_k on \partial\Omega
+
+  // define the Jacobian
+  prob.addMatrixOperator(sot(1.0));
+  prob.addMatrixOperator(zot(4*pow<3>(prob.solution())));
+
+  // define objective function
+  prob.addVectorOperator(makeOperator(tag::gradtest{}, -gradientOf(prob.solution())));
+  prob.addVectorOperator(zot(-pow<4>(prob.solution())));
+  prob.addVectorOperator(zot(f, 6));
+
+  // set boundary condition
+  prob.addDirichletBC([](auto){return true;}, -prob.solution());
+
+  AdaptInfo adaptInfo{"adapt"};
+  AdaptStationary adaptStat{"adapt", newton, adaptInfo};
+  adaptStat.adapt();
+
+  return 0;
+}