diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index dd4b7e1340857f788438626cc16ebe3cda01b3ce..bb2490892ac8a5931059ee27ef31b106d093ced0 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/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"))
{
}
};
diff --git a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
index d6f9656ff41726467b8084de1b1cc74040a0a6c9..b5e86111ff2b494054aed7cf354db86eb5894431 100644
--- a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
+++ b/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]
diff --git a/dune/phasefield/porenetwork/fem_base_porenetwork.hh b/dune/phasefield/porenetwork/fem_base_porenetwork.hh
index 601757c01d5dad6338178ae3e25fd0d5c0a0810f..ceb0b62d253c5ef93ce2b22dd24630039f8cf301 100644
--- a/dune/phasefield/porenetwork/fem_base_porenetwork.hh
+++ b/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),
diff --git a/dune/phasefield/porenetwork/pn_1pflow.hh b/dune/phasefield/porenetwork/pn_1pflow.hh
index 71f32f12c91fd5091c3a063a4b55fb963847c2af..544cab16a36245344840f23bfb4540d075a35a92 100644
--- a/dune/phasefield/porenetwork/pn_1pflow.hh
+++ b/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);
- }
-
-}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 2fce0700447e78837b550cfd06e1a5912c2ef776..986d7a311fd1ff764814cfc56d26daa3b9ca4ccd 100644
--- a/src/CMakeLists.txt
+++ b/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)
diff --git a/src/paper_thinstrip_mixed/tmixed_parameters.hh b/src/paper_thinstrip_mixed/tmixed_parameters.hh
index 535347e8fbf3d89b3a45786f24cb4a62fadba95b..3fdeef307dfc3a681e836bf63000709df124951e 100644
--- a/src/paper_thinstrip_mixed/tmixed_parameters.hh
+++ b/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;
diff --git a/src/pntest/CMakeLists.txt b/src/pntest/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..baf654051d9d06d2d04855eca95213b15748abd1
--- /dev/null
+++ b/src/pntest/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(pntest pntest.cc)
+dune_symlink_to_source_files(FILES grids pn1p.ini)
diff --git a/src/pntest/chns.ini b/src/pntest/chns.ini
new file mode 100644
index 0000000000000000000000000000000000000000..c7addb30770b280cf2f3b29c269fe671b245b851
--- /dev/null
+++ b/src/pntest/chns.ini
@@ -0,0 +1,66 @@
+[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
diff --git a/src/pntest/pntest.cc b/src/pntest/pntest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2e6aab1bd65a923c0e7fef2444eb936ac9cedb94
--- /dev/null
+++ b/src/pntest/pntest.cc
@@ -0,0 +1,66 @@
+#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;
+ }
+}