Added new fluid-solid example (non-monolithic)

dune/phasefield/localoperator/fem_fs_chns_cahnhilliard.hh

+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.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/gridfunctionspace/gridfunctionspace.hh>
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
+//Parameter class (for functions), FEM space and RangeField
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+class FEM_FS_CHNS_CahnHilliard :
+  public Dune::PDELab::FullVolumePattern,
+  public Dune::PDELab::LocalOperatorDefaultFlags,
+  public Dune::PDELab::NumericalJacobianVolume<FEM_FS_CHNS_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+  public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+  public Dune::PDELab::NumericalJacobianBoundary<FEM_FS_CHNS_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+  public Dune::PDELab::NumericalJacobianApplyBoundary<FEM_FS_CHNS_CahnHilliard<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+{
+private:
+  typedef typename CHContainer::field_type RF;
+  typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity   ENTITY;
+
+  CHLocalBasisCache<CHGFS> chCache;
+  VLocalBasisCache<NSGFS> vCache;
+
+  typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+  mutable CHLView chOldLocal;
+
+  typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+  mutable NSLView nsLocal;
+
+
+  Param& param; // parameter functions
+
+public:
+  enum {doPatternVolume=true};
+  enum {doAlphaVolume=true};
+
+  //constructor
+  FEM_FS_CHNS_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, NSGFS& nsGfs_, NSContainer& nsV) :
+    chOldLocal(chGfs_,chVOld_),
+    nsLocal(nsGfs_,nsV),
+    param(param_)
+  {
+  }
+
+
+  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
+  {
+    //Quadrature Rule
+    const int dim = EG::Entity::dimension;
+    const int order = 5;
+    auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+    chOldLocal.DoEvaluation(eg);
+    //chOldLocal.DoEvaluation(eg);
+    nsLocal.DoSubEvaluation(eg,child(nsLocal.lfs,Dune::Indices::_0));
+
+    for(const auto& ip : rule)
+    {
+      const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+      auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+      CHVariables_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(),chOldLocal.lv, S);
+      NSVariables<RF,dim,decltype(nsLocal.lfs)> ns(vCache, ip.position(), nsLocal.lfs, nsLocal.lv, S);
+
+      RF pf_f = ch.pf + 2*param.delta*(1-ch.pf);
+      Dune::FieldVector<RF,dim> gradpf_f(0.0);
+      gradpf_f.axpy((1-2*param.delta),ch.gradpf);
+
+      for (size_t i=0; i<ch.pfspace.size(); i++)
+      {
+        RF J1 = - param.mobility.eval(ch.pf)*param.eps * (ch.gradxi*ch.basis.gradphi[i][0]);
+
+        r.accumulate(ch.pfspace,i,factor*(
+          (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+          -J1
+          -ch.pf*(ns.v*ch.basis.gradphi[i][0])
+        ));
+
+        r.accumulate(ch.xispace,i,factor*(
+          ch.xi*ch.basis.phi[i]
+          - 1/param.eps* param.dw.eval_d(ch.pf) *ch.basis.phi[i]
+          - param.eps * (ch.gradpf*ch.basis.gradphi[i][0])
+        ));
+      }
+    }
+  }
+};

dune/phasefield/localoperator/fem_fs_chns_navierstokes.hh

+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+//#include<dune/pdelab/localoperator/idefault.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/gridfunctionspace/gridfunctionspace.hh>
+
+#include<dune/phasefield/localoperator/fem_base_chns.hh>
+
+
+
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+    class FEM_FS_CHNS_NavierStokes :
+      public Dune::PDELab::FullVolumePattern,
+      public Dune::PDELab::LocalOperatorDefaultFlags,
+      public Dune::PDELab::NumericalJacobianVolume<FEM_FS_CHNS_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+      public Dune::PDELab::NumericalJacobianApplyVolume<FEM_FS_CHNS_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
+    {
+    private:
+      typedef typename CHContainer::field_type RF;
+      typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity   ENTITY;
+
+      CHLocalBasisCache<CHGFS> chCache;
+      VLocalBasisCache<NSGFS> vCache;
+      PLocalBasisCache<NSGFS> pCache;
+
+      typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+      mutable CHLView chLocal;
+      //mutable CHLView chOldLocal;
+
+      typedef GFS_LocalViewer<NSGFS, NSContainer, ENTITY, RF> NSLView;
+      mutable NSLView nsOldLocal;
+
+      Param& param; // parameter functions
+
+
+    public:
+
+      // pattern assembly flags
+      enum { doPatternVolume = true };
+
+      // residual assembly flags
+      enum { doAlphaVolume = true };
+
+
+      FEM_FS_CHNS_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
+                         NSGFS& nsGfs_, NSContainer& nsVOld_) :
+        chLocal(chGfs_,chV_),
+        //chOldLocal(chGfs_,chVOld_),
+        nsOldLocal(nsGfs_,nsVOld_),
+        param(param_)
+      {
+      }
+
+
+      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
+      {
+        //Quadrature Rule
+        const int dim = EG::Entity::dimension;
+        const int order = 5; //TODO: Choose more educated
+        auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+        chLocal.DoEvaluation(eg);
+        //chOldLocal.DoEvaluation(eg);
+        nsOldLocal.DoSubEvaluation(eg,child(nsOldLocal.lfs,Dune::Indices::_0));
+
+        for(const auto& ip : rule)
+        {
+          const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+          auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+          CHVariables<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
+          NSVariables_P_VOld_Extension<RF,dim,LFSV> ns(vCache, pCache, ip.position(), lfsv, x.base(), nsOldLocal.lv, S);
+
+          RF rho_f=param.phys.rho_1*ch.pf;
+          RF rho_f_reg= rho_f + param.delta;
+
+          RF pf_f = ch.pf + 2*param.delta*(1-ch.pf);
+          Dune::FieldVector<RF,dim> gradpf_f(0.0);
+          gradpf_f.axpy((1-2*param.delta),ch.gradpf);
+
+          for (size_t i=0; i<ns.pspace.size(); i++)
+          {
+            //r.accumulate(pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor);
+            r.accumulate(ns.pspace,i, -1.0 * (pf_f*ns.div_v  + (gradpf_f*ns.v) )* ns.pbasis.phi[i] * factor);
+          }
+
+          for (size_t i=0; i<ns.vfemspace.size(); i++)
+          {
+            for(int k=0; k<dim; k++)
+            {
+              const auto v_nabla_v = ns.v * ns.jacv[k];
+              const auto Jf_nabla_v = - param.phys.rho_1 * param.mobility.eval(ch.pf)*param.eps * (ch.gradxi* ns.jacv[k]);
+
+              r.accumulate(child(ns.vspace,k),i, (
+                  rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt                 //TimeEvolution
+                  + ( rho_f* v_nabla_v  + Jf_nabla_v  )* ns.vbasis.phi[i]                       //Momentum Transport
+                  //- ns.p * (pf_f*ns.vbasis.gradphi[i][0][k] + ns.vbasis.phi[i]*gradpf_f[k])   //Pressure in weak form
+                  + pf_f * (ns.gradp[k] * ns.vbasis.phi[i])                                     //Pressure in strong form
+                  + param.phys.mu * (ns.jacv[k]*ns.vbasis.gradphi[i][0])                        //Viscosity
+                  + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i]                              //Velocity Dissipation in Solid
+                                              )*factor);
+            // integrate mu * (grad u)^T * grad phi_i                                                  //Transpose for Viscosity
+            for(int dd=0; dd<dim; ++dd)
+                r.accumulate(child(ns.vspace,k),i, param.phys.mu * (ns.jacv[dd][k] * (ns.vbasis.gradphi[i][0][dd])) * factor);
+            }
+          }
+        }
+      }
+
+    };

src/CMakeLists.txt

@@ -3,3 +3,4 @@ target_link_dune_default_libraries("dune-phasefield")
 add_subdirectory(example_2p_cahnhilliard)
 add_subdirectory(example_2p_cahnhilliard_navierstokes)
 add_subdirectory(example_3p_ffs_monolithic)
+add_subdirectory(example_2p_fs)

src/example_2p_fs/2p_driver.hh

+#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_cahnhilliard.hh>
+#include <dune/phasefield/localoperator/fem_fs_chns_navierstokes.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/ff_chns.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+
+
+//_________________________________________________________________________________________________________________________
+
+template<typename Grid, typename GV, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
+void driver_2p(Grid& grid, GV& gv, ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+{
+  const int dim = GV::dimension;
+  typedef double RF;
+
+  typedef Params_fs< RF,
+                     DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,
+                     Mobility_Const<RF>,
+                     VelDissipation_Quadratic_Shifted<RF>
+                   > Parameters;
+  Parameters param(ptree);
+
+  TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
+
+  //_________________________________________________________________________________________________________________________
+  // Make grid function spaces
+
+  typedef GFS_ff_chns<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
+  GFS gfs(es,vFem,pFem,phiFem);
+  typedef typename GFS::NS NS_GFS;
+  typedef typename GFS::CH CH_GFS;
+
+
+  //_________________________________________________________________________________________________________________________
+  //Constraints
+
+  typedef Bdry_ff_chns<2, Dune::PDELab::DirichletConstraintsParameters,         // v_x, v_y
+              Dune::PDELab::NoDirichletConstraintsParameters,                   // p
+              Dune::PDELab::NoDirichletConstraintsParameters,                   // phi
+              Dune::PDELab::NoDirichletConstraintsParameters>                   // mu
+              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_ff_chns< VelocityInitial_Cavity<GV,RF,Parameters,dim>,       //v
+                      ZeroInitial<GV,RF,Parameters>,                       //p
+                      PhiInitial_Sphere<GV,RF,Parameters>,                 //phi
+                      ZeroInitial<GV,RF,Parameters>>  Initial;             //mu
+  Initial initial(gv,param);
+  Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+  Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+  CH_V chVOld = chV;
+  NS_V nsVOld = nsV;
+
+  //_________________________________________________________________________________________________________________________
+  // LocalOperators
+
+
+  typedef FEM_FS_CHNS_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V > CH_LOP;
+  CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
+
+  typedef FEM_FS_CHNS_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V > NS_LOP;
+  NS_LOP nsLop(param, gfs.ch, chV, gfs.ns, nsVOld);
+
+  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_ff_chns<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);
+
+
+  // Refine Intial Values
+  for(int i=0; i<ptree.get("RefinementInitial.maxRefineLevel",2); i++)
+  {
+        // mark elements for refinement
+        typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND0;
+        RIND0 rInd0(gv,gfs.ch,ptree.sub("RefinementInitial"));
+        rInd0.calculateIndicator(gv,chV,false);
+        rInd0.markGrid(grid,2);
+        rInd0.refineGrid(grid,chV,nsV,gfs.ch,gfs.ns,3);
+        //rInd0.myRefineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1);   // if used with PeriodicGridView
+        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);
+
+
+
+
+  //_________________________________________________________________________________________________________________________
+  printVector(chV,10,"chV");
+  printVector(nsV,10,"nsV");
+  vtkwriter.write(param.time.t,Dune::VTK::ascii);
+
+  //____________________________________________________________________________________________________________________
+  chVOld = chV;
+  nsVOld = nsV;
+  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,600,200,3,1);
+
+      typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
+                  <NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,1000,200,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);
+      int result = timestepManager.doTimestepToConvergence(ns_newton,ch_newton,vtkwriter,nsV,chV);
+
+      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);
+      rInd.refineGrid(grid,chV,chVOld,nsV,nsVOld,gfs.ch,gfs.ns,3);
+      //rInd.myRefineGrid(grid,gv,chV,chVOld,nsV,nsVOld,gfs.ch,gfs.ns,1);  //For PeriodicGridView
+
+      // recompute constraints
+      gfs.updateConstraintsContainer(bdry);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+      interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
+    }
+}

src/example_2p_fs/2p_fs.cc

+#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_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)
+{
+  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("chns2p.ini",ptree);
+    ptreeparser.readOptions(argc,argv,ptree);
+
+    const int refinement = ptree.get<int>("domain.level");
+
+    typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+    typedef Grid::ctype DF;
+
+    std::vector<int> boundary_index_map;
+    std::vector<int> element_index_map;
+
+    std::string filename = ptree.get("domain.filename", "square.msh");
+
+    Dune::GridFactory<Grid> factory;
+    Dune::GmshReader<Grid>::read(factory,filename,true,true);
+    std::shared_ptr<Grid> gridp(factory.createGrid());
+
+    typedef Grid::LeafGridView GV;
+    GV gv = gridp->leafGridView();
+
+    using ES = Dune::PDELab::AllEntitySet<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 = 2;
+    typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+
+    V_FEM vFem(es);
+    P_FEM pFem(es);
+    PHI_FEM phiFem(es);
+
+    driver_2p(*gridp, gv, 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;
+  }
+}

src/example_2p_fs/CMakeLists.txt

+add_executable(2p_fs 2p_fs.cc)
+dune_symlink_to_source_files(FILES grids chns2p.ini)

src/example_2p_fs/chns2p.ini

+[Physics]
+D = 0.02
+rho1 = 2.0
+rho2 = 1.0
+mu=0.01
+
+[domain]
+level = 0
+filename = grids/rectangle_periodic_coarse.msh
+
+[output]
+filename = output1
+
+[Phasefield]
+eps=0.03
+delta=0.03
+DoubleWellDelta=0.02
+DoubleWellGamma=0.015
+Mobility=0.04
+
+[Time]
+tMax = 20
+dt = 0.02
+dtMax = 0.1
+dtMin = 0.0025
+#saveintervall = 0.5
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+
+[RefinementInitial]
+etaRefine  = 0.2
+etaCoarsen = 0.18
+maxRefineLevel = 4
+minRefineLevel = -1
+
+[Refinement]
+etaRefine  = 0.2
+etaCoarsen = 0.1
+maxRefineLevel = 4
+minRefineLevel = -1
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 20
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 2
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 100
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 2