diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 33d0dbe9b5e16d160fee706482d5792e555f2f8e..ed136822072a542abe33e8043b1ef2f1c5b0b316 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -54,13 +54,46 @@ public:
}
};
+template<typename Number>
+class Mobility_pf2p
+{
+public:
+ Number mobility;
+ Number mobilityPressureSaturation;
+ Mobility_pf2p(const Dune::ParameterTree & param) :
+ mobility(param.get("Mobility",(Number)0.1)),
+ mobilityPressureSaturation(param.get("MobilityPressureSaturation",(Number)0.1))
+ {
+ }
+
+ Number eval( const Number& pf) const
+ {
+ return mobility;
+ }
+
+ Number eval( const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return mobility;
+ }
+
+ Number evalPressureSaturation( const Number& pf) const
+ {
+ return mobilityPressureSaturation;
+ }
+
+ Number evalPressureSaturation( const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return mobilityPressureSaturation;
+ }
+};
+
//______________________________________________________________________________________________________
template<typename Number>
class Reaction_Const
{
public:
- const Number reactionRate;
+ Number reactionRate;
Reaction_Const(const Dune::ParameterTree & param) :
reactionRate(param.get("ReactionRate",(Number)1))
{
@@ -397,6 +430,34 @@ public:
{
}
};
+
+
+template <typename Number, typename TW, typename HTW, typename Mobility, typename Reaction, typename VelDissipation>
+class Params_pn_2pflow : public Params_base<Number, Mobility>
+{
+public:
+ Number delta;
+
+ TW tw;
+ HTW htw;
+ Reaction reaction;
+ VelDissipation vDis;
+
+ Number SurfaceTensionScaling;
+
+
+ Params_pn_2pflow(const Dune::ParameterTree& ptree) :
+ Params_base<Number, Mobility>(ptree),
+ delta(ptree.get("Phasefield.delta",this->eps)),
+ tw(ptree.sub("Phasefield")),
+ htw(ptree.sub("Phasefield")),
+ reaction(ptree.sub("Physics")),
+ vDis(ptree.sub("Physics")),
+ SurfaceTensionScaling(ptree.get("Physics.SurfaceTensionScaling",(Number)1))
+ {
+ }
+};
+
//
// template<typename Number>
// class Problem
diff --git a/dune/phasefield/base_potentials.hh b/dune/phasefield/base_potentials.hh
index dd199383184392109f798e20d70a8f0a4647ed20..b9f7d8667500c2759ed070473ec25a579e658225 100644
--- a/dune/phasefield/base_potentials.hh
+++ b/dune/phasefield/base_potentials.hh
@@ -364,3 +364,69 @@ public:
return dw.eval_q(pf1/(1-pf2))*(1-pf2)*(1-pf2); //A function being q on the 13 interface and 0 on the others
}
};
+
+template<typename Number>
+class TPPentalty_6thOrder
+{
+public:
+ TPPentalty_6thOrder(const Dune::ParameterTree & ptree)
+ {
+ }
+
+ Number eval(const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return pf1*pf1*pf2*pf2*pf3*pf3;
+ }
+
+ Number eval_dpf1(const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return 2*pf1*pf2*pf2*pf3*pf3;
+ }
+
+ Number eval_dpf2(const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return 2*pf1*pf1*pf2*pf3*pf3;
+ }
+
+ Number eval_dpf3(const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return 2*pf1*pf1*pf2*pf2*pf3;
+ }
+};
+
+template<typename Number, typename DW, typename TPPenalty>
+class TripleWellHierarchicalV1
+{
+public:
+ const Number Sigma1;
+ const Number Sigma2;
+ const Number Sigma3;
+ const Number SigmaT;
+ const Number Lambda;
+
+ DW dw;
+ TPPenalty tppen;
+
+ TripleWellHierarchicalV1(const Dune::ParameterTree & param) :
+ Sigma1(param.get("Sigma1",(Number)1)),
+ Sigma2(param.get("Sigma2",(Number)1)),
+ Sigma3(param.get("Sigma3",(Number)1)),
+ SigmaT(1/(1/Sigma1 + 1/Sigma2)),
+ Lambda(param.get("Lambda",(Number)1)),
+ dw(param),
+ tppen(param)
+ {
+ }
+
+ Number eval_dpf1 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (1 - SigmaT/Sigma1)*(Sigma1*dw.eval_d(pf1) + Lambda*tppen.eval_dpf1(pf1,pf2,pf3))
+ -SigmaT/Sigma2*(Sigma2*dw.eval_d(pf2) + Lambda*tppen.eval_dpf2(pf1,pf2,pf3));
+ }
+
+ Number eval_dpf2 (const Number& pf1, const Number& pf2, const Number& pf3) const
+ {
+ return (1 - SigmaT/Sigma2)*(Sigma2*dw.eval_d(pf2) + Lambda*tppen.eval_dpf2(pf1,pf2,pf3))
+ -SigmaT/Sigma1*(Sigma1*dw.eval_d(pf1) + Lambda*tppen.eval_dpf1(pf1,pf2,pf3));
+ }
+};
diff --git a/dune/phasefield/boundary_rectangular.hh b/dune/phasefield/boundary_rectangular.hh
index 3a7728ce8b2b4250298bba0f0da40143afa17ec5..4c6fa95528f00f8d3dba3bdc30f4f6dcd0233e28 100644
--- a/dune/phasefield/boundary_rectangular.hh
+++ b/dune/phasefield/boundary_rectangular.hh
@@ -146,3 +146,50 @@ public:
return false;
}
};
+
+
+template<typename DomainSize>
+class TopRightCorner : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ TopRightCorner(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (( xg[0] > ds.right - (ds.right-ds.left)/20) && ( xg[1]>ds.top - ds.domainEps))
+ {
+ return true;
+ }
+ return false;
+ }
+};
+
+template<typename DomainSize>
+class LeftBottomCorner : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ LeftBottomCorner(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (( xg[0] < ds.left + (ds.right-ds.left)/20) && ( xg[1]<ds.bottom + ds.domainEps))
+ {
+ return true;
+ }
+ return false;
+ }
+};
diff --git a/dune/phasefield/fff_chns.hh b/dune/phasefield/fff_chns.hh
index 332ff5d979ad3bd17bdc944fc9b3a1970455fb9a..9ab1c30a1ac915bd66fa17ae9e4d28e220030591 100644
--- a/dune/phasefield/fff_chns.hh
+++ b/dune/phasefield/fff_chns.hh
@@ -214,3 +214,38 @@ public:
{
}
};
+
+template<typename BdryV1, typename BdryV2, typename BdryP,
+ typename BdryPhi, typename BdryMu, typename BdryPhi2, typename BdryMu2>
+class Bdry_fff_chns_compositeV_2d_Explicit_Domain
+{
+public:
+ BdryV1 bdryV1;
+ BdryV2 bdryV2;
+ BdryP bdryP;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ BdryPhi2 bdryPhi2;
+ BdryMu2 bdryMu2;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV1,BdryV2> BdryV;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryV, BdryP> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu, BdryPhi2, BdryMu2> BdryCH;
+ BdryV bdryV;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ template<typename DomainSize>
+ Bdry_fff_chns_compositeV_2d_Explicit_Domain(DomainSize ds) :
+ bdryV1(ds),
+ bdryV2(ds),
+ bdryP(ds),
+ bdryPhi(ds),
+ bdryMu(ds),
+ bdryPhi2(ds),
+ bdryMu2(ds),
+ bdryV(bdryV1,bdryV2),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu,bdryPhi2,bdryMu2)
+ {
+ }
+};
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index 8b98df41364bf4f2a4fb045d6f3e651003977690..e9df86c4f249d4f8730bf77becc4fccc3e2f43c6 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -657,6 +657,86 @@ public:
};
+template<typename GV, typename RF, typename Params, const int Phase>
+class PhiInitial_Bubble_Spinodal :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Bubble_Spinodal<GV,RF,Params,Phase> >
+{
+ RF eps;
+ RF r;
+ RF mean;
+ RF deviation;
+ int seed;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_Bubble_Spinodal(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_Bubble_Spinodal<GV,RF,Params,Phase> > (gv)
+ {
+ mean = params.initial.ptree.get("Mean",0.5);
+ deviation = params.initial.ptree.get("Deviation",0.1);
+ seed = params.initial.ptree.get("Seed",1);
+ eps = params.eps;
+ r = params.initial.ptree.get("Radius",0.2);
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ RF dist = std::sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5));
+ RF rf = 0.5*(1+tanh(1/eps/sqrt(2)*(dist-r)));
+
+ int xint = (int) (x[0]*10000);
+ int yint = (int) (x[1]*10000);
+ std::srand(xint + yint*10061 + seed);
+ RF noise = mean + deviation*(((RF)std::rand())/RAND_MAX);
+ if(Phase == 0)
+ {
+ y[0] = (1-rf)*noise;
+ }
+ else
+ {
+ y[0] = (1-rf)*(1-noise);
+ }
+ }
+};
+
+
+template<typename GV, typename RF, typename Params, const int Phase>
+class PhiInitial_Bubble_Interior :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Bubble_Interior<GV,RF,Params,Phase> >
+{
+ RF eps;
+ RF r;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_Bubble_Interior(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_Bubble_Interior<GV,RF,Params,Phase> > (gv)
+ {
+ eps = params.eps;
+ r = params.initial.ptree.get("Radius",0.2);
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ RF dist = std::sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5));
+ RF rf = 0.5*(1+tanh(1/eps/sqrt(2)*(dist-r)));
+ RF fluid = 0.5*(1+tanh(1/eps/sqrt(2)*(x[1]-0.5)));
+
+ if(Phase == 0)
+ {
+ y[0] = (1-rf)*fluid;
+ }
+ else
+ {
+ y[0] = (1-rf)*(1-fluid);
+ }
+ }
+};
+
//______________________________________________________________________________________________________________________
#include <dune/phasefield/boundary_fn.hh>
diff --git a/dune/phasefield/localoperator/fem_3p_cahnhilliard_hierarchical_top.hh b/dune/phasefield/localoperator/fem_3p_cahnhilliard_hierarchical_top.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ff4e1944f06ec25c3480d56dc41a66df3158daa2
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_3p_cahnhilliard_hierarchical_top.hh
@@ -0,0 +1,97 @@
+#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>
+class FEM_3P_CahnHilliard_HierarchicalTop :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_3P_CahnHilliard_HierarchicalTop<Param, CHGFS, CHContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_3P_CahnHilliard_HierarchicalTop<Param, CHGFS, CHContainer> >
+{
+private:
+ typedef typename CHContainer::field_type RF;
+ typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ CHLocalBasisCache<CHGFS> chCache;
+
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chOldLocal;
+
+ Param& param; // parameter functions
+
+public:
+ enum {doPatternVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ FEM_3P_CahnHilliard_HierarchicalTop(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_) :
+ chOldLocal(chGfs_,chVOld_),
+ 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 = 4; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ chOldLocal.DoEvaluation(eg);
+
+ for(const auto& ip : rule)
+ {
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+ CHVariables_Threephase_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(), chOldLocal.lv, S);
+
+ RF pf3 = 1-ch.pf-ch.pf2;
+
+ for (size_t i=0; i<ch.pfspace.size(); i++)
+ {
+ RF gradxi3TimesGradphi = ch.gradxi*ch.basis.gradphi[i][0];
+ RF J3 = param.mobility.eval(ch.pf,ch.pf2,pf3) * param.eps/param.tw.Sigma3 * gradxi3TimesGradphi;
+ RF J1 = (-(ch.pf +param.delta/2)/(ch.pf+ch.pf2+param.delta))* J3 ;
+ RF J2 = (-(ch.pf2+param.delta/2)/(ch.pf+ch.pf2+param.delta))* J3;
+ // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
+ r.accumulate(ch.pfspace,i,factor*(
+ (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+ +J1
+ ));
+
+ r.accumulate(ch.pf2space,i,factor*(
+ (ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+ +J2
+ ));
+
+ r.accumulate(ch.xispace,i,factor*(
+ ch.xi*ch.basis.phi[i]
+ - 1/param.eps* param.tw.dw.eval_d(pf3) *ch.basis.phi[i]
+ + param.eps * param.tw.Sigma3* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.tw.Sigma3* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+
+ r.accumulate(ch.xi2space,i,factor*(
+ ch.xi2*ch.basis.phi[i]
+ //- 1/param.eps* param.tw.eval_dpf2(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+ //- param.eps * param.tw.Sigma2* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+ }
+ }
+ }
+
+};
diff --git a/dune/phasefield/localoperator/fem_3p_cahnhilliard_hierarchicalv1.hh b/dune/phasefield/localoperator/fem_3p_cahnhilliard_hierarchicalv1.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3d8c92fc8e2de010dd82237d15f6b27b9adb2972
--- /dev/null
+++ b/dune/phasefield/localoperator/fem_3p_cahnhilliard_hierarchicalv1.hh
@@ -0,0 +1,98 @@
+#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>
+class FEM_3P_CahnHilliard_HierarchicalV1 :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_3P_CahnHilliard_HierarchicalV1<Param, CHGFS, CHContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_3P_CahnHilliard_HierarchicalV1<Param, CHGFS, CHContainer> >
+{
+private:
+ typedef typename CHContainer::field_type RF;
+ typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ CHLocalBasisCache<CHGFS> chCache;
+
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chOldLocal;
+
+ Param& param; // parameter functions
+
+public:
+ enum {doPatternVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ FEM_3P_CahnHilliard_HierarchicalV1(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_) :
+ chOldLocal(chGfs_,chVOld_),
+ 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 = 4; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ chOldLocal.DoEvaluation(eg);
+
+ for(const auto& ip : rule)
+ {
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+ CHVariables_Threephase_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(), chOldLocal.lv, S);
+
+ RF pf3 = 1-ch.pf-ch.pf2;
+
+ for (size_t i=0; i<ch.pfspace.size(); i++)
+ {
+ RF J1 = param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.htw.Sigma1 * (ch.gradxi*ch.basis.gradphi[i][0]);
+
+ // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
+ r.accumulate(ch.pfspace,i,factor*(
+ (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+ +J1
+ ));
+
+ r.accumulate(ch.pf2space,i,factor*(
+ (ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+ +param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.htw.Sigma2 * (ch.gradxi2*ch.basis.gradphi[i][0])
+ ));
+
+ r.accumulate(ch.xispace,i,factor*(
+ ch.xi*ch.basis.phi[i]
+ - 1/param.eps* param.htw.eval_dpf1(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+ - param.eps * param.htw.Sigma1* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+
+ r.accumulate(ch.xi2space,i,factor*(
+ ch.xi2*ch.basis.phi[i]
+ - 1/param.eps* param.htw.eval_dpf2(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+ - param.eps * param.htw.Sigma2* (ch.gradpf2*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+ }
+ }
+ }
+
+};
diff --git a/dune/phasefield/porenetwork/fem_2p_cahnhilliard_precipitation.hh b/dune/phasefield/porenetwork/fem_2p_cahnhilliard_precipitation.hh
new file mode 100644
index 0000000000000000000000000000000000000000..54e23701d036e91e7840f0c1c72d917b5a2b608d
--- /dev/null
+++ b/dune/phasefield/porenetwork/fem_2p_cahnhilliard_precipitation.hh
@@ -0,0 +1,102 @@
+#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>
+class FEM_2P_CahnHilliard_Precipitation :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_2P_CahnHilliard_Precipitation<Param, CHGFS, CHContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_2P_CahnHilliard_Precipitation<Param, CHGFS, CHContainer> >
+{
+private:
+ typedef typename CHContainer::field_type RF;
+ typedef typename CHGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ CHLocalBasisCache<CHGFS> chCache;
+
+ typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
+ mutable CHLView chOldLocal;
+
+ Param& param; // parameter functions
+
+public:
+ enum {doPatternVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ FEM_2P_CahnHilliard_Precipitation(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_) :
+ chOldLocal(chGfs_,chVOld_),
+ 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 = 4; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ chOldLocal.DoEvaluation(eg);
+
+ for(const auto& ip : rule)
+ {
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+
+ CHVariables_Threephase_PFOld_Extension<RF,dim,LFSV> ch(chCache, ip.position(), lfsv, x.base(), chOldLocal.lv, S);
+
+ RF pf3 = 1-ch.pf-ch.pf2;
+ RF pf3Old = 1-ch.pfOld-ch.pf2Old;
+
+ for (size_t i=0; i<ch.pfspace.size(); i++)
+ {
+
+ RF Rspeed = ch.pfOld < 0.95 ? param.tw.eval_q(ch.pfOld,ch.pf2Old,pf3Old) : 0.0f;
+ RF R1 = -Rspeed/param.eps*param.reaction.eval(1)*ch.basis.phi[i];
+
+ RF gradxi3TimesGradphi = ch.gradxi*ch.basis.gradphi[i][0];
+ RF J3 = param.mobility.eval(ch.pf,ch.pf2,pf3) * param.eps/param.tw.Sigma3 * gradxi3TimesGradphi;
+ RF J1 = (-(ch.pf +param.delta/2)/(ch.pf+ch.pf2+param.delta))* J3 ;
+ RF J2 = (-(ch.pf2+param.delta/2)/(ch.pf+ch.pf2+param.delta))* J3;
+ // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
+ r.accumulate(ch.pfspace,i,factor*(
+ (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+ +J1-R1
+ ));
+
+ r.accumulate(ch.pf2space,i,factor*(
+ (ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+ +J2
+ ));
+
+ r.accumulate(ch.xispace,i,factor*(
+ ch.xi*ch.basis.phi[i]
+ - 1/param.eps* param.tw.dw.eval_d(pf3) *ch.basis.phi[i]
+ + param.eps * param.tw.Sigma3* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.tw.Sigma3* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+
+ r.accumulate(ch.xi2space,i,factor*(
+ ch.xi2*ch.basis.phi[i]
+ //- 1/param.eps* param.tw.eval_dpf2(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+ //- param.eps * param.tw.Sigma2* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+ }
+ }
+ }
+
+};
diff --git a/dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh b/dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d36cd10feb59412775ba3c174097131067d106f9
--- /dev/null
+++ b/dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh
@@ -0,0 +1,121 @@
+#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_2P_CahnHilliard_PressureSaturation :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_2P_CahnHilliard_PressureSaturation<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_2P_CahnHilliard_PressureSaturation<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 {doLambdaVolume=true};
+ enum {doAlphaVolume=true};
+
+ //constructor
+ FEM_2P_CahnHilliard_PressureSaturation(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 = 4; //TODO: Choose more educated
+ 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_Threephase_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 pf3 = 1-ch.pf-ch.pf2;
+ RF pf3Old = 1-ch.pfOld-ch.pf2Old;
+ //RF xi3 = -param.tw.Sigma3*(ch.xi/param.tw.Sigma1 + ch.xi2/param.tw.Sigma2);
+
+
+ for (size_t i=0; i<ch.pfspace.size(); i++)
+ {
+
+ RF Rspeed = ch.pfOld < 0.95 ? param.tw.eval_q(ch.pfOld,pf3Old,ch.pf2Old) : 0.0f;
+ RF R1 = -Rspeed/param.eps*param.reaction.eval(1)*ch.basis.phi[i];
+ RF J1 = param.mobility.evalPressureSaturation(ch.pf,ch.pf2,pf3)*param.eps/param.htw.Sigma1 * (ch.gradxi*ch.basis.gradphi[i][0]);
+ RF J2 = param.mobility.evalPressureSaturation(ch.pf,ch.pf2,pf3)*param.eps/param.htw.Sigma2 * (ch.gradxi2*ch.basis.gradphi[i][0]);
+
+ // integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
+ r.accumulate(ch.pfspace,i,factor*(
+ (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+ +J1
+ //+(ch.pf*ns.div_v + (ch.gradpf*ns.v))*ch.basis.phi[i]
+ -ch.pf*(ns.v*ch.basis.gradphi[i][0]) //weak transport
+ +R1
+ ));
+
+ r.accumulate(ch.pf2space,i,factor*(
+ (ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+ +J2
+ //+(ch.pf2*ns.div_v + (ch.gradpf2*ns.v))*ch.basis.phi[i]
+ -ch.pf2*(ns.v*ch.basis.gradphi[i][0])
+ -R1
+ ));
+
+ r.accumulate(ch.xispace,i,factor*(
+ ch.xi*ch.basis.phi[i]
+ - 1/param.eps* param.htw.eval_dpf1(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+ - param.eps * param.htw.Sigma1* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+
+ r.accumulate(ch.xi2space,i,factor*(
+ ch.xi2*ch.basis.phi[i]
+ - 1/param.eps* param.htw.eval_dpf2(ch.pf,ch.pf2,pf3) *ch.basis.phi[i]
+ - param.eps * param.htw.Sigma2* (ch.gradpf2*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf*ch.basis.gradphi[i][0])
+ + param.eps * param.htw.SigmaT* (ch.gradpf2*ch.basis.gradphi[i][0])
+ ));
+
+
+ }
+ }
+ }
+
+};
diff --git a/dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh b/dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7dda41197e79436c1a7b596fa0b05d10ade646fa
--- /dev/null
+++ b/dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh
@@ -0,0 +1,143 @@
+#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>
+
+
+
+
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+ class FEM_2P_NavierStokes_PressureSaturation :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_2P_NavierStokes_PressureSaturation<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_2P_NavierStokes_PressureSaturation<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
+
+ //FOR DETECTING ENTITY changes
+ typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
+ mutable ENT OldEnt;
+
+ public:
+
+
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+ //enum { doLambdaVolume = true };
+
+
+
+ FEM_2P_NavierStokes_PressureSaturation (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
+ 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_Threephase_Extension<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 pf3=1-ch.pf-ch.pf2;
+ RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;
+ Dune::FieldVector<RF,dim> gradpf_f(0.0);
+ gradpf_f.axpy((1-2*param.delta),ch.gradpf+ch.gradpf2);
+
+ RF rho_f=param.phys.rho_f(ch.pf,ch.pf2);
+ //RF rho_f_reg = rho_f + (param.phys.rho_1+param.phys.rho_2)*param.delta;
+
+ for (size_t i=0; i<ns.pspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ // integrate: (div v)*phi
+ //r.accumulate(ns.pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor); Weak form
+ 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++)
+ {
+ // integrate: ( (rho+rhoOld)/2 * v_n+1 - rhoOld * v_n)/dt * phi + nu * Dv_n+1 : D phi - p_n+1 div phi
+
+ //std::cout << (param.nu * (jacv[k]*gradvphi[i][0]) - p * gradvphi[i][0][k]
+ //- param.rho*g[k]*vphi[i] )*factor << std::endl;
+
+ //RF lambda = 3;
+ //const auto v_nabla_v = v * jacv[k];
+ RF gradxi3k = -param.tw.Sigma3*(ch.gradxi[k]/param.tw.Sigma1 + ch.gradxi2[k]/param.tw.Sigma2);
+ RF S = 0;
+ S += - ch.xi2*(ch.gradpf[k]*pf_f - ch.pf * gradpf_f[k])/(pf_f*pf_f); //TODO!!!!!!! Other pf_f scaling!!!!!
+ S += - ch.xi *(ch.gradpf2[k]*pf_f - ch.pf2 * gradpf_f[k])/(pf_f*pf_f);
+ S += - 2*param.delta*pf3*(gradxi3k-ch.gradxi2[k]-ch.gradxi[k]);
+
+ //Stronger form - flow in interfaces possible
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k] - 2*param.delta*pf3*gradxi3k;
+
+ //Works but is not the one in the theory (same as above, w/o regularization)
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k];
+
+ S = S*param.SurfaceTensionScaling*ns.vbasis.phi[i];
+
+ r.accumulate(child(ns.vspace,k),i, (// rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt //TimeEvolution
+ //- 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.eval_mu(pf_f) * (ns.jacv[k]*ns.vbasis.gradphi[i][0]) //+vphi[i]*(jacv[k]*gradpf) ) //Viscosity
+ + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i] //Velocity Dissipation in Solid
+ //- rho*g[k]*vphi[i] //Gravitation
+ - S //Surface Tension
+ )*factor);
+ } //for k
+ } // for i
+ } // for ip
+ } //alpha volume
+
+ };
diff --git a/dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh b/dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh
new file mode 100644
index 0000000000000000000000000000000000000000..eaa3364d03d74a11cf1b5f9c6a765464f5d6150f
--- /dev/null
+++ b/dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh
@@ -0,0 +1,148 @@
+#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>
+
+
+
+
+template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
+ class FEM_2P_NavierStokes_Transmissibility :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<FEM_2P_NavierStokes_Transmissibility<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<FEM_2P_NavierStokes_Transmissibility<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
+
+ //FOR DETECTING ENTITY changes
+ typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENT;
+ mutable ENT OldEnt;
+
+ public:
+
+
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+ //enum { doLambdaVolume = true };
+
+
+
+ FEM_2P_NavierStokes_Transmissibility (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
+ 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_Threephase_Extension<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 pf3=1-ch.pf-ch.pf2;
+ RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;
+ Dune::FieldVector<RF,dim> gradpf_f(0.0);
+ gradpf_f.axpy((1-2*param.delta),ch.gradpf+ch.gradpf2);
+
+ RF rho_f=param.phys.rho_f(ch.pf,ch.pf2);
+ RF rho_f_reg = rho_f + (param.phys.rho_1+param.phys.rho_2)*param.delta;
+
+ for (size_t i=0; i<ns.pspace.size(); i++)
+ {
+ for(int k=0; k<dim; k++)
+ {
+ // integrate: (div v)*phi
+ //r.accumulate(ns.pspace,i, -1.0* pf_f*(v*gradphi[i][0])*factor); Weak form
+ 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++)
+ {
+ // integrate: ( (rho+rhoOld)/2 * v_n+1 - rhoOld * v_n)/dt * phi + nu * Dv_n+1 : D phi - p_n+1 div phi
+
+ //std::cout << (param.nu * (jacv[k]*gradvphi[i][0]) - p * gradvphi[i][0][k]
+ //- param.rho*g[k]*vphi[i] )*factor << std::endl;
+
+ //RF lambda = 3;
+ //const auto v_nabla_v = v * jacv[k];
+ RF gradxi3k = -param.tw.Sigma3*(ch.gradxi[k]/param.tw.Sigma1 + ch.gradxi2[k]/param.tw.Sigma2);
+ RF S = 0;
+ S += - ch.xi2*(ch.gradpf[k]*pf_f - ch.pf * gradpf_f[k])/(pf_f*pf_f); //TODO!!!!!!! Other pf_f scaling!!!!!
+ S += - ch.xi *(ch.gradpf2[k]*pf_f - ch.pf2 * gradpf_f[k])/(pf_f*pf_f);
+ S += - 2*param.delta*pf3*(gradxi3k-ch.gradxi2[k]-ch.gradxi[k]);
+
+ //Stronger form - flow in interfaces possible
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k] - 2*param.delta*pf3*gradxi3k;
+
+ //Works but is not the one in the theory (same as above, w/o regularization)
+ //S += -ch.pf*ch.gradxi[k] - ch.pf2*ch.gradxi2[k];
+
+ S = S*param.SurfaceTensionScaling*ns.vbasis.phi[i];
+
+ r.accumulate(child(ns.vspace,k),i, ( rho_f_reg*(ns.v[k]-ns.vOld[k])*ns.vbasis.phi[i]/param.time.dt //TimeEvolution
+ //- 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.eval_mu(pf_f) * (ns.jacv[k]*ns.vbasis.gradphi[i][0]) //+vphi[i]*(jacv[k]*gradpf) ) //Viscosity
+ + param.vDis.eval(pf_f)*ns.v[k]*ns.vbasis.phi[i] //Velocity Dissipation in Solid
+ //- rho*g[k]*vphi[i] //Gravitation
+ - S //Surface Tension
+ )*factor);
+ if(k==0)
+ {
+ r.accumulate(child(ns.vspace,k),i, ( pf_f * (param.phys.outerMomentumXForce * ns.vbasis.phi[i]) // Momentum Forcing Term
+ )*factor);
+ }
+ } //for k
+ } // for i
+ } // for ip
+ } //alpha volume
+
+ };
diff --git a/dune/phasefield/porenetwork/pn_2pflow.hh b/dune/phasefield/porenetwork/pn_2pflow.hh
new file mode 100644
index 0000000000000000000000000000000000000000..02fc246a89ff3b6e3765a994fea8f22d00634680
--- /dev/null
+++ b/dune/phasefield/porenetwork/pn_2pflow.hh
@@ -0,0 +1,449 @@
+#pragma once
+
+// C includes
+#include<sys/stat.h>
+// C++ includes
+#include<iostream>
+#include<vector>
+// Dune includes
+#include<dune/grid/io/file/vtk.hh>
+#include<dune/pdelab/newton/newton.hh>
+#include<dune/pdelab/common/functionutilities.hh>
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+#include<dune/pdelab/function/callableadapter.hh>
+
+#include <dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh>
+#include <dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh>
+#include <dune/phasefield/porenetwork/fem_2p_cahnhilliard_precipitation.hh>
+#include <dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh>
+#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
+
+#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
+
+#include <dune/phasefield/boundary_rectangular.hh>
+#include <dune/phasefield/porenetwork/pn_initial.hh>
+#include <dune/phasefield/chns_base.hh>
+#include <dune/phasefield/fff_chns.hh>
+#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/debug/vector_debug.hh>
+
+template< typename Grid>
+std::shared_ptr<Grid> initPnGrid(std::string filename)
+{
+ Dune::GridFactory<Grid> factory;
+ Dune::GmshReader<Grid>::read(factory,filename,true,true);
+ return std::shared_ptr<Grid>(factory.createGrid());
+}
+
+
+template<typename Parameters>
+class Pn_2PFlow{
+private:
+
+ typedef double RF;
+ static const int dim = 2;
+
+ typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
+ typedef Grid::ctype DF;
+ typedef Grid::LeafGridView GV;
+ typedef Dune::PDELab::AllEntitySet<GV> ES;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,2> V_FEM;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> P_FEM;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> PHI_FEM;
+
+ typedef GFS_fff_chns<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
+
+ typedef typename GFS::NS NS_GFS;
+ typedef typename GFS::CH CH_GFS;
+ typedef Dune::PDELab::Backend::Vector<CH_GFS, RF> CH_V;
+ typedef Dune::PDELab::Backend::Vector<NS_GFS, RF> NS_V;
+
+ typedef Ini_fff_chns< ZeroInitialVector<GV,RF, dim, Parameters>, //V
+ ZeroInitial<GV,RF,Parameters>, //P
+ PhiInitial_2p_Quartersquare<GV,RF,Parameters,0>,
+ ZeroInitial<GV,RF,Parameters>,
+ PhiInitial_2p_Quartersquare<GV,RF,Parameters,1>,
+ ZeroInitial<GV,RF,Parameters>> Initial;
+ typedef DomainSize<RF> DS;
+ typedef Bdry_fff_chns_compositeV_2d_Explicit_Domain<AllBoundary<DS>, //V
+ AllBoundary<DS>, //V
+ TopRightCorner<DS>, //P
+ NoBoundary<DS>,
+ NoBoundary<DS>,
+ NoBoundary<DS>,
+ NoBoundary<DS>> Bdry;
+
+ typedef FEM_2P_NavierStokes_PressureSaturation<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_PS_LOP;
+ typedef FEM_2P_CahnHilliard_PressureSaturation<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_PS_LOP;
+ typedef FEM_2P_CahnHilliard_Precipitation<Parameters, CH_GFS, CH_V> CH_R_LOP;
+ typedef FEM_2P_NavierStokes_Transmissibility<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_T_LOP;
+
+ typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+
+ typedef DGF_fff_chns<GFS, NS_V, CH_V> DGF;
+
+ typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_PS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_PS_GO;
+ typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_PS_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_PS_GO;
+ typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_R_LOP, MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_R_GO;
+ typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_T_LOP, MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_T_GO;
+
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<NS_PS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_PS_LS;
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<CH_PS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_PS_LS;
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<CH_R_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_R_LS;
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<NS_T_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_T_LS;
+
+ typedef Dune::PDELab::Newton<NS_PS_GO,NS_PS_LS,NS_V> NS_PS_Newton;
+ typedef Dune::PDELab::Newton<CH_PS_GO,CH_PS_LS,CH_V> CH_PS_Newton;
+ typedef Dune::PDELab::Newton<CH_R_GO,CH_R_LS,CH_V> CH_R_Newton;
+ typedef Dune::PDELab::Newton<NS_T_GO,NS_T_LS,NS_V> NS_T_Newton;
+
+ Parameters param;
+ std::shared_ptr<Grid> gridp;
+ GV gv;
+ ES es;
+ V_FEM vFem;
+ P_FEM pFem;
+ PHI_FEM phiFem;
+
+ GFS gfs;
+
+ CH_V chV, chVOld;
+ NS_V nsV, nsVOld;
+ Initial initial;
+ DS ds;
+ Bdry bdry;
+
+ MBE mbe;
+
+ RF kf = 0;
+ RF phiSolid = 0;
+
+ RF PoreReferenceLength;
+ RF initialPhiSolid = 0;
+
+ std::unique_ptr<NS_PS_LOP> nspsLop;
+ std::unique_ptr<CH_PS_LOP> chpsLop;
+ std::unique_ptr<CH_R_LOP> chrLop;
+ std::unique_ptr<NS_T_LOP> nstLop;
+ std::unique_ptr<DGF> dgf;
+ std::unique_ptr<NS_PS_GO> nspsGo;
+ std::unique_ptr<CH_PS_GO> chpsGo;
+ std::unique_ptr<CH_R_GO> chrGo;
+ std::unique_ptr<NS_T_GO> nstGo;
+
+ std::unique_ptr<NS_PS_LS> nspsLs;
+ std::unique_ptr<CH_PS_LS> chpsLs;
+ std::unique_ptr<CH_R_LS> chrLs;
+ std::unique_ptr<NS_T_LS> nstLs;
+
+ std::unique_ptr<NS_PS_Newton> nspsNewton;
+ std::unique_ptr<CH_PS_Newton> chpsNewton;
+ std::unique_ptr<CH_R_Newton> chrNewton;
+ std::unique_ptr<NS_T_Newton> nstNewton;
+ std::unique_ptr<Dune::VTKSequenceWriter<GV>> vtkwriter;
+
+ //------------------------------------------------------------------------------------
+
+
+ void initializeGridOperators()
+ {
+ //Grid Operator
+ nspsGo = std::make_unique<NS_PS_GO>(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,*nspsLop,mbe);
+ chpsGo = std::make_unique<CH_PS_GO>(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,*chpsLop,mbe);
+ chrGo = std::make_unique<CH_R_GO> (gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,*chrLop,mbe);
+ nstGo = std::make_unique<NS_T_GO> (gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,*nstLop,mbe);
+ // Linear solver
+ nspsLs = std::make_unique<NS_PS_LS>(*nspsGo,1000,200,3,0);
+ chpsLs = std::make_unique<CH_PS_LS>(*chpsGo,500,100,3,0);
+ chrLs = std::make_unique<CH_R_LS>(*chrGo,500,100,3,0);
+ nstLs = std::make_unique<NS_T_LS>(*nstGo,1000,200,3,0);
+ // Nonlinear solver
+ nspsNewton = std::make_unique<NS_PS_Newton>(*nspsGo,nsV,*nspsLs);
+ nspsNewton->setParameters(param.nsNewtonTree);
+
+ chpsNewton = std::make_unique<CH_PS_Newton>(*chpsGo,chV,*chpsLs);
+ chpsNewton->setParameters(param.chNewtonTree);
+
+ chrNewton = std::make_unique<CH_R_Newton>(*chrGo,chV,*chrLs);
+ chrNewton->setParameters(param.chNewtonTree);
+
+ nstNewton = std::make_unique<NS_T_Newton>(*nstGo,nsV,*nstLs);
+ nstNewton->setParameters(param.nsNewtonTree);
+
+ //printVector(chV,10,"chV");
+ //printVector(nsV,10,"nsV");
+ }
+
+ void remesh()
+ {
+ // mark elements for refinement
+ typedef RefinementIndicator<GV,RF,PF_3P_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+ RIND rInd(gv,gfs.ch,param.pTree.sub("Refinement"));
+ rInd.calculateIndicator(gv,chV,false);
+ rInd.markGrid(*gridp,2);
+ rInd.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,0);
+ // recompute constraints
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+ initializeGridOperators();
+ }
+
+ void applyNsPressureSaturation (int verbosity = 0)
+ {
+ nsVOld = nsV;
+ if (verbosity > 0) std::cout << "= Begin newtonapply NS_PS:" << std::endl;
+ nspsNewton->apply();
+ }
+
+ void applyChPressureSaturation(int verbosity = 0)
+ {
+ chVOld = chV;
+ if (verbosity > 0) std::cout << "= Begin newtonapply CH_PS:" << std::endl;
+ chpsNewton->apply();
+ //chVOld = chV;
+ remesh();
+ }
+
+ void applyNsTransmissibility(int verbosity = 0)
+ {
+ if (verbosity > 0) std::cout << "= Begin newtonapply NS_T:" << std::endl;
+ nstNewton->apply();
+ }
+
+ void applyChReaction(int verbosity = 0)
+ {
+ chVOld = chV;
+ if (verbosity > 0) std::cout << "= Begin newtonapply CH_R:" << std::endl;
+ chrNewton->apply();
+ //chVOld = chV;
+ remesh();
+ }
+
+/* void updateKf()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> w;
+ dgf->vxDgf.evaluate(element, xlocal, w);
+ Dune::FieldVector<double,1> phi1;
+ dgf->phiDgf.evaluate(element, xlocal, phi1);
+ return Dune::FieldVector<double,1>(phi1*w);
+ //Dune::FieldVector<double,1> phi2;
+ //dgf->phi2Dgf.evaluate(element, xlocal, phi2);
+ //return Dune::FieldVector<double,1>((phi1+phi2)*w);
+ });
+ Dune::FieldVector<double,1> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ // Calculation of kf: 4* because of Quarter-Circle,
+ // Scaling with x_ref^4
+ kf = 4 * integral[0] * PoreReferenceLength*PoreReferenceLength*PoreReferenceLength*PoreReferenceLength;
+ }
+
+ void updatePhiSolid()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> phi1;
+ dgf->phiDgf.evaluate(element, xlocal, phi1);
+ return Dune::FieldVector<double,1>(1-phi1);
+ //Dune::FieldVector<double,1> phi2;
+ //dgf->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];
+ }
+
+ RF calculateGrowthSpeed()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> phi1;
+ dgf->phiDgf.evaluate(element, xlocal, phi1);
+
+ RF Rspeed = phi1[0] < 0.95 ? param.dw.eval_q(phi1[0]) : 0.0f;
+ RF R1 = -Rspeed/param.eps*param.reaction.eval(1);
+ return Dune::FieldVector<double,1>(-R1);
+ });
+ Dune::FieldVector<double,1> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ return integral[0];
+ }
+
+*/
+
+//------------------------------------------------------------------------------------
+
+public:
+
+ //Constructor
+ Pn_2PFlow(std::string GridFilename, std::string OutputFilename, Parameters param_, RF PoreBaseRadius_ = 10e-6) :
+ param(param_),
+ gridp(initPnGrid<Grid>(GridFilename)),
+ gv(gridp->leafGridView()),
+ es(gv),
+ vFem(es), pFem(es), phiFem(es),
+ gfs(es,vFem,pFem,phiFem),
+ chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch), nsVOld(gfs.ns),
+ initial(gv,param), ds(0,1,0,1,1e-6), bdry(ds), mbe(10),
+ PoreReferenceLength(PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) ) )
+ {
+ Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+ Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+ // Refine Intial Values
+ for(int i=0; i<param.pTree.get("Refinement.maxRefineLevel",2); i++)
+ {
+ // mark elements for refinement
+ typedef RefinementIndicator<GV,RF,PF_3P_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND0;
+ RIND0 rInd0(gv,gfs.ch,param.pTree.sub("Refinement"));
+ rInd0.calculateIndicator(gv,chV,false);
+ rInd0.markGrid(*gridp,2);
+ rInd0.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,1);
+ //rInd0.myRefineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1);
+ 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);
+
+ //chVOld = chV;
+
+ //Local Operator
+ nspsLop = std::make_unique<NS_PS_LOP>(param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
+ chpsLop = std::make_unique<CH_PS_LOP>(param, gfs.ch, chVOld, gfs.ns, nsV);
+ chrLop = std::make_unique<CH_R_LOP> (param, gfs.ch, chVOld);
+ nstLop = std::make_unique<NS_T_LOP> (param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
+
+ dgf = std::make_unique<DGF>(gfs,nsV,chV);
+
+ struct stat st;
+ if( stat( OutputFilename.c_str(), &st ) != 0 )
+ {
+ int stat = 0;
+ stat = mkdir(OutputFilename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+ if( stat != 0 && stat != -1)
+ std::cout << "Error: Cannot create directory " << OutputFilename << std::endl;
+ }
+ auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
+ vtkwriter = std::make_unique<Dune::VTKSequenceWriter<GV> >(stationaryvtkwriter,OutputFilename.c_str(),OutputFilename.c_str(),"");
+ dgf->addToVTKWriter(*vtkwriter);
+
+
+ initializeGridOperators();
+
+/* applyNs();
+ updatePhiSolid();
+ initialPhiSolid = phiSolid;
+ updateKf();
+*/
+
+ RF dtmax = param.pTree.get("Time.dtMax",0.01);
+ param.time.dt = dtmax;
+ for(int i=0; i<1000; i++)
+ {
+ writeVTK(i);
+ if(i < 10)
+ {
+ applyChPressureSaturation(1);
+ }
+ else if(i<40)
+ {
+ param.reaction.reactionRate = 0;
+ applyChPressureSaturation(1);
+ applyNsPressureSaturation(1);
+ }
+ else
+ {
+ param.reaction.reactionRate = 1;
+ applyChReaction(1);
+ }
+ }
+ }
+
+ void writeVTK(RF time)
+ {
+ vtkwriter->write(time,Dune::VTK::ascii);
+ }
+
+/*
+
+ void grow(RF length) // How much the the solid-fluid interface should grow in normal direction (in m)
+ {
+ length = length/PoreReferenceLength;
+ if(length <= 0)
+ {
+ std::cout << "Need positive length (dissolution not yet implemented)" << std::endl;
+ return;
+ }
+ RF dtmax = param.pTree.get("Time.dtMax",0.01);
+ while(length > dtmax)
+ {
+ //writeVTK(1-length);
+ length -= dtmax;
+ param.time.dt = dtmax;
+ applyCh(0);
+ }
+ param.time.dt = length;
+ applyCh(0);
+ applyNs(0);
+ updateKf();
+ updatePhiSolid();
+ }
+
+ void growArea(RF area) // How much Area the solid phase should grow (in m^2)
+ {
+ RF phiSolid_beforeGrowth = phiSolid;
+ area = area/(PoreReferenceLength*PoreReferenceLength*4); //4 because of Symmetry
+ if(area <= 0)
+ {
+ std::cout << "Need positive area (dissolution not yet implemented)" << std::endl;
+ return;
+ }
+ RF dtmax = param.pTree.get("Time.dtMax",0.01);
+
+ while(area - (phiSolid - phiSolid_beforeGrowth) > dtmax*calculateGrowthSpeed())
+ {
+ //writeVTK(1-length);
+ param.time.dt = dtmax;
+ applyCh(0);
+ updatePhiSolid();
+ }
+ param.time.dt = (area - (phiSolid - phiSolid_beforeGrowth)) / calculateGrowthSpeed();
+ applyCh(0);
+ applyNs(0);
+ updateKf();
+ updatePhiSolid();
+ }
+
+ RF getKf()
+ {
+ return kf;
+ }
+
+ RF getPhiSolid()
+ {
+ return phiSolid;
+ }
+
+ RF getPoreArea()
+ {
+ return 4*PoreReferenceLength*PoreReferenceLength*(1-phiSolid); //4 because of Symmetry
+ }
+
+ RF getPoreAreaFraction()
+ {
+ return (1-phiSolid)/(1-initialPhiSolid);
+ }
+
+ void setPoreBaseRadius(RF PoreBaseRadius_)
+ {
+ PoreReferenceLength = PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) );
+ }
+
+ RF getPoreReferenceLength()
+ {
+ return PoreReferenceLength;
+ }
+
+*/
+};
diff --git a/dune/phasefield/porenetwork/pn_initial.hh b/dune/phasefield/porenetwork/pn_initial.hh
index 2f8c77d89249944a0fa5a93bf23abbd09f960eb1..be66d9ed1dfeeb6afb767b9205ef84dede559ede 100644
--- a/dune/phasefield/porenetwork/pn_initial.hh
+++ b/dune/phasefield/porenetwork/pn_initial.hh
@@ -47,3 +47,61 @@ public:
}
}
};
+
+
+template<typename GV, typename RF, typename Params, const int Phase>
+class PhiInitial_2p_Quartersquare :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_2p_Quartersquare<GV,RF,Params,Phase> >
+{
+ const Params& param;
+ RF eps;
+ RF r;
+ int myshape;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_2p_Quartersquare(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_2p_Quartersquare<GV,RF,Params,Phase> > (gv),
+ param(params)
+ {
+ eps = params.eps;
+ r = params.initial.ptree.get("Radius",0.2);
+ std::cout << "EPS" << eps;
+ myshape = params.initial.shape;
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ RF Pore = 0;
+ if(myshape == INITIAL_SHAPE_SQUARE)
+ {
+ RF rf1 = param.tw.dw.eval_shape((x[0]-r)*1/eps);
+ RF rf2 = param.tw.dw.eval_shape((x[1]-r)*1/eps);
+ //RF rf1 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-r)));
+ //RF rf2 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[1]-r)));
+ Pore=rf1*rf2;
+ }
+ else if(myshape == INITIAL_SHAPE_CIRCLE)
+ {
+ RF rad = sqrt((x[0]-1)*(x[0]-1) + (x[1]-1)*(x[1]-1));
+ Pore = param.tw.dw.eval_shape(((1-r)-rad)*1/eps);
+ }
+ else
+ {
+ std::cout << "WARNING: pn_initial.hh cannot find shape" << myshape << std::endl;
+ }
+
+ RF Cornerflow = param.tw.dw.eval_shape((x[0]+x[1]-4*r)*1/eps);
+
+ if(Phase == 0)
+ {
+ y = Pore*(1-Cornerflow);
+ }
+ if(Phase == 1)
+ {
+ y = Pore*Cornerflow;
+ }
+ }
+};
diff --git a/src/pntest/CMakeLists.txt b/src/pntest/CMakeLists.txt
index baf654051d9d06d2d04855eca95213b15748abd1..05b01592c924f23d1f16e2d14b31b887687cc08b 100644
--- a/src/pntest/CMakeLists.txt
+++ b/src/pntest/CMakeLists.txt
@@ -1,2 +1,3 @@
-add_executable(pntest pntest.cc)
-dune_symlink_to_source_files(FILES grids pn1p.ini)
+#add_executable(pntest pntest.cc)
+add_executable(pn2ptest pn2ptest.cc)
+dune_symlink_to_source_files(FILES grids pn1p.ini pn2p.ini)
diff --git a/src/pntest/pn2p.ini b/src/pntest/pn2p.ini
new file mode 100644
index 0000000000000000000000000000000000000000..4a774f424c8c7972d56f99ff4346dd6a7200b230
--- /dev/null
+++ b/src/pntest/pn2p.ini
@@ -0,0 +1,58 @@
+[Physics]
+uAst = 1
+D = 0.02
+rho1 = 1.0
+rho2 = 1.0
+mu=0.05
+ReactionRate=1#1.5
+#ReactionLimiter = 1000
+SurfaceTensionScaling=1
+VelDissipationRate = 4000
+CurvatureAlpha=0#0.001
+
+[domain]
+filename = grids/square_periodic.msh
+
+[output]
+filename = outp3
+
+[Phasefield]
+eps=0.02
+delta=0.01#0.03
+DoubleWellDelta=0.05
+Mobility=0.0005
+MobilityPressureSaturation=0.5
+Sigma1=1
+Sigma2=20
+Sigma3=20
+Lambda=0
+
+[Time]
+dtMax = 0.002
+
+[Initial]
+shape = 1
+# 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/pn2ptest.cc b/src/pntest/pn2ptest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f37a01d0b5f9194488930ed6eb62852484d76444
--- /dev/null
+++ b/src/pntest/pn2ptest.cc
@@ -0,0 +1,94 @@
+#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>
+#endif
+// phasefield includes
+#include<dune/phasefield/porenetwork/pn_2pflow.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("pn2p.ini",ptree);
+ ptreeparser.readOptions(argc,argv,ptree);
+
+ // Define Template Options
+ typedef double RF;
+ typedef Params_pn_2pflow< RF,
+ TripleWell<RF,
+ DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >
+ >,
+ TripleWellHierarchicalV1<RF,
+ DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,
+ TPPentalty_6thOrder<RF>
+ >,
+ Mobility_pf2p<RF>,
+ Reaction_Const<RF>,
+ VelDissipation_Quadratic<RF>
+ //VelDissipation_Quadratic_Shifted<RF>
+ > Parameters;
+ Parameters param(ptree);
+
+
+ //Set input/output filenamens
+ std::string GridFilename = ptree.get("domain.filename","square.msh");
+ std::string OutputFilename = ptree.get("output.filename","output");
+
+ //Set initial Shape
+ param.initial.shape = INITIAL_SHAPE_SQUARE;
+ //param.initial.shape = INITIAL_SHAPE_CIRCLE;
+
+ //Create Simulation
+ Pn_2PFlow<Parameters> pn_2pFlow(GridFilename,OutputFilename,param, 10e-6); //Last Parameter is initial pore Radius in m
+
+/*
+ //Example of usecase for Pn_1PFlow
+ std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
+ std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
+ std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
+
+ pn_1pFlow.writeVTK(0.0);
+ //pn_1pFlow.grow(2e-6); //Growth is measured in m
+ pn_1pFlow.growArea((10e-6)*(10e-6)*0.1); // Growth of Area in m^2
+ pn_1pFlow.writeVTK(1.0);
+ std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
+ std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl;
+ std::cout << "Kf has unit m^4" << std::endl;
+
+ std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
+ std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
+*/
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 1;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ return 1;
+ }
+}