Made Input/Output dimensional (Added new Parameter for constructor to define...

Made Input/Output dimensional (Added new Parameter for constructor to define initial pore radius, Kf has now dimension m^4). Paramerter classes now used by value instead of by reference

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)),

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

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;
+  }
 };

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

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;