Skip to content
Snippets Groups Projects
Commit 5d4ba304 authored by Lars von Wolff's avatar Lars von Wolff
Browse files

Added support for ff, ffs, fs. Next up: Proper inflow/outflow bc

parent 1767541e
No related branches found
No related tags found
No related merge requests found
#pragma once
// C includes
#include<sys/stat.h>
// C++ includes
#include<iostream>
// Dune includes
#include<dune/grid/io/file/vtk.hh>
#include<dune/pdelab/newton/newton.hh>
#include <dune/phasefield/localoperator/fem_fs_chns_r_cahnhilliard.hh>
#include <dune/phasefield/localoperator/fem_fs_chns_r_instatstokes.hh>
#include <dune/phasefield/localoperator/pf_diff_estimator.hh>
#include <dune/phasefield/utility/ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES.hh>
#include <dune/phasefield/chns_base.hh>
#include <dune/phasefield/fs_chns_r.hh>
#include <dune/phasefield/initial_fn.hh>
#include <dune/phasefield/boundary_fn.hh>
#include <dune/phasefield/debug/vector_debug.hh>
//_________________________________________________________________________________________________________________________
template<typename Grid, typename ES, typename V_FEM, typename P_FEM, typename PHI_FEM>
void driver_fs (Grid& grid, const ES& es, const V_FEM& vFem, const P_FEM pFem, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
{
typedef typename Grid::LeafGridView GV;
//std::bitset<2> periodic_directions(std::string("01"));
//typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
GV gv = grid.leafGridView();
const int dim = GV::dimension;
typedef double RF;
typedef Params_fs_r< RF,
DoubleWell_poly<RF>,
Mobility_Const<RF>,
Reaction_Quadratic_Limited<RF>,
VelDissipation_Quadratic<RF> > Parameters;
Parameters param(ptree);
TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
//_________________________________________________________________________________________________________________________
// Make grid function spaces
typedef GFS_fs_chns_r<RF, ES, V_FEM, P_FEM, PHI_FEM, PHI_FEM> GFS;
GFS gfs(es,vFem,pFem,phiFem,phiFem);
typedef typename GFS::NS NS_GFS;
typedef typename GFS::CH CH_GFS;
//_________________________________________________________________________________________________________________________
//Constraints
typedef Bdry_fs_chns_r<dim, InflowOutflow_v,
InflowOutflow_p,
InflowOutflow_phi,
Dune::PDELab::NoDirichletConstraintsParameters,
InflowOutflow_u> Bdry;
Bdry bdry;
gfs.updateConstraintsContainer(bdry);
//_________________________________________________________________________________________________________________________
//Coefficient vectors
using CH_V = Dune::PDELab::Backend::Vector<CH_GFS, RF>;
using NS_V = Dune::PDELab::Backend::Vector<NS_GFS, RF>;
CH_V chV(gfs.ch);
NS_V nsV(gfs.ns);
//_________________________________________________________________________________________________________________________
//Initial Values
typedef Ini_fs_chns_r< VelocityInitial_Inflow<GV,RF,Parameters,dim>,
ZeroInitial<GV,RF,Parameters>,
PhiInitial_Sphere<GV,RF,Parameters>,
ZeroInitial<GV,RF,Parameters>,
UInitial_Inflow<GV,RF,Parameters> > Initial;
Initial initial(gv,param);
Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
Dune::PDELab::interpolate(initial.iniNS , gfs.ns, nsV);
//____________________________________________________________________________________________________________________
// Refine Initial Data
for(int i=0; i<4; i++)
{
// mark elements for refinement
typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
rInd.calculateIndicator(gv,chV,false);
rInd.markGrid(grid,2);
// do refinement
auto chT = transferSolutions(gfs.ch,4,chV);
auto nsT = transferSolutions(gfs.ns,4,nsV);
Dune::PDELab::adaptGrid(grid, chT,nsT);
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);
CH_V chVOld = chV;
NS_V nsVOld = nsV;
//_________________________________________________________________________________________________________________________
// LocalOperators
typedef FEM_FS_CHNS_R_InstatStokes<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> NS_LOP;
NS_LOP nsLop(param, gfs.ch, chV, chVOld, gfs.ns, nsVOld);
typedef FEM_FS_CHNS_R_CahnHilliard<Parameters, CH_GFS, CH_V, NS_GFS, NS_V> CH_LOP;
CH_LOP chLop(param, gfs.ch, chVOld, gfs.ns, nsV);
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
MBE mbe(10);// Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
//_________________________________________________________________________________________________________________________
// Construct discrete Grid Functions
typedef DGF_fs_chns_r<GFS, NS_V, CH_V> DGF;
DGF dgf(gfs, nsV, chV);
//_________________________________________________________________________________________________________________________
// prepare VTK writer
std::string filename=ptree.get("output.filename","output");
struct stat st;
if( stat( filename.c_str(), &st ) != 0 )
{
int stat = 0;
stat = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
if( stat != 0 && stat != -1)
std::cout << "Error: Cannot create directory " << filename << std::endl;
}
auto stationaryvtkwriter = std::make_shared<Dune::SubsamplingVTKWriter<GV> >(gv,0);
Dune::VTKSequenceWriter<GV> vtkwriter(stationaryvtkwriter,filename.c_str(),filename.c_str(),"");
dgf.addToVTKWriter(vtkwriter);
vtkwriter.write(param.time.t,Dune::VTK::ascii);
//_________________________________________________________________________________________________________________________
printVector(chV,10,"chV");
printVector(nsV,10,"nsV");
//____________________________________________________________________________________________________________________
//int i_it=0;
while (!timestepManager.isFinished())
{
// Grid Operators
typedef Dune::PDELab::GridOperator<NS_GFS,NS_GFS,NS_LOP,MBE,RF,RF,RF,typename GFS::NS_CC,typename GFS::NS_CC> NS_GO;
NS_GO nsGo(gfs.ns,gfs.nsCC,gfs.ns,gfs.nsCC,nsLop,mbe);
typedef Dune::PDELab::GridOperator<CH_GFS,CH_GFS,CH_LOP,MBE,RF,RF,RF,typename GFS::CH_CC,typename GFS::CH_CC> CH_GO;
CH_GO chGo(gfs.ch,gfs.chCC,gfs.ch,gfs.chCC,chLop,mbe);
// Linear solver
typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
<CH_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> CH_LS; CH_LS ls_ch(chGo,500,3,1);
typedef ISTLBackend_NOVLP_BASE_PREC_FOR_GMRES
<NS_GO, Dune::SeqILU0, Dune::RestartedGMResSolver> NS_LS; NS_LS ls_ns(nsGo,500,3,1);
// Nonlinear solver
typedef Dune::PDELab::Newton<NS_GO,NS_LS,NS_V> PDESolver1;
PDESolver1 ns_newton(nsGo,nsV,ls_ns);
ns_newton.setParameters(ptree.sub("NS_Newton"));
typedef Dune::PDELab::Newton<CH_GO,CH_LS,CH_V> PDESolver2;
PDESolver2 ch_newton(chGo,chV,ls_ch);
ch_newton.setParameters(ptree.sub("Phasefield_Newton"));
//_________________________________________________________________________________________________________________________
// do time step
interpolateToBdry(gfs,initial.iniNS,initial.iniCH, nsV, chV);
int result = timestepManager.doTimestep(ns_newton,ch_newton,vtkwriter);
if(result <0)
{
return;
}
else if(result == 0)
{
nsV = nsVOld;
chV = chVOld;
continue;
}
chVOld = chV;
nsVOld = nsV;
//_________________________________________________________________________________________________________________________
// mark elements for refinement
typedef RefinementIndicator<GV,RF,PF_DIFF_ESTIMATOR<CH_GFS>,CH_GFS> RIND;
RIND rInd(gv,gfs.ch,ptree.sub("Refinement"));
rInd.calculateIndicator(gv,chV,false);
rInd.markGrid(grid,2);
printVector(chV,3,"chV");
printVector(nsV,3,"nsV");
// do refinement
std::cout << "adapt grid and solution" << std::endl;
auto chTransfer = transferSolutions(gfs.ch,4,chV,chVOld);
auto nsTransfer = transferSolutions(gfs.ns,4,nsV,nsVOld);
Dune::PDELab::adaptGrid(grid, chTransfer,nsTransfer);
std::cout << "adapt grid and solution ended" << std::endl;
printVector(chV,3,"chV");
printVector(nsV,3,"nsV");
//vtkwriter.write(timestepManager.t+1e-6,Dune::VTK::ascii);
// recompute constraints
gfs.updateConstraintsContainer(bdry);
interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsVOld,chVOld);
//vtkwriter.write(timestepManager.t+2e-6,Dune::VTK::ascii);
}
}
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include<iostream>
// dune includes
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/common/timer.hh>
//#if HAVE_UG
//#include<dune/grid/uggrid.hh>
//#endif
#if HAVE_DUNE_ALUGRID
#include<dune/alugrid/grid.hh>
#include<dune/alugrid/dgf.hh>
#include<dune/grid/io/file/dgfparser/dgfparser.hh>
#endif
// pdelab includes
#include<dune/pdelab/finiteelementmap/pkfem.hh>
#include "2p_driver.hh"
//===============================================================
// Main program with grid setup
//===============================================================
int main(int argc, char** argv)
{
#if !HAVE_SUPERLU
std::cout << "Error: These examples work only if SuperLU is available." << std::endl;
exit(1);
#endif
#if HAVE_UG
try{
// Maybe initialize Mpi
Dune::MPIHelper&
helper = Dune::MPIHelper::instance(argc, argv);
if(Dune::MPIHelper::isFake)
std::cout<< "This is a sequential program." << std::endl;
else
std::cout << "Parallel code run on "
<< helper.size() << " process(es)" << std::endl;
// open ini file
Dune::ParameterTree ptree;
Dune::ParameterTreeParser ptreeparser;
ptreeparser.readINITree("2pchns_r.ini",ptree);
ptreeparser.readOptions(argc,argv,ptree);
// read ini file
// const int dim = ptree.get<int>("grid.dim");
const int refinement = ptree.get<int>("domain.level");
if (true/*dim==2*/)
{
typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
//typedef Dune::UGGrid<2> Grid;
typedef Grid::ctype DF;
std::vector<int> boundary_index_map;
std::vector<int> element_index_map;
std::string filename = ptree.get("domain.filename",
"turbtube2d.msh");
Dune::GridFactory<Grid> factory;
if (helper.rank()==0)
Dune::GmshReader<Grid>::read(factory,filename,boundary_index_map,element_index_map,true,false);
std::shared_ptr<Grid> gridp(factory.createGrid());
//gridp->globalRefine(refinement);
//gridp->loadBalance();
//std::cout << " after load balance /" << helper.rank() << "/ " << gridp->size(0) << std::endl;
//typedef Grid::LeafGridView LGV;
//std::bitset<2> periodic_directions(std::string("01"));
//typedef Dune::XT::Grid::PeriodicGridView<LGV> GV;
//GV gv(gridp->leafGridView(),periodic_directions);
typedef Grid::LeafGridView GV;
GV gv = gridp->leafGridView();
std::cout << "OK" << std::endl;
using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
ES es(gv);
const int k = 2;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k> V_FEM;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,k-1> P_FEM;
const int l = 1;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
V_FEM vFem(es);
P_FEM pFem(es);
PHI_FEM phiFem(es);
driver_fs(*gridp, es,vFem,pFem,phiFem, ptree, helper);
}
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;
return 1;
}
catch (...){
std::cerr << "Unknown exception thrown!" << std::endl;
return 1;
}
#endif
}
[physics]
uAst = 2
D = 0.005
rho1 = 1.0
rho2 = 1.0
mu=0.01
ReactionRate=0 #=3
#ReactionLimiter = 1000
#VelDissipationRate = 10
CurvatureAlpha=0.03
[domain]
level = 0
#filename = grids/unitsquare.msh
#filename = grids/square_27k_del.msh
#filename = grids/square_7k_del.msh
filename = grids/rectangle_3k_per.msh
#filename = grids/rectangle_53k_per.msh
#filename = grids/square_2k_del_per.msh
[Phasefield]
eps=0.005
delta=0.02
Mobility=0.005
[Time]
tMax = 30
dt = 0.015
dtMax = 0.05
dtMin = 0.0015
#saveintervall = 0.5
#dtIncreaseFactor=1.2
#dtDecreaseFactor=0.5
[Initial]
MeanInflow=1
uConst=1.0
uInflow=1.3
Radius=0.15
[Refinement]
etaRefine = 0.2
etaCoarsen = 0.1
maxRefineLevel = 7
minRefineLevel = 0
[NS_Newton]
ReassembleThreshold = 0.0
LineSearchMaxIterations = 30
MaxIterations = 30
AbsoluteLimit = 1e-9
Reduction = 1e-7
VerbosityLevel = 1
[Phasefield_Newton]
ReassembleThreshold = 0.0
LineSearchMaxIterations = 30
MaxIterations = 30
VerbosityLevel = 1
add_executable(2p_equation 2p_equation.cc)
dune_symlink_to_source_files(FILES grids 2pchns_r.ini)
#pragma once
// C includes
#include<sys/stat.h>
// C++ includes
#include<math.h>
#include<iostream>
// dune-common includes
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/common/timer.hh>
// dune-geometry includes
#include<dune/geometry/referenceelements.hh>
#include<dune/geometry/quadraturerules.hh>
// dune-grid includes
#include<dune/grid/onedgrid.hh>
#include<dune/grid/yaspgrid.hh>
#include<dune/grid/utility/structuredgridfactory.hh>
#include<dune/grid/io/file/vtk.hh>
#include<dune/grid/io/file/gmshreader.hh>
#if HAVE_UG
#include<dune/grid/uggrid.hh>
#endif
#if HAVE_DUNE_ALUGRID
#include<dune/alugrid/grid.hh>
#include<dune/alugrid/dgf.hh>
#include<dune/grid/io/file/dgfparser/dgfparser.hh>
#endif
// dune-istl included by pdelab
// dune-pdelab includes
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/common/vtkexport.hh>
#include<dune/pdelab/finiteelementmap/pkfem.hh>
#include<dune/pdelab/finiteelementmap/qkfem.hh>
#include<dune/pdelab/finiteelementmap/p0fem.hh>
#include<dune/pdelab/constraints/common/constraints.hh>
#include<dune/pdelab/constraints/common/constraintsparameters.hh>
#include<dune/pdelab/constraints/conforming.hh>
#include<dune/pdelab/function/callableadapter.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/subspace.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include<dune/pdelab/gridfunctionspace/interpolate.hh>
#include<dune/pdelab/gridfunctionspace/vtk.hh>
#include<dune/pdelab/gridoperator/gridoperator.hh>
#include<dune/pdelab/gridoperator/onestep.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/backend/istl.hh>
#include<dune/pdelab/stationary/linearproblem.hh>
#include<dune/pdelab/instationary/onestep.hh>
#include<dune/pdelab/newton/newton.hh>
#include<dune/pdelab/adaptivity/adaptivity.hh>
//#include <dune/xt/grid/grids.hh>
//#include <dune/xt/grid/view/periodic.hh>
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment