diff --git a/src/test_2pchns_r/2p_driver.hh b/src/test_2pchns_r/2p_driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f949feb8c7297f5421c930b251c4a26ba3777350
--- /dev/null
+++ b/src/test_2pchns_r/2p_driver.hh
@@ -0,0 +1,240 @@
+#pragma once
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<iostream>
+// Dune includes
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/pdelab/newton/newton.hh>
+
+#include <dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh>
+#include <dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh>
+#include <dune/phasefield/localoperator/pf_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/fs_chns_r.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/boundary_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename Grid, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+
+  typedef typename Grid::LeafGridView GV;
+  //std::bitset<2> periodic_directions(std::string("01"));
+  //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+  GV gv = grid.leafGridView();
+  const int dim = GV::dimension;
+  typedef double RF;
+
+
+  typedef Params_fs_r< RF,
+                        DoubleWell_poly<RF>,
+                        Mobility_Const<RF>,
+                        Reaction_Quadratic_Limited<RF>,
+                        VelDissipation_Quadratic<RF> > Parameters;
+  Parameters param(ptree);
+  TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+
+//_________________________________________________________________________________________________________________________
+// Make grid function spaces
+
+
+  typedef GFS_fs_chns_r<RF, ES, V_FEM, P_FEM, PHI_FEM, PHI_FEM> GFS;
+  GFS gfs(es,vFem,pFem,phiFem,phiFem);
+  typedef typename GFS::NS NS_GFS;
+  typedef typename GFS::CH CH_GFS;
+
+
+  //_________________________________________________________________________________________________________________________
+  //Constraints
+
+  typedef Bdry_fs_chns_r<dim, InflowOutflow_v,
+              InflowOutflow_p,
+              InflowOutflow_phi,
+              Dune::PDELab::NoDirichletConstraintsParameters,
+              InflowOutflow_u> Bdry;
+  Bdry bdry;
+  gfs.updateConstraintsContainer(bdry);
+
+  //_________________________________________________________________________________________________________________________
+  //Coefficient vectors
+
+  using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
+  using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
+  CH_V chV(gfs.ch);
+  NS_V nsV(gfs.ns);
+
+  //_________________________________________________________________________________________________________________________
+  //Initial Values
+
+    typedef Ini_fs_chns_r< VelocityInitial_Inflow<GV,RF,Parameters,dim>,
+                      ZeroInitial<GV,RF,Parameters>,
+                      PhiInitial_Sphere<GV,RF,Parameters>,
+                      ZeroInitial<GV,RF,Parameters>,
+                      UInitial_Inflow<GV,RF,Parameters> >  Initial;
+    Initial initial(gv,param);
+    Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+    Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+  //____________________________________________________________________________________________________________________
+  // Refine Initial Data
+
+  for(int i=0; i<4; i++)
+  {
+    // mark elements for refinement
+    typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+    RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+    rInd.calculateIndicator(gv,chV,false);
+    rInd.markGrid(grid,2);
+
+    // do refinement
+    auto chT = transferSolutions(gfs.ch,4,chV);
+    auto nsT = transferSolutions(gfs.ns,4,nsV);
+    Dune::PDELab::adaptGrid(grid, chT,nsT);
+
+    Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+    Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+  }
+
+
+  gfs.updateConstraintsContainer(bdry);
+  interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+
+  CH_V chVOld = chV;
+  NS_V nsVOld = nsV;
+
+  //_________________________________________________________________________________________________________________________
+  // LocalOperators
+
+  typedef FEM_FS_CHNS_R_InstatStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+  NS_LOP nsLop(param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
+
+  typedef FEM_FS_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
+  CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
+
+  typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+  MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+
+
+  //_________________________________________________________________________________________________________________________
+  // Construct discrete Grid Functions
+  typedef DGF_fs_chns_r<GFS, NS_V, CH_V> DGF;
+  DGF dgf(gfs, nsV, chV);
+
+  //_________________________________________________________________________________________________________________________
+  // prepare VTK writer
+  std::string filename=ptree.get("output.filename","output");
+  struct stat st;
+  if( stat( filename.c_str(), &st ) != 0 )
+  {
+    int stat = 0;
+    stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+    if( stat != 0 && stat != -1)
+      std::cout << "Error: Cannot create directory "  << filename << std::endl;
+  }
+  auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
+  Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
+  dgf.addToVTKWriter(vtkwriter);
+
+
+
+  vtkwriter.write(param.time.t,Dune::VTK::ascii);
+
+  //_________________________________________________________________________________________________________________________
+  printVector(chV,10,"chV");
+  printVector(nsV,10,"nsV");
+
+
+  //____________________________________________________________________________________________________________________
+  //int i_it=0;
+
+
+  while (!timestepManager.isFinished())
+    {
+
+      // Grid Operators
+      typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
+      NS_GO nsGo(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,nsLop,mbe);
+
+      typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
+      CH_GO chGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,chLop,mbe);
+
+      // Linear solver
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,3,1);
+
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,500,3,1);
+
+      // Nonlinear solver
+      typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
+      PDESolver1 ns_newton(nsGo,nsV,ls_ns);
+      ns_newton.setParameters(ptree.sub("NS_Newton"));
+
+      typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
+      PDESolver2 ch_newton(chGo,chV,ls_ch);
+      ch_newton.setParameters(ptree.sub("Phasefield_Newton"));
+
+      //_________________________________________________________________________________________________________________________
+      // do time step
+
+
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
+
+      int result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
+      if(result <0)
+      {
+        return;
+      }
+      else if(result == 0)
+      {
+        nsV = nsVOld;
+        chV = chVOld;
+        continue;
+      }
+
+      chVOld = chV;
+      nsVOld = nsV;
+
+      //_________________________________________________________________________________________________________________________
+
+      // mark elements for refinement
+      typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+      RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
+      rInd.calculateIndicator(gv,chV,false);
+      rInd.markGrid(grid,2);
+
+      printVector(chV,3,"chV");
+      printVector(nsV,3,"nsV");
+
+      // do refinement
+      std::cout << "adapt grid and solution" << std::endl;
+      auto chTransfer = transferSolutions(gfs.ch,4,chV,chVOld);
+      auto nsTransfer = transferSolutions(gfs.ns,4,nsV,nsVOld);
+      Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
+      std::cout << "adapt grid and solution ended" << std::endl;
+
+      printVector(chV,3,"chV");
+      printVector(nsV,3,"nsV");
+
+      //vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
+
+      // recompute constraints
+      gfs.updateConstraintsContainer(bdry);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
+
+
+      //vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
+    }
+}
diff --git a/src/test_2pchns_r/2p_equation.cc b/src/test_2pchns_r/2p_equation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..00e5afdf6c2f17211d2404ee708a5a872014b87c
--- /dev/null
+++ b/src/test_2pchns_r/2p_equation.cc
@@ -0,0 +1,116 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include<iostream>
+// dune includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.hh>
+//#if HAVE_UG
+//#include<dune/grid/uggrid.hh>
+//#endif
+#if HAVE_DUNE_ALUGRID
+#include<dune/alugrid/grid.hh>
+#include<dune/alugrid/dgf.hh>
+#include<dune/grid/io/file/dgfparser/dgfparser.hh>
+#endif
+// pdelab includes
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+
+#include "2p_driver.hh"
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+
+
+
+
+int main(int argc, char** argv)
+{
+  #if !HAVE_SUPERLU
+    std::cout << "Error: These examples work only if SuperLU is available." << std::endl;
+    exit(1);
+  #endif
+
+  #if HAVE_UG
+  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 << "Parallel code run on "
+                << helper.size() << " process(es)" << std::endl;
+
+    // open ini file
+    Dune::ParameterTree ptree;
+    Dune::ParameterTreeParser ptreeparser;
+    ptreeparser.readINITree("2pchns_r.ini",ptree);
+    ptreeparser.readOptions(argc,argv,ptree);
+
+    // read ini file
+//    const int dim = ptree.get<int>("grid.dim");
+    const int refinement = ptree.get<int>("domain.level");
+
+
+    if (true/*dim==2*/)
+        {
+
+          typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+          //typedef Dune::UGGrid<2> Grid;
+          typedef Grid::ctype DF;
+
+          std::vector<int> boundary_index_map;
+          std::vector<int> element_index_map;
+
+          std::string filename = ptree.get("domain.filename",
+                                           "turbtube2d.msh");
+          Dune::GridFactory<Grid> factory;
+          if (helper.rank()==0)
+            Dune::GmshReader<Grid>::read(factory,filename,boundary_index_map,element_index_map,true,false);
+          std::shared_ptr<Grid> gridp(factory.createGrid());
+          //gridp->globalRefine(refinement);
+          //gridp->loadBalance();
+          //std::cout << " after load balance /" << helper.rank() << "/ " << gridp->size(0) << std::endl;
+
+          //typedef Grid::LeafGridView LGV;
+
+          //std::bitset<2> periodic_directions(std::string("01"));
+          //typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
+          //GV gv(gridp->leafGridView(),periodic_directions);
+
+          typedef Grid::LeafGridView GV;
+          GV gv = gridp->leafGridView();
+
+          std::cout << "OK" << std::endl;
+
+          using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
+          ES es(gv);
+
+          const int k = 2;
+          typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
+          typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k-1> P_FEM;
+          const int l = 1;
+          typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+          V_FEM vFem(es);
+          P_FEM pFem(es);
+          PHI_FEM phiFem(es);
+
+          driver_fs(*gridp, es,vFem,pFem,phiFem, ptree, helper);
+          }
+
+        }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+  #endif
+}
diff --git a/src/test_2pchns_r/2pchns_r.ini b/src/test_2pchns_r/2pchns_r.ini
new file mode 100644
index 0000000000000000000000000000000000000000..fcdb12312ce985341d944de91a1d63ba7744131d
--- /dev/null
+++ b/src/test_2pchns_r/2pchns_r.ini
@@ -0,0 +1,59 @@
+[physics]
+uAst = 2
+D = 0.005
+rho1 = 1.0
+rho2 = 1.0
+mu=0.01
+ReactionRate=0 #=3
+#ReactionLimiter = 1000
+#VelDissipationRate = 10
+CurvatureAlpha=0.03
+
+[domain]
+level = 0
+#filename = grids/unitsquare.msh
+#filename = grids/square_27k_del.msh
+#filename = grids/square_7k_del.msh
+filename = grids/rectangle_3k_per.msh
+#filename = grids/rectangle_53k_per.msh
+#filename = grids/square_2k_del_per.msh
+
+[Phasefield]
+eps=0.005
+delta=0.02
+Mobility=0.005
+
+[Time]
+tMax = 30
+dt = 0.015
+dtMax = 0.05
+dtMin = 0.0015
+#saveintervall = 0.5
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+MeanInflow=1
+uConst=1.0
+uInflow=1.3
+Radius=0.15
+
+[Refinement]
+etaRefine  = 0.2
+etaCoarsen = 0.1
+maxRefineLevel = 7
+minRefineLevel = 0
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 1
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+VerbosityLevel = 1
diff --git a/src/test_2pchns_r/CMakeLists.txt b/src/test_2pchns_r/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4fa37dd41326644365959ae4b097b9dd1d09db08
--- /dev/null
+++ b/src/test_2pchns_r/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(2p_equation 2p_equation.cc)
+dune_symlink_to_source_files(FILES grids 2pchns_r.ini)
diff --git a/src/test_2pchns_r/myincludes.hh b/src/test_2pchns_r/myincludes.hh
new file mode 100644
index 0000000000000000000000000000000000000000..86fe4121616bbfa13ebf9ecde4742476c0dfff24
--- /dev/null
+++ b/src/test_2pchns_r/myincludes.hh
@@ -0,0 +1,59 @@
+#pragma once
+
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<math.h>
+#include<iostream>
+// dune-common includes
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/timer.hh>
+// dune-geometry includes
+#include<dune/geometry/referenceelements.hh>
+#include<dune/geometry/quadraturerules.hh>
+// dune-grid includes
+#include<dune/grid/onedgrid.hh>
+#include<dune/grid/yaspgrid.hh>
+#include<dune/grid/utility/structuredgridfactory.hh>
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/grid/io/file/gmshreader.hh>
+#if HAVE_UG
+#include<dune/grid/uggrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
+#include<dune/alugrid/grid.hh>
+#include<dune/alugrid/dgf.hh>
+#include<dune/grid/io/file/dgfparser/dgfparser.hh>
+#endif
+// dune-istl included by pdelab
+// dune-pdelab includes
+#include<dune/pdelab/common/function.hh>
+#include<dune/pdelab/common/vtkexport.hh>
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+#include<dune/pdelab/finiteelementmap/qkfem.hh>
+#include<dune/pdelab/finiteelementmap/p0fem.hh>
+#include<dune/pdelab/constraints/common/constraints.hh>
+#include<dune/pdelab/constraints/common/constraintsparameters.hh>
+#include<dune/pdelab/constraints/conforming.hh>
+#include<dune/pdelab/function/callableadapter.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include<dune/pdelab/gridfunctionspace/subspace.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
+#include<dune/pdelab/gridfunctionspace/interpolate.hh>
+#include<dune/pdelab/gridfunctionspace/vtk.hh>
+#include<dune/pdelab/gridoperator/gridoperator.hh>
+#include<dune/pdelab/gridoperator/onestep.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/backend/istl.hh>
+#include<dune/pdelab/stationary/linearproblem.hh>
+#include<dune/pdelab/instationary/onestep.hh>
+#include<dune/pdelab/newton/newton.hh>
+#include<dune/pdelab/adaptivity/adaptivity.hh>
+
+//#include <dune/xt/grid/grids.hh>
+//#include <dune/xt/grid/view/periodic.hh>