From 6005339fdbce858e700ac0e3e96efb7d64794cb1 Mon Sep 17 00:00:00 2001
From: Lars von Wolff <lars.von-wolff@mathematik.uni-stuttgart.de>
Date: Tue, 14 Jul 2020 09:35:12 +0200
Subject: [PATCH] Final Summerschool Simulation and Mixed Thinstrip working
---
dune/phasefield/base_parameters.hh | 2 +-
dune/phasefield/initial_fn.hh | 39 +++-
src/paper_thinstrip_mixed/chns.ini | 20 +-
src/paper_thinstrip_mixed/driver_mixed.hh | 210 ++++++++++++++++--
src/paper_thinstrip_mixed/micro_space.hh | 156 ++++++++-----
src/paper_thinstrip_mixed/tmixed_boundary.hh | 10 +-
.../tmixed_chns_r_cahnhilliard.hh | 67 +++---
src/paper_thinstrip_mixed/tmixed_equation.cc | 6 +-
src/paper_thinstrip_mixed/tmixed_upscaledp.hh | 89 ++++++++
src/rectangle_periodic.geo | 4 +-
src/summerschool/ss_calcite.ini | 8 +-
src/summerschool_inflow/inflow.cc | 2 +-
src/summerschool_inflow/ss_calcite.ini | 4 +-
13 files changed, 485 insertions(+), 132 deletions(-)
create mode 100644 src/paper_thinstrip_mixed/tmixed_upscaledp.hh
diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index 9cef08f..f1b033c 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -162,7 +162,7 @@ template<typename Number>
class VelDissipation_Quadratic_Shifted
{
public:
- const Number velDissipationRate;
+ Number velDissipationRate;
const Number dMax;
VelDissipation_Quadratic_Shifted(const Dune::ParameterTree & param) :
velDissipationRate(param.get("VelDissipationRate",(Number)10)),
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index be95900..a1f760f 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -425,6 +425,14 @@ private:
return (749 - px_y)/((RF)224/1.5);
}
+ inline static RF PxToCoord_x4(const RF px_x)
+ {
+ return (px_x - 1081)/((RF)189/0.75);
+ }
+ inline static RF PxToCoord_y4(const RF px_y)
+ {
+ return (2037 - px_y)/((RF)189/0.75);
+ }
public:
typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
@@ -435,7 +443,7 @@ public:
radius = params.initial.ptree.get("Radius",0.2);
-
+/*
//Experiment 1
xValues.push_back(PxToCoord_x(154)); yValues.push_back(PxToCoord_y(75));
xValues.push_back(PxToCoord_x(317)); yValues.push_back(PxToCoord_y(85));
@@ -458,6 +466,8 @@ public:
//xValues.push_back(PxToCoord_x(1432)); yValues.push_back(PxToCoord_y(101));
xValues.push_back(PxToCoord_x(1731)); yValues.push_back(PxToCoord_y(115));
xValues.push_back(PxToCoord_x(1915)); yValues.push_back(PxToCoord_y(77));
+*/
+
/*
//Experiment 7
@@ -480,6 +490,33 @@ public:
xValues.push_back(PxToCoord_x3(548)); yValues.push_back(PxToCoord_y3(687));
xValues.push_back(PxToCoord_x3(594)); yValues.push_back(PxToCoord_y3(680));
*/
+
+
+//Inflow S1
+xValues.push_back(PxToCoord_x4(1090)); yValues.push_back(PxToCoord_y4(1858));
+xValues.push_back(PxToCoord_x4(1102)); yValues.push_back(PxToCoord_y4(1885));
+xValues.push_back(PxToCoord_x4(1126)); yValues.push_back(PxToCoord_y4(1892));
+xValues.push_back(PxToCoord_x4(1134)); yValues.push_back(PxToCoord_y4(1905));
+xValues.push_back(PxToCoord_x4(1186)); yValues.push_back(PxToCoord_y4(1901));
+xValues.push_back(PxToCoord_x4(1248)); yValues.push_back(PxToCoord_y4(1899));
+xValues.push_back(PxToCoord_x4(1296)); yValues.push_back(PxToCoord_y4(1908));
+xValues.push_back(PxToCoord_x4(1093)); yValues.push_back(PxToCoord_y4(1951));
+xValues.push_back(PxToCoord_x4(1122)); yValues.push_back(PxToCoord_y4(1936));
+xValues.push_back(PxToCoord_x4(1143)); yValues.push_back(PxToCoord_y4(1942));
+xValues.push_back(PxToCoord_x4(1095)); yValues.push_back(PxToCoord_y4(1975));
+xValues.push_back(PxToCoord_x4(1141)); yValues.push_back(PxToCoord_y4(1970));
+xValues.push_back(PxToCoord_x4(1155)); yValues.push_back(PxToCoord_y4(1985));
+xValues.push_back(PxToCoord_x4(1171)); yValues.push_back(PxToCoord_y4(1967));
+xValues.push_back(PxToCoord_x4(1230)); yValues.push_back(PxToCoord_y4(1949));
+xValues.push_back(PxToCoord_x4(1249)); yValues.push_back(PxToCoord_y4(1957));
+xValues.push_back(PxToCoord_x4(1273)); yValues.push_back(PxToCoord_y4(1963));
+xValues.push_back(PxToCoord_x4(1122)); yValues.push_back(PxToCoord_y4(2018));
+xValues.push_back(PxToCoord_x4(1134)); yValues.push_back(PxToCoord_y4(2010));
+xValues.push_back(PxToCoord_x4(1184)); yValues.push_back(PxToCoord_y4(2021));
+xValues.push_back(PxToCoord_x4(1199)); yValues.push_back(PxToCoord_y4(1998));
+xValues.push_back(PxToCoord_x4(1283)); yValues.push_back(PxToCoord_y4(2029));
+
+
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
diff --git a/src/paper_thinstrip_mixed/chns.ini b/src/paper_thinstrip_mixed/chns.ini
index b46a925..d20b964 100644
--- a/src/paper_thinstrip_mixed/chns.ini
+++ b/src/paper_thinstrip_mixed/chns.ini
@@ -7,15 +7,15 @@ mu=0.005
ReactionRate=0.5
#ReactionLimiter = 1000
SurfaceTensionScaling=0.1
-VelDissipationRate = 1000
-CurvatureAlpha=0.001
+VelDissipationRate = 100
+CurvatureAlpha=0#0.001
[domain]
level = 0
filename = grids/rectangle_periodic_coarse.msh
#YScalingFactor = 0.5
-YScaling = 0.5
-YResolution = 100
+YScaling = 1
+YResolution = 200
XResolution = 100
[output]
@@ -23,18 +23,18 @@ filename = output2
[Phasefield]
eps=0.03
-delta=0.03
+delta=0.00001#0.03
DoubleWellDelta=0.05
-Mobility=0.01
+Mobility=1#0.01
Sigma1=1
Sigma2=1
Sigma3=1
[Time]
-tMax = 1
-dt = 0.015
-dtMax = 0.15
-dtMin = 0.0015
+tMax = 100
+dt = 0.0001
+dtMax = 0.0001
+dtMin = 0.000015
#saveintervall = 0.5
#dtIncreaseFactor=1.2
#dtDecreaseFactor=0.5
diff --git a/src/paper_thinstrip_mixed/driver_mixed.hh b/src/paper_thinstrip_mixed/driver_mixed.hh
index 4c2c685..1434230 100644
--- a/src/paper_thinstrip_mixed/driver_mixed.hh
+++ b/src/paper_thinstrip_mixed/driver_mixed.hh
@@ -11,6 +11,7 @@
#include "tmixed_chns_r_cahnhilliard.hh"
#include "tmixed_chns_navierstokes.hh"
+#include "tmixed_upscaledp.hh"
#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
@@ -41,6 +42,13 @@ public:
PhiInitial_Thinstrip_Layers<GV,RF,Parameters,1>,
ZeroInitial<GV,RF,Parameters>> Initial;
+ typedef Bdry_micro<Dune::PDELab::NoDirichletConstraintsParameters, //VX=W
+ Dune::PDELab::DirichletConstraintsParameters, //VY
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters,
+ Dune::PDELab::NoDirichletConstraintsParameters> Bdry;
+
typedef TMIXED_CHNS_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
typedef TMIXED_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
@@ -56,15 +64,21 @@ 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;
+ GV& gv;
CH_V chV, chVOld, chVUpwind;
NS_V nsV, nsVOld, nsVUpwind;
Initial initial;
+ Bdry bdry;
RF xpos;
RF dxp = 1;
RF dxpUpwind = 1;
+ RF c = 1;
+ RF dx;
MBE mbe;
+ RF kf = 0;
+
std::unique_ptr<NS_LOP> nsLop;
std::unique_ptr<CH_LOP> chLop;
std::unique_ptr<DgfMicro> dgfMicro;
@@ -75,13 +89,15 @@ public:
std::unique_ptr<NS_Newton> nsNewton;
std::unique_ptr<CH_Newton> chNewton;
- YSlice(GV& gv, GFS& gfs, Parameters& param, RF xpos_) :
- 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)
+ 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),
+ initial(gv,param), xpos(xpos_), mbe(10), dx(dx_)
{
param.TempXpos = xpos;
Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
chVOld = chV;
nsVOld = nsV;
@@ -90,7 +106,7 @@ public:
//Local Operator
nsLop = std::make_unique<NS_LOP>(param, gfs.ch, chV, gfs.ns);
- chLop = std::make_unique<CH_LOP>(param, gfs.ch, chVOld, chVUpwind, gfs.ns, nsV, nsVUpwind, dxp, dxpUpwind);
+ chLop = std::make_unique<CH_LOP>(param, gfs.ch, chVOld, chVUpwind, gfs.ns, nsV, nsVUpwind, dxp, dxpUpwind, c, dx);
dgfMicro = std::make_unique<DgfMicro>(gfs,nsV,chV);
@@ -113,25 +129,83 @@ public:
printVector(nsV,10,"nsV");
}
- void applyTimestep()
+ void applyNs()
{
int verb = 2;
if (verb > 1) std::cout << "= Begin newtonapply NS:" << std::endl;
nsNewton->apply();
+ nsVOld = nsV;
+ }
- if (verb > 1) std::cout << "= Begin newtonapply CH:" << std::endl;
+ void applyCh()
+ {
+ int verb=2;
+ if (verb > 1) std::cout << "= Begin newtonapply CH:" << " with dxp= " << dxp << " dxpUpwind= " << dxpUpwind << " c= " << c << std::endl;
chNewton->apply();
-
chVOld = chV;
- nsVOld = nsV;
+ }
+
+ void updateKf()
+ {
+ auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+ Dune::FieldVector<double,1> w;
+ dgfMicro->vxDgf.evaluate(element, xlocal, w);
+ Dune::FieldVector<double,1> phi1;
+ 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> integral;
+ Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+ kf = integral[0];
}
};
+template<typename RF, typename YSimS>
+class P_ParamAdaptor
+{
+public:
+ const int xs;
+ YSimS& ySimS;
+ P_ParamAdaptor(YSimS& ySimS_) : ySimS(ySimS_), xs(ySimS_.size())
+ {
+ }
+
+ RF evalKf(const auto& xg)
+ {
+ int i = std::max(0, std::min(xs-1, (int)(xg[0]*xs + 1e-6) ));
+ return ySimS[i].kf;
+ }
+};
+
+//_________________________________________________________________________________________________________________________
+
+template< typename XGV, typename PGrad_DGF, typename SimS>
+void evaluateDxP(XGV& xgv, PGrad_DGF& pGradDgf, SimS& simS)
+{
+ Dune::FieldVector<double,1> pgrad(0);
+ Dune::FieldVector<double,1> pgradUpwind(0);
+ int xs = simS.size();
+ auto it = elements(xgv).begin();
+ for(int i=0; i<xs; i++)
+ {
+ pGradDgf.evaluate(*it,it->geometry().local(it->geometry().center()), pgrad );
+ simS[i].dxp = pgrad[0];
+ simS[i].dxpUpwind = pgradUpwind[0];
+ it++;
+ pgradUpwind = pgrad;
+ }
+ //Periodic:
+ simS[0].dxpUpwind = simS[xs-1].dxp;
+}
//_________________________________________________________________________________________________________________________
template<typename GV, typename ES, typename TGV, typename T_FEM, typename V_FEM, typename P_FEM, typename PHI_FEM>
-void driver_mixed(GV& xgv, int xsize, GV& ygv, ES& y_es, TGV& tgv, const T_FEM& tFem, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
+void driver_mixed(GV& xgv, const P_FEM pFem, int xsize,
+ GV& ygv, ES& y_es, TGV& tgv, const T_FEM& tFem, const V_FEM& vFem, const PHI_FEM& phiFem,
+ Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
{
const int dim = GV::dimension;
@@ -143,13 +217,15 @@ void driver_mixed(GV& xgv, int xsize, GV& ygv, ES& y_es, TGV& tgv, const T_FEM&
TripleWell<RF,DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> > >,
Mobility_Const<RF>,
Reaction_Thinstrip<RF>,
- VelDissipation_Quadratic<RF> > Parameters;
+ VelDissipation_Quadratic<RF>
+ //VelDissipation_Quadratic_Shifted<RF>
+ > Parameters;
Parameters param(ptree);
TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
//_________________________________________________________________________________________________________________________
-// Make grid function spaces
+// Make micro problems
typedef GFS_micro<RF, ES, V_FEM, PHI_FEM> GFS;
@@ -161,13 +237,69 @@ void driver_mixed(GV& xgv, int xsize, GV& ygv, ES& y_es, TGV& tgv, const T_FEM&
ySimS.reserve(xsize);
for(int i=0; i<xsize; i++)
{
- ySimS.emplace_back(ygv, gfs, param, ((RF)i)/xsize);
+ ySimS.emplace_back(ygv, gfs, param, ((RF)i)/((RF)xsize), ((RF)1)/((RF)xsize) );
}
+ //_________________________________________________________________________________________________________________________
+ // Make macro problems
+ using XES = Dune::PDELab::NonOverlappingEntitySet<GV>;
+ XES x_es(xgv);
+
+ typedef Dune::PDELab::ISTL::VectorBackend<> VectorBackend;
+ typedef Dune::PDELab::ConformingDirichletConstraints ConstraintsAssembler;
+ ConstraintsAssembler conP();
+ typedef Dune::PDELab::GridFunctionSpace<XES, P_FEM, ConstraintsAssembler, VectorBackend> P_GFS;
+ P_GFS pGfs(x_es,pFem);
+
+ auto blambda = [](const auto& x){return true;};
+ auto b = Dune::PDELab::
+ makeBoundaryConditionFromCallable(xgv,blambda);
+ //typedef Dune::PDELab::DirichletConstraintsParameters BdryP;
+ //BdryP pBdry();
+ typedef typename P_GFS::template ConstraintsContainer<RF>::Type P_CC;
+ P_CC pCC;
+ //pCC.clear();
+ Dune::PDELab::constraints(b,pGfs,pCC);
+
+ typedef Dune::PDELab::Backend::Vector<P_GFS, RF> P_V;
+ P_V pV(pGfs);
+
+ auto p0lambda = [&](const auto& xg)
+ {return (1-xg[0])*(1-xg[0]);};
+ auto p0 = Dune::PDELab::
+ makeGridFunctionFromCallable(xgv,p0lambda);
+ Dune::PDELab::interpolate(p0,pGfs,pV);
+
+
+
+
+ typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
+ MBE mbe(10);
+ typedef P_ParamAdaptor<RF, YSimS> P_PAdaptor;
+ P_PAdaptor pPAdaptor(ySimS);
+ typedef TMIXED_UpscaledP<P_PAdaptor,P_FEM,RF> P_LOP;
+ P_LOP pLop(pPAdaptor);
+
+ typedef Dune::PDELab::GridOperator<P_GFS, P_GFS, P_LOP, MBE,
+ RF, RF, RF, P_CC, P_CC> P_GO;
+ P_GO pGo(pGfs, pCC, pGfs, pCC, pLop, mbe);
+ typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<P_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> P_LS;
+ P_LS pLs(pGo,500,100,3,1);
+ typedef Dune::PDELab::Newton<P_GO,P_LS,P_V> P_Newton;
+ P_Newton pNewton(pGo,pV,pLs);
+ pNewton.setReassembleThreshold(0.0); // always reassemble J
+ pNewton.setVerbosityLevel(3); // be verbose
+ pNewton.setReduction(1e-10); // total reduction
+ pNewton.setMinLinearReduction(1e-4); // min. red. in lin. solve
+ pNewton.setMaxIterations(25); // limit number of its
+ pNewton.setLineSearchMaxIterations(10); // limit line search
+ pNewton.setAbsoluteLimit(1e-10);
+
//_________________________________________________________________________________________________________________________
// prepare VTK writer
std::string filename=ptree.get("output.filename","output");
+ std::string x_filename = "xoutput";
struct stat st;
if( stat( filename.c_str(), &st ) != 0 )
{
@@ -176,27 +308,69 @@ void driver_mixed(GV& xgv, int xsize, GV& ygv, ES& y_es, TGV& tgv, const T_FEM&
if( stat != 0 && stat != -1)
std::cout << "Error: Cannot create directory " << filename << std::endl;
}
+ if( stat( x_filename.c_str(), &st ) != 0 )
+ {
+ int stat = 0;
+ stat = mkdir(x_filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+ if( stat != 0 && stat != -1)
+ std::cout << "Error: Cannot create directory " << x_filename << std::endl;
+ }
auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<TGV> >(tgv,0);
Dune::VTKSequenceWriter<TGV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
- DGF_macro<RF, GV, TGV, YSimS, Dune::VTKSequenceWriter<TGV>> dgfMacro(ygv, xsize,tgv,ySimS,vtkwriter);
-
+ DGF_macro<RF, GV, TGV, YSimS> dgfMacro(ygv, tgv,ySimS);
+ dgfMacro.addToVTKWriter(vtkwriter);
+ vtkwriter.write(param.time.t,Dune::VTK::base64);
+ auto x_stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(xgv,0);
+ Dune::VTKSequenceWriter<GV> x_vtkwriter(x_stationaryvtkwriter,x_filename.c_str(),x_filename.c_str(),"");
+ typedef Dune::PDELab::DiscreteGridFunction<P_GFS,P_V> P_DGF;
+ P_DGF pDgf(pGfs,pV);
+ typedef Dune::PDELab::DiscreteGridFunctionGradient<P_GFS,P_V> PGrad_DGF;
+ PGrad_DGF pGradDgf(pGfs,pV);
+ typedef Dune::PDELab::VTKGridFunctionAdapter<P_DGF> P_VTKF;
+ x_vtkwriter.addVertexData(std::shared_ptr<P_VTKF>(
+ new P_VTKF(pDgf,"p")));
+ x_vtkwriter.write(param.time.t,Dune::VTK::base64);
+//_________________________________________________________________________________________________________________________
while (!timestepManager.isFinished())
{
- for(int i=1; i<xsize; i++)
+ for(int i=1; i<xsize; i++)
{
ySimS[i].chVUpwind = ySimS[i-1].chV;
+ }
+ ySimS[0].chVUpwind = ySimS[xsize-1].chV; //Periodic
+
+ for(int i=0; i<xsize; i++)
+ {
+ ySimS[i].applyNs();
+ ySimS[i].updateKf();
+ }
+
+ for(int i=1; i<xsize; i++)
+ {
ySimS[i].nsVUpwind = ySimS[i-1].nsV;
}
- for(int i=0; i<3; i++)
+ ySimS[0].nsVUpwind = ySimS[xsize-1].nsV; //Periodic
+
+ //For Debug
+ //timestepManager.applyTimestep(vtkwriter);
+ //x_vtkwriter.write(param.time.t);
+
+
+ std::cout << "Begin newtonapply Pressure:" << std::endl;
+ pNewton.apply();
+ evaluateDxP(xgv, pGradDgf,ySimS);
+ for(int i=0; i<xsize; i++)
{
- ySimS[i].applyTimestep();
+ ySimS[i].applyCh();
}
- timestepManager.applyTimestep();
+
+ timestepManager.applyTimestep(vtkwriter);
+ x_vtkwriter.write(param.time.t);
}
//_________________________________________________________________________________________________________________________
diff --git a/src/paper_thinstrip_mixed/micro_space.hh b/src/paper_thinstrip_mixed/micro_space.hh
index a0dca7f..d9f5a21 100644
--- a/src/paper_thinstrip_mixed/micro_space.hh
+++ b/src/paper_thinstrip_mixed/micro_space.hh
@@ -113,8 +113,8 @@ public:
typedef Dune::PDELab::DiscreteGridFunction<VYSubGFS,CH_V> VYDGF;
VYDGF vyDgf;
- typedef PF_3P_Color_DGF<typename GFS::CH,CH_V,3> ColorDGF;
- ColorDGF colorDgf;
+ //typedef PF_3P_Color_DGF<typename GFS::CH,CH_V,3> ColorDGF;
+ //ColorDGF colorDgf;
DGF_micro(const GFS& gfs, const NS_V& nsV, const CH_V& chV) :
@@ -131,96 +131,138 @@ public:
mu2Dgf(mu2SubGfs,chV),
vySubGfs(gfs.ch),
- vyDgf(vySubGfs,chV),
+ vyDgf(vySubGfs,chV)
- colorDgf(gfs.ch,chV)
+ //colorDgf(gfs.ch,chV)
{
}
-/* 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<UDGF> >(uDgf,"u"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ColorDGF> >(colorDgf,"_coloring_"));
- }*/
};
-template<typename RF, typename GV, typename TGV, typename YSimS, typename VTKWRITER>
-class DGF_macro
+
+
+template<typename RF, typename GV, typename YSimS>
+class DGF_evaluator
{
public:
- //const int xsize;
- //int getSimSIndex(RF xglobal)
- //{
- // return std::max(0, std::min(1, (int)(xsize*xglobal) ));
- //}
const GV& ygv;
const int xs;
YSimS& ySimS;
- const TGV& tgv;
+ const int ys;
+ typedef Dune::FieldVector<RF,1> ScalarField;
+ typedef Dune::FieldVector<RF,1> LocalC;
+ typedef Dune::FieldVector<RF,2> GlobalC;
struct Pos{
int i;
- Dune::FieldVector<RF,1> localCoord;
+ LocalC localCoord;
typename GV::Traits::template Codim<0> :: Iterator it;
};
Pos FindPosition(const auto& xyglobal)
{
- int i = std::max(0, std::min(1, (int)(xs*xyglobal[0]) ));
- int j = (int)(xs*xyglobal[1]);
+ int i = std::max(0, std::min(xs-1, (int)(xs*xyglobal[0] + 0.5) ));
+ int j = std::max(0, std::min(ys-1, (int)(ys*xyglobal[1]) ));
+ //std::cout << "i" << i <<"j" << j << "xs" << xs << "ys" << ys << "bouns" << ySimS.size() << std::endl;
auto it = elements(ygv).begin();
for(int k=0; k<j; k++){ ++it; }
- Dune::FieldVector<RF,1> yglob(xyglobal[1]);
- Dune::FieldVector<RF,1> localCoord = it->geometry().local(yglob);
+ LocalC yglob(xyglobal[1]);
+ LocalC localCoord = it->geometry().local(yglob);
return Pos{i,localCoord,it};
}
- Dune::FieldVector<RF,1> EvalPhi(const auto& xyglobal)
+ /*ScalarField EvalPhi(const GlobalC& xyglobal)
{
Pos p = FindPosition(xyglobal);
- Dune::FieldVector<RF,1> result(0);
+ ScalarField result(0);
ySimS[p.i].dgfMicro->phiDgf.evaluate(*(p.it),p.localCoord,result);
return result;
+ }*/
+
+ DGF_evaluator(const GV& ygv_, YSimS& ySimS_)
+ : ygv(ygv_), ySimS(ySimS_), xs(ySimS_.size()), ys(ygv.size(0))
+ {
+ }
+};
+
+
+template<typename GV, typename RF, typename Evl, int Variable>
+class MacroToMicroGridFun
+ : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<
+ GV, RF, 1, Dune::FieldVector<RF,1> >,
+ MacroToMicroGridFun<GV,RF,Evl,Variable>
+ >
+{
+public:
+ Evl& evaluator;
+ GV gv;
+
+ typedef Dune::PDELab::
+ GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> > Traits;
+ MacroToMicroGridFun(GV gv_, Evl& evl_): gv(gv_), evaluator(evl_)
+ {
}
- DGF_macro(const GV& ygv_, const int xs_, const TGV& tgv_, YSimS& ySimS_, VTKWRITER& vtkwriter)
- : ygv(ygv_), xs(xs_), ySimS(ySimS_), tgv(tgv_)
+ inline void evaluate (const typename Traits::ElementType& e,
+ const typename Traits::DomainType& x,
+ typename Traits::RangeType& y) const
+ {
+ auto p = evaluator.FindPosition(e.geometry().global(x));
+ switch(Variable)
+ {
+ case 0: evaluator.ySimS[p.i].dgfMicro->phiDgf.evaluate(*(p.it),p.localCoord,y);break;
+ case 1: evaluator.ySimS[p.i].dgfMicro->muDgf.evaluate(*(p.it),p.localCoord,y);break;
+ case 2: evaluator.ySimS[p.i].dgfMicro->phi2Dgf.evaluate(*(p.it),p.localCoord,y);break;
+ case 3: evaluator.ySimS[p.i].dgfMicro->mu2Dgf.evaluate(*(p.it),p.localCoord,y);break;
+ case 4: evaluator.ySimS[p.i].dgfMicro->vyDgf.evaluate(*(p.it),p.localCoord,y);break;
+ case 5: evaluator.ySimS[p.i].dgfMicro->vxDgf.evaluate(*(p.it),p.localCoord,y);break;
+ }
+ }
+
+ inline const GV& getGridView () const {return gv;}
+};
+
+
+template<typename RF, typename GV, typename TGV, typename YSimS>
+class DGF_macro
+{
+public:
+ typedef DGF_evaluator<RF,GV,YSimS> Evl;
+ Evl evaluator;
+ MacroToMicroGridFun<TGV,RF,Evl,0> phiGridFun;
+ MacroToMicroGridFun<TGV,RF,Evl,1> muGridFun;
+ MacroToMicroGridFun<TGV,RF,Evl,2> phi2GridFun;
+ MacroToMicroGridFun<TGV,RF,Evl,3> mu2GridFun;
+ MacroToMicroGridFun<TGV,RF,Evl,4> vyGridFun;
+ MacroToMicroGridFun<TGV,RF,Evl,5> vxGridFun;
+
+ DGF_macro(const GV& ygv_, const TGV& tgv_, YSimS& ySimS_)
+ : evaluator(ygv_,ySimS_),
+ phiGridFun(tgv_,evaluator),
+ muGridFun(tgv_,evaluator),
+ phi2GridFun(tgv_,evaluator),
+ mu2GridFun(tgv_,evaluator),
+ vyGridFun(tgv_,evaluator),
+ vxGridFun(tgv_,evaluator)
{
- auto PhiFunction = Dune::PDELab::makeGridFunctionFromCallable (tgv, [&](const auto& xyglobal){
- int i = std::max(0, std::min(1, (int)(xs*xyglobal[0]) ));
- int j = (int)(xs*xyglobal[1]);
- auto it = elements(ygv).begin();
- for(int k=0; k<j; k++){ ++it; }
- Dune::FieldVector<RF,1> yglob(xyglobal[1]);
- Dune::FieldVector<RF,1> result(0);
- Dune::FieldVector<RF,1> localCoord = it->geometry().local(yglob);
- ySimS[i].dgfMicro->phiDgf.evaluate(*it,localCoord,result);
- return result;
- });
- vtkwriter.addVertexData(
- std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<decltype(PhiFunction)> >(PhiFunction,"v"));
- vtkwriter.write(0,Dune::VTK::ascii);
}
-/* template<typename VTKWRITER>
+ 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<UDGF> >(uDgf,"u"));
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ColorDGF> >(colorDgf,"_coloring_"));
- }*/
+ vtkwriter.addVertexData(std::make_shared<
+ Dune::PDELab::VTKGridFunctionAdapter<MacroToMicroGridFun<TGV,RF,Evl,0>> >(phiGridFun,"phi"));
+ vtkwriter.addVertexData(std::make_shared<
+ Dune::PDELab::VTKGridFunctionAdapter<MacroToMicroGridFun<TGV,RF,Evl,1>> >(muGridFun,"mu"));
+ vtkwriter.addVertexData(std::make_shared<
+ Dune::PDELab::VTKGridFunctionAdapter<MacroToMicroGridFun<TGV,RF,Evl,2>> >(phi2GridFun,"phi2"));
+ vtkwriter.addVertexData(std::make_shared<
+ Dune::PDELab::VTKGridFunctionAdapter<MacroToMicroGridFun<TGV,RF,Evl,3>> >(mu2GridFun,"mu2"));
+ vtkwriter.addVertexData(std::make_shared<
+ Dune::PDELab::VTKGridFunctionAdapter<MacroToMicroGridFun<TGV,RF,Evl,4>> >(vyGridFun,"vy"));
+ vtkwriter.addVertexData(std::make_shared<
+ Dune::PDELab::VTKGridFunctionAdapter<MacroToMicroGridFun<TGV,RF,Evl,5>> >(vxGridFun,"vx"));
+ }
};
//_________________________________________________________________________________________________________________________
diff --git a/src/paper_thinstrip_mixed/tmixed_boundary.hh b/src/paper_thinstrip_mixed/tmixed_boundary.hh
index 4681339..6c81e30 100644
--- a/src/paper_thinstrip_mixed/tmixed_boundary.hh
+++ b/src/paper_thinstrip_mixed/tmixed_boundary.hh
@@ -169,16 +169,18 @@ public:
eps = params.eps;
dSolid = params.initial.ptree.get("dSolid",0.3);
d1 = params.initial.ptree.get("d1",0.3);
- yS = params.YScaling;
}
inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
{
- RF rf = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid*yS));
- RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[1] - dSolid*yS - d1*yS));
+ 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;
+ //std::cout << "TempXPos" << pref.TempXpos;
+ RF rf = pref.tw.dw.eval_shape(1/eps*(x[0] - dSolid -mod));
+ RF rf2 = pref.tw.dw.eval_shape(1/eps*(x[0] - dSolid - d1 ));
if(Phase == 1)
{
- y = rf2;
+ y = rf*rf2;
return;
}
y=rf*(1-rf2);
diff --git a/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh b/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh
index 5dcf3e0..c7ccbe4 100644
--- a/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh
+++ b/src/paper_thinstrip_mixed/tmixed_chns_r_cahnhilliard.hh
@@ -36,6 +36,8 @@ private:
Param& param; // parameter functions
RF& dxp;
RF& dxpUpwind;
+ RF& c;
+ const RF dx;
public:
enum {doPatternVolume=true};
@@ -44,14 +46,16 @@ public:
//constructor
TMIXED_CHNS_R_CahnHilliard(Param& param_, CHGFS& chGfs_, CHContainer& chVOld_, CHContainer& chVUpwind_,
- NSGFS& nsGfs_, NSContainer& nsV, NSContainer& nsVUpwind_, RF& dxp_, RF& dxpUpwind_) :
+ NSGFS& nsGfs_, NSContainer& nsV, NSContainer& nsVUpwind_, RF& dxp_, RF& dxpUpwind_, RF& c_, RF& dx_) :
chOldLocal(chGfs_,chVOld_),
chUpwindLocal(chGfs_,chVUpwind_),
nsLocal(nsGfs_,nsV),
nsUpwindLocal(nsGfs_,nsVUpwind_),
param(param_),
dxp(dxp_),
- dxpUpwind(dxpUpwind_)
+ dxpUpwind(dxpUpwind_),
+ c(c_),
+ dx(dx_)
{
}
@@ -89,36 +93,52 @@ public:
//TODO
-/*
+
RF pf3 = 1-ch.pf-ch.pf2;
+ RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;
+ RF pf3Upwind = 1-chUpwind.pf-chUpwind.pf2;
+ RF pf_fUpwind = chUpwind.pf+chUpwind.pf2+2*param.delta*pf3Upwind;
+ RF pf3Old = 1-chOld.pf-chOld.pf2;
+ RF pf_fOld=chOld.pf+chOld.pf2+2*param.delta*pf3Old;
+
+ auto gradpf3 = -ch.gradpf - ch.gradpf2;
+ auto gradpf_f = ch.gradpf + ch.gradpf2+2*param.delta*gradpf3;
+
+ auto gradpf3Old = -chOld.gradpf - chOld.gradpf2;
+ auto gradpf_fOld = chOld.gradpf + chOld.gradpf2+2*param.delta*gradpf3Old;
+
RF xi3 = -param.tw.Sigma3*(ch.xi/param.tw.Sigma1 + ch.xi2/param.tw.Sigma2);
- RF pf_c = ch.pf;
- RF pf_c_Old = ch.pfOld;
- RF pf_c_reg = pf_c+param.delta;
- RF pf_c_reg_Old = pf_c_Old + param.delta;
+ RF velocityX = -dxp*ns.vx;
+ RF velocityXUpwind = -dxpUpwind*nsUpwind.vx;
for (size_t i=0; i<ch.pfspace.size(); i++)
{
//RF Rspeed = std::max(param.tw.eval_q(ch.pf,ch.pf2,pf3) - 0.1, 0.0);
- RF Rspeed = param.tw.eval_q(ch.pf,ch.pf2,pf3);
- RF R1 = -Rspeed/param.eps*
- (param.reaction.eval(ch.u)+param.phys.curvatureAlpha*(ch.xi-xi3))*ch.basis.phi[i];
+
+ //RF Rspeed = param.tw.eval_q(ch.pf,ch.pf2,pf3);
+ //RF R1 = -Rspeed/param.eps*
+ // (param.reaction.eval(c)+param.phys.curvatureAlpha*(ch.xi-xi3))*ch.basis.phi[i];
+
+ RF R1 = 0;
+
RF J1 = param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma1 * (ch.gradxi*ch.basis.gradphi[i][0]);
// integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
r.accumulate(ch.pfspace,i,factor*(
- (ch.pf-ch.pfOld)*ch.basis.phi[i]/param.time.dt
+ (ch.pf-chOld.pf)*ch.basis.phi[i]/param.time.dt
+J1
- +(ch.pf*ns.div_v + (ch.gradpf*ns.v))*ch.basis.phi[i]
+ + ( chOld.pf * velocityX - chUpwind.pf * velocityXUpwind) /dx * ch.basis.phi[i]
+ + ( (ch.pf * ch.gradvy[0]) + (ch.gradpf[0]*ch.vy) )*ch.basis.phi[i]
//-pf1*(v*gradphi[i][0]) weak transport
-R1
));
r.accumulate(ch.pf2space,i,factor*(
- (ch.pf2-ch.pf2Old)*ch.basis.phi[i]/param.time.dt
+ (ch.pf2-chOld.pf2)*ch.basis.phi[i]/param.time.dt
+param.mobility.eval(ch.pf,ch.pf2,pf3)*param.eps/param.tw.Sigma2 * (ch.gradxi2*ch.basis.gradphi[i][0])
- +(ch.pf2*ns.div_v + (ch.gradpf2*ns.v))*ch.basis.phi[i]
+ + ( chOld.pf2 * velocityX - chUpwind.pf2 * velocityXUpwind) /dx * ch.basis.phi[i]
+ + ( (ch.pf2 * ch.gradvy[0]) + (ch.gradpf2[0]*ch.vy) )*ch.basis.phi[i]
//-pf2*(v*gradphi[i][0])
));
@@ -134,22 +154,13 @@ public:
- param.eps * param.tw.Sigma2* (ch.gradpf2*ch.basis.gradphi[i][0])
));
- // integrate (pf-pfold)*(u-uast)*phi_i/delta_t + pfold*(u-uold)*phi_i/delta_t + D*pf*gradu*grad phi_i = 0
-
//Dune::FieldVector<RF, dim> xg = eg.geometry().global( ip.position() );
- r.accumulate(ch.uspace,i,factor*(
-
- //(u-param.uAst)*phi[i]/param.dt //For Testing
-
- (pf_c_reg-pf_c_reg_Old)*(ch.u-param.phys.uAst)*ch.basis.phi[i]/param.time.dt
- + pf_c_reg_Old*(ch.u-ch.uOld)*ch.basis.phi[i]/param.time.dt
- + param.phys.D*pf_c_reg_Old*(ch.gradu*ch.basis.gradphi[i][0])
- + J1*(ch.u-param.phys.uAst)
- //- pf_c*(u-param.phys.uAst)*(v*gradphi[i][0])
- + (pf_c*(ch.u-param.phys.uAst)*ns.div_v + (ch.u-param.phys.uAst)*(ch.gradpf*ns.v) + pf_c*(ch.gradu*ns.v) )*ch.basis.phi[i]
- ));
+ r.accumulate(ch.vyspace,i,factor*(
+ ( (pf_f * ch.gradvy[0]) + (gradpf_f*ch.vy) )*ch.basis.phi[i]
+ + ( pf_fOld * velocityX - pf_fUpwind * velocityXUpwind) /dx * ch.basis.phi[i]
+ ));
}
-*/
+
}
}
diff --git a/src/paper_thinstrip_mixed/tmixed_equation.cc b/src/paper_thinstrip_mixed/tmixed_equation.cc
index 2c38b53..2ad9de3 100644
--- a/src/paper_thinstrip_mixed/tmixed_equation.cc
+++ b/src/paper_thinstrip_mixed/tmixed_equation.cc
@@ -83,7 +83,6 @@ int main(int argc, char** argv)
int A = ptree.get("domain.XResolution", 10);
int B = ptree.get("domain.YResolution", 10);
- A=B; //Same resolution for now is easier
Dune::FieldVector<DF,1> L(1);
Dune::FieldVector<DF,2> LL(1);
std::shared_ptr<Grid> xgridp(new Grid(L,std::array<int,1>{A}));
@@ -103,17 +102,16 @@ int main(int argc, char** argv)
const int k = 2;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> V_FEM;
- typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> P_FEM;
const int l = 2;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> PHI_FEM;
-
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> P_FEM;
typedef Dune::PDELab::PkLocalFiniteElementMap<TGV,DF,double,1> Tensor_FEM;
V_FEM vFem(y_es);
P_FEM pFem(y_es);
PHI_FEM phiFem(y_es);
Tensor_FEM tensorFem(tgv);
- driver_mixed(xgv, A, ygv, y_es, tgv, tensorFem,vFem,pFem,phiFem, ptree, helper);
+ driver_mixed(xgv, pFem, A, ygv, y_es, tgv, tensorFem,vFem,phiFem, ptree, helper);
}
}
catch (Dune::Exception &e){
diff --git a/src/paper_thinstrip_mixed/tmixed_upscaledp.hh b/src/paper_thinstrip_mixed/tmixed_upscaledp.hh
new file mode 100644
index 0000000..300c180
--- /dev/null
+++ b/src/paper_thinstrip_mixed/tmixed_upscaledp.hh
@@ -0,0 +1,89 @@
+#pragma once
+
+#include<dune/pdelab/common/quadraturerules.hh>
+#include<dune/pdelab/finiteelement/localbasiscache.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include<dune/pdelab/localoperator/variablefactories.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+
+
+
+
+
+template<typename Param, typename FEM, typename RF>
+ class TMIXED_UpscaledP :
+ public Dune::PDELab::FullVolumePattern,
+ public Dune::PDELab::LocalOperatorDefaultFlags,
+ public Dune::PDELab::NumericalJacobianVolume<TMIXED_UpscaledP<Param, FEM, RF> >,
+ public Dune::PDELab::NumericalJacobianApplyVolume<TMIXED_UpscaledP<Param, FEM, RF> >
+ {
+ private:
+ //typedef typename CHContainer::field_type RF;
+ //typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity ENTITY;
+
+ typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
+ Dune::PDELab::LocalBasisCache<LocalBasis> cache;
+
+ Param& param; // parameter functions
+
+
+ public:
+ // pattern assembly flags
+ enum { doPatternVolume = true };
+
+ // residual assembly flags
+ enum { doAlphaVolume = true };
+ //enum { doLambdaVolume = true };
+
+
+
+ TMIXED_UpscaledP (Param& param_) :
+ param(param_)
+ {
+ }
+
+
+ template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+ void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+ {
+// std::cout << "hi_alpha" << std::endl;
+ //Quadrature Rule
+ const int dim = EG::Entity::dimension;
+ const int order = 2; //TODO: Choose more educated
+ auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
+
+ for(const auto& ip : rule)
+ {
+ // evaluate shape functions and gradient of shape functions
+ auto& phi = cache.evaluateFunction(ip.position(), lfsv.finiteElement().localBasis());
+ auto& gradphihat = cache.evaluateJacobian(ip.position(), lfsv.finiteElement().localBasis());
+
+ // transform gradients of shape functions to real element
+ const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
+ auto gradphi = makeJacobianContainer(lfsv);
+ for (size_t i=0; i<lfsv.size(); i++)
+ {
+ S.mv(gradphihat[i][0],gradphi[i][0]);
+ }
+ RF p=0.0;
+ Dune::FieldVector<RF,dim> gradp(0.0);
+ for (size_t i=0; i<lfsu.size(); i++)
+ {
+ gradp.axpy(x(lfsv,i),gradphi[i][0]);
+ p += x(lfsv,i)*phi[i];
+ }
+ auto xg = eg.geometry().global(ip.position());
+ //auto xg = eg.geometry().center();
+ //std::cout << "pos: " << xg[0] << " kf: " << kfresult-kfold << std::endl;
+ auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());
+ for (size_t i=0; i<lfsv.size(); i++)
+ {
+ r.accumulate(lfsv,i,( (gradp*gradphi[i][0])*param.evalKf(xg) )*factor);
+ } // for i
+ } // for ip
+// std::cout << "bye_alpha" << std::endl;
+ } //alpha volume
+
+ };
diff --git a/src/rectangle_periodic.geo b/src/rectangle_periodic.geo
index 9e9fc8d..434fe3b 100644
--- a/src/rectangle_periodic.geo
+++ b/src/rectangle_periodic.geo
@@ -2,8 +2,8 @@
lc = 0.015;
Point(1) = {0, 0, 0, lc};
-Point(2) = {1.5, 0, 0, lc};
-Point(3) = {1.5, 0.75, 0, lc};
+Point(2) = {0.9375, 0, 0, lc};
+Point(3) = {0.9375, 0.75, 0, lc};
Point(4) = {0, 0.75, 0, lc};
Line(1) = {2,1};
diff --git a/src/summerschool/ss_calcite.ini b/src/summerschool/ss_calcite.ini
index 0440d73..98622da 100644
--- a/src/summerschool/ss_calcite.ini
+++ b/src/summerschool/ss_calcite.ini
@@ -2,7 +2,7 @@
uAst = 27.1 #Molar density of CaCO3, in mol/l
uEquilibrium = 0.00014 #Solubility of CaCO3 in water, in mol/l
ReactionRate = 10
-BulkRate = 0.00001141 #0.00005716 #Urease reaction rate, depending un urea concentration and experimentally matched, in mol/(l*s) # Was 0.00001796
+BulkRate = 0.00005716#0.00001141 #0.00005716 #Urease reaction rate, depending un urea concentration and experimentally matched, in mol/(l*s) # Was 0.00001796
D = 0.000804 #of CO3^2- in Water see Zeebe,2011, in mm^2/s
rho1 = 0.997 #Water, in 10^-6 kg/mm^3
rho2 = 2.71 #Kalcite, in 10^-6 kg/mm^3
@@ -15,7 +15,7 @@ filename = grids/summerschool_full_v2.msh
CellHeight=0.085
[output]
-filename = centertest
+filename = finalresult
[Phasefield]
@@ -39,8 +39,8 @@ reactionOnlyStart = 20 #NEW
[Initial]
MeanInflow=1.882 #mm/s (in mean-> *1.5 for max velocity)
-uConst=0.000367667 #FROM INFLOW SIMULATION
-uInflow=0.000367667 #FROM INFLOW SIMULATION
+uConst=0.000315003 #FROM INFLOW SIMULATION
+uInflow=0.000315003 #FROM INFLOW SIMULATION
Radius=0.02
Radius1=0.015
diff --git a/src/summerschool_inflow/inflow.cc b/src/summerschool_inflow/inflow.cc
index 476b4ba..1029df3 100644
--- a/src/summerschool_inflow/inflow.cc
+++ b/src/summerschool_inflow/inflow.cc
@@ -68,7 +68,7 @@ int main(int argc, char** argv)
matrix[1][0] = 0;
std::cout << "b" << std::endl;
Dune::GridFactory<Grid>::WorldVector shift;
- shift[0] = 1;
+ shift[0] = 0.9375;
shift[1] = 0;
std::cout << "c" << std::endl;
factory.insertFaceTransformation(matrix,shift );
diff --git a/src/summerschool_inflow/ss_calcite.ini b/src/summerschool_inflow/ss_calcite.ini
index d6d5aad..30cbe3c 100644
--- a/src/summerschool_inflow/ss_calcite.ini
+++ b/src/summerschool_inflow/ss_calcite.ini
@@ -2,7 +2,7 @@
uAst = 27.1 #Molar density of CaCO3, in mol/l
uEquilibrium = 0.00014 #Solubility of CaCO3 in water, in mol/l
ReactionRate = 10
-BulkRate = 0.00001141 # 0.00001796 #Urease reaction rate, depending un urea concentration and experimentally matched, in mol/(l*s)
+BulkRate = 0.00005716 #0.00001141 # 0.00001796 #Urease reaction rate, depending un urea concentration and experimentally matched, in mol/(l*s)
D = 0.000804 #of CO3^2- in Water see Zeebe,2011, in mm^2/s
rho1 = 0.997 #Water, in 10^-6 kg/mm^3
rho2 = 2.71 #Kalcite, in 10^-6 kg/mm^3
@@ -11,7 +11,7 @@ VelDissipationRate = 1000000
CurvatureAlpha=0
[domain]
-filename = grids/inflow_periodic.msh
+filename = grids/inflow_periodic_s1.msh
CellHeight=0.085
periodic=01
--
GitLab