Select Git revision
pn_1pflow.hh 8.88 KiB
#pragma once
// C includes
#include<sys/stat.h>
// C++ includes
#include<iostream>
#include<vector>
// 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 <dune/phasefield/localoperator/pf_diff_estimator.hh>
#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.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>
#include <dune/phasefield/debug/vector_debug.hh>
template< typename Grid>
std::shared_ptr<Grid> initPnGrid(std::string filename)
{
Dune::GridFactory<Grid> factory;
Dune::GmshReader<Grid>::read(factory,filename,true,true);
return std::shared_ptr<Grid>(factory.createGrid());
}
template<typename Parameters>
class Pn_1PFlow{
private:
typedef double RF;
const int dim = 2;
typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
typedef Grid::ctype DF;
typedef Grid::LeafGridView GV;
typedef Dune::PDELab::AllEntitySet<GV> ES;
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;
typedef typename GFS::NS NS_GFS;
typedef typename GFS::CH CH_GFS;
typedef Dune::PDELab::Backend::Vector<CH_GFS, RF> CH_V;
typedef Dune::PDELab::Backend::Vector<NS_GFS, RF> NS_V;
typedef Ini_1p_chns< ZeroInitial<GV,RF,Parameters>, //V
PhiInitial_Quartersquare<GV,RF,Parameters>,
ZeroInitial<GV,RF,Parameters>> Initial;
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;
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
typedef DGF_1p_chns<GFS, NS_V, CH_V> DGF;
typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS;
typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES<NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS;
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;
std::shared_ptr<Grid> gridp;
GV gv;
ES es;
V_FEM vFem;
PHI_FEM phiFem;
GFS gfs;
CH_V chV, chVOld;
NS_V nsV;
Initial initial;
DS ds;
Bdry bdry;
MBE mbe;
RF kf = 0;
RF phiSolid = 0;
std::unique_ptr<NS_LOP> nsLop;
std::unique_ptr<CH_LOP> chLop;
std::unique_ptr<DGF> dgf;
std::unique_ptr<NS_GO> nsGo;
std::unique_ptr<CH_GO> chGo;
std::unique_ptr<NS_LS> nsLs;
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;
//------------------------------------------------------------------------------------
void initializeGridOperators()
{
//Grid Operator
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);
// 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");
}
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(int verbosity = 0)
{
if (verbosity > 0) std::cout << "= Begin newtonapply NS:" << std::endl;
nsNewton->apply();
}
void applyCh(int verbosity = 0)
{
chVOld = chV;
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;
dgf->vxDgf.evaluate(element, xlocal, w);
Dune::FieldVector<double,1> phi1;
dgf->phiDgf.evaluate(element, xlocal, phi1);
return Dune::FieldVector<double,1>(phi1*w);
//Dune::FieldVector<double,1> phi2;
//dgf->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];
}
void updatePhiSolid()
{
auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
Dune::FieldVector<double,1> phi1;
dgf->phiDgf.evaluate(element, xlocal, phi1);
return Dune::FieldVector<double,1>(1-phi1);
//Dune::FieldVector<double,1> phi2;
//dgf->phi2Dgf.evaluate(element, xlocal, phi2);
//return Dune::FieldVector<double,1>(1-phi1-phi2);
});
Dune::FieldVector<double,1> integral;
Dune::PDELab::integrateGridFunction(productFunction,integral,2);
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;
}
};