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

Added first example

parent f7b9e9d4
No related branches found
No related tags found
No related merge requests found
......@@ -15,3 +15,4 @@ add_subdirectory(paper_thinstrip_full)
add_subdirectory(paper_thinstrip_mixed)
add_subdirectory(summerschool_inflow)
add_subdirectory(pntest)
add_subdirectory(example_2p_cahnhilliard)
#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_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)
{
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("ch2p.ini",ptree);
ptreeparser.readOptions(argc,argv,ptree);
const int refinement = ptree.get<int>("domain.level");
typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
typedef Grid::ctype DF;
std::vector<int> boundary_index_map;
std::vector<int> element_index_map;
std::string filename = ptree.get("domain.filename", "square.msh");
Dune::GridFactory<Grid> factory;
Dune::GmshReader<Grid>::read(factory,filename,true,true);
std::shared_ptr<Grid> gridp(factory.createGrid());
typedef Grid::LeafGridView GV;
GV gv = gridp->leafGridView();
using ES = Dune::PDELab::AllEntitySet<GV>;
ES es(gv);
const int l = 2;
typedef Dune::PDELab::PkLocalFiniteElementMap<ES,DF,double,l> PHI_FEM;
PHI_FEM phiFem(es);
driver_2p(*gridp, gv, es, 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;
}
}
#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_2p_cahnhilliard.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/ff_chns.hh>
#include <dune/phasefield/initial_fn.hh>
#include <dune/phasefield/debug/vector_debug.hh>
//_________________________________________________________________________________________________________________________
template<typename Grid, typename GV, typename ES, typename PHI_FEM>
void driver_2p(Grid& grid, GV& gv, ES& es, const PHI_FEM& phiFem, Dune::ParameterTree& ptree, Dune::MPIHelper& helper)
{
const int dim = GV::dimension;
typedef double RF;
// We use the Paramter File for the Cahn--Hilliard Navier--Stokes Model,
//as it includes all relevant parameters of the Cahn--Hilliard Model
typedef Params_ff< RF, DoubleWell_poly<RF> , Mobility_Const<RF> > Parameters;
Parameters param(ptree);
TimestepManager timestepManager(ptree.sub("Time"),param.time,2);
//_________________________________________________________________________________________________________________________
// Make grid function spaces
// Again, we use the Cahn--Hilliard Navier--Stokes Model, in operator splitting form the Cahn--Hilliard grid function space
// is exactly the space we need.
typedef GFS_ff_chns<RF, ES, PHI_FEM, PHI_FEM, PHI_FEM> GFS;
GFS gfs(es,phiFem,phiFem,phiFem);
typedef typename GFS::NS NS_GFS;
typedef typename GFS::CH CH_GFS;
//_________________________________________________________________________________________________________________________
//Constraints
typedef Bdry_ff_chns<2, Dune::PDELab::DirichletConstraintsParameters, // v_x, v_y
Dune::PDELab::NoDirichletConstraintsParameters, // p
Dune::PDELab::NoDirichletConstraintsParameters, // phi
Dune::PDELab::NoDirichletConstraintsParameters> // mu
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_ff_chns< ZeroInitial<GV,RF,Parameters>, //v
ZeroInitial<GV,RF,Parameters>, //p
PhiInitial_SpinodalDecomposition<GV,RF,Parameters>, //phi
ZeroInitial<GV,RF,Parameters>> Initial; //mu
Initial initial(gv,param);
Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
CH_V chVOld = chV;
//_________________________________________________________________________________________________________________________
// LocalOperators
typedef FEM_2P_CahnHilliard<Parameters, CH_GFS, CH_V > CH_LOP;
CH_LOP chLop(param, gfs.ch, chVOld);
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_ff_chns<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);
// Refine Intial Values
for(int i=0; i<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,ptree.sub("Refinement"));
rInd0.calculateIndicator(gv,chV,false);
rInd0.markGrid(grid,2);
rInd0.refineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1);
//rInd0.myRefineGrid(grid,gv,chV,nsV,gfs.ch,gfs.ns,1); // if used with PeriodicGridView
Dune::PDELab::interpolate(initial.iniCH , gfs.ch, chV);
}
gfs.updateConstraintsContainer(bdry);
interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
//_________________________________________________________________________________________________________________________
printVector(chV,10,"chV");
printVector(nsV,10,"nsV");
//____________________________________________________________________________________________________________________
chVOld = chV;
nsVOld = nsV;
while (!timestepManager.isFinished())
{
// Grid Operators
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,600,200,3,1);
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.doOnlyCHTimestep(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);
rInd.refineGrid(grid,chV,chVOld,gfs.ch,2);
//rInd.myRefineGrid(grid,gv,chV,chVOld,nsV,nsVOld,gfs.ch,gfs.ns,1); //For PeriodicGridView
// recompute constraints
gfs.updateConstraintsContainer(bdry);
interpolateToBdry(gfs,initial.iniCH,chV);
interpolateToBdry(gfs,initial.iniCH,chVOld);
}
}
add_executable(2p_cahnhilliard_equation 2p_cahnhilliard_equation.cc)
dune_symlink_to_source_files(FILES grids ch2p.ini)
[Physics]
uAst = 1
D = 0.02
rho1 = 1.0
rho2 = 1.0
mu=0.005
ReactionRate=0.5
#ReactionLimiter = 1000
SurfaceTensionScaling=0.1
VelDissipationRate = 1000
#VelDissipationMax = 0.9
CurvatureAlpha=0.001
[domain]
level = 0
filename = grids/tper1.msh
len = 1
periodic = 01
#YScalingFactor = 0.5
YScaling = 1 #0.125
[output]
filename = output1
[Phasefield]
eps=0.03 #0.03
delta=0.03 #0.03
DoubleWellDelta=0.05 #0.05
Mobility=0.04
Sigma1=1
Sigma2=1
Sigma3=1
[Time]
tMax = 1
dt = 0.0015
dtMax = 0.01
dtMin = 0.00015
#saveintervall = 0.5
#dtIncreaseFactor=1.2
#dtDecreaseFactor=0.5
[Initial]
# FlowVyy= 10# pConst/(Domain_width * mu)
pConst=0.05 #0.1
uConst=0.5 #0.5
uInflow=0.3 #0.3
[RefinementInitial]
etaRefine = 0.2
etaCoarsen = 0.18
maxRefineLevel = 4
minRefineLevel = -1
[Refinement]
etaRefine = 0.2
etaCoarsen = 0.1
maxRefineLevel = 4
minRefineLevel = -1
[NS_Newton]
ReassembleThreshold = 0.0
LineSearchMaxIterations = 30
MaxIterations = 20
AbsoluteLimit = 1e-9
Reduction = 1e-7
VerbosityLevel = 2
[Phasefield_Newton]
ReassembleThreshold = 0.0
LineSearchMaxIterations = 30
MaxIterations = 100
AbsoluteLimit = 1e-9
Reduction = 1e-7
VerbosityLevel = 2
DGF
INTERVAL
0 0 % lower left corner
1 1 % upper right corner
8 8 % number of cells in each direction
#
PERIODICFACETRANSFORMATION
1 0, 0 1 + 1 0 % make periodic in x-direction
1 0, 0 1 + 0 1 % make periodic in y-direction
#
BOUNDARYDOMAIN
DEFAULT 1
#
GRIDPARAMETER
OVERLAP 0
#
#
DGF
INTERVAL
0 0 % lower left corner
1 1 % upper right corner
8 8 % number of cells in each direction
#
PERIODICFACETRANSFORMATION
1 0, 0 1 + 1 0 % make periodic in x-direction
#
BOUNDARYDOMAIN
DEFAULT 1
#
GRIDPARAMETER
OVERLAP 0
#
#
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment