diff --git a/dune.module b/dune.module
index 0d03071577edd79376d29371d4f919b1c1e598ce..c2db464d3fdf65bd783732a858ef4fd2fc58acba 100644
--- a/dune.module
+++ b/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
diff --git a/src/paper_thinstrip_full/thinstrip_parameters.hh b/src/paper_thinstrip_full/thinstrip_parameters.hh
index 8a41e71c4749d2b09212f5233f6c89a039815263..1c2232f84f7639d1938786cbd9c6ee3602b9034d 100644
--- a/src/paper_thinstrip_full/thinstrip_parameters.hh
+++ b/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;
diff --git a/src/paper_thinstrip_mixed/chns.ini b/src/paper_thinstrip_mixed/chns.ini
index 1d14a993120d0753f0053d22175ac302280ad3ff..7723e3c1b3c675901a00b3e7c303ffca63b912db 100644
--- a/src/paper_thinstrip_mixed/chns.ini
+++ b/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
diff --git a/src/paper_thinstrip_mixed/chns_reactive.ini b/src/paper_thinstrip_mixed/chns_reactive.ini
new file mode 100644
index 0000000000000000000000000000000000000000..5acf69a20caffa1d4869ab7994fbdc39bd624a80
--- /dev/null
+++ b/src/paper_thinstrip_mixed/chns_reactive.ini
@@ -0,0 +1,68 @@
+[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
diff --git a/src/paper_thinstrip_mixed/driver_mixed.hh b/src/paper_thinstrip_mixed/driver_mixed.hh
index 69c830d28659839523879c0ec1b1757f41214fd4..5b201f03f66eabe0ba62bf23bc2616eceb95c113 100644
--- a/src/paper_thinstrip_mixed/driver_mixed.hh
+++ b/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);
diff --git a/src/paper_thinstrip_mixed/tmixed_boundary.hh b/src/paper_thinstrip_mixed/tmixed_boundary.hh
index 60af9cd2f3aa4e0841e51de1f9c1fd8e8304f45a..bc7c4efb3ebed18440163473194b2b2bb0976765 100644
--- a/src/paper_thinstrip_mixed/tmixed_boundary.hh
+++ b/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)));
diff --git a/src/paper_thinstrip_mixed/tmixed_parameters.hh b/src/paper_thinstrip_mixed/tmixed_parameters.hh
index 3fdeef307dfc3a681e836bf63000709df124951e..0c48e1e04c3bd8d6482b44223f345366104c46b8 100644
--- a/src/paper_thinstrip_mixed/tmixed_parameters.hh
+++ b/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;
diff --git a/src/paper_thinstrip_mixed/tmixed_upscaledc.hh b/src/paper_thinstrip_mixed/tmixed_upscaledc.hh
index bea9c727208e7948e46f74b911128d83e2453ced..a7620a26fa635df401d22467e2f813900d210e40 100644
--- a/src/paper_thinstrip_mixed/tmixed_upscaledc.hh
+++ b/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;