diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index ed136822072a542abe33e8043b1ef2f1c5b0b316..47e93008ed5af59af34f174cad9380c8d72b2ad0 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -438,6 +438,8 @@ class Params_pn_2pflow : public Params_base<Number, Mobility>
public:
Number delta;
+ Number mu_ratio;
+
TW tw;
HTW htw;
Reaction reaction;
@@ -445,15 +447,20 @@ public:
Number SurfaceTensionScaling;
+ Number pf1Forcing;
+ Number pf2Forcing;
Params_pn_2pflow(const Dune::ParameterTree& ptree) :
Params_base<Number, Mobility>(ptree),
delta(ptree.get("Phasefield.delta",this->eps)),
+ mu_ratio(ptree.get("Pnysics.mu_ratio",1)),
tw(ptree.sub("Phasefield")),
htw(ptree.sub("Phasefield")),
reaction(ptree.sub("Physics")),
vDis(ptree.sub("Physics")),
- SurfaceTensionScaling(ptree.get("Physics.SurfaceTensionScaling",(Number)1))
+ SurfaceTensionScaling(ptree.get("Physics.SurfaceTensionScaling",(Number)1)),
+ pf1Forcing(0),
+ pf2Forcing(0)
{
}
};
diff --git a/dune/phasefield/chns_base.hh b/dune/phasefield/chns_base.hh
index fbd474aaf9841439791be727412848e90d7195a4..817ad4a122671988d816810b69453103fdfa4e15 100644
--- a/dune/phasefield/chns_base.hh
+++ b/dune/phasefield/chns_base.hh
@@ -242,6 +242,23 @@ public:
//_________________________________________________________________________________________________________________________
+template<typename GFS, typename NS_BDRY_FN, typename CH_BDRY_FN, typename NSX_BDRY_FN,
+ typename NS_V, typename CH_V, typename NSX_V>
+void interpolateToBdry(const GFS& gfs,
+ const NS_BDRY_FN& nsBdryFn, const CH_BDRY_FN& chBdryFn, const NSX_BDRY_FN& nsxBdryFn,
+ NS_V& nsV, CH_V& chV, NSX_V& nsxV)
+{
+ NS_V nsInterp(gfs.ns);
+ CH_V chInterp(gfs.ch);
+ NSX_V nsxInterp(gfs.nsx);
+ Dune::PDELab::interpolate(nsBdryFn, gfs.ns, nsInterp);
+ Dune::PDELab::interpolate(chBdryFn, gfs.ch, chInterp);
+ Dune::PDELab::interpolate(nsxBdryFn, gfs.nsx, nsxInterp);
+ Dune::PDELab::copy_constrained_dofs(gfs.nsCC, nsInterp, nsV);
+ Dune::PDELab::copy_constrained_dofs(gfs.chCC, chInterp, chV);
+ Dune::PDELab::copy_constrained_dofs(gfs.nsxCC, nsxInterp, nsxV);
+}
+
template<typename GFS, typename NS_BDRY_FN, typename CH_BDRY_FN, typename NS_V, typename CH_V>
void interpolateToBdry(const GFS& gfs, const NS_BDRY_FN& nsBdryFn, const CH_BDRY_FN& chBdryFn, NS_V& nsV, CH_V& chV)
{
@@ -543,6 +560,34 @@ public:
}
}
+ template<typename Grid, typename CH_V, typename NS_V, typename NSX_V, typename CH_GFS, typename NS_GFS, typename NSX_GFS>
+ void refineGrid_NSX(Grid &grid, CH_V &chV, NS_V &nsV, NSX_V &nsxV, CH_GFS chgfs, NS_GFS nsgfs, NSX_GFS nsxgfs, int verbosity = 0)
+ {
+ if(verbosity > 1)
+ {
+ printVector(chV,verbosity-2,"chV");
+ printVector(nsV,verbosity-2,"nsV");
+ printVector(nsxV,verbosity-2,"nsxV");
+ }
+
+ // do refinement
+ if(verbosity > 0) std::cout << "adapt grid and solution... ";
+ auto chTransfer = transferSolutions(chgfs,4,chV);
+ auto nsTransfer = transferSolutions(nsgfs,4,nsV);
+ auto nsxTransfer = transferSolutions(nsxgfs,4,nsxV);
+ Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer,nsxTransfer);
+
+ if(verbosity > 0) std::cout << " ...adapt ended" << std::endl;
+
+ if(verbosity > 1)
+ {
+
+ printVector(chV,verbosity-2,"chV");
+ printVector(nsV,verbosity-2,"nsV");
+ printVector(nsxV,verbosity-2,"nsxV");
+ }
+ }
+
template<typename Grid, typename CH_V, typename NS_V, typename CH_GFS, typename NS_GFS>
void myRefineGrid(Grid &grid, GV& gv, CH_V &chV, NS_V &nsV, CH_GFS chgfs, NS_GFS nsgfs, int verbosity = 0)
{
diff --git a/dune/phasefield/porenetwork/2p_chns.hh b/dune/phasefield/porenetwork/2p_chns.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ab593161cf67c1af0c1bcac02a9e301a8d35afb6
--- /dev/null
+++ b/dune/phasefield/porenetwork/2p_chns.hh
@@ -0,0 +1,293 @@
+#pragma once
+
+#include<dune/pdelab/constraints/common/constraints.hh>
+#include<dune/pdelab/constraints/conforming.hh>
+
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include<dune/pdelab/gridfunctionspace/subspace.hh>
+#include<dune/pdelab/gridfunctionspace/vectorgridfunctionspace.hh>
+#include<dune/phasefield/utility/pf_3p_color_dgf.hh>
+
+
+//_________________________________________________________________________________________________________________________
+
+
+template<typename RF, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM, typename VX_FEM>
+class GFS_2p_chns
+{
+public:
+ typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
+ typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
+ typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> VB_Block;
+ typedef Dune::PDELab::LexicographicOrderingTag OrderingTagL;
+ typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTagB;
+
+ ConstraintsAssembler conp, conphi, conmu;
+
+ typedef Dune::PDELab::VectorGridFunctionSpace<ES,V_FEM,ES::GridView::dimension,
+ VectorBackend,VectorBackend,ConstraintsAssembler> V_GFS;
+ V_GFS vGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ P_GFS pGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> PHI_GFS;
+ PHI_GFS phiGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, PHI_FEM, ConstraintsAssembler, VectorBackend> MU_GFS;
+ MU_GFS muGfs;
+
+ typedef Dune::PDELab::GridFunctionSpace<ES, VX_FEM, ConstraintsAssembler, VectorBackend> VX_GFS;
+ VX_GFS vxGfs;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, VX_GFS> NSX;
+ NSX nsx;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VectorBackend, OrderingTagL, V_GFS, P_GFS> NS;
+ NS ns;
+
+ typedef Dune::PDELab::CompositeGridFunctionSpace<VB_Block, OrderingTagB, PHI_GFS, MU_GFS, PHI_GFS, MU_GFS> CH;
+ CH ch;
+
+ typedef typename NS::template ConstraintsContainer<RF>::Type NS_CC;
+ NS_CC nsCC;
+
+ typedef typename CH::template ConstraintsContainer<RF>::Type CH_CC;
+ CH_CC chCC;
+
+ typedef typename NSX::template ConstraintsContainer<RF>::Type NSX_CC;
+ NSX_CC nsxCC;
+
+ GFS_2p_chns(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem, const VX_FEM& vxFem) :
+ conp(),conphi(),conmu(),
+ vGfs(es,vFem),
+ pGfs(es,pFem, conp),
+ phiGfs(es,phiFem, conphi),
+ muGfs(es,phiFem, conmu),
+ vxGfs(es,vxFem),
+ nsx(vxGfs),
+ ns(vGfs, pGfs),
+ ch(phiGfs, muGfs, phiGfs, muGfs),
+ nsCC(),
+ chCC(),
+ nsxCC()
+ {
+ vxGfs.name("Xvelocity");
+ vGfs.name("velocity");
+ pGfs.name("pressure");
+ phiGfs.name("phasefield");
+ muGfs.name("potential");
+ }
+
+ template<typename Bdry>
+ void updateConstraintsContainer(const Bdry& bdry)
+ {
+ Dune::PDELab::constraints(bdry.bdryNS,ns,nsCC, 0);
+ Dune::PDELab::constraints(bdry.bdryCH,ch,chCC, 0);
+ Dune::PDELab::constraints(bdry.bdryNSX,nsx,nsxCC, 0);
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+
+template<typename GFS, typename NS_V, typename CH_V, typename NSX_V>
+class DGF_2p_chns
+{
+public:
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<0> > VSubGFS;
+ VSubGFS vSubGfs;
+ typedef Dune::PDELab::VectorDiscreteGridFunction<VSubGFS,NS_V> VDGF;
+ VDGF vDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NS,Dune::TypeTree::TreePath<1> > PSubGFS;
+ PSubGFS pSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PSubGFS,NS_V> PDGF;
+ PDGF pDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<0> > PhiSubGFS;
+ PhiSubGFS phiSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<PhiSubGFS,CH_V> PhiDGF;
+ PhiDGF phiDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<1> > MuSubGFS;
+ MuSubGFS muSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<MuSubGFS,CH_V> MuDGF;
+ MuDGF muDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<2> > Phi2SubGFS;
+ Phi2SubGFS phi2SubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<Phi2SubGFS,CH_V> Phi2DGF;
+ Phi2DGF phi2Dgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::CH,Dune::TypeTree::TreePath<3> > Mu2SubGFS;
+ Mu2SubGFS mu2SubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<Mu2SubGFS,CH_V> Mu2DGF;
+ Mu2DGF mu2Dgf;
+
+ typedef PF_3P_Color_DGF<typename GFS::CH,CH_V,3> ColorDGF;
+ ColorDGF colorDgf;
+
+ typedef Dune::PDELab::GridFunctionSubSpace<typename GFS::NSX,Dune::TypeTree::TreePath<0> > VXSubGFS;
+ VXSubGFS vxSubGfs;
+ typedef Dune::PDELab::DiscreteGridFunction<VXSubGFS,NSX_V> VXDGF;
+ VXDGF vxDgf;
+
+
+ DGF_2p_chns(const GFS& gfs, const NS_V& nsV, const CH_V& chV, const NSX_V& nsxV) :
+ vSubGfs(gfs.ns),
+ vDgf(vSubGfs,nsV),
+ pSubGfs(gfs.ns),
+ pDgf(pSubGfs,nsV),
+
+ phiSubGfs(gfs.ch),
+ phiDgf(phiSubGfs,chV),
+ muSubGfs(gfs.ch),
+ muDgf(muSubGfs,chV),
+ phi2SubGfs(gfs.ch),
+ phi2Dgf(phi2SubGfs,chV),
+ mu2SubGfs(gfs.ch),
+ mu2Dgf(mu2SubGfs,chV),
+
+ colorDgf(gfs.ch,chV),
+
+ vxSubGfs(gfs.nsx),
+ vxDgf(vxSubGfs,nsxV)
+ {
+ }
+
+ template<typename VTKWRITER>
+ void addToVTKWriter(VTKWRITER& vtkwriter)
+ {
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vDgf,"v"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PDGF> >(pDgf,"p"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<Phi2DGF> >(phi2Dgf,"phi2"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<Mu2DGF> >(mu2Dgf,"mu2"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ColorDGF> >(colorDgf,"_coloring_"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VXDGF> >(vxDgf,"vx"));
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename IniV, typename IniP, typename IniPhi, typename IniMu, typename IniPhi2, typename IniMu2, typename IniVX>
+class Ini_2p_chns
+{
+public:
+ IniV iniV;
+ IniP iniP;
+ IniPhi iniPhi;
+ IniMu iniMu;
+ IniPhi2 iniPhi2;
+ IniMu2 iniMu2;
+ typedef Dune::PDELab::CompositeGridFunction<IniV,IniP> IniNS;
+ typedef Dune::PDELab::CompositeGridFunction<IniPhi,IniMu,IniPhi2,IniMu2> IniCH;
+ IniNS iniNS;
+ IniCH iniCH;
+
+ IniVX iniVx;
+ typedef Dune::PDELab::CompositeGridFunction<IniVX> IniNSX;
+ IniNSX iniNSX;
+
+ template<typename GV, typename Params>
+ Ini_2p_chns(const GV& gv, const Params& params) :
+ iniV(gv, params),
+ iniP(gv, params),
+ iniPhi(gv, params),
+ iniMu(gv, params),
+ iniPhi2(gv, params),
+ iniMu2(gv, params),
+ iniNS(iniV,iniP),
+ iniCH(iniPhi,iniMu,iniPhi2,iniMu2),
+ iniVx(gv,params),
+ iniNSX(iniVx)
+ {
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template<typename BdryV1, typename BdryV2, typename BdryP,
+ typename BdryPhi, typename BdryMu, typename BdryPhi2, typename BdryMu2, typename BdryVX>
+class Bdry_2p_chns_compositeV_2d
+{
+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;
+
+ BdryVX bdryVx;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryVX> BdryNSX;
+ BdryNS bdryNSX;
+
+ Bdry_2p_chns_compositeV_2d() :
+ bdryV1(),
+ bdryV2(),
+ bdryP(),
+ bdryPhi(),
+ bdryMu(),
+ bdryPhi2(),
+ bdryMu2(),
+ bdryV(bdryV1,bdryV2),
+ bdryNS(bdryV,bdryP),
+ bdryCH(bdryPhi,bdryMu,bdryPhi2,bdryMu2),
+ bdryVx(),
+ bdryNSX(bdryVx)
+ {
+ }
+};
+
+template<typename BdryV1, typename BdryV2, typename BdryP,
+ typename BdryPhi, typename BdryMu, typename BdryPhi2, typename BdryMu2, typename BdryVX>
+class Bdry_2p_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;
+
+ BdryVX bdryVx;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryVX> BdryNSX;
+ BdryNSX bdryNSX;
+
+
+ template<typename DomainSize>
+ Bdry_2p_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),
+ bdryVx(ds),
+ bdryNSX(bdryVx)
+ {
+ }
+};
diff --git a/dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh b/dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh
index d36cd10feb59412775ba3c174097131067d106f9..c9d4759b26f10b44ff3a5720d8822375c4fe442d 100644
--- a/dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh
+++ b/dune/phasefield/porenetwork/fem_2p_cahnhilliard_pressuresaturation.hh
@@ -85,7 +85,7 @@ public:
(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
+ -(ch.pf+pf3*param.delta)*(ns.v*ch.basis.gradphi[i][0]) //weak transport, modified by \delta \phi_3 to ensure d_t \phi_3 = 0
+R1
));
@@ -93,7 +93,7 @@ public:
(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])
+ -(ch.pf2+pf3*param.delta)*(ns.v*ch.basis.gradphi[i][0])
-R1
));
diff --git a/dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh b/dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh
index 7dda41197e79436c1a7b596fa0b05d10ade646fa..bd177638398cd9bdd4b1788306575b33cbd8fed0 100644
--- a/dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh
+++ b/dune/phasefield/porenetwork/fem_2p_navierstokes_pressuresaturation.hh
@@ -90,7 +90,7 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
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=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++)
diff --git a/dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh b/dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh
index eaa3364d03d74a11cf1b5f9c6a765464f5d6150f..0177e54693664b77b5254dfa8efca0e23f54b209 100644
--- a/dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh
+++ b/dune/phasefield/porenetwork/fem_2p_navierstokes_transmissibility.hh
@@ -9,7 +9,7 @@
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/phasefield/localoperator/fem_base_chns.hh>
-
+#include<dune/phasefield/porenetwork/fem_base_porenetwork.hh>
@@ -25,15 +25,10 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
CHLocalBasisCache<CHGFS> chCache;
- VLocalBasisCache<NSGFS> vCache;
- PLocalBasisCache<NSGFS> pCache;
+ Micro_VXLocalBasisCache<NSGFS> vCache;
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
@@ -53,11 +48,8 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
- FEM_2P_NavierStokes_Transmissibility (Param& param_, CHGFS& chGfs_, CHContainer& chV_, CHContainer& chVOld_,
- NSGFS& nsGfs_, NSContainer& nsVOld_) :
+ FEM_2P_NavierStokes_Transmissibility (Param& param_, CHGFS& chGfs_, CHContainer& chV_) :
chLocal(chGfs_,chV_),
- //chOldLocal(chGfs_,chVOld_),
- nsOldLocal(nsGfs_,nsVOld_),
param(param_)
{
}
@@ -73,8 +65,6 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
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)
{
@@ -82,65 +72,28 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
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);
+ NSPorenetworkVariables<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;
- 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;
+ RF mu_min = std::max(1.,param.mu_ratio);
+ RF mu = std::min(mu_min, (ch.pf + pf3)*1 + ch.pf2*param.mu_ratio );
+ //Dune::FieldVector<RF,dim> gradpf_f(0.0);
+ //gradpf_f.axpy((1-2*param.delta),ch.gradpf+ch.gradpf2);
- 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++)
+ //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.vxspace.size(); i++) //TODO Change for 3p
{
- 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
+ r.accumulate(ns.vxspace, i, factor*(
+ param.vDis.eval(pf_f)*ns.vx*ns.vbasis.phi[i]
+ + mu*(ns.gradvx * ns.vbasis.gradphi[i][0])
+ - param.pf1Forcing * ch.pf * ns.vbasis.phi[i]
+ - param.pf2Forcing * ch.pf2* ns.vbasis.phi[i]
+ ));
} // for i
} // for ip
} //alpha volume
diff --git a/dune/phasefield/porenetwork/pn_2pflow.hh b/dune/phasefield/porenetwork/pn_2pflow.hh
index 02fc246a89ff3b6e3765a994fea8f22d00634680..a2bcca9a144132849ac9c5dc611c53707e19b8be 100644
--- a/dune/phasefield/porenetwork/pn_2pflow.hh
+++ b/dune/phasefield/porenetwork/pn_2pflow.hh
@@ -23,7 +23,7 @@
#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/porenetwork/2p_chns.hh>
#include <dune/phasefield/initial_fn.hh>
#include <dune/phasefield/debug/vector_debug.hh>
@@ -50,52 +50,58 @@ private:
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 Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,2> VX_FEM;
- typedef GFS_fff_chns<RF, ES, V_FEM, P_FEM, PHI_FEM> GFS;
+ typedef GFS_2p_chns<RF, ES, V_FEM, P_FEM, PHI_FEM, VX_FEM> GFS;
typedef typename GFS::NS NS_GFS;
typedef typename GFS::CH CH_GFS;
+ typedef typename GFS::NSX NSX_GFS;
typedef Dune::PDELab::Backend::Vector<CH_GFS, RF> CH_V;
typedef Dune::PDELab::Backend::Vector<NS_GFS, RF> NS_V;
+ typedef Dune::PDELab::Backend::Vector<NSX_GFS, RF> NSX_V;
- typedef Ini_fff_chns< ZeroInitialVector<GV,RF, dim, Parameters>, //V
+ typedef Ini_2p_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>,
ZeroInitial<GV,RF,Parameters>> Initial;
typedef DomainSize<RF> DS;
- typedef Bdry_fff_chns_compositeV_2d_Explicit_Domain<AllBoundary<DS>, //V
+ typedef Bdry_2p_chns_compositeV_2d_Explicit_Domain<AllBoundary<DS>, //V
AllBoundary<DS>, //V
TopRightCorner<DS>, //P
NoBoundary<DS>,
NoBoundary<DS>,
NoBoundary<DS>,
- NoBoundary<DS>> Bdry;
+ NoBoundary<DS>,
+ LeftBottomBoundary<DS> //VX
+ > 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 FEM_2P_NavierStokes_Transmissibility<Parameters, CH_GFS, CH_V, NSX_GFS, NSX_V> NSX_T_LOP;
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
- typedef DGF_fff_chns<GFS, NS_V, CH_V> DGF;
+ typedef DGF_2p_chns<GFS, NS_V, CH_V, NSX_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 Dune::PDELab::GridOperator<NSX_GFS,NSX_GFS,NSX_T_LOP, MBE,RF,RF,RF,typename GFS::NSX_CC,typename GFS::NSX_CC> NSX_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 ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<NSX_T_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NSX_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;
+ typedef Dune::PDELab::Newton<NSX_T_GO,NSX_T_LS,NSX_V> NSX_T_Newton;
Parameters param;
std::shared_ptr<Grid> gridp;
@@ -104,42 +110,55 @@ private:
V_FEM vFem;
P_FEM pFem;
PHI_FEM phiFem;
+ VX_FEM vxFem;
GFS gfs;
CH_V chV, chVOld;
NS_V nsV, nsVOld;
+ NSX_V nsxV;
Initial initial;
DS ds;
Bdry bdry;
MBE mbe;
- RF kf = 0;
- RF phiSolid = 0;
+ RF kf_11 = 0;
+ RF kf_12 = 0;
+ RF kf_21 = 0;
+ RF kf_22 = 0;
+
+ RF phiFluid1 = 0; // = Integral(phi1)
+ RF phiFluid2 = 0; // = Integral(phi2)
+ RF phiSolid = 0; // = Integral(phi3)
+ RF saturation = 0; // = phiFluid1/(phiFluid1+phiFluid2)
+
+ RF pressureFluid1 = 0;
+ RF pressureFluid2 = 0;
+ RF capillaryPressure = 0;
- RF PoreReferenceLength;
+ 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<NSX_T_LOP> nsxtLop;
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<NSX_T_GO> nsxtGo;
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<NSX_T_LS> nsxtLs;
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<NSX_T_Newton> nsxtNewton;
std::unique_ptr<Dune::VTKSequenceWriter<GV>> vtkwriter;
//------------------------------------------------------------------------------------
@@ -151,12 +170,12 @@ private:
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);
+ nsxtGo = std::make_unique<NSX_T_GO> (gfs.nsx,gfs.nsxCC,gfs.nsx,gfs.nsxCC,*nsxtLop,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);
+ nsxtLs = std::make_unique<NSX_T_LS>(*nsxtGo,1000,200,3,0);
// Nonlinear solver
nspsNewton = std::make_unique<NS_PS_Newton>(*nspsGo,nsV,*nspsLs);
nspsNewton->setParameters(param.nsNewtonTree);
@@ -167,8 +186,8 @@ private:
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);
+ nsxtNewton = std::make_unique<NSX_T_Newton>(*nsxtGo,nsxV,*nsxtLs);
+ nsxtNewton->setParameters(param.nsNewtonTree);
//printVector(chV,10,"chV");
//printVector(nsV,10,"nsV");
@@ -181,10 +200,10 @@ private:
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);
+ rInd.refineGrid_NSX(*gridp,chV,nsV,nsxV,gfs.ch,gfs.ns,gfs.nsx,0);
// recompute constraints
gfs.updateConstraintsContainer(bdry);
- interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,initial.iniNSX,nsV,chV,nsxV);
initializeGridOperators();
}
@@ -206,8 +225,8 @@ private:
void applyNsTransmissibility(int verbosity = 0)
{
- if (verbosity > 0) std::cout << "= Begin newtonapply NS_T:" << std::endl;
- nstNewton->apply();
+ if (verbosity > 0) std::cout << "= Begin newtonapply NSX_T:" << std::endl;
+ nsxtNewton->apply();
}
void applyChReaction(int verbosity = 0)
@@ -219,41 +238,87 @@ private:
remesh();
}
-/* void updateKf()
+ void updateKf()
+ {
+ param.pf1Forcing = 1;
+ param.pf2Forcing = 0;
+ applyNsTransmissibility();
+ updateKf_integration(true);
+ param.pf1Forcing = 0;
+ param.pf2Forcing = 1;
+ applyNsTransmissibility();
+ updateKf_integration(false);
+ }
+
+ void updateKf_integration(bool phase1)
{
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> phi2;
+ dgf->phi2Dgf.evaluate(element, xlocal, phi2);
+ return Dune::FieldVector<double,2>{phi1*w,phi2*w};
});
- Dune::FieldVector<double,1> integral;
+ Dune::FieldVector<double,2> 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;
+ RF scaling_factor = 4*poreReferenceLength*poreReferenceLength*poreReferenceLength*poreReferenceLength;
+ if(phase1)
+ {
+ kf_11 = scaling_factor*integral[0];
+ kf_21 = scaling_factor*integral[1];
+ }
+ else
+ {
+ kf_12 = scaling_factor*integral[0];
+ kf_22 = scaling_factor*integral[1];
+ }
}
- void updatePhiSolid()
+
+ void updateSaturation()
+ {
+ 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,3>{phi1[0], phi2[0], 1-phi1[0]-phi2[0]};
+ });
+ Dune::FieldVector<double,3> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ phiFluid1 = integral[0];
+ phiFluid2 = integral[1];
+ phiSolid = integral[2];
+ saturation = phiFluid1/(phiFluid1+phiFluid2); //TODO Check div by 0
+ }
+
+ void updateCapillaryPressure()
{
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> phi2;
+ dgf->phi2Dgf.evaluate(element, xlocal, phi2);
+ Dune::FieldVector<double,1> p;
+ dgf->pDgf.evaluate(element, xlocal, p);
+ return Dune::FieldVector<double,4>{ (phi1[0]>0.75)*p[0], (phi1[0]>0.75), (phi2[0]>0.75)*p[0], (phi2[0]>0.75)};
});
- Dune::FieldVector<double,1> integral;
+ Dune::FieldVector<double,4> integral;
Dune::PDELab::integrateGridFunction(productFunction,integral,2);
- phiSolid = integral[0];
+ //phiSolid = integral[0];
+
+ pressureFluid1 = integral[0]/integral[1];
+ pressureFluid2 = integral[2]/integral[3];
+ capillaryPressure = pressureFluid2-pressureFluid1;
}
- RF calculateGrowthSpeed()
+/* RF calculateGrowthSpeed()
{
auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
Dune::FieldVector<double,1> phi1;
@@ -280,14 +345,15 @@ public:
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),
+ vFem(es), pFem(es), phiFem(es),vxFem(es),
+ gfs(es,vFem,pFem,phiFem,vxFem),
+ chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch), nsVOld(gfs.ns), nsxV(gfs.nsx),
initial(gv,param), ds(0,1,0,1,1e-6), bdry(ds), mbe(10),
- PoreReferenceLength(PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) ) )
+ 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);
+ Dune::PDELab::interpolate(initial.iniNSX , gfs.nsx, nsxV);
// Refine Intial Values
for(int i=0; i<param.pTree.get("Refinement.maxRefineLevel",2); i++)
@@ -297,14 +363,15 @@ public:
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.refineGrid_NSX(*gridp,chV,nsV,nsxV,gfs.ch,gfs.ns,gfs.nsx,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);
+ Dune::PDELab::interpolate(initial.iniNSX , gfs.nsx, nsxV);
}
gfs.updateConstraintsContainer(bdry);
- interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,initial.iniNSX, nsV,chV,nsxV);
//chVOld = chV;
@@ -312,9 +379,9 @@ public:
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);
+ nsxtLop = std::make_unique<NSX_T_LOP> (param, gfs.ch, chV);
- dgf = std::make_unique<DGF>(gfs,nsV,chV);
+ dgf = std::make_unique<DGF>(gfs,nsV,chV,nsxV);
struct stat st;
if( stat( OutputFilename.c_str(), &st ) != 0 )
@@ -339,10 +406,68 @@ public:
RF dtmax = param.pTree.get("Time.dtMax",0.01);
param.time.dt = dtmax;
- for(int i=0; i<1000; i++)
+
+ int InitialGrowthTime = param.pTree.get("test.InitialGrowthTime",0);
+ int FillTime = param.pTree.get("test.FillTime",100);
+ int EmptyTime = param.pTree.get("test.EmptyTime",200);
+
+ param.reaction.reactionRate = 0;
+ for(int j=0; j<100; j++)
{
- writeVTK(i);
- if(i < 10)
+ applyChPressureSaturation(1);
+ applyNsPressureSaturation(1);
+ }
+ for(int j=0; j<InitialGrowthTime; j++)
+ {
+ param.reaction.reactionRate = 1;
+ applyChReaction(1);
+ }
+ param.reaction.reactionRate = 0;
+ for(int j=0; j<100; j++)
+ {
+ applyChPressureSaturation(1);
+ applyNsPressureSaturation(1);
+ }
+
+ for(int i=0; i<FillTime+EmptyTime; i++)
+ {
+ if(i<FillTime)
+ {
+ param.reaction.reactionRate = 1;
+ }
+ else
+ {
+ param.reaction.reactionRate = -1;
+ }
+ applyChPressureSaturation(1);
+ param.reaction.reactionRate = 0;
+ for(int j=0; j<20; j++)
+ {
+ applyChPressureSaturation(1);
+ applyNsPressureSaturation(1);
+ //writeVTK(i+((RF)j+1)/21);
+ }
+ writeVTK(i+1);
+ updateSaturation();
+ updateCapillaryPressure();
+ updateKf();
+
+ try
+ {
+ std::ofstream myfile;
+ myfile.open("pressuresaturation.csv",std::ofstream::app);
+ myfile << i << ", " << phiFluid1 << ", " << phiFluid2 << ", " << phiSolid << ", " << saturation << ", "
+ << pressureFluid1 << ", " << pressureFluid2 << ", " << capillaryPressure << ", "
+ << kf_11 << ", " << kf_12 << ", " << kf_21 << ", " << kf_22 << ", " << std::endl;
+ myfile.close();
+ }
+ catch(...)
+ {
+ std::cout << "Error in integration or CSV Files" << std::endl;
+ }
+
+
+ /*if(i < 10)
{
applyChPressureSaturation(1);
}
@@ -356,7 +481,7 @@ public:
{
param.reaction.reactionRate = 1;
applyChReaction(1);
- }
+ }*/
}
}
@@ -369,7 +494,7 @@ public:
void grow(RF length) // How much the the solid-fluid interface should grow in normal direction (in m)
{
- length = length/PoreReferenceLength;
+ length = length/poreReferenceLength;
if(length <= 0)
{
std::cout << "Need positive length (dissolution not yet implemented)" << std::endl;
@@ -393,7 +518,7 @@ public:
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
+ area = area/(poreReferenceLength*poreReferenceLength*4); //4 because of Symmetry
if(area <= 0)
{
std::cout << "Need positive area (dissolution not yet implemented)" << std::endl;
@@ -427,7 +552,7 @@ public:
RF getPoreArea()
{
- return 4*PoreReferenceLength*PoreReferenceLength*(1-phiSolid); //4 because of Symmetry
+ return 4*poreReferenceLength*poreReferenceLength*(1-phiSolid); //4 because of Symmetry
}
RF getPoreAreaFraction()
@@ -437,12 +562,12 @@ public:
void setPoreBaseRadius(RF PoreBaseRadius_)
{
- PoreReferenceLength = PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) );
+ poreReferenceLength = PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) );
}
RF getPoreReferenceLength()
{
- return PoreReferenceLength;
+ return poreReferenceLength;
}
*/
diff --git a/src/pntest/pn2p.ini b/src/pntest/pn2p.ini
index 4a774f424c8c7972d56f99ff4346dd6a7200b230..b955c96f961bb4b367f99d61e7e29b231489afac 100644
--- a/src/pntest/pn2p.ini
+++ b/src/pntest/pn2p.ini
@@ -3,13 +3,20 @@ uAst = 1
D = 0.02
rho1 = 1.0
rho2 = 1.0
-mu=0.05
+mu=0.05 # For the Stokes to Equilibrium, tuned to the timestep dtMax
+mu_ratio=0.018 # 0.00018 / 0.01 (Viscosity Air / Viscosity Water)
ReactionRate=1#1.5
#ReactionLimiter = 1000
SurfaceTensionScaling=1
VelDissipationRate = 4000
CurvatureAlpha=0#0.001
+[test]
+InitialGrowthTime = 10
+FillTime = 150
+EmptyTime = 400
+
+
[domain]
filename = grids/square_periodic.msh
@@ -18,8 +25,9 @@ filename = outp3
[Phasefield]
eps=0.02
-delta=0.01#0.03
-DoubleWellDelta=0.05
+delta=0.05
+DoubleWellDelta=0.025
+DoubleWellGamma=0.02
Mobility=0.0005
MobilityPressureSaturation=0.5
Sigma1=1