Fixed Kc, Kf, phiC in thinstrip_mixed, added ion source to thinstrip_full

dune.module

@@ -7,4 +7,4 @@ Module: dune-phasefield
 Version: 0.1
 Maintainer: lars.von-wolff@mathematik.uni-stuttgart.de
 #depending on #dune-xt-grid
-Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-typetree dune-istl dune-functions dune-alugrid dune-pdelab dune-xt
+Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-typetree dune-istl dune-functions dune-alugrid dune-pdelab dune-xt-grid dune-xt-common

src/paper_thinstrip_full/thinstrip_parameters.hh

@@ -11,12 +11,14 @@ public:
   const Number YScaling;
   const Number len;
   Number gradP;
+  const Number ionSource;
 
   Params_thinstrip(const Dune::ParameterTree& ptree) :
     Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>(ptree),
     YScaling(ptree.get("domain.YScaling",1.0)),
     len(ptree.get("domain.len",1.0)),
-    gradP(ptree.get("Initial.pConst",1.0))
+    gradP(ptree.get("Initial.pConst",1.0)),
+    ionSource(ptree.get("Initial.ionSource",1.0))
   {
     ModifyParameters();
     std::cout << "YScaling " << YScaling << " Eps " << this->eps << std::endl;

src/paper_thinstrip_mixed/chns.ini

@@ -20,7 +20,7 @@ YResolution = 200
 XResolution = 400
 
 [output]
-filename = comparison_hyperbolic
+filename = comparison_hyperbolic_bugfix
 
 [Phasefield]
 eps=0.03
@@ -45,6 +45,7 @@ saveintervall = 0.005
 pConst=0.05
 uConst=0.5
 uInflow=0.3
+ionSource=0.0
 
 [Refinement]
 etaRefine  = 0.2

src/paper_thinstrip_mixed/chns_reactive.ini

+[Physics]
+uAst = 1
+D = 0.02
+rho1 = 1.0
+rho2 = 1.0
+mu=0.005
+ReactionRate=2#0.5
+#ReactionLimiter = 1000
+SurfaceTensionScaling=0.1
+VelDissipationRate = 1000
+#VelDissipationMax = 0.9
+CurvatureAlpha=0.001 #0
+
+[domain]
+level = 0
+filename = grids/rectangle_periodic_coarse.msh
+#YScalingFactor = 0.5
+YScaling = 1
+YResolution = 100#200
+XResolution = 100#400
+
+[output]
+filename = comparison_diffusion
+
+[Phasefield]
+eps=0.03
+delta= 0.03 #0.00001
+DoubleWellDelta=0.05
+Mobility=0.04 #1 #0.01
+Sigma1=1
+Sigma2=1
+Sigma3=1
+
+[Time]
+tMax = 4
+dt = 0.002#0.0001
+dtMax = 0.002#0.0001
+dtMin = 0.000015
+saveintervall = 0.005
+#dtIncreaseFactor=1.2
+#dtDecreaseFactor=0.5
+
+[Initial]
+# FlowVyy= 10# pConst/(Domain_width * mu)
+pConst=0.05
+uConst=0.5
+uInflow=0.3
+ionSource=1.0
+
+[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

src/paper_thinstrip_mixed/driver_mixed.hh

@@ -65,6 +65,7 @@ public:
   typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> NS_Newton;
   typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> CH_Newton;
 
+  Parameters& parameters;
   GV& gv;
   CH_V chV, chVOld, chVUpwind;
   NS_V nsV, nsVOld, nsVUpwind;
@@ -96,7 +97,7 @@ public:
   std::unique_ptr<CH_Newton>  chNewton;
 
   YSlice(GV& gv_, GFS& gfs, Parameters& param, RF xpos_, RF dx_) :
-  gv(gv_), chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch), nsVOld(gfs.ns), chVUpwind(gfs.ch), nsVUpwind(gfs.ns),
+  parameters(param), gv(gv_), chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch), nsVOld(gfs.ns), chVUpwind(gfs.ch), nsVUpwind(gfs.ns),
   initial(gv,param), xpos(xpos_), mbe(10), dx(dx_)
   {
     param.TempXpos = xpos;
@@ -169,7 +170,8 @@ public:
         dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
         Dune::FieldVector<double,1> phi2;
         dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
-        return Dune::FieldVector<double,1>((phi1+phi2)*w);
+        Dune::FieldVector<double,1> phi3 = 1-phi1-phi2;
+        return Dune::FieldVector<double,1>((phi1+phi2+2*parameters.delta*phi3)*w);
     });
     Dune::FieldVector<double,1> integral;
     Dune::PDELab::integrateGridFunction(productFunction,integral,2);
@@ -182,7 +184,7 @@ public:
     auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
         Dune::FieldVector<double,1> phi1;
         dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
-        return Dune::FieldVector<double,1>(phi1);
+        return Dune::FieldVector<double,1>(phi1+parameters.delta);
     });
     Dune::FieldVector<double,1> integral;
     Dune::PDELab::integrateGridFunction(productFunction,integral,2);
@@ -210,7 +212,7 @@ public:
         dgfMicro->vxDgf.evaluate(element, xlocal, w);
         Dune::FieldVector<double,1> phi1;
         dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
-        return Dune::FieldVector<double,1>(phi1*w);
+        return Dune::FieldVector<double,1>((phi1+parameters.delta)*w);
     });
     Dune::FieldVector<double,1> integral;
     Dune::PDELab::integrateGridFunction(productFunction,integral,2);

src/paper_thinstrip_mixed/tmixed_boundary.hh

@@ -175,7 +175,9 @@ public:
   {
     //bool inMiddle = (pref.TempXpos > 0.4) && (pref.TempXpos < 0.6);
     //RF mod = 0.1*std::sin(2*3.141592*pref.TempXpos);//inMiddle?0.3:0;
-    RF d_sine = 0.15*sin(pref.TempXpos*2*M_PI);
+
+    RF d_sine = /*0;*/0.15*sin(pref.TempXpos*2*M_PI);
+
     //std::cout << "TempXPos" << pref.TempXpos;
     RF rf = pref.tw.dw.eval_shape(1/eps*(x[0] - dSolid));
     RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[0] - (dSolid + d1 + d_sine)));

src/paper_thinstrip_mixed/tmixed_parameters.hh

@@ -9,11 +9,13 @@ class Params_tmixed : public Params_ffs_r<Number, TW, Mobility, Reaction, VelDis
 {
 public:
   const Number YScaling;
+  const Number ionSource;
   Number TempXpos = 0;
 
   Params_tmixed(const Dune::ParameterTree& ptree) :
     Params_ffs_r<Number, TW, Mobility, Reaction, VelDissipation>(ptree),
-    YScaling(ptree.get("domain.YScaling",1.0))
+    YScaling(ptree.get("domain.YScaling",1.0)),
+    ionSource(ptree.get("Initial.ionSource",1.0))
   {
     ModifyParameters();
     std::cout << "YScaling " << YScaling << " Eps " << this->eps << std::endl;

src/paper_thinstrip_mixed/tmixed_upscaledc.hh

@@ -81,11 +81,11 @@ template<typename Param, typename PAdaptor, typename CGFS, typename CContainer,
         RF source = 0;
         if(xg > 0.1 && xg < 0.3)
         {
-          source = (xg-0.1)*(0.3-xg);
+          source = (xg-0.1)*(0.3-xg)/0.04;
         }
         r.accumulate(lfsv,0, eg.geometry().volume()*
           ((phic*c - phicOld*cOld + param.phys.uAst*(phiSolid-phiSolidOld))/param.time.dt
-          + source));
+          - phicOld*param.ionSource*source));
         //std::cout << "bye alpha" << std::endl;
       } //alpha volume
 
@@ -152,8 +152,8 @@ template<typename Param, typename PAdaptor, typename CGFS, typename CContainer,
 
         // contribution to residual on inside and outside elements
         auto dudn = (outsideC-insideC)/distance;
-        RF peclet = 5;  //TODO THIS HAS TO BE FIXED IN OPTIONS!!!!!!!!!!
-        auto DiffusionCoeff = 1/peclet * 0.5*(pAdapt.evalPhiC(inside_global) + pAdapt.evalPhiC(outside_global));
+        //RF peclet = 100;  //TODO THIS HAS TO BE FIXED IN OPTIONS!!!!!!!!!!
+        auto DiffusionCoeff = param.phys.D * 0.5*(pAdapt.evalPhiC(inside_global) + pAdapt.evalPhiC(outside_global));
         r_i.accumulate(lfsv_i,0,-dudn*DiffusionCoeff + flux);
         r_o.accumulate(lfsv_o,0, dudn*DiffusionCoeff - flux);
         //std::cout << "bye skeleton" << std::endl;