From b43ec9f2f519ac615ac9397d4d324726738b9a39 Mon Sep 17 00:00:00 2001
From: Lars von Wolff <lars.von-wolff@ians.uni-stuttgart.de>
Date: Tue, 19 Jan 2021 13:48:49 +0100
Subject: [PATCH] Porenetwork 1p first mockup finished
---
dune/phasefield/base_parameters.hh | 18 ++
dune/phasefield/boundary_rectangular.hh | 148 +++++++++++++
dune/phasefield/ffs_chns_r.hh | 16 --
dune/phasefield/initial_fn.hh | 2 +-
dune/phasefield/porenetwork/1p_chns.hh | 29 ++-
.../porenetwork/fem_1p_cahnhilliard_r.hh | 6 +-
.../porenetwork/fem_1p_navierstokes.hh | 2 +-
dune/phasefield/porenetwork/pn_1pflow.hh | 204 ++++++++++++------
dune/phasefield/porenetwork/pn_initial.hh | 33 +++
.../thinstrip_boundary.hh | 102 +--------
src/pntest/pn1p.ini | 19 +-
src/pntest/pntest.cc | 21 +-
src/summerschool/2p_driver.hh | 17 +-
src/summerschool/ss_calcite.ini | 4 +-
14 files changed, 404 insertions(+), 217 deletions(-)
create mode 100644 dune/phasefield/boundary_rectangular.hh
create mode 100644 dune/phasefield/porenetwork/pn_initial.hh
diff --git a/dune/phasefield/base_parameters.hh b/dune/phasefield/base_parameters.hh
index bb24908..86f1a1d 100644
--- a/dune/phasefield/base_parameters.hh
+++ b/dune/phasefield/base_parameters.hh
@@ -56,6 +56,22 @@ public:
//______________________________________________________________________________________________________
+template<typename Number>
+class Reaction_Const
+{
+public:
+ const Number reactionRate;
+ Reaction_Const(const Dune::ParameterTree & param) :
+ reactionRate(param.get("ReactionRate",(Number)1))
+ {
+ }
+
+ Number eval( const Number& u) const
+ {
+ return reactionRate;
+ }
+};
+
template<typename Number>
class Reaction_Quadratic_Limited
{
@@ -260,6 +276,7 @@ public:
Params_Initial initial;
+ const Dune::ParameterTree& pTree;
const Dune::ParameterTree& nsNewtonTree;
const Dune::ParameterTree& chNewtonTree;
@@ -269,6 +286,7 @@ public:
phys(ptree.sub("Physics"),eps),
time(0,0),
initial(ptree.sub("Initial")),
+ pTree(ptree),
nsNewtonTree(ptree.sub("NS_Newton")),
chNewtonTree(ptree.sub("Phasefield_Newton"))
{
diff --git a/dune/phasefield/boundary_rectangular.hh b/dune/phasefield/boundary_rectangular.hh
new file mode 100644
index 0000000..3a7728c
--- /dev/null
+++ b/dune/phasefield/boundary_rectangular.hh
@@ -0,0 +1,148 @@
+
+#pragma once
+#include <cmath>
+#include <dune/pdelab/common/function.hh>
+#include <dune/pdelab/constraints/common/constraintsparameters.hh>
+
+template<typename RF>
+class DomainSize
+{
+public:
+ const RF bottom;
+ const RF top;
+ const RF left;
+ const RF right;
+ const RF domainEps;
+
+ DomainSize(RF b, RF t, RF l, RF r, RF eps):
+ bottom(b), top(t), left(l), right(r), domainEps(eps)
+ {
+ }
+};
+
+template<typename DomainSize>
+class LeftBottomBoundary : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ LeftBottomBoundary(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (( xg[0] > ds.left + ds.domainEps) || ( xg[1]>ds.bottom + ds.domainEps))
+ {
+ return false;
+ }
+ return true;
+ }
+};
+
+
+template<typename DomainSize>
+class LeftBoundary : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ LeftBoundary(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if ( xg[0] > ds.domainEps)
+ {
+ return false;
+ }
+ return true;
+ }
+};
+
+template<typename DomainSize>
+class BottomBoundary : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ BottomBoundary(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if ( xg[1]>ds.bottom + ds.domainEps)
+ {
+ return false;
+ }
+ return true;
+ }
+};
+
+template<typename DomainSize>
+class LeftRightBoundary : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ LeftRightBoundary(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ Dune::FieldVector<typename I::ctype, I::coorddimension>
+ xg = intersection.geometry().global( coord );
+ if (xg[1]>ds.bottom + ds.domainEps && xg[1]<ds.top - ds.domainEps)
+ {
+ return true;
+ }
+ return false;
+ }
+};
+
+template<typename DomainSize>
+class AllBoundary : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ AllBoundary(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ return true;
+ }
+};
+
+template<typename DomainSize>
+class NoBoundary : public Dune::PDELab::DirichletConstraintsParameters
+{
+private:
+ const DomainSize ds;
+public:
+ NoBoundary(const DomainSize d) : ds(d)
+ {
+ }
+
+ template<typename I>
+ bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
+ {
+ return false;
+ }
+};
diff --git a/dune/phasefield/ffs_chns_r.hh b/dune/phasefield/ffs_chns_r.hh
index 3e403cd..0f8a811 100644
--- a/dune/phasefield/ffs_chns_r.hh
+++ b/dune/phasefield/ffs_chns_r.hh
@@ -492,19 +492,3 @@ public:
{
}
};
-
-template<typename RF>
-class DomainSize
-{
-public:
- const RF bottom;
- const RF top;
- const RF left;
- const RF right;
- const RF domainEps;
-
- DomainSize(RF b, RF t, RF l, RF r, RF eps):
- bottom(b), top(t), left(l), right(r), domainEps(eps)
- {
- }
-};
diff --git a/dune/phasefield/initial_fn.hh b/dune/phasefield/initial_fn.hh
index 601bf05..8b98df4 100644
--- a/dune/phasefield/initial_fn.hh
+++ b/dune/phasefield/initial_fn.hh
@@ -446,7 +446,7 @@ public:
//Experiment 1
- xValues.push_back(PxToCoord_x(114)); yValues.push_back(PxToCoord_y(81)); //Extra Nucleus, to check the influence on Nucleus 1
+ //xValues.push_back(PxToCoord_x(114)); yValues.push_back(PxToCoord_y(81)); //Extra Nucleus, to check the influence on Nucleus 1
xValues.push_back(PxToCoord_x(154)); yValues.push_back(PxToCoord_y(75)); //Nucleus 1
xValues.push_back(PxToCoord_x(317)); yValues.push_back(PxToCoord_y(85));
diff --git a/dune/phasefield/porenetwork/1p_chns.hh b/dune/phasefield/porenetwork/1p_chns.hh
index 843dc73..d05c13d 100644
--- a/dune/phasefield/porenetwork/1p_chns.hh
+++ b/dune/phasefield/porenetwork/1p_chns.hh
@@ -44,7 +44,7 @@ public:
typedef typename CH::template ConstraintsContainer<RF>::Type CH_CC;
CH_CC chCC;
- GFS_1p_chns(const ES& es, const V_FEM& vFem, const P_FEM& pFem, const PHI_FEM& phiFem) :
+ GFS_1p_chns(const ES& es, const V_FEM& vFem, const PHI_FEM& phiFem) :
conphi(),conmu(),
vxGfs(es,vFem),
phiGfs(es,phiFem, conphi),
@@ -104,7 +104,7 @@ public:
template<typename VTKWRITER>
void addToVTKWriter(VTKWRITER& vtkwriter)
{
- vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VDGF> >(vxDgf,"vx"));
+ vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<VXDGF> >(vxDgf,"vx"));
vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<PhiDGF> >(phiDgf,"phi"));
vtkwriter.addVertexData(std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<MuDGF> >(muDgf,"mu"));
}
@@ -149,7 +149,7 @@ public:
BdryNS bdryNS;
BdryCH bdryCH;
- Bdry_ff_chns() :
+ Bdry_1p_chns() :
bdryVx(),
bdryPhi(),
bdryMu(),
@@ -158,3 +158,26 @@ public:
{
}
};
+
+template<typename BdryVX, typename BdryPhi, typename BdryMu>
+class Bdry_1p_chns_Explicit_Domain
+{
+public:
+ BdryVX bdryVx;
+ BdryPhi bdryPhi;
+ BdryMu bdryMu;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryVX> BdryNS;
+ typedef Dune::PDELab::CompositeConstraintsParameters<BdryPhi, BdryMu> BdryCH;
+ BdryNS bdryNS;
+ BdryCH bdryCH;
+
+ template<typename DomainSize>
+ Bdry_1p_chns_Explicit_Domain(DomainSize ds) :
+ bdryVx(ds),
+ bdryPhi(ds),
+ bdryMu(ds),
+ bdryNS(bdryVx),
+ bdryCH(bdryPhi,bdryMu)
+ {
+ }
+};
diff --git a/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh b/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh
index 4815840..b34ec3e 100644
--- a/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh
+++ b/dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh
@@ -47,7 +47,7 @@ public:
{
//Quadrature Rule
const int dim = EG::Entity::dimension;
- const int order = 4; //TODO: Choose more educated
+ const int order = 5; //TODO: Choose more educated
auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);
chOldLocal.DoEvaluation(eg);
@@ -65,9 +65,9 @@ public:
{
//RF Rspeed = param.dw.eval_q((ch.pfOld-0.05)/0.7)/0.7;
- RF Rspeed = param.dw.eval_q(ch.pfOld);
+ RF Rspeed = ch.pfOld < 0.95 ? param.dw.eval_q(ch.pfOld) : 0.0f;
RF R1 = -Rspeed/param.eps*param.reaction.eval(1)*ch.basis.phi[i];
-
+ //std::cout << "mobility" << param.mobility.eval(ch.pf) << "R1" << R1 << "dt" << param.time.dt << std::endl;
RF J1 = param.mobility.eval(ch.pf) * (ch.gradxi*ch.basis.gradphi[i][0]);
// integrate: (pf-pfold)*phi_i/delta_t + gradpf * grad phi_i + ...
diff --git a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
index b5e8611..aa3294a 100644
--- a/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
+++ b/dune/phasefield/porenetwork/fem_1p_navierstokes.hh
@@ -48,7 +48,7 @@ template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, t
- TMIXED_CHNS_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
+ FEM_1P_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
NSGFS& nsGfs_) :
chLocal(chGfs_,chV_),
param(param_)
diff --git a/dune/phasefield/porenetwork/pn_1pflow.hh b/dune/phasefield/porenetwork/pn_1pflow.hh
index 544cab1..fa1e81f 100644
--- a/dune/phasefield/porenetwork/pn_1pflow.hh
+++ b/dune/phasefield/porenetwork/pn_1pflow.hh
@@ -8,17 +8,17 @@
// Dune includes
#include<dune/grid/io/file/vtk.hh>
#include<dune/pdelab/newton/newton.hh>
+#include<dune/pdelab/common/functionutilities.hh>
+#include<dune/pdelab/function/callableadapter.hh>
#include <dune/phasefield/porenetwork/fem_1p_cahnhilliard_r.hh>
#include <dune/phasefield/porenetwork/fem_1p_navierstokes.hh>
-#include "tmixed_upscaledp.hh"
-#include "tmixed_upscaledc.hh"
-#include <dune/phasefield/localoperator/pf_3p_diff_estimator.hh>
+#include <dune/phasefield/localoperator/pf_diff_estimator.hh>
#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
-#include "tmixed_parameters.hh"
-#include "tmixed_boundary.hh"
+#include <dune/phasefield/boundary_rectangular.hh>
+#include <dune/phasefield/porenetwork/pn_initial.hh>
#include <dune/phasefield/chns_base.hh>
#include <dune/phasefield/porenetwork/1p_chns.hh>
#include <dune/phasefield/initial_fn.hh>
@@ -35,7 +35,7 @@ std::shared_ptr<Grid> initPnGrid(std::string filename)
template<typename Parameters>
class Pn_1PFlow{
-public:
+private:
typedef double RF;
const int dim = 2;
@@ -44,10 +44,8 @@ public:
typedef Grid::ctype DF;
typedef Grid::LeafGridView GV;
typedef Dune::PDELab::AllEntitySet<GV> ES;
- const int k = 2;
- typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
- const int l = 1;
- typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,2> V_FEM;
+ typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,1> PHI_FEM;
typedef GFS_1p_chns<RF, ES, V_FEM, PHI_FEM> GFS;
@@ -57,12 +55,12 @@ public:
typedef Dune::PDELab::Backend::Vector<NS_GFS, RF> NS_V;
typedef Ini_1p_chns< ZeroInitial<GV,RF,Parameters>, //V
- PhiInitial_Thinstrip_Layers<GV,RF,Parameters,0>, //TODO
+ PhiInitial_Quartersquare<GV,RF,Parameters>,
ZeroInitial<GV,RF,Parameters>> Initial;
-
- typedef Bdry_1p_chns<Dune::PDELab::DirichletConstraintsParameters, //V //TODO
- Dune::PDELab::NoDirichletConstraintsParameters,
- Dune::PDELab::NoDirichletConstraintsParameters> Bdry;
+ typedef DomainSize<RF> DS;
+ typedef Bdry_1p_chns_Explicit_Domain<LeftBottomBoundary<DS>, //V
+ NoBoundary<DS>,
+ NoBoundary<DS>> Bdry;
typedef FEM_1P_NavierStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
typedef FEM_1P_CahnHilliard_R<Parameters, CH_GFS, CH_V> CH_LOP;
@@ -84,7 +82,6 @@ public:
GV gv;
ES es;
V_FEM vFem;
- P_Fem pFem;
PHI_FEM phiFem;
GFS gfs;
@@ -92,6 +89,7 @@ public:
CH_V chV, chVOld;
NS_V nsV;
Initial initial;
+ DS ds;
Bdry bdry;
MBE mbe;
@@ -108,31 +106,10 @@ public:
std::unique_ptr<CH_LS> chLs;
std::unique_ptr<NS_Newton> nsNewton;
std::unique_ptr<CH_Newton> chNewton;
+ std::unique_ptr<Dune::VTKSequenceWriter<GV>> vtkwriter;
- Pn_1PFlow(std::string filename, Parameters& param_) :
- param(param_),
- gridp(initPnGrid<Grid>(filename)),
- gv(gridp->leafGridView()),
- es(gv),
- vFem(es), phiFem(es),
- gfs(es,vFem,phiFem),
- chV(gfs.ch), nsV(gfs.ns), chVOld(gfs.ch),
- initial(gv,param), mbe(10),
- {
- 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;
+ //------------------------------------------------------------------------------------
- //Local Operator
- nsLop = std::make_unique<NS_LOP>(param, gfs.ch, chV, gfs.ns);
- chLop = std::make_unique<CH_LOP>(param, gfs.ch, chVOld);
-
- dgf = std::make_unique<DGF>(gfs,nsV,chV);
- initializeGridOperators();
- }
void initializeGridOperators()
{
@@ -143,52 +120,55 @@ public:
chLs = std::make_unique<CH_LS>(*chGo,500,100,3,1);
nsLs = std::make_unique<NS_LS>(*nsGo,1000,200,3,1);
// Nonlinear solver
-
nsNewton = std::make_unique<NS_Newton>(*nsGo,nsV,*nsLs);
nsNewton->setParameters(param.nsNewtonTree);
chNewton = std::make_unique<CH_Newton>(*chGo,chV,*chLs);
chNewton->setParameters(param.chNewtonTree);
- printVector(chV,10,"chV");
- printVector(nsV,10,"nsV");
+ //printVector(chV,10,"chV");
+ //printVector(nsV,10,"nsV");
+ }
- updatePhiC();
- updatePhiSolid();
+ void remesh()
+ {
+ // mark elements for refinement
+ typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
+ RIND rInd(gv,gfs.ch,param.pTree.sub("Refinement"));
+ rInd.calculateIndicator(gv,chV,false);
+ rInd.markGrid(*gridp,2);
+ rInd.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,1);
+ // recompute constraints
+ gfs.updateConstraintsContainer(bdry);
+ interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
+ initializeGridOperators();
}
- void applyNs()
+ void applyNs(int verbosity = 0)
{
- int verb = 0;
- if (verb > 1) std::cout << "= Begin newtonapply NS:" << std::endl;
+ if (verbosity > 0) std::cout << "= Begin newtonapply NS:" << std::endl;
nsNewton->apply();
- nsVOld = nsV;
- updateKf();
- updateKc();
}
- void applyCh()
+ void applyCh(int verbosity = 0)
{
- phicOld = phic;
- phiSolidOld = phiSolid;
- int verb=0;
- if (verb > 1) std::cout << "= Begin newtonapply CH:" << std::endl;
- chNewton->apply();
chVOld = chV;
- updatePhiC();
- updatePhiSolid();
+ if (verbosity > 0) std::cout << "= Begin newtonapply CH:" << std::endl;
+ chNewton->apply();
+ //chVOld = chV;
+ remesh();
}
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);
+ dgf->vxDgf.evaluate(element, xlocal, w);
Dune::FieldVector<double,1> phi1;
- dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
+ dgf->phiDgf.evaluate(element, xlocal, phi1);
return Dune::FieldVector<double,1>(phi1*w);
//Dune::FieldVector<double,1> phi2;
- //dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
+ //dgf->phi2Dgf.evaluate(element, xlocal, phi2);
//return Dune::FieldVector<double,1>((phi1+phi2)*w);
});
Dune::FieldVector<double,1> integral;
@@ -200,10 +180,10 @@ public:
{
auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
Dune::FieldVector<double,1> phi1;
- dgfMicro->phiDgf.evaluate(element, xlocal, phi1);
+ dgf->phiDgf.evaluate(element, xlocal, phi1);
return Dune::FieldVector<double,1>(1-phi1);
//Dune::FieldVector<double,1> phi2;
- //dgfMicro->phi2Dgf.evaluate(element, xlocal, phi2);
+ //dgf->phi2Dgf.evaluate(element, xlocal, phi2);
//return Dune::FieldVector<double,1>(1-phi1-phi2);
});
Dune::FieldVector<double,1> integral;
@@ -211,4 +191,102 @@ public:
phiSolid = integral[0];
}
+//------------------------------------------------------------------------------------
+
+public:
+
+ //Constructor
+ Pn_1PFlow(std::string GridFilename, std::string OutputFilename, Parameters& param_) :
+ param(param_),
+ gridp(initPnGrid<Grid>(GridFilename)),
+ gv(gridp->leafGridView()),
+ es(gv),
+ 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)
+ {
+ Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
+ Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
+
+ // Refine Intial Values
+ for(int i=0; i<param.pTree.get("Refinement.maxRefineLevel",2); i++)
+ {
+ // mark elements for refinement
+ typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND0;
+ RIND0 rInd0(gv,gfs.ch,param.pTree.sub("Refinement"));
+ rInd0.calculateIndicator(gv,chV,false);
+ rInd0.markGrid(*gridp,2);
+ rInd0.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,1);
+ //rInd0.myRefineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1);
+ 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;
+
+ //Local Operator
+ nsLop = std::make_unique<NS_LOP>(param, gfs.ch, chV, gfs.ns);
+ chLop = std::make_unique<CH_LOP>(param, gfs.ch, chVOld);
+
+ dgf = std::make_unique<DGF>(gfs,nsV,chV);
+
+ struct stat st;
+ if( stat( OutputFilename.c_str(), &st ) != 0 )
+ {
+ int stat = 0;
+ stat = mkdir(OutputFilename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
+ if( stat != 0 && stat != -1)
+ std::cout << "Error: Cannot create directory " << OutputFilename << std::endl;
+ }
+ auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
+ vtkwriter = std::make_unique<Dune::VTKSequenceWriter<GV> >(stationaryvtkwriter,OutputFilename.c_str(),OutputFilename.c_str(),"");
+ dgf->addToVTKWriter(*vtkwriter);
+
+
+ initializeGridOperators();
+ applyNs();
+ updatePhiSolid();
+ updateKf();
+ }
+
+ void writeVTK(RF time)
+ {
+ vtkwriter->write(time,Dune::VTK::ascii);
+ }
+
+ void grow(RF length)
+ {
+ if(length <= 0)
+ {
+ std::cout << "Need positive length" << std::endl;
+ return;
+ }
+ RF dtmax = param.pTree.get("Time.dtMax",0.01);
+ while(length > dtmax)
+ {
+ writeVTK(1-length);
+ length -= dtmax;
+ param.time.dt = dtmax;
+ applyCh(0);
+ }
+ param.time.dt = length;
+ applyCh(0);
+ applyNs(0);
+ updateKf();
+ updatePhiSolid();
+ }
+
+ RF getKf()
+ {
+ return kf;
+ }
+
+ RF getPhiSolid()
+ {
+ return phiSolid;
+ }
};
diff --git a/dune/phasefield/porenetwork/pn_initial.hh b/dune/phasefield/porenetwork/pn_initial.hh
new file mode 100644
index 0000000..df0c621
--- /dev/null
+++ b/dune/phasefield/porenetwork/pn_initial.hh
@@ -0,0 +1,33 @@
+#pragma once
+
+template<typename GV, typename RF, typename Params>
+class PhiInitial_Quartersquare :
+ public Dune::PDELab::AnalyticGridFunctionBase<
+ Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>, PhiInitial_Quartersquare<GV,RF,Params> >
+{
+ const Params& param;
+ RF eps;
+ RF r;
+
+public:
+ typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+
+ PhiInitial_Quartersquare(const GV & gv, const Params & params) :
+ Dune::PDELab::AnalyticGridFunctionBase<Traits,PhiInitial_Quartersquare<GV,RF,Params> > (gv),
+ param(params)
+ {
+ eps = params.eps;
+ r = params.initial.ptree.get("Radius",0.2);
+ std::cout << "EPS" << eps;
+ }
+
+ inline void evaluateGlobal(const typename Traits::DomainType & x, typename Traits::RangeType & y) const
+ {
+ RF rf1 = param.dw.eval_shape((x[0]-r)*1/eps);
+ RF rf2 = param.dw.eval_shape((x[1]-r)*1/eps);
+ //RF rf1 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[0]-r)));
+ //RF rf2 = 0.5*(1+tanh(1/eps/sqrt(2)*(x[1]-r)));
+
+ y=rf1*rf2;
+ }
+};
diff --git a/src/paper_thinstrip_full/thinstrip_boundary.hh b/src/paper_thinstrip_full/thinstrip_boundary.hh
index 4ae0251..6b90749 100644
--- a/src/paper_thinstrip_full/thinstrip_boundary.hh
+++ b/src/paper_thinstrip_full/thinstrip_boundary.hh
@@ -5,112 +5,12 @@
#include <dune/pdelab/constraints/common/constraintsparameters.hh>
#include <dune/phasefield/initial_fn.hh>
+#include <dune/phasefield/boundary_rectangular.hh>
#define GRID_TOP_THINSTRIP 1
-template<typename DomainSize>
-class LeftBoundary : public Dune::PDELab::DirichletConstraintsParameters
-{
-private:
- const DomainSize ds;
-public:
- LeftBoundary(const DomainSize d) : ds(d)
- {
- }
-
- template<typename I>
- bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
- {
- Dune::FieldVector<typename I::ctype, I::coorddimension>
- xg = intersection.geometry().global( coord );
- if ( xg[0] > ds.domainEps)
- {
- return false;
- }
- return true;
- }
-};
-template<typename DomainSize>
-class BottomBoundary : public Dune::PDELab::DirichletConstraintsParameters
-{
-private:
- const DomainSize ds;
-public:
- BottomBoundary(const DomainSize d) : ds(d)
- {
- }
-
- template<typename I>
- bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
- {
- Dune::FieldVector<typename I::ctype, I::coorddimension>
- xg = intersection.geometry().global( coord );
- if ( xg[1]>ds.bottom + ds.domainEps)
- {
- return false;
- }
- return true;
- }
-};
-
-template<typename DomainSize>
-class LeftRightBoundary : public Dune::PDELab::DirichletConstraintsParameters
-{
-private:
- const DomainSize ds;
-public:
- LeftRightBoundary(const DomainSize d) : ds(d)
- {
- }
-
- template<typename I>
- bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
- {
- Dune::FieldVector<typename I::ctype, I::coorddimension>
- xg = intersection.geometry().global( coord );
- if (xg[1]>ds.bottom + ds.domainEps && xg[1]<ds.top - ds.domainEps)
- {
- return true;
- }
- return false;
- }
-};
-
-template<typename DomainSize>
-class AllBoundary : public Dune::PDELab::DirichletConstraintsParameters
-{
-private:
- const DomainSize ds;
-public:
- AllBoundary(const DomainSize d) : ds(d)
- {
- }
-
- template<typename I>
- bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
- {
- return true;
- }
-};
-
-template<typename DomainSize>
-class NoBoundary : public Dune::PDELab::DirichletConstraintsParameters
-{
-private:
- const DomainSize ds;
-public:
- NoBoundary(const DomainSize d) : ds(d)
- {
- }
-
- template<typename I>
- bool isDirichlet(const I & intersection, const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
- {
- return false;
- }
-};
template<typename GV, typename RF, typename Params, int dim>
class VelocityInitial_InOutFlow :
diff --git a/src/pntest/pn1p.ini b/src/pntest/pn1p.ini
index c7addb3..4a4a76d 100644
--- a/src/pntest/pn1p.ini
+++ b/src/pntest/pn1p.ini
@@ -4,40 +4,29 @@ D = 0.02
rho1 = 1.0
rho2 = 1.0
mu=0.005
-ReactionRate=1.5
+ReactionRate=1#1.5
#ReactionLimiter = 1000
SurfaceTensionScaling=0.1
VelDissipationRate = 100
CurvatureAlpha=0#0.001
[domain]
-level = 0
filename = grids/square_periodic.msh
-#YScalingFactor = 0.5
-YScaling = 1
-YResolution = 200
-XResolution = 100
[output]
filename = output2
[Phasefield]
eps=0.03
-delta=0.00001#0.03
+delta=0.001#0.03
DoubleWellDelta=0.05
-Mobility=1#0.01
+Mobility=0.0005
Sigma1=1
Sigma2=1
Sigma3=1
[Time]
-tMax = 100
-dt = 0.0001
-dtMax = 0.0001
-dtMin = 0.000015
-saveintervall = 0.005
-#dtIncreaseFactor=1.2
-#dtDecreaseFactor=0.5
+dtMax = 0.002
[Initial]
# FlowVyy= 10# pConst/(Domain_width * mu)
diff --git a/src/pntest/pntest.cc b/src/pntest/pntest.cc
index 2e6aab1..2900352 100644
--- a/src/pntest/pntest.cc
+++ b/src/pntest/pntest.cc
@@ -15,7 +15,8 @@
// pdelab includes
#include<dune/pdelab/finiteelementmap/pkfem.hh>
#include<dune/pdelab/finiteelementmap/qkfem.hh>
-#include<dune/pdelab/porenetwork/1p_chns.hh>
+// phasefield includes
+#include<dune/phasefield/porenetwork/pn_1pflow.hh>
//===============================================================
// Main program
@@ -40,20 +41,26 @@ int main(int argc, char** argv)
ptreeparser.readINITree("pn1p.ini",ptree);
ptreeparser.readOptions(argc,argv,ptree);
+ typedef double RF;
typedef Params_fs_r< RF,
//DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,
DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> >,
Mobility_Const<RF>,
- Reaction_Thinstrip<RF>,
+ Reaction_Const<RF>,
VelDissipation_Quadratic<RF>
//VelDissipation_Quadratic_Shifted<RF>
> Parameters;
Parameters param(ptree);
- std::string filename = ptree.get("domain.filename","square.msh");
- Pn_1PFlow<Parameters> pn_1pFlow(filename,param);
- pn_1pFlow.applyCh();
- pn_1pFlow.applyNs();
- pn_1pFlow.updateKf();
+ std::string GridFilename = ptree.get("domain.filename","square.msh");
+ std::string OutputFilename = ptree.get("output.filename","output");
+ Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param);
+
+ std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
+ pn_1pFlow.writeVTK(0.0);
+ pn_1pFlow.grow(0.2);
+ pn_1pFlow.writeVTK(1.0);
+ std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
+
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;
diff --git a/src/summerschool/2p_driver.hh b/src/summerschool/2p_driver.hh
index b1a6e70..9245d9e 100644
--- a/src/summerschool/2p_driver.hh
+++ b/src/summerschool/2p_driver.hh
@@ -246,7 +246,7 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
//_________________________________________________________________________________________________________________________
// WRITE CSV OUTPUT
std::string csvname=ptree.get("output.csvname","centers");
- for( int radiusFactor = 3; radiusFactor<9; radiusFactor++)
+ for( int radiusFactor = 10; radiusFactor<19; radiusFactor++)
{
try
{
@@ -257,7 +257,7 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
auto analyticFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& x){
return ( (x[0]-initial.iniPhi.xValues[i])*(x[0]-initial.iniPhi.xValues[i])
+ (x[1]-initial.iniPhi.yValues[i])*(x[1]-initial.iniPhi.yValues[i]) )
- < radiusFactor/2.0*radiusFactor/2.0*initial.iniPhi.radius*initial.iniPhi.radius?
+ < radiusFactor/10.0*radiusFactor/10.0*initial.iniPhi.radius*initial.iniPhi.radius?
1.0:0.0;
});
auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
@@ -265,18 +265,23 @@ void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, c
analyticFunction.evaluate(element, xlocal, y1);
Dune::FieldVector<double,1> y2;
dgf.phiDgf.evaluate(element, xlocal, y2);
+ Dune::FieldVector<double,2> vv;
+ dgf.vDgf.evaluate(element, xlocal, vv);
Dune::FieldVector<double, dim> xg = element.geometry().global( xlocal );
- Dune::FieldVector<double,3> result;
+ Dune::FieldVector<double,5> result;
result[0] = (1-y2[0])*y1[0];
result[1] = result[0]*xg[0];
result[2] = result[0]*xg[1];
+ result[3] = y2[0]*y1[0];
+ result[4] = result[3]*std::sqrt(vv[0]*vv[0]+vv[1]*vv[1]);
return result;
});
- Dune::FieldVector<double,3> integral;
+ Dune::FieldVector<double,5> integral;
Dune::PDELab::integrateGridFunction(productFunction,integral,2);
- std::cout << "Integral: " << integral[0] << ", " << integral[1]/integral[0] << ", " << integral[2]/integral[0]<< std::endl;
- myfile << integral[0] << ", " << integral[1]/integral[0] << ", " << integral[2]/integral[0] << ", ";
+ std::cout << "Integral: " << integral[0] << ", " << integral[1]/integral[0] << ", " << integral[2]/integral[0]<< " -- " << integral[3] << " , " << integral[4] << " , " << integral[4]/integral[3] << std::endl;
+ //myfile << integral[0] << ", " << integral[1]/integral[0] << ", " << integral[2]/integral[0] << ", ";
+ myfile << integral[3] << ", " <<integral[4]/integral[3] << ", ";
}
myfile << param.time.t << std::endl;
myfile.close();
diff --git a/src/summerschool/ss_calcite.ini b/src/summerschool/ss_calcite.ini
index 3665b2b..cb67663 100644
--- a/src/summerschool/ss_calcite.ini
+++ b/src/summerschool/ss_calcite.ini
@@ -20,7 +20,9 @@ CellHeight=0.085
[output]
#filename = precip_christian_v2
#filename = tpexperimental
-filename = extra_nucleus
+#filename = extra_nucleus
+filename = flowaround
+csvname = test_flow
[Phasefield]
eps=0.002 #0.002
--
GitLab