From b00c3fe44df1a9a514e2006f164dbff4619d5c40 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Tue, 9 Aug 2016 13:03:33 +0200
Subject: [PATCH] Solve nonlinear problem
---
dune/perftool/vtkpredicate.hh | 4 +-
python/dune/perftool/pdelab/driver.py | 56 +++++-
test/CMakeLists.txt | 1 +
test/nonlinear/CMakeLists.txt | 15 ++
test/nonlinear/nonlinear.ufl | 8 +
test/nonlinear/nonlinear_ref.vtu | 3 +
test/nonlinear/nonlinear_symdiff.mini | 11 ++
test/nonlinear/reference_driver.hh | 124 +++++++++++++
test/nonlinear/reference_main.cc | 33 ++++
.../reference_nonlinearpoissonfem.hh | 172 ++++++++++++++++++
test/nonlinear/reference_problem.hh | 56 ++++++
11 files changed, 474 insertions(+), 9 deletions(-)
create mode 100644 test/nonlinear/CMakeLists.txt
create mode 100644 test/nonlinear/nonlinear.ufl
create mode 100644 test/nonlinear/nonlinear_ref.vtu
create mode 100644 test/nonlinear/nonlinear_symdiff.mini
create mode 100644 test/nonlinear/reference_driver.hh
create mode 100644 test/nonlinear/reference_main.cc
create mode 100644 test/nonlinear/reference_nonlinearpoissonfem.hh
create mode 100644 test/nonlinear/reference_problem.hh
diff --git a/dune/perftool/vtkpredicate.hh b/dune/perftool/vtkpredicate.hh
index fa77e2a3..d567e01e 100644
--- a/dune/perftool/vtkpredicate.hh
+++ b/dune/perftool/vtkpredicate.hh
@@ -1,6 +1,8 @@
#ifndef DUNE_PERFTOOL_VTKPREDICATE_HH
#define DUNE_PERFTOOL_VTKPREDICATE_HH
+#include <dune/typetree/nodeinterface.hh>
+
/** A predicate for vtk output
* that cuts all vector grid function spaces with a length greater than 4.
* This is necessary because VTK cannot handle them.
@@ -14,4 +16,4 @@ struct CuttingPredicate
}
};
-#endif
\ No newline at end of file
+#endif
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 4fe77cdf..3db5c3c7 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/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
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 8b9c34a9..76611b31 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,3 +1,4 @@
add_subdirectory(laplace)
add_subdirectory(poisson)
add_subdirectory(stokes)
+add_subdirectory(nonlinear)
diff --git a/test/nonlinear/CMakeLists.txt b/test/nonlinear/CMakeLists.txt
new file mode 100644
index 00000000..ff8e1624
--- /dev/null
+++ b/test/nonlinear/CMakeLists.txt
@@ -0,0 +1,15 @@
+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)
diff --git a/test/nonlinear/nonlinear.ufl b/test/nonlinear/nonlinear.ufl
new file mode 100644
index 00000000..080010cb
--- /dev/null
+++ b/test/nonlinear/nonlinear.ufl
@@ -0,0 +1,8 @@
+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]
diff --git a/test/nonlinear/nonlinear_ref.vtu b/test/nonlinear/nonlinear_ref.vtu
new file mode 100644
index 00000000..eb718c0a
--- /dev/null
+++ b/test/nonlinear/nonlinear_ref.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:58c404919a3b6e9fea3a0609fc9ae56a9c5b98f5f1de0e86d133b8c8fe34895c
+size 225446
diff --git a/test/nonlinear/nonlinear_symdiff.mini b/test/nonlinear/nonlinear_symdiff.mini
new file mode 100644
index 00000000..31e2aba0
--- /dev/null
+++ b/test/nonlinear/nonlinear_symdiff.mini
@@ -0,0 +1,11 @@
+__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
diff --git a/test/nonlinear/reference_driver.hh b/test/nonlinear/reference_driver.hh
new file mode 100644
index 00000000..db0b072f
--- /dev/null
+++ b/test/nonlinear/reference_driver.hh
@@ -0,0 +1,124 @@
+#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
diff --git a/test/nonlinear/reference_main.cc b/test/nonlinear/reference_main.cc
new file mode 100644
index 00000000..4304d9b0
--- /dev/null
+++ b/test/nonlinear/reference_main.cc
@@ -0,0 +1,33 @@
+#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;
+ }
+}
diff --git a/test/nonlinear/reference_nonlinearpoissonfem.hh b/test/nonlinear/reference_nonlinearpoissonfem.hh
new file mode 100644
index 00000000..556d72a6
--- /dev/null
+++ b/test/nonlinear/reference_nonlinearpoissonfem.hh
@@ -0,0 +1,172 @@
+#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);
+ }
+ }
+};
diff --git a/test/nonlinear/reference_problem.hh b/test/nonlinear/reference_problem.hh
new file mode 100644
index 00000000..df35410c
--- /dev/null
+++ b/test/nonlinear/reference_problem.hh
@@ -0,0 +1,56 @@
+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;
+ }
+};
--
GitLab