Solve nonlinear problem

python/dune/perftool/pdelab/driver.py

@@ -50,6 +50,18 @@ def set_form(f):
     _form = f
 
 
+def is_linear(form):
+    '''Test if form is linear in test function'''
+    from ufl import derivative
+    from ufl.algorithms import expand_derivatives
+    jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
+    for coeff in jacform.coefficients():
+        if 0 == coeff.count():
+            return False
+
+    return True
+
+
 def isLagrange(fem):
     return fem._short_name is 'CG'
 
@@ -802,12 +814,27 @@ def typedef_stationarylinearproblemsolver(name):
     return "typedef Dune::PDELab::StationaryLinearProblemSolver<{}, {}, {}> {};".format(gotype, lstype, xtype, name)
 
 
+@preamble
+def typedef_stationarynonlinearproblemsolver(name):
+    include_file("dune/pdelab/newton/newton.hh", filetag="driver")
+    gotype = type_gridoperator()
+    lstype = type_linearsolver()
+    xtype = type_vector()
+    return "typedef Dune::PDELab::Newton<{}, {}, {}> {};".format(gotype, lstype, xtype, name)
+
+
 @symbol
 def type_stationarylinearproblemsolver():
     typedef_stationarylinearproblemsolver("SLP")
     return "SLP"
 
 
+@symbol
+def type_stationarynonlinearproblemssolver():
+    typedef_stationarynonlinearproblemsolver("SNP")
+    return "SNP"
+
+
 @preamble
 def define_stationarylinearproblemsolver(name):
     slptype = type_stationarylinearproblemsolver()
@@ -818,22 +845,35 @@ def define_stationarylinearproblemsolver(name):
     return "{} {}({}, {}, {}, {});".format(slptype, name, go, ls, x, red)
 
 
+@preamble
+def define_stationarynonlinearproblemsolver(name):
+    snptype = type_stationarynonlinearproblemssolver()
+    go = name_gridoperator()
+    x = name_vector()
+    ls = name_linearsolver()
+    return "{}  {}({}, {}, {});".format(snptype, name, go, x, ls)
+
 @symbol
 def name_stationarylinearproblemsolver():
     define_stationarylinearproblemsolver("slp")
     return "slp"
 
 
+@symbol
+def name_stationarynonlinearproblemsolver():
+    define_stationarynonlinearproblemsolver("snp")
+    return "snp"
+
+
 @preamble
 def dune_solve():
-    from ufl.algorithms.predicates import is_multilinear
-    # This is crap as it does check for linearity of the rank 1 form,
-    # which is by definition True.
-#     if is_multilinear(_form):
-    slp = name_stationarylinearproblemsolver()
-    return "{}.apply();".format(slp)
-#     else:
-#         raise NotImplementedError
+    # Test if form is linear in ansatzfunction
+    if is_linear(_form):
+        slp = name_stationarylinearproblemsolver()
+        return "{}.apply();".format(slp)
+    else:
+        snp = name_stationarynonlinearproblemsolver()
+        return "{}.apply();".format(snp)
 
 
 @preamble

test/CMakeLists.txt

 add_subdirectory(laplace)
 add_subdirectory(poisson)
 add_subdirectory(stokes)
+add_subdirectory(nonlinear)

test/nonlinear/CMakeLists.txt

+add_generated_executable(UFLFILE nonlinear.ufl
+                         TARGET nonlinear_symdiff
+                         FORM_COMPILER_ARGS
+                         )
+
+
+dune_add_system_test(TARGET nonlinear_symdiff
+                     INIFILE nonlinear_symdiff.mini
+                     SCRIPT dune_vtkcompare.py
+                     )
+
+# Add the reference code with the PDELab localoperator that produced
+# the reference vtk file
+add_executable(nonlinear_ref reference_main.cc)
+set_target_properties(nonlinear_ref PROPERTIES EXCLUDE_FROM_ALL 1)

test/nonlinear/nonlinear.ufl

+f = Expression("return -4;")
+g = Expression("return x.two_norm2();")
+
+V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+forms = [(inner(grad(u), grad(v)) + u*u*v - f*v)*dx]

test/nonlinear/nonlinear_symdiff.mini

+__name = nonlinear_symdiff
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = nonlinear_ref
+extension = vtu

test/nonlinear/reference_driver.hh

+#ifndef DUNE_PERFTOOL_TEST_NONLINEAR_REFERENCE_DRIVER_HH
+#define DUNE_PERFTOOL_TEST_NONLINEAR_REFERENCE_DRIVER_HH
+
+#include <dune/common/parametertreeparser.hh>
+#include <dune/pdelab/constraints/conforming.hh>
+#include <dune/pdelab/finiteelementmap/pkfem.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/alugrid/grid.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/function/callableadapter.hh>
+#include <dune/testtools/gridconstruction.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/pdelab/gridoperator/gridoperator.hh>
+#include <dune/pdelab/gridfunctionspace/vtk.hh>
+#include <dune/pdelab/newton/newton.hh>
+#include <dune/perftool/vtkpredicate.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <string>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/finiteelementmap/opbfem.hh>
+
+#include "reference_problem.hh"
+#include "reference_nonlinearpoissonfem.hh"
+
+void driver (int argc, char** argv){
+  // Dimension and important types
+  const int dim = 2;
+  typedef double RF;
+
+  // Open ini file
+  Dune::ParameterTree initree;
+  Dune::ParameterTreeParser::readINITree(argv[1], initree);
+
+  // Create grid
+  typedef Dune::ALUGrid<dim, dim, Dune::simplex, Dune::conforming> Grid;
+  IniGridFactory<Grid> factory(initree);
+  std::shared_ptr<Grid> grid = factory.getGrid();
+  typedef Grid::LeafGridView GV;
+  GV gv = grid->leafGridView();
+
+  // Make PDE parameter class
+  RF eta = 1.0;
+  Problem<RF> problem(eta);
+  auto glambda = [&](const auto& e, const auto& x)
+    {return problem.g(e,x);};
+  auto g = Dune::PDELab::
+    makeGridFunctionFromCallable(gv,glambda);
+  auto blambda = [&](const auto& i, const auto& x)
+    {return problem.b(i,x);};
+  auto b = Dune::PDELab::
+    makeBoundaryConditionFromCallable(gv,blambda);
+
+  // Make grid function space
+  typedef Grid::ctype DF;
+
+  typedef Dune::PDELab::PkLocalFiniteElementMap<GV,DF,double,1> FEM;
+  FEM fem(gv);
+  typedef Dune::PDELab::ConformingDirichletConstraints CON;
+  typedef Dune::PDELab::istl::VectorBackend<> VBE;
+  typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE> GFS;
+  GFS gfs(gv,fem);
+  gfs.name("Vh");
+
+  // Assemble constraints
+  typedef typename GFS::template
+    ConstraintsContainer<RF>::Type CC;
+  CC cc;
+  Dune::PDELab::constraints(b,gfs,cc); // assemble constraints
+  std::cout << "constrained dofs=" << cc.size() << " of "
+            << gfs.globalSize() << std::endl;
+
+  // A coefficient vector
+  using Z = Dune::PDELab::Backend::Vector<GFS,RF>;
+  Z z(gfs); // initial value
+
+  // Make a grid function out of it
+  typedef Dune::PDELab::DiscreteGridFunction<GFS,Z> ZDGF;
+  ZDGF zdgf(gfs,z);
+
+  // Fill the coefficient vector
+  Dune::PDELab::interpolate(g,gfs,z);
+
+  // Make a local operator
+  typedef NonlinearPoissonFEM<Problem<RF>,FEM> LOP;
+  LOP lop(problem);
+
+  // Make a global operator
+  typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
+  int degree = initree.get("fem.degree",(int)1);
+  MBE mbe((int)pow(1+2*degree,dim));
+  typedef Dune::PDELab::GridOperator<
+    GFS,GFS,  /* ansatz and test space */
+    LOP,      /* local operator */
+    MBE,      /* matrix backend */
+    RF,RF,RF, /* domain, range, jacobian field type*/
+    CC,CC     /* constraints for ansatz and test space */
+    > GO;
+  GO go(gfs,cc,gfs,cc,lop,mbe);
+
+  // Select a linear solver backend
+  typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LS;
+  LS ls(false);
+
+  // Set up nonlinear solver
+  Dune::PDELab::Newton<GO,LS,Z> newton(go,z,ls);
+  // newton.setReassembleThreshold(0.0); // always reassemble J
+  // newton.setVerbosityLevel(3);        // be verbose
+  // newton.setReduction(1e-10);         // total reduction
+  // newton.setMinLinearReduction(1e-4); // min. red. in lin. solve
+  // newton.setMaxIterations(25);        // limit number of its
+  // newton.setLineSearchMaxIterations(10); // limit line search
+
+  // Solve nonlinear problem
+  newton.apply();
+
+  // Write VTK output file
+  int subsampling = initree.get("vtk.subsamplinglevel",(int)0);
+  Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,subsampling);
+  Dune::PDELab::addSolutionToVTKWriter(vtkwriter, gfs, z);
+  std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output");
+  vtkwriter.write(vtkfile, Dune::VTK::ascii);
+}
+
+#endif

test/nonlinear/reference_main.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+
+#include "reference_driver.hh"
+
+int main(int argc, char** argv)
+{
+  try{
+    //Maybe initialize Mpi
+    Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
+    if(Dune::MPIHelper::isFake)
+      std::cout<< "This is a sequential program." << std::endl;
+    else
+      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
+        <<" processes!"<<std::endl;
+
+    driver(argc, argv);
+
+    return 0;
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+}

test/nonlinear/reference_nonlinearpoissonfem.hh

+#include<vector>
+
+#include<dune/common/exceptions.hh>
+#include<dune/common/fvector.hh>
+#include<dune/geometry/type.hh>
+
+#include<dune/pdelab/common/referenceelements.hh>
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+
+/** a local operator for solving the nonlinear Poisson equation with conforming FEM
+ *
+ * \f{align*}{
+ *   -\Delta u(x) + q(u(x)) &=& f(x) x\in\Omega,  \\
+ *                     u(x) &=& g(x) x\in\partial\Omega_D \\
+ *  -\nabla u(x) \cdot n(x) &=& j(x) x\in\partial\Omega_N \\
+ * \f}
+ *
+ */
+template<typename Param, typename FEM>
+class NonlinearPoissonFEM :
+  public Dune::PDELab::
+    NumericalJacobianVolume<NonlinearPoissonFEM<Param,FEM> >,
+  public Dune::PDELab::
+    NumericalJacobianApplyVolume<NonlinearPoissonFEM<Param,FEM> >,
+  public Dune::PDELab::FullVolumePattern,
+  public Dune::PDELab::LocalOperatorDefaultFlags
+{
+  typedef typename FEM::Traits::FiniteElementType::
+     Traits::LocalBasisType LocalBasis;
+  Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+  Param& param; // parameter functions
+  int incrementorder; // increase of integration order
+
+public:
+  // pattern assembly flags
+  enum { doPatternVolume = true };
+
+  // residual assembly flags
+  enum { doLambdaVolume = true };
+  enum { doLambdaBoundary = true };
+  enum { doAlphaVolume = true };
+
+  //! constructor stores a copy of the parameter object
+  NonlinearPoissonFEM (Param& param_, int incrementorder_=0)
+    : param(param_), incrementorder(incrementorder_)
+  {}
+
+  //! right hand side integral
+  template<typename EG, typename LFSV, typename R>
+  void lambda_volume (const EG& eg, const LFSV& lfsv,
+                      R& r) const
+  {
+    // select quadrature rule
+    auto geo = eg.geometry();
+    const int order = incrementorder+
+      2*lfsv.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(geo,order);
+
+    // loop over quadrature points
+    for (const auto& ip : rule)
+      {
+        // evaluate basis functions
+        auto& phihat = cache.evaluateFunction(ip.position(),
+                             lfsv.finiteElement().localBasis());
+
+        // integrate -f*phi_i
+        decltype(ip.weight()) factor = ip.weight()*
+          geo.integrationElement(ip.position());
+        auto f=param.f(eg.entity(),ip.position());
+        for (size_t i=0; i<lfsv.size(); i++)
+          r.accumulate(lfsv,i,-f*phihat[i]*factor);
+      }
+  }
+
+  // Neumann boundary integral
+  template<typename IG, typename LFSV, typename R>
+  void lambda_boundary (const IG& ig, const LFSV& lfsv,
+                        R& r) const
+  {
+    // evaluate boundary condition type
+    auto localgeo = ig.geometryInInside();
+    auto facecenterlocal = Dune::PDELab::
+      referenceElement(localgeo).position(0,0);
+    bool isdirichlet=param.b(ig.intersection(),facecenterlocal);
+
+    // skip rest if we are on Dirichlet boundary
+    if (isdirichlet) return;
+
+    // select quadrature rule
+    auto globalgeo = ig.geometry();
+    const int order = incrementorder+
+      2*lfsv.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(globalgeo,order);
+
+    // loop over quadrature points and integrate normal flux
+    for (const auto& ip : rule)
+      {
+        // quadrature point in local coordinates of element
+        auto local = localgeo.global(ip.position());
+
+        // evaluate shape functions (assume Galerkin method)
+        auto& phihat = cache.evaluateFunction(local,
+                             lfsv.finiteElement().localBasis());
+
+        // integrate j
+        decltype(ip.weight()) factor = ip.weight()*
+          globalgeo.integrationElement(ip.position());
+        auto j = param.j(ig.intersection(),ip.position());
+        for (size_t i=0; i<lfsv.size(); i++)
+          r.accumulate(lfsv,i,j*phihat[i]*factor);
+      }
+  }
+
+  //! volume integral depending on test and ansatz functions
+  template<typename EG, typename LFSU, typename X,
+           typename LFSV, typename R>
+  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
+                     const LFSV& lfsv, R& r) const
+  {
+    // types & dimension
+    const int dim = EG::Entity::dimension;
+    typedef decltype(Dune::PDELab::
+                     makeZeroBasisFieldValue(lfsu)) RF;
+
+    // select quadrature rule
+    auto geo = eg.geometry();
+    const int order = incrementorder+
+      2*lfsu.finiteElement().localBasis().order();
+    auto rule = Dune::PDELab::quadratureRule(geo,order);
+
+    // loop over quadrature points
+    for (const auto& ip : rule)
+      {
+        // evaluate basis functions
+        auto& phihat = cache.evaluateFunction(ip.position(),
+                             lfsu.finiteElement().localBasis());
+
+        // evaluate u
+        RF u=0.0;
+        for (size_t i=0; i<lfsu.size(); i++)
+          u += x(lfsu,i)*phihat[i];
+
+        // evaluate gradient of shape functions
+        auto& gradphihat = cache.evaluateJacobian(ip.position(),
+                             lfsu.finiteElement().localBasis());
+
+        // transform gradients of shape functions to real element
+        const auto S = geo.jacobianInverseTransposed(ip.position());
+        auto gradphi = makeJacobianContainer(lfsu);
+        for (size_t i=0; i<lfsu.size(); i++)
+          S.mv(gradphihat[i][0],gradphi[i][0]);
+
+        // compute gradient of u
+        Dune::FieldVector<RF,dim> gradu(0.0);
+        for (size_t i=0; i<lfsu.size(); i++)
+          gradu.axpy(x(lfsu,i),gradphi[i][0]);
+
+        // integrate (grad u)*grad phi_i + q(u)*phi_i
+        auto factor = ip.weight()*
+          geo.integrationElement(ip.position());
+        auto q = param.q(u);
+        for (size_t i=0; i<lfsu.size(); i++)
+          r.accumulate(lfsu,i,(gradu*gradphi[i][0]+
+                               q*phihat[i])*factor);
+      }
+  }
+};

test/nonlinear/reference_problem.hh

+template<typename Number>
+class Problem
+{
+  Number eta;
+public:
+  typedef Number value_type;
+
+  //! Constructor without arg sets nonlinear term to zero
+  Problem () : eta(0.0) {}
+
+  //! Constructor takes eta parameter
+  Problem (const Number& eta_) : eta(eta_) {}
+
+  //! nonlinearity
+  Number q (Number u) const
+  {
+    return eta*u*u;
+  }
+
+  //! derivative of nonlinearity
+  Number qprime (Number u) const
+  {
+    return 2*eta*u;
+  }
+
+  //! right hand side
+  template<typename E, typename X>
+  Number f (const E& e, const X& x) const
+  {
+    return -2.0*x.size();
+  }
+
+  //! boundary condition type function (true = Dirichlet)
+  template<typename I, typename X>
+  bool b (const I& i, const X& x) const
+  {
+    return true;
+  }
+
+  //! Dirichlet extension
+  template<typename E, typename X>
+  Number g (const E& e, const X& x) const
+  {
+    auto global = e.geometry().global(x);
+    Number s=0.0;
+    for (std::size_t i=0; i<global.size(); i++) s+=global[i]*global[i];
+    return s;
+  }
+
+  //! Neumann boundary condition
+  template<typename I, typename X>
+  Number j (const I& i, const X& x) const
+  {
+    return 0.0;
+  }
+};