diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 86f1a1d6cc0dbe3b3c9a677e8ae05062c3d8b0cf..950a9dae589269e9c79f1de18e12cb7dcef0d6d3 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -255,7 +255,7 @@ public:
class Params_Initial
{
public:
- const Dune::ParameterTree& ptree;
+ const Dune::ParameterTree ptree;
Params_Initial(const Dune::ParameterTree& tree) :
ptree(tree)
{
@@ -276,9 +276,9 @@ public:
Params_Initial initial;
- const Dune::ParameterTree& pTree;
- const Dune::ParameterTree& nsNewtonTree;
- const Dune::ParameterTree& chNewtonTree;
+ const Dune::ParameterTree pTree;
+ const Dune::ParameterTree nsNewtonTree;
+ const Dune::ParameterTree chNewtonTree;
Params_base(const Dune::ParameterTree& ptree) :
eps(ptree.get("Phasefield.eps",(Number)0.1)),
diff --git a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
index 1f6c2c763c6ba67e1198ff025befea1f6f7cb29c..26cec5c6a6891e3a1af6de41b0c6442ec65619c0 100644
--- a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
+++ b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
@@ -80,7 +80,7 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
{
r.accumulate(ns.vxspace, i, factor*(
param.vDis.eval(pf_f)*ns.vx*ns.vbasis.phi[i]
- + param.phys.eval_mu(pf_f) * (ns.gradvx * ns.vbasis.gradphi[i][0])
+ + (ns.gradvx * ns.vbasis.gradphi[i][0])
- pf_f * ns.vbasis.phi[i]
));
} // for i
diff --git a/dune/phasefield/porenetwork/pn_1pflow.hh b/dune/phasefield/porenetwork/pn_1pflow.hh
index cac91eac0da8ea0af65d9a37788b4a48141403fa..cc60fce716c9776e5afe4bb4c4928f0f84a13859 100644
--- a/dune/phasefield/porenetwork/pn_1pflow.hh
+++ b/dune/phasefield/porenetwork/pn_1pflow.hh
@@ -78,7 +78,7 @@ private:
typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> NS_Newton;
typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> CH_Newton;
- Parameters& param;
+ Parameters param;
std::shared_ptr<Grid> gridp;
GV gv;
ES es;
@@ -98,6 +98,8 @@ private:
RF kf = 0;
RF phiSolid = 0;
+ RF PoreReferenceLength;
+
std::unique_ptr<NS_LOP> nsLop;
std::unique_ptr<CH_LOP> chLop;
std::unique_ptr<DGF> dgf;
@@ -118,8 +120,8 @@ private:
nsGo = std::make_unique<NS_GO>(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,*nsLop,mbe);
chGo = std::make_unique<CH_GO>(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,*chLop,mbe);
// Linear solver
- chLs = std::make_unique<CH_LS>(*chGo,500,100,3,1);
- nsLs = std::make_unique<NS_LS>(*nsGo,1000,200,3,1);
+ chLs = std::make_unique<CH_LS>(*chGo,500,100,3,0);
+ nsLs = std::make_unique<NS_LS>(*nsGo,1000,200,3,0);
// Nonlinear solver
nsNewton = std::make_unique<NS_Newton>(*nsGo,nsV,*nsLs);
nsNewton->setParameters(param.nsNewtonTree);
@@ -174,7 +176,9 @@ private:
});
Dune::FieldVector<double,1> integral;
Dune::PDELab::integrateGridFunction(productFunction,integral,2);
- kf = integral[0];
+ // Calculation of kf: 4* because of Quarter-Circle,
+ // Scaling with x_ref^4
+ kf = 4 * integral[0] * PoreReferenceLength*PoreReferenceLength*PoreReferenceLength*PoreReferenceLength;
}
void updatePhiSolid()
@@ -197,7 +201,7 @@ private:
public:
//Constructor
- Pn_1PFlow(std::string GridFilename, std::string OutputFilename, Parameters& param_) :
+ Pn_1PFlow(std::string GridFilename, std::string OutputFilename, Parameters param_, RF PoreBaseRadius_ = 10e-6) :
param(param_),
gridp(initPnGrid<Grid>(GridFilename)),
gv(gridp->leafGridView()),
@@ -205,7 +209,8 @@ public:
vFem(es), phiFem(es),
gfs(es,vFem,phiFem),
chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch),
- initial(gv,param), ds(0,1,0,1,1e-6), bdry(ds),mbe(10)
+ 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);
@@ -261,6 +266,7 @@ public:
void grow(RF length)
{
+ length = length/PoreReferenceLength;
if(length <= 0)
{
std::cout << "Need positive length" << std::endl;
@@ -269,7 +275,7 @@ public:
RF dtmax = param.pTree.get("Time.dtMax",0.01);
while(length > dtmax)
{
- writeVTK(1-length);
+ //writeVTK(1-length);
length -= dtmax;
param.time.dt = dtmax;
applyCh(0);
@@ -290,4 +296,14 @@ public:
{
return phiSolid;
}
+
+ void setPoreBaseRadius(RF PoreBaseRadius_)
+ {
+ PoreReferenceLength = PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) );
+ }
+
+ RF getPoreReferenceLength()
+ {
+ return PoreReferenceLength;
+ }
};
diff --git a/src/pntest/pn1p.ini b/src/pntest/pn1p.ini
index 4a4a76dbba81e6a517bfc7acb82db4b45e962a95..e6954242fccec148a58c598e32ffde9ce9f9ef45 100644
--- a/src/pntest/pn1p.ini
+++ b/src/pntest/pn1p.ini
@@ -3,7 +3,7 @@ uAst = 1
D = 0.02
rho1 = 1.0
rho2 = 1.0
-mu=0.005
+#mu=0.005
ReactionRate=1#1.5
#ReactionLimiter = 1000
SurfaceTensionScaling=0.1
diff --git a/src/pntest/pntest.cc b/src/pntest/pntest.cc
index d1a8d99647101b497d5d8a9e934e836fd92f017b..1d66b2e542a9eae1ad047b9e63aea27184d65fe4 100644
--- a/src/pntest/pntest.cc
+++ b/src/pntest/pntest.cc
@@ -53,16 +53,17 @@ int main(int argc, char** argv)
std::string OutputFilename = ptree.get("output.filename","output");
//Create Simulation
- Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param);
+ Pn_1PFlow<Parameters> pn_1pFlow(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;
pn_1pFlow.writeVTK(0.0);
- pn_1pFlow.grow(0.2);
+ pn_1pFlow.grow(2e-6); //Growth is measured in m
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;
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;