diff --git a/test/navier-stokes/howto/CMakeLists.txt b/test/navier-stokes/howto/CMakeLists.txt
index a81818a7e75a3ddb67d292b8cd1efcdac54de992..27fd5d441df271aee459d1a62401458e3b4dd258 100644
--- a/test/navier-stokes/howto/CMakeLists.txt
+++ b/test/navier-stokes/howto/CMakeLists.txt
@@ -1,3 +1,4 @@
add_executable(dgstokes dgstokes.cc)
+add_executable(taylor-green taylor-green.cc)
-dune_symlink_to_source_files(FILES dgstokes.ini)
+dune_symlink_to_source_files(FILES dgstokes.ini taylor-green.ini)
diff --git a/test/navier-stokes/howto/taylor-green.cc b/test/navier-stokes/howto/taylor-green.cc
new file mode 100644
index 0000000000000000000000000000000000000000..00178bbaa4c8c9314d8baa495e83935de7ecee4f
--- /dev/null
+++ b/test/navier-stokes/howto/taylor-green.cc
@@ -0,0 +1,351 @@
+// -*- tab-width: 2; indent-tabs-mode: nil -*-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <iostream>
+#include <vector>
+#include <map>
+#include <string>
+#include <random>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/io.hh>
+
+#include "dune/pdelab/finiteelementmap/qkdg.hh"
+#include <dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
+#include <dune/pdelab/gridfunctionspace/subspace.hh>
+#include <dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
+#include <dune/pdelab/gridfunctionspace/vtk.hh>
+#include <dune/pdelab/gridoperator/gridoperator.hh>
+#include <dune/pdelab/gridfunctionspace/interpolate.hh>
+#include <dune/pdelab/localoperator/dgnavierstokes.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/finiteelementmap/monomfem.hh>
+#include <dune/pdelab/common/function.hh>
+#include <dune/pdelab/common/vtkexport.hh>
+#include <dune/pdelab/constraints/p0.hh>
+#include<dune/pdelab/gridoperator/onestep.hh>
+#include<dune/pdelab/newton/newton.hh>
+#include "dune/perftool/vtkpredicate.hh"
+#include "dune/grid/io/file/vtk/vtksequencewriter.hh"
+
+#include "taylor-green.hh"
+
+
+//===============================================================
+// Problem setup and solution
+//===============================================================
+template<typename GV, typename RF, int vOrder, int pOrder>
+void taylor_green (const GV& gv, const Dune::ParameterTree& configuration, std::string filename)
+{
+ // Some types
+ using ES = Dune::PDELab::AllEntitySet<GV>;
+ ES es(gv);
+ typedef typename ES::Grid::ctype DF;
+ static const unsigned int dim = ES::dimension;
+ Dune::Timer watch;
+ std::cout << "=== Initialize" << std::endl;
+
+ // Create grid function spaces
+ watch.reset();
+ typedef Dune::PDELab::QkDGLocalFiniteElementMap<DF, RF, 2, 2> FEM_V;
+ typedef Dune::PDELab::QkDGLocalFiniteElementMap<DF, RF, 1, 2> FEM_P;
+ FEM_V fem_v;
+ FEM_P fem_p;
+
+ static const unsigned int vBlockSize = Dune::QkStuff::QkSize<vOrder,dim>::value;
+ static const unsigned int pBlockSize = Dune::QkStuff::QkSize<pOrder,dim>::value;
+ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none,vBlockSize> VVectorBackend;
+ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none,pBlockSize> PVectorBackend;
+ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none> VectorBackend;
+ // typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VectorBackend;
+
+ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none> VVectorBackend;
+ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none> PVectorBackend;
+ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none> VectorBackend;
+
+ typedef Dune::PDELab::istl::VectorBackend<> VelocityVectorBackend;
+
+
+ using Con = Dune::PDELab::P0ParallelConstraints;
+ Con con;
+
+ // velocity
+ Dune::dinfo << "--- v^dim" << std::endl;
+ // typedef Dune::PDELab::EntityBlockedOrderingTag VelocityOrderingTag;
+ typedef Dune::PDELab::LexicographicOrderingTag VelocityOrderingTag;
+ typedef Dune::PDELab::VectorGridFunctionSpace<
+ ES,vFEM,dim,
+ VelocityVectorBackend,
+ VVectorBackend,
+ Con,
+ VelocityOrderingTag
+ > velocityGFS;
+ velocityGFS velocityGfs(es,vFem);
+ velocityGfs.name("v");
+ // p
+ Dune::dinfo << "--- p" << std::endl;
+ typedef Dune::PDELab::GridFunctionSpace<
+ ES,
+ pFEM,
+ Con,
+ PVectorBackend> pGFS;
+ pGFS pGfs(es,pFem);
+ pGfs.name("p");
+ // GFS
+ Dune::dinfo << "--- v^dim,p" << std::endl;
+
+
+ // typedef Dune::PDELab::EntityBlockedOrderingTag StokesOrderingTag;
+ typedef Dune::PDELab::LexicographicOrderingTag StokesOrderingTag;
+
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<
+ VectorBackend, StokesOrderingTag,
+ velocityGFS, pGFS> GFS;
+ GFS gfs(velocityGfs, pGfs);
+ using namespace Dune::Indices;
+ velocityGfs.child(_0).name("velocity_0");
+ velocityGfs.child(_1).name("velocity_1");
+ pGfs.name("pressure");
+ gfs.name("test");
+ gfs.update();
+ using CC = typename GFS::template ConstraintsContainer<double>::Type;
+ CC cc;
+ cc.clear();
+ std::cout << "gfs with " << gfs.size() << " dofs generated "<< std::endl;
+ std::cout << "cc with " << cc.size() << " dofs generated "<< std::endl;
+ std::cout << "=== function space setup " << watch.elapsed() << " s" << std::endl;
+
+ // <<<3>>> Make coefficient Vector and initialize it from a function
+ using BType = BCTypeParamGlobalDirichlet;
+ BType b;
+ using FType = ZeroVectorFunction<ES,RF,dim>;
+ FType f(es);
+ using VType = TaylorGreenVelocity<ES,RF,dim>;
+ VType v(es);
+ using PType = TaylorGreenPressure<ES,RF>;
+ PType p(es);
+
+ // <<<4>>> Make grid Function operator
+ watch.reset();
+ typedef Dune::PDELab::DefaultInteriorPenalty<RF> PenaltyTerm;
+
+ typedef Dune::PDELab::DGNavierStokesParameters<ES,RF,FType,BType,VType,PType,true,false,PenaltyTerm>
+ LocalDGOperatorParameters;
+ LocalDGOperatorParameters lop_params(configuration.sub("parameters"),f,b,v,p);
+ typedef Dune::PDELab::DGNavierStokes<LocalDGOperatorParameters> LOP;
+ const int superintegration_order = 0;
+ LOP lop(lop_params,superintegration_order);
+
+ typedef Dune::PDELab::NavierStokesMass<LocalDGOperatorParameters> MLOP;
+ MLOP mlop(lop_params,1);
+
+
+ typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
+ MBE mbe(75); // Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+ typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,CC,CC> GO_R;
+ GO_R go_r(gfs,cc,gfs,cc,lop,mbe);
+ typedef Dune::PDELab::GridOperator<GFS,GFS,MLOP,MBE,RF,RF,RF,CC,CC> GO_M;
+ GO_M go_m(gfs,cc,gfs,cc,mlop,mbe);
+
+ typedef Dune::PDELab::OneStepGridOperator<GO_R,GO_M> IGO;
+ IGO igo(go_r,go_m);
+ using V = typename IGO::Traits::Domain;
+ std::cout << "=== grid operator space setup " << watch.elapsed() << " s" << std::endl;
+
+
+ // Linear solver
+ using LinearSolver = Dune::PDELab::ISTLBackend_OVLP_BCGS_ILU0<GFS,CC>;
+ LinearSolver ls(gfs,cc);
+
+ // Solve (possibly) nonlinear problem
+ std::cout << "=== Setup Newton:" << std::endl;
+ watch.reset();
+ typedef Dune::PDELab::Newton<IGO,LinearSolver,V> PDESolver;
+ PDESolver newton(igo,ls);
+ // newton.setReassembleThreshold(0.0);
+ // newton.setVerbosityLevel(2);
+ // newton.setMaxIterations(50);
+ // newton.setLineSearchMaxIterations(30);
+
+ // Dune::PDELab::FractionalStepParameter<RF> method;
+ using TSM = Dune::PDELab::OneStepThetaParameter<RF>;
+ TSM tsm(1.0);
+ Dune::PDELab::OneStepMethod<RF,IGO,PDESolver,V,V> osm(tsm,igo,newton);
+ // osm.setVerbosityLevel(2);
+
+ // Create initial solution
+ using InitialVelocity = TaylorGreenVelocity<GV,RF,2>;
+ InitialVelocity initial_velocity(gv);
+ using InitialPressure = TaylorGreenPressure<GV,RF>;
+ InitialPressure initial_pressure(gv);
+ typedef Dune::PDELab::CompositeGridFunction<InitialVelocity,InitialPressure> InitialSolution;
+ InitialSolution initial_solution(initial_velocity, initial_pressure);
+
+ // Make coefficent vector and initialize it from a function
+ V xold(gfs);
+ xold = 0.0;
+ Dune::PDELab::interpolate(initial_solution,gfs,xold);
+ std::cout << "=== Finished interpolation:" << watch.elapsed() << std::endl;
+ watch.reset();
+
+ watch.reset();
+ RF time = 0.0;
+ RF time_end = 1.0;
+ RF dt = 1e-2;
+ RF dt_min = 1e-6;
+
+ // Visualize initial condition
+ using VTKSW = Dune::VTKSequenceWriter<GV>;
+ using VTKWriter = Dune::SubsamplingVTKWriter<GV>;
+ VTKWriter vtkwriter(gv, 1);
+ VTKSW vtkSequenceWriter(std::make_shared<VTKWriter>(vtkwriter), filename);
+ CuttingPredicate predicate;
+ Dune::PDELab::addSolutionToVTKWriter(vtkSequenceWriter, gfs, xold, Dune::PDELab::vtk::defaultNameScheme(), predicate);
+ vtkSequenceWriter.write(time, Dune::VTK::appendedraw);
+
+
+ V x(gfs,0.0);
+ while (time < time_end - dt_min*0.5){
+ osm.apply(time,dt,xold,x);
+
+ xold = x;
+ time += dt;
+
+ vtkSequenceWriter.write(time, Dune::VTK::appendedraw);
+ }
+
+ std::cout << "=== Total time:" << watch.elapsed() << std::endl;
+
+
+// typedef typename GOS::Jacobian M;
+// watch.reset();
+// M m(gos);
+// std::cout << m.patternStatistics() << std::endl;
+// std::cout << "=== matrix setup " << watch.elapsed() << " s" << std::endl;
+// m = 0.0;
+// watch.reset();
+// gos.jacobian(x,m);
+// std::cout << "=== jacobian assembly " << watch.elapsed() << " s" << std::endl;
+
+// // std::ofstream matrix("Matrix");
+// // Dune::printmatrix(matrix, m.base(), "M", "r", 6, 3);
+
+// // evaluate residual w.r.t initial guess
+// V r(gfs);
+// r = 0.0;
+// watch.reset();
+// x = 0.0;
+// gos.residual(x,r);
+// std::cout << "=== residual evaluation " << watch.elapsed() << " s" << std::endl;
+
+// bool verbose = true;
+
+// using Dune::PDELab::Backend::native;
+// typedef Dune::PDELab::Backend::Native<M> ISTLM;
+// typedef Dune::PDELab::Backend::Native<V> ISTLV;
+// #ifdef USE_SUPER_LU // use lu decomposition as solver
+// #if HAVE_SUPERLU
+// // make ISTL solver
+// Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(native(m));
+// Dune::SuperLU<ISTLM> solver(native(m), verbose?1:0);
+// Dune::InverseOperatorResult stat;
+// #else
+// #error No superLU support, please install and configure it.
+// #endif
+// #else // Use iterative solver
+// // make ISTL solver
+// Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(native(m));
+// Dune::SeqILU0<ISTLM,ISTLV,ISTLV> ilu0(native(m),1.0);
+// Dune::BiCGSTABSolver<ISTLV> solver(opa,ilu0,1E-10,20000, verbose?2:1);
+// Dune::InverseOperatorResult stat;
+// #endif
+
+// // solve the jacobian system
+// r *= -1.0; // need -residual
+// x = r;
+// solver.apply(native(x),native(r),stat);
+
+
+ // if (false) {
+ // using Dune::PDELab::Backend::native;
+ // V r(x);
+ // // Setup random input
+ // std::size_t seed = 0;
+ // auto rng = std::mt19937_64(seed);
+ // auto dist = std::uniform_real_distribution<>(-1., 1.);
+ // for (auto& v : x)
+ // v = dist(rng);
+ // r=0.0;
+ // gos.residual(x, r);
+ // Dune::printvector(std::cout, native(r), "residual vector", "row");
+ // }
+ // if (false) {
+ // // Setup random input
+ // std::size_t seed = 0;
+ // auto rng = std::mt19937_64(seed);
+ // auto dist = std::uniform_real_distribution<>(1., 10.);
+ // for (auto& v : x)
+ // v = dist(rng);
+ // //using M = typename GO::Traits::Jacobian;
+ // M m(gos);
+ // gos.jacobian(x, m);
+ // using Dune::PDELab::Backend::native;
+ // Dune::printmatrix(std::cout, native(m),"global stiffness matrix","row",9,1);
+ // }
+}
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+
+int main(int argc, char** argv)
+{
+ try{
+ //Maybe initialize Mpi
+ Dune::MPIHelper::instance(argc, argv);
+
+ Dune::ParameterTree configuration;
+ const std::string config_filename("taylor-green.ini");
+ std::cout << "Reading ini-file \""<< config_filename
+ << "\"" << std::endl;
+
+ Dune::ParameterTreeParser::readINITree(config_filename, configuration);
+
+ // make grid
+ const int dim = 2;
+ Dune::FieldVector<double,dim> lowerleft(-1.0);
+ Dune::FieldVector<double,dim> upperright(1.0);
+ std::array<int, dim> cells(Dune::fill_array<int, dim>(8));
+ std::bitset<dim> periodic(true);
+ using Grid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+ Grid grid(lowerleft, upperright, cells, periodic, 1);
+
+ // get view
+ typedef Grid::LeafGridView GV;
+ const GV gv=grid.leafGridView();
+
+ // solve problem
+ Dune::dinfo.push(false);
+ taylor_green<GV,double,2,1>(gv,configuration,"taylor-green");
+
+ // test passed
+ 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/navier-stokes/howto/taylor-green.hh b/test/navier-stokes/howto/taylor-green.hh
new file mode 100644
index 0000000000000000000000000000000000000000..44cc59a201c7dd4aa320c92df8cc06ccbed9fe17
--- /dev/null
+++ b/test/navier-stokes/howto/taylor-green.hh
@@ -0,0 +1,145 @@
+#ifndef TAYLOR_GREEN_HH
+#define TAYLOR_GREEN_HH
+
+//===============================================================
+// Define parameter functions f,g,j and \partial\Omega_D/N
+//===============================================================
+
+// constraints parameter class for selecting boundary condition type
+class BCTypeParamGlobalDirichlet
+{
+public :
+ typedef Dune::PDELab::StokesBoundaryCondition BC;
+
+ struct Traits
+ {
+ typedef BC::Type RangeType;
+ };
+
+ BCTypeParamGlobalDirichlet() {}
+
+ template<typename I>
+ inline void evaluate (const I & intersection, /*@\label{bcp:name}@*/
+ const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord,
+ BC::Type& y) const
+ {
+ y = BC::VelocityDirichlet;
+ }
+
+ template<typename T>
+ void setTime(T t){
+ }
+};
+
+
+template<typename GV, typename RF, int dim>
+class TaylorGreenVelocity :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim>,
+ TaylorGreenVelocity<GV,RF,dim> >
+{
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim> Traits;
+ typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,TaylorGreenVelocity<GV,RF,dim> > BaseT;
+
+ typedef typename Traits::DomainType DomainType;
+ typedef typename Traits::RangeType RangeType;
+
+ TaylorGreenVelocity(const GV & gv) : BaseT(gv)
+ {
+ time = 0.0;
+ }
+
+ inline void evaluateGlobal(const DomainType & x, RangeType & y) const
+ {
+ // TODO Get mu and rho from somewhere else!
+ RF pi = M_PI;
+ RF mu = 1.0/100;
+ RF rho = 1.0;
+ RF nu = mu/rho;
+ y[0] = -exp(-2.0*pi*pi*nu*time)*cos(pi*x[0])*sin(pi*x[1]);
+ y[1] = exp(-2.0*pi*pi*nu*time)*sin(pi*x[0])*cos(pi*x[1]);
+ }
+
+ template <typename T>
+ void setTime(T t){
+ time = t;
+ }
+
+private:
+ RF time;
+};
+
+
+template<typename GV, typename RF>
+class TaylorGreenPressure
+ : public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
+ TaylorGreenPressure<GV,RF> >
+{
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+ typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,TaylorGreenPressure<GV,RF> > BaseT;
+
+ typedef typename Traits::DomainType DomainType;
+ typedef typename Traits::RangeType RangeType;
+
+ TaylorGreenPressure (const GV& gv) : BaseT(gv)
+ {
+ time = 0.0;
+ }
+
+ inline void evaluateGlobal (const typename Traits::DomainType& x,
+ typename Traits::RangeType& y) const
+ {
+ RF pi = M_PI;
+ RF mu = 1.0/100;
+ RF rho = 1.0;
+ RF nu = mu/rho;
+ y = -0.25*rho*exp(-4.0*pi*pi*nu*time)*(cos(2.0*pi*x[0])+cos(2.0*pi*x[1]));
+ }
+
+ template<typename T>
+ void setTime(T t){
+ time = t;
+ }
+
+private:
+ RF time;
+};
+
+
+
+template<typename GV, typename RF, std::size_t dim_range>
+class ZeroVectorFunction :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range>,
+ ZeroVectorFunction<GV,RF,dim_range> >,
+ public Dune::PDELab::InstationaryFunctionDefaults
+{
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,dim_range> Traits;
+ typedef Dune::PDELab::AnalyticGridFunctionBase<Traits, ZeroVectorFunction> BaseT;
+
+ typedef typename Traits::DomainType DomainType;
+ typedef typename Traits::RangeType RangeType;
+
+ ZeroVectorFunction(const GV & gv) : BaseT(gv) {}
+
+ inline void evaluateGlobal(const DomainType & x, RangeType & y) const
+ {
+ y=0;
+ }
+};
+
+template<typename GV, typename RF>
+class ZeroScalarFunction
+ : public ZeroVectorFunction<GV,RF,1>
+{
+public:
+
+ ZeroScalarFunction(const GV & gv) : ZeroVectorFunction<GV,RF,1>(gv) {}
+
+};
+
+#endif
diff --git a/test/navier-stokes/howto/taylor-green.ini b/test/navier-stokes/howto/taylor-green.ini
new file mode 100644
index 0000000000000000000000000000000000000000..e88e6a0674b0d79d03afcbb1c807fcf5bafb1655
--- /dev/null
+++ b/test/navier-stokes/howto/taylor-green.ini
@@ -0,0 +1,9 @@
+# Configuration file for the example programs - dgstokes, dgnavierstokesvelvecfem -
+
+[parameters]
+rho = 1.0
+mu = 1.0
+[parameters.dg]
+epsilon = -1
+sigma = 1.0
+beta = 1.0