Added pore network test

dune/phasefield/base_parameters.hh

@@ -260,13 +260,17 @@ public:
 
   Params_Initial initial;
 
+  const Dune::ParameterTree& nsNewtonTree;
+  const Dune::ParameterTree& chNewtonTree;
 
   Params_base(const Dune::ParameterTree& ptree) :
     eps(ptree.get("Phasefield.eps",(Number)0.1)),
     mobility(ptree.sub("Phasefield")),
     phys(ptree.sub("Physics"),eps),
     time(0,0),
-    initial(ptree.sub("Initial"))
+    initial(ptree.sub("Initial")),
+    nsNewtonTree(ptree.sub("NS_Newton")),
+    chNewtonTree(ptree.sub("Phasefield_Newton"))
   {
   }
 };

dune/phasefield/porenetwork/fem_1p_navierstokes.hh

@@ -72,12 +72,11 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
           const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
           auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
 
-          CHVariables_Threephase_Micro_Extension<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
-          NSMicroVariables<RF,dim,LFSV> ns(vCache, ip.position(), lfsv, x.base(), S);
-          RF pf3 = 1-ch.pf - ch.pf2;
-          RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;
+          CHVariables<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
+          NSPorenetworkVariables<RF,dim,LFSV> ns(vCache, ip.position(), lfsv, x.base(), S);
+          RF pf_f=ch.pf+param.delta;
 
-          for (size_t i=0; i<ns.vxspace.size(); i++)
+          for (size_t i=0; i<ns.vxspace.size(); i++) //TODO Check
           {
             r.accumulate(ns.vxspace, i, factor*(
               param.vDis.eval(pf_f)*ns.vx*ns.vbasis.phi[i]

dune/phasefield/porenetwork/fem_base_porenetwork.hh

@@ -4,7 +4,7 @@
 
 
 template<typename RF, const int dim, typename NSLFS>
-class NSPorenetworkFlow
+class NSPorenetworkVariables
 {
 public:
   const typename NSLFS::template Child<0>::Type& vxspace;
@@ -14,7 +14,7 @@ public:
   Dune::FieldVector<RF,dim> gradvx;
 
   template<typename Pos, typename vCache, typename LV, typename S>
-  NSPorenetworkFlow(const vCache& vcache, const Pos& position, const NSLFS& nsLfs, const LV& lv, S& s) :
+  NSPorenetworkVariables(const vCache& vcache, const Pos& position, const NSLFS& nsLfs, const LV& lv, S& s) :
     vxspace(child(nsLfs,Dune::Indices::_0)),
     vbasis(vcache, position, vxspace, s),
     vx(0.0),

dune/phasefield/porenetwork/pn_1pflow.hh

@@ -10,7 +10,7 @@
 #include<dune/pdelab/newton/newton.hh>
 
 #include <dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh>
-#include "tmixed_chns_navierstokes.hh"
+#include <dune/phasefield/porenetwork/fem_1p_navierstokes.hh>
 #include "tmixed_upscaledp.hh"
 #include "tmixed_upscaledc.hh"
 #include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
@@ -64,7 +64,7 @@ public:
                     Dune::PDELab::NoDirichletConstraintsParameters,
                     Dune::PDELab::NoDirichletConstraintsParameters> Bdry;
 
-  typedef TMIXED_CHNS_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
+  typedef FEM_1P_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
   typedef FEM_1P_CahnHilliard_R<Parameters, CH_GFS, CH_V> CH_LOP;
   typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
 
@@ -79,7 +79,7 @@ public:
   typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> NS_Newton;
   typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> CH_Newton;
 
-
+  Parameters& param;
   std::shared_ptr<Grid> gridp;
   GV gv;
   ES es;
@@ -94,17 +94,10 @@ public:
   Initial initial;
   Bdry bdry;
 
-
-  RF dxp = 1;
-  RF c = 1;
   MBE mbe;
 
   RF kf = 0;
-  RF phic = 0;
-  RF phicOld = 0;
   RF phiSolid = 0;
-  RF phiSolidOld = 0;
-  RF kc = 0;
 
   std::unique_ptr<NS_LOP> nsLop;
   std::unique_ptr<CH_LOP> chLop;
@@ -116,7 +109,8 @@ public:
   std::unique_ptr<NS_Newton>  nsNewton;
   std::unique_ptr<CH_Newton>  chNewton;
 
-  Pn_1PFlow(std::string filename, Parameters& param) :
+  Pn_1PFlow(std::string filename, Parameters& param_) :
+  param(param_),
   gridp(initPnGrid<Grid>(filename)),
   gv(gridp->leafGridView()),
   es(gv),
@@ -137,9 +131,10 @@ public:
     chLop = std::make_unique<CH_LOP>(param, gfs.ch, chVOld);
 
     dgf = std::make_unique<DGF>(gfs,nsV,chV);
+    initializeGridOperators();
   }
 
-  void initializeGridOperators(Parameters& param)
+  void initializeGridOperators()
   {
     //Grid Operator
     nsGo = std::make_unique<NS_GO>(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,*nsLop,mbe);
@@ -177,7 +172,7 @@ public:
     phicOld = phic;
     phiSolidOld = phiSolid;
     int verb=0;
-    if (verb > 1)   std::cout << "= Begin newtonapply CH:" << " with dxp= " << dxp << " dxpUpwind= " << dxpUpwind << " c= " << c << std::endl;
+    if (verb > 1)   std::cout << "= Begin newtonapply CH:" << std::endl;
     chNewton->apply();
     chVOld = chV;
     updatePhiC();
@@ -191,377 +186,29 @@ public:
         dgfMicro->vxDgf.evaluate(element, xlocal, w);
         Dune::FieldVector<double,1> phi1;
         dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
-        Dune::FieldVector<double,1> phi2;
-        dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
-        return Dune::FieldVector<double,1>((phi1+phi2)*w);
+        return Dune::FieldVector<double,1>(phi1*w);
+        //Dune::FieldVector<double,1> phi2;
+        //dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
+        //return Dune::FieldVector<double,1>((phi1+phi2)*w);
     });
     Dune::FieldVector<double,1> integral;
     Dune::PDELab::integrateGridFunction(productFunction,integral,2);
     kf = integral[0];
   }
 
-  void updatePhiC()
-  {
-    auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
-        Dune::FieldVector<double,1> phi1;
-        dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
-        return Dune::FieldVector<double,1>(phi1);
-    });
-    Dune::FieldVector<double,1> integral;
-    Dune::PDELab::integrateGridFunction(productFunction,integral,2);
-    phic = integral[0];
-  }
-
   void updatePhiSolid()
   {
     auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
       Dune::FieldVector<double,1> phi1;
       dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
-      Dune::FieldVector<double,1> phi2;
-      dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
-      return Dune::FieldVector<double,1>(1-phi1-phi2);
+      return Dune::FieldVector<double,1>(1-phi1);
+      //Dune::FieldVector<double,1> phi2;
+      //dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
+      //return Dune::FieldVector<double,1>(1-phi1-phi2);
     });
     Dune::FieldVector<double,1> integral;
     Dune::PDELab::integrateGridFunction(productFunction,integral,2);
     phiSolid = integral[0];
   }
 
-  void updateKc()
-  {
-    auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
-        Dune::FieldVector<double,1> w;
-        dgfMicro->vxDgf.evaluate(element, xlocal, w);
-        Dune::FieldVector<double,1> phi1;
-        dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
-        return Dune::FieldVector<double,1>(phi1*w);
-    });
-    Dune::FieldVector<double,1> integral;
-    Dune::PDELab::integrateGridFunction(productFunction,integral,2);
-    kc = integral[0];
-  }
-};
-
-template<typename RF, typename YSimS>
-class P_ParamAdaptor
-{
-public:
-  const int xs;
-  YSimS& ySimS;
-  P_ParamAdaptor(YSimS& ySimS_) : ySimS(ySimS_), xs(ySimS_.size())
-  {
-  }
-
-  RF evalKf(const auto& xg)
-  {
-    int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
-    return ySimS[i].kf;
-  }
-
-  RF evalKc(const auto& xg)
-  {
-    int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
-    return ySimS[i].kc;
-  }
-
-  RF evalDxP(const auto& xg)
-  {
-    int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
-    return ySimS[i].dxp;
-  }
-
-
-  RF evalPhiC(const auto& xg)
-  {
-    int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
-    return ySimS[i].phic;
-  }
-
-  RF evalPhiCOld(const auto& xg)
-  {
-    int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
-    return ySimS[i].phicOld;
-  }
-
-  RF evalPhiSolid(const auto& xg)
-  {
-    int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
-    return ySimS[i].phiSolid;
-  }
-
-  RF evalPhiSolidOld(const auto& xg)
-  {
-    int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
-    return ySimS[i].phiSolidOld;
-  }
 };
-
-//_________________________________________________________________________________________________________________________
-
-template< typename XGV, typename PGrad_DGF, typename SimS>
-void evaluateDxP(XGV& xgv, PGrad_DGF& pGradDgf, SimS& simS)
-{
-  Dune::FieldVector<double,1> pgrad(0);
-  Dune::FieldVector<double,1> pgradUpwind(0);
-  int xs = simS.size();
-  auto it = elements(xgv).begin();
-  for(int i=0; i<xs; i++)
-  {
-    pGradDgf.evaluate(*it,it->geometry().local(it->geometry().center()), pgrad );
-    pgrad[0] = pgrad[0] - 1; //TODO make more general with pressure drop/distance
-    simS[i].dxp = pgrad[0];
-    simS[i].dxpUpwind = pgradUpwind[0];
-    it++;
-    pgradUpwind = pgrad;
-  }
-  //Periodic:
-  simS[0].dxpUpwind = simS[xs-1].dxp;
-}
-
-template< typename XGV, typename C_DGF, typename SimS>
-void evaluateC(XGV& xgv, C_DGF& cDgf, SimS& simS)
-{
-  Dune::FieldVector<double,1> cTmp(0);
-  int xs = simS.size();
-  auto it = elements(xgv).begin();
-  for(int i=0; i<xs; i++)
-  {
-    cDgf.evaluate(*it,it->geometry().local(it->geometry().center()), cTmp );
-    simS[i].c = cTmp[0];
-    it++;
-  }
-}
-
-//_________________________________________________________________________________________________________________________
-
-template<typename XGV, typename YGV, typename X_ES, typename Y_ES, typename TGV, typename T_FEM, typename V_FEM, typename P_FEM, typename PHI_FEM, typename C_FEM>
-void driver_mixed(XGV& xgv, X_ES& x_es, const P_FEM& pFem, const C_FEM& cFem, const int xsize,
-                  YGV& ygv, Y_ES& y_es, TGV& tgv, const T_FEM& tFem, const V_FEM& vFem, const PHI_FEM& phiFem,
-                  Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
-{
-
-  const int dim = 1;
-  typedef double RF;
-
-
-  typedef Params_tmixed< RF,
-                        //TripleWell<RF,DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> > >,
-                        TripleWell<RF,DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> > >,
-                        Mobility_Const<RF>,
-                        Reaction_Thinstrip<RF>,
-                        VelDissipation_Quadratic<RF>
-                        //VelDissipation_Quadratic_Shifted<RF>
-                        > Parameters;
-  Parameters param(ptree);
-  TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
-
-
-//_________________________________________________________________________________________________________________________
-// Make micro problems
-
-
-  typedef GFS_micro<RF, Y_ES, V_FEM, PHI_FEM> GFS;
-  GFS gfs(y_es,vFem,phiFem);
-
-  typedef YSlice<YGV, GFS, Parameters> YSim;
-  typedef std::vector<YSim> YSimS;
-  YSimS ySimS;
-  ySimS.reserve(xsize);
-  for(int i=0; i<xsize; i++)
-  {
-    ySimS.emplace_back(ygv, gfs, param, ((RF)i)/((RF)xsize), ((RF)1)/((RF)xsize) );
-  }
-
-
-  //_________________________________________________________________________________________________________________________
-  // Make macro problems
-  typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
-  typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
-  typedef Dune::PDELab::NoConstraints NCon;
-  ConstraintsAssembler conP();
-  typedef Dune::PDELab::GridFunctionSpace<X_ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
-  P_GFS pGfs(x_es,pFem);
-
-  auto blambda = [](const auto& x){return true;};
-  auto b = Dune::PDELab::
-    makeBoundaryConditionFromCallable(xgv,blambda);
-  //typedef Dune::PDELab::DirichletConstraintsParameters BdryP;
-  //BdryP pBdry();
-  typedef typename P_GFS::template ConstraintsContainer<RF>::Type P_CC;
-  P_CC pCC;
-  //pCC.clear();
-  Dune::PDELab::constraints(b,pGfs,pCC);
-
-  typedef Dune::PDELab::Backend::Vector<P_GFS, RF> P_V;
-  P_V pV(pGfs);
-
-  auto p0lambda = [&](const auto& xg)
-    {return (1-xg[0])*(1-xg[0]);};
-  auto p0 = Dune::PDELab::
-      makeGridFunctionFromCallable(xgv,p0lambda);
-  Dune::PDELab::interpolate(p0,pGfs,pV);
-
-
-
-
-  typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
-  MBE mbe(10);
-  typedef P_ParamAdaptor<RF, YSimS> P_PAdaptor;
-  P_PAdaptor pPAdaptor(ySimS);
-  typedef TMIXED_UpscaledP<P_PAdaptor,P_FEM,RF> P_LOP;
-  P_LOP pLop(pPAdaptor);
-
-  typedef Dune::PDELab::GridOperator<P_GFS, P_GFS, P_LOP, MBE,
-                                     RF, RF, RF, P_CC, P_CC> P_GO;
-  P_GO pGo(pGfs, pCC, pGfs, pCC, pLop, mbe);
-  typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<P_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> P_LS;
-  P_LS pLs(pGo,500,100,3,1);
-  typedef Dune::PDELab::Newton<P_GO,P_LS,P_V> P_Newton;
-  P_Newton pNewton(pGo,pV,pLs);
-  pNewton.setReassembleThreshold(0.0); // always reassemble J
-  pNewton.setVerbosityLevel(1);        // be verbose
-  pNewton.setReduction(1e-10);         // total reduction
-  pNewton.setMinLinearReduction(1e-4); // min. red. in lin. solve
-  pNewton.setMaxIterations(25);        // limit number of its
-  pNewton.setLineSearchMaxIterations(10); // limit line search
-  pNewton.setAbsoluteLimit(1e-10);
-
-
-
-
-
-
-  typedef Dune::PDELab::GridFunctionSpace<X_ES, C_FEM,Dune::PDELab::NoConstraints,VectorBackend> C_GFS;
-  C_GFS cGfs(x_es,cFem);
-  typedef Dune::PDELab::Backend::Vector<C_GFS, RF> C_V;
-  C_V cV(cGfs);
-  C_V cVOld(cGfs);
-  auto c0lambda = [&](const auto& xg)
-    //{return 0.5+0.4*std::sin(4*3.14152*xg[0]);};
-    {return 0.5;};
-  auto c0 = Dune::PDELab::
-      makeGridFunctionFromCallable(xgv,c0lambda);
-  Dune::PDELab::interpolate(c0,cGfs,cV);
-  cVOld = cV;
-  //typedef P_ParamAdaptor<RF, YSimS> P_PAdaptor;
-  //P_PAdaptor pPAdaptor(ySimS);
-  typedef TMIXED_UpscaledC<Parameters,P_PAdaptor,C_GFS,C_V,RF> C_LOP;
-  C_LOP cLop(param,pPAdaptor,cGfs,cVOld);
-
-  typedef Dune::PDELab::EmptyTransformation NoCC;
-  //typedef typename C_GFS::template ConstraintsContainer<RF>::Type C_CC;
-  //C_CC cCC;
-
-  typedef Dune::PDELab::GridOperator<C_GFS, C_GFS, C_LOP, MBE,
-                                     RF, RF, RF, NoCC, NoCC> C_GO;
-  MBE mbe2(10);
-  C_GO cGo(cGfs, cGfs, cLop, mbe);
-  typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<C_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> C_LS;
-  C_LS cLs(cGo,500,100,3,1);
-  //typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<C_GO> C_LS;
-  //C_LS cLs(100,2);
-
-  typedef Dune::PDELab::Newton<C_GO,C_LS,C_V> C_Newton;
-  C_Newton cNewton(cGo,cV,cLs);
-  cNewton.setReassembleThreshold(0.0); // always reassemble J
-  cNewton.setVerbosityLevel(1);        // be verbose
-  cNewton.setReduction(1e-10);         // total reduction
-  cNewton.setMinLinearReduction(1e-4); // min. red. in lin. solve
-  cNewton.setMaxIterations(25);        // limit number of its
-  cNewton.setLineSearchMaxIterations(10); // limit line search
-  cNewton.setAbsoluteLimit(1e-10);
-
-
-
-
-
-
-  //_________________________________________________________________________________________________________________________
-  // prepare VTK writer
-  std::string filename=ptree.get("output.filename","output");
-  std::string x_filename = "xoutput";
-  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;
-  }
-  if( stat( x_filename.c_str(), &st ) != 0 )
-  {
-    int stat = 0;
-    stat = mkdir(x_filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
-    if( stat != 0 && stat != -1)
-      std::cout << "Error: Cannot create directory "  << x_filename << std::endl;
-  }
-  auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<TGV> >(tgv,0);
-  Dune::VTKSequenceWriter<TGV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
-
-  DGF_macro<RF, YGV, TGV, YSimS> dgfMacro(ygv, tgv,ySimS);
-  dgfMacro.addToVTKWriter(vtkwriter);
-  vtkwriter.write(param.time.t,Dune::VTK::base64);
-
-  auto x_stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<XGV> >(xgv,0);
-  Dune::VTKSequenceWriter<XGV> x_vtkwriter(x_stationaryvtkwriter,x_filename.c_str(),x_filename.c_str(),"");
-
-  typedef Dune::PDELab::DiscreteGridFunction<P_GFS,P_V> P_DGF;
-  P_DGF pDgf(pGfs,pV);
-  typedef Dune::PDELab::DiscreteGridFunctionGradient<P_GFS,P_V> PGrad_DGF;
-  PGrad_DGF pGradDgf(pGfs,pV);
-  typedef Dune::PDELab::DiscreteGridFunction<C_GFS,C_V> C_DGF;
-  C_DGF cDgf(cGfs,cV);
-
-  typedef Dune::PDELab::VTKGridFunctionAdapter<P_DGF> P_VTKF;
-  typedef Dune::PDELab::VTKGridFunctionAdapter<C_DGF> C_VTKF;
-  x_vtkwriter.addVertexData(std::shared_ptr<P_VTKF>(
-                             new P_VTKF(pDgf,"p")));
-  x_vtkwriter.addVertexData(std::shared_ptr<C_VTKF>(
-                             new C_VTKF(cDgf,"c")));
-  x_vtkwriter.write(param.time.t,Dune::VTK::base64);
-
-//_________________________________________________________________________________________________________________________
-while (!timestepManager.isFinished())
-{
-  for(int i=1; i<xsize; i++)
-    {
-      ySimS[i].chVUpwind = ySimS[i-1].chV;
-    }
-    ySimS[0].chVUpwind = ySimS[xsize-1].chV;  //Periodic
-
-    for(int i=0; i<xsize; i++)
-    {
-      ySimS[i].applyNs();
-    }
-
-    for(int i=1; i<xsize; i++)
-    {
-      ySimS[i].nsVUpwind = ySimS[i-1].nsV;
-    }
-    ySimS[0].nsVUpwind = ySimS[xsize-1].nsV;  //Periodic
-
-
-    //For Debug
-    //timestepManager.applyTimestep(vtkwriter);
-    //x_vtkwriter.write(param.time.t);
-
-
-    std::cout << "Begin newtonapply Pressure:" << std::endl;
-    pNewton.apply();
-    evaluateDxP(xgv, pGradDgf,ySimS);
-
-    for(int i=0; i<xsize; i++)
-    {
-      ySimS[i].applyCh();
-    }
-  //while (!timestepManager.isFinished())
-  //{
-    std::cout << "Begin newtonapply Concentration:" << std::endl;
-    cNewton.apply();
-    cVOld =cV;
-    evaluateC(xgv,cDgf,ySimS);
-
-
-    timestepManager.applyTimestep(vtkwriter, x_vtkwriter);
-  }
-
-}

src/CMakeLists.txt

@@ -14,3 +14,4 @@ add_subdirectory(example_nucleus)
 add_subdirectory(paper_thinstrip_full)
 add_subdirectory(paper_thinstrip_mixed)
 add_subdirectory(summerschool_inflow)
+add_subdirectory(pntest)

src/paper_thinstrip_mixed/tmixed_parameters.hh

@@ -9,15 +9,11 @@ class Params_tmixed : public Params_ffs_r<Number, TW, Mobility, Reaction, VelDis
 {
 public:
   const Number YScaling;
-  const Dune::ParameterTree& nsNewtonTree;
-  const Dune::ParameterTree& chNewtonTree;
   Number TempXpos = 0;
 
   Params_tmixed(const Dune::ParameterTree& ptree) :
     Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>(ptree),
-    YScaling(ptree.get("domain.YScaling",1.0)),
-    nsNewtonTree(ptree.sub("NS_Newton")),
-    chNewtonTree(ptree.sub("Phasefield_Newton"))
+    YScaling(ptree.get("domain.YScaling",1.0))
   {
     ModifyParameters();
     std::cout << "YScaling " << YScaling << " Eps " << this->eps << std::endl;

src/pntest/CMakeLists.txt

+add_executable(pntest pntest.cc)
+dune_symlink_to_source_files(FILES grids pn1p.ini)

src/pntest/chns.ini

+[Physics]
+uAst = 1
+D = 0.02
+rho1 = 1.0
+rho2 = 1.0
+mu=0.005
+ReactionRate=1.5
+#ReactionLimiter = 1000
+SurfaceTensionScaling=0.1
+VelDissipationRate = 100
+CurvatureAlpha=0#0.001
+
+[domain]
+level = 0
+filename = grids/square_periodic.msh
+#YScalingFactor = 0.5
+YScaling = 1
+YResolution = 200
+XResolution = 100
+
+[output]
+filename = output2
+
+[Phasefield]
+eps=0.03
+delta=0.00001#0.03
+DoubleWellDelta=0.05
+Mobility=1#0.01
+Sigma1=1
+Sigma2=1
+Sigma3=1
+
+[Time]
+tMax = 100
+dt = 0.0001
+dtMax = 0.0001
+dtMin = 0.000015
+saveintervall = 0.005
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+# FlowVyy= 10# pConst/(Domain_width * mu)
+pConst=0.1
+uConst=0.5
+uInflow=0.3
+
+[Refinement]
+etaRefine  = 0.2
+etaCoarsen = 0.1
+maxRefineLevel = 6
+minRefineLevel = -1
+
+[NS_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+AbsoluteLimit = 1e-9
+Reduction = 1e-7
+VerbosityLevel = 0
+
+[Phasefield_Newton]
+ReassembleThreshold = 0.0
+LineSearchMaxIterations = 30
+MaxIterations = 30
+VerbosityLevel = 0

src/pntest/pntest.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<dune/pdelab/finiteelementmap/qkfem.hh>
+#include<dune/pdelab/porenetwork/1p_chns.hh>
+
+//===============================================================
+// Main program
+//===============================================================
+
+
+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("pn1p.ini",ptree);
+    ptreeparser.readOptions(argc,argv,ptree);
+
+    typedef Params_fs_r< RF,
+                          //DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,
+                          DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> >,
+                          Mobility_Const<RF>,
+                          Reaction_Thinstrip<RF>,
+                          VelDissipation_Quadratic<RF>
+                          //VelDissipation_Quadratic_Shifted<RF>
+                          > Parameters;
+    Parameters param(ptree);
+    std::string filename = ptree.get("domain.filename","square.msh");
+    Pn_1PFlow<Parameters> pn_1pFlow(filename,param);
+    pn_1pFlow.applyCh();
+    pn_1pFlow.applyNs();
+    pn_1pFlow.updateKf();
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+}