DUNE PDELab (2.8)
recipe-linear-system-solution-istl.cc
See explanation at Solving a linear system using ISTL
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/filledarray.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/pdelab.hh>
template<typename GV, typename RF>
{
typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
public:
{
return true;
}
{
typename Traits::PermTensorType I;
return I;
}
{
typename Traits::RangeType v(0.0);
return v;
}
{
return 0.0;
}
typename Traits::RangeFieldType
{
return 0.0;
}
/* return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet for Dirichlet boundary conditions
* return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann for flux boundary conditions
* return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Outflow for outflow boundary conditions
*/
BCType
bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
}
typename Traits::RangeFieldType
{
typename Traits::DomainType x = e.geometry().global(xlocal);
if (x[0] < 0.5)
return 1.0;
return 0.0;
}
typename Traits::RangeFieldType
j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return 0.0;
}
typename Traits::RangeFieldType
o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return 0.0;
}
};
int main(int argc, char **argv)
{
// initialize MPI, finalize is done automatically on exit
Dune::MPIHelper::instance(argc,argv);
// define parameters
typedef double NumberType;
// need a grid in order to test grid functions
constexpr unsigned int dim = 1;
std::array<int,dim> N(Dune::filledArray<dim,int>(5));
Grid grid(L,N);
// make grid
typedef Dune::YaspGrid<dim> GM;
// make problem parameters
typedef GenericEllipticProblem<typename GM::LeafGridView,NumberType> Problem;
Problem problem;
BCType bctype(grid.leafGridView(),problem);
typedef typename GM::ctype DF;
FEM fem(grid.leafGridView());
Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed,1>> GFS;
GFS gfs(grid.leafGridView(),fem);
typedef typename GFS::template ConstraintsContainer<NumberType>::Type CC;
CC cc;
Dune::PDELab::constraints(bctype,gfs,cc);
// local operator for finite elemenent problem
LOP lop(problem);
// [Grid operator]
auto go = GO(gfs,cc,gfs,cc,lop,MBE(nonzeros));
// [Make degree of freedom vector]
typedef Dune::PDELab::Backend::Vector<GFS,NumberType> X;
X x(gfs,0.0);
// [Set it to match boundary conditions]
G g(grid.leafGridView(),problem);
Dune::PDELab::interpolate(g,gfs,x);
// [Assemble residual]
X d(gfs,0.0);
go.residual(x,d);
// [Assemble matrix]
typedef GO::Jacobian M;
M A(go);
go.jacobian(x,A);
// [Extract ISTL types]
typedef Dune::PDELab::Backend::Native<M> ISTLM;
typedef Dune::PDELab::Backend::Native<X> ISTLX;
// [Access to ISTL objects]
using Dune::PDELab::Backend::native;
// [Define operator]
Operator matrixop(native(A));
// [Define preconditioner and solver]
Dune::SeqJac<ISTLM,ISTLX,ISTLX> preconditioner(native(A), 1, .5);
Dune::CGSolver<ISTLX> solver(matrixop, preconditioner, 1e-3, 10, 2);
// [Run solver]
X v(gfs,0.0);
solver.apply(native(v), native(d), res);
// [Subtract defect correction]
x -= v;
// [Solution output]
}
Dirichlet Constraints construction.
Definition: conforming.hh:38
Definition: convectiondiffusionparameter.hh:221
Definition: convectiondiffusionparameter.hh:325
Definition: convectiondiffusionfem.hh:47
Definition: qkfem.hh:20
Definition: recipe-geometry-grid.cc:44
Traits::RangeFieldType o(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
outflow boundary condition
Definition: recipe-geometry-grid.cc:116
Traits::PermTensorType A(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
tensor diffusion coefficient
Definition: recipe-geometry-grid.cc:58
BCType bctype(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
boundary condition type function
Definition: recipe-geometry-grid.cc:95
static constexpr bool permeabilityIsConstantPerCell()
tensor diffusion constant per cell? return false if you want more than one evaluation of A per cell.
Definition: recipe-geometry-grid.cc:51
Traits::RangeType b(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
velocity field
Definition: recipe-geometry-grid.cc:69
Traits::RangeFieldType c(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
sink term
Definition: recipe-geometry-grid.cc:77
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
Dirichlet boundary condition value.
Definition: recipe-geometry-grid.cc:102
Traits::RangeFieldType j(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
flux boundary condition
Definition: recipe-geometry-grid.cc:109
Traits::RangeFieldType f(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
source term
Definition: recipe-geometry-grid.cc:84
Utility to generate an array with a certain value.
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:749
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
void printvector(std::ostream &s, const V &v, std::string title, std::string rowtext, int columns=1, int width=10, int precision=2)
Print an ISTL vector.
Definition: io.hh:86
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:76
Helpers for dealing with MPI.
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
Statistics about the application of an inverse operator.
Definition: solver.hh:46
Traits class for convection diffusion parameters.
Definition: convectiondiffusionparameter.hh:76
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: convectiondiffusionparameter.hh:93
@ dimDomain
dimension of the domain
Definition: convectiondiffusionparameter.hh:83
Dune::FieldMatrix< RangeFieldType, dimDomain, dimDomain > PermTensorType
permeability tensor type
Definition: convectiondiffusionparameter.hh:102
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: convectiondiffusionparameter.hh:99
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: convectiondiffusionparameter.hh:105
RF RangeFieldType
Export type for range field.
Definition: convectiondiffusionparameter.hh:96
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: convectiondiffusionparameter.hh:90
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:188
![Creative Commons License](https://i.creativecommons.org/l/by-sa/4.0/80x15.png)