...
 
Commits (7)
......@@ -33,9 +33,7 @@ public:
lop.jacobian_apply_volume(eg, lfsu, x_view, lfsv, y_view);
const auto alpha = 1.0;
y.axpy(alpha * 1e-2, x);
if(p_decomp->intersects_exterior_boundary(eg)){
p_decomp->copy_exterior_boundary(eg, x, y);
}
p_decomp->copy_exterior_boundary(eg, x, y);
}
}
......
......@@ -2,7 +2,7 @@
#define BLOCKSTRUCTURED_PRECONDITIONER_DECOMPOSITION_HH
#include <dune/pdelab/localoperator/flags.hh>
#include <dune/pdelab/gridoperator/blockstructured.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/backend/istl.hh>
......@@ -12,7 +12,7 @@ class SchurDecomposition{
using RT = double;
using MB = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
using GOP = Dune::PDELab::BlockstructuredGridOperator<GFS, GFS, DomainCount, MB, DF, RT, RT, CC, CC>;
using GOP = Dune::PDELab::GridOperator<GFS, GFS, DomainCount, MB, DF, RT, RT, CC, CC>;
using X = typename GOP::Domain;
using V = typename Dune::PDELab::Backend::Native<X>;
......
......@@ -2,7 +2,7 @@
#define BLOCKSTRUCTURED_PRECONDITIONER_INTERPOLATION_HH
#include "dune/localfunctions/lagrange/lagrangecube.hh"
#include <dune/pdelab/gridoperator/blockstructured.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/finiteelement/localbasiscache.hh>
......@@ -15,7 +15,7 @@ class InterpolationOperator{
constexpr static int k = F_GFS::Traits::FiniteElementMap::blocks;
using MB = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
using GOP = Dune::PDELab::BlockstructuredGridOperator<C_GFS, F_GFS, LOP, MB, DF, RT, RT, C_CC, F_CC>;
using GOP = Dune::PDELab::GridOperator<C_GFS, F_GFS, LOP, MB, DF, RT, RT, C_CC, F_CC>;
public:
using C_X = typename GOP::Domain;
......
......@@ -2,7 +2,7 @@
#define BLOCKSTRUCTURED_PRECONDITIONER_RESTRICTION_HH
#include "dune/localfunctions/lagrange/lagrangecube.hh"
#include <dune/pdelab/gridoperator/blockstructured.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/finiteelement/localbasiscache.hh>
......@@ -14,7 +14,7 @@ class RestrictionOperator{
constexpr static int k = F_GFS::Traits::FiniteElementMap::blocks;
using MB = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
using GOP = Dune::PDELab::BlockstructuredGridOperator<F_GFS, C_GFS, LOP, MB, DF, RT, RT, F_CC, C_CC>;
using GOP = Dune::PDELab::GridOperator<F_GFS, C_GFS, LOP, MB, DF, RT, RT, F_CC, C_CC>;
public:
using F_X = typename GOP::Domain;
......
......@@ -11,12 +11,11 @@
namespace Dune{
namespace PDELab{
template <typename GO, typename V>
void solveMatrixFree(GO& go, V& x ){
template <typename GO, typename V, typename Preconditioner = Dune::Richardson<V,V>>
void solveMatrixFree(GO& go, V& x, Preconditioner& pre = Preconditioner{1.0}){
using ISTLOnTheFlyOperator = Dune::PDELab::OnTheFlyOperator<V,V,GO>;
ISTLOnTheFlyOperator opb(go);
Dune::Richardson<V,V> richardson(1.0);
Dune::BiCGSTABSolver<V> solverb(opb,richardson,1E-10,5000,2);
Dune::BiCGSTABSolver<V> solverb(opb,pre,1E-10,5000,2);
Dune::InverseOperatorResult stat;
// evaluate residual w.r.t initial guess
using TrialGridFunctionSpace = typename GO::Traits::TrialGridFunctionSpace;
......
......@@ -54,111 +54,6 @@ namespace Dune {
const X* u_;
};
template<template<class> class Solver>
class ISTLBackend_SEQ_Richardson
: public SequentialNorm, public LinearResultStorage
{
public:
/*! \brief make a linear solver object
\param[in] maxiter_ maximum number of iterations to do
\param[in] verbose_ print messages if true
*/
explicit ISTLBackend_SEQ_Richardson(unsigned maxiter_=5000, int verbose_=1)
: maxiter(maxiter_), verbose(verbose_)
{}
/*! \brief solve the given linear system
\param[in] A the given matrix
\param[out] z the solution vector to be computed
\param[in] r right hand side
\param[in] reduction to be achieved
*/
template<class M, class V, class W>
void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
{
using Backend::Native;
using Backend::native;
Dune::MatrixAdapter<Native<M>,
Native<V>,
Native<W>> opa(native(A));
Dune::Richardson<Native<V>,Native<W> > prec(0.7);
Solver<Native<V> > solver(opa, prec, reduction, maxiter, verbose);
Dune::InverseOperatorResult stat;
solver.apply(native(z), native(r), stat);
res.converged = stat.converged;
res.iterations = stat.iterations;
res.elapsed = stat.elapsed;
res.reduction = stat.reduction;
res.conv_rate = stat.conv_rate;
}
private:
unsigned maxiter;
int verbose;
};
template<class Operator, template<class> class Solver>
class ISTLBackend_SEQ_MatrixFree_Richardson
: public SequentialNorm, public LinearResultStorage
{
public:
/*! \brief make a linear solver object
\param[in] maxiter_ maximum number of iterations to do
\param[in] verbose_ print messages if true
*/
explicit ISTLBackend_SEQ_MatrixFree_Richardson(Operator& op, unsigned maxiter=5000, int verbose=1)
: opa_(op), u_(static_cast<typename Operator::domain_type*>(0))
, maxiter_(maxiter)
, verbose_(verbose)
{}
/*! \brief solve the given linear system
\param[in] A the given matrix
\param[out] z the solution vector to be computed
\param[in] r right hand side
\param[in] reduction to be achieved
*/
template<class V, class W>
void apply(V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
{
Dune::Richardson<V,W> prec(0.7);
Solver<V> solver(opa_, prec, reduction, maxiter_, verbose_);
Dune::InverseOperatorResult stat;
solver.apply(z, r, stat);
res.converged = stat.converged;
res.iterations = stat.iterations;
res.elapsed = stat.elapsed;
res.reduction = stat.reduction;
res.conv_rate = stat.conv_rate;
}
//! Set position of jacobian.
//! Must be called before apply().
void setLinearizationPoint(const typename Operator::domain_type& u)
{
u_ = &u;
opa_.setLinearizationPoint(u);
}
private:
Operator& opa_;
const typename Operator::domain_type* u_;
unsigned maxiter_;
int verbose_;
};
} // end namespace PDELab
} // end namespace Dune
#endif
......@@ -82,7 +82,7 @@ def is_linear(form=None):
def isLagrange(fem):
return fem._short_name is 'CG'
return fem._short_name in ('CG', 'P', 'Q')
def isSimplical(cell):
......
......@@ -2,18 +2,22 @@
#include <string>
#include "dune/common/parallel/mpihelper.hh"
#include "dune/pdelab/stationary/linearproblem.hh"
#include "dune/common/parametertree.hh"
#include "dune/common/parametertreeparser.hh"
#include "dune/grid/yaspgrid.hh"
#include "dune/grid/io/file/vtk/subsamplingvtkwriter.hh"
#include "dune/pdelab/constraints/conforming.hh"
#include "dune/pdelab/gridoperator/gridoperator.hh"
#include "dune/common/parametertree.hh"
#include "dune/common/parametertreeparser.hh"
#include "dune/pdelab/function/callableadapter.hh"
#include "dune/pdelab/common/functionutilities.hh"
#include "dune/pdelab/stationary/linearproblem.hh"
#include "dune/pdelab/gridfunctionspace/vtk.hh"
#include "dune/pdelab/gridfunctionspace/gridfunctionadapter.hh"
#include "dune/testtools/gridconstruction.hh"
#include "dune/codegen/matrixfree.hh"
#include "dune/codegen/blockstructured/blockstructuredqkfem.hh"
#include "dune/codegen/blockstructured/preconditioner/preconditioner.hh"
#include "dune/codegen/vtkpredicate.hh"
#include "blockstructured_poisson_preconditioner_rOperator_file.hh"
#include "rCoarseOperator_file.hh"
......@@ -114,37 +118,61 @@ int main(int argc, char** argv){
InterpolationOperator interpolation(i_lop, gfs, cc, c_gfs, c_cc, decomp);
using I = decltype(interpolation);
using LS = IterativeMatrixFreeLocalSolver<GV, LOP, LocalDecompositionOperator<GFS, GFS>, NBLOCKS>;
LS ls(gv, gfs, lop);
using LOP_NN = NeumannNeumannOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_NN = Dune::PDELab::GridOperator<GFS, GFS, LOP_NN, MB, DF, RangeType, RangeType, CC, CC>;
LOP_NN lop_nn(gv, ls);
GOP_NN gop_nn(gfs, cc, gfs, cc, lop_nn, mb);
using LOP_SCHUR_APPLY = SchurApplyOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_SCHUR_APPLY = Dune::PDELab::GridOperator<GFS, GFS, LOP_SCHUR_APPLY, MB, DF, RangeType, RangeType, CC, CC>;
LOP_SCHUR_APPLY lop_schur_apply(gv, ls);
GOP_SCHUR_APPLY gop_schur_apply(gfs, cc, gfs, cc, lop_schur_apply, mb);
using LOP_SCHUR_RHS = SchurRHSOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_SCHUR_RHS = Dune::PDELab::GridOperator<GFS, GFS, LOP_SCHUR_RHS, MB, DF, RangeType, RangeType, CC, CC>;
LOP_SCHUR_RHS lop_schur_rhs(gv, ls);
GOP_SCHUR_RHS gop_schur_rhs(gfs, cc, gfs, cc, lop_schur_rhs, mb);
using LOP_SCHUR_BACK_TRAFO = SchurBackTrafoOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_SCHUR_BACK_TRAFO = Dune::PDELab::GridOperator<GFS, GFS, LOP_SCHUR_BACK_TRAFO, MB, DF, RangeType, RangeType, CC, CC>;
LOP_SCHUR_BACK_TRAFO lop_schur_back_trafo(gv, ls);
GOP_SCHUR_BACK_TRAFO gop_schur_back_trafo(gfs, cc, gfs, cc, lop_schur_back_trafo, mb);
Dune::Richardson<X, X> no_pre;
SchurOperatorMatrixFree<GOP_SCHUR_APPLY, GOP_SCHUR_RHS, GOP_SCHUR_BACK_TRAFO, Decomposition, X>
schurOperatorMatrixFree(gop_schur_apply, gop_schur_rhs, gop_schur_back_trafo, decomp);
NeumannNeumann<GOP_NN, Decomposition, X> pre_nn(gop_nn, decomp);
CoarseGridCorrection<I, R, Decomposition, X> pre_coarse(interpolation, restriction, decomp, c_gfs, c_cc, c_lop);
AdditiveTwoLevel<decltype(pre_nn), decltype(pre_coarse), X> pre_twoLevel(pre_nn, pre_coarse);
SchurComplement schurComplement(schurOperatorMatrixFree, pre_twoLevel, initree.get<bool>("only_preconditioning", true));
Dune::PDELab::solveMatrixFree(gop, x_r, schurComplement);
}
\ No newline at end of file
using LS = IterativeMatrixFreeLocalSolver<GV, LOP, LocalDecompositionOperator<GFS, GFS>, NBLOCKS>;
LS ls(gv, gfs, lop);
using LOP_NN = NeumannNeumannOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_NN = Dune::PDELab::GridOperator<GFS, GFS, LOP_NN, MB, DF, RangeType, RangeType, CC, CC>;
LOP_NN lop_nn(gv, ls);
GOP_NN gop_nn(gfs, cc, gfs, cc, lop_nn, mb);
using LOP_SCHUR_APPLY = SchurApplyOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_SCHUR_APPLY = Dune::PDELab::GridOperator<GFS, GFS, LOP_SCHUR_APPLY, MB, DF, RangeType, RangeType, CC, CC>;
LOP_SCHUR_APPLY lop_schur_apply(gv, ls);
GOP_SCHUR_APPLY gop_schur_apply(gfs, cc, gfs, cc, lop_schur_apply, mb);
using LOP_SCHUR_RHS = SchurRHSOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_SCHUR_RHS = Dune::PDELab::GridOperator<GFS, GFS, LOP_SCHUR_RHS, MB, DF, RangeType, RangeType, CC, CC>;
LOP_SCHUR_RHS lop_schur_rhs(gv, ls);
GOP_SCHUR_RHS gop_schur_rhs(gfs, cc, gfs, cc, lop_schur_rhs, mb);
using LOP_SCHUR_BACK_TRAFO = SchurBackTrafoOperator<GV, LS, RangeType, NBLOCKS>;
using GOP_SCHUR_BACK_TRAFO = Dune::PDELab::GridOperator<GFS, GFS, LOP_SCHUR_BACK_TRAFO, MB, DF, RangeType, RangeType, CC, CC>;
LOP_SCHUR_BACK_TRAFO lop_schur_back_trafo(gv, ls);
GOP_SCHUR_BACK_TRAFO gop_schur_back_trafo(gfs, cc, gfs, cc, lop_schur_back_trafo, mb);
Dune::Richardson<X, X> no_pre;
SchurOperatorMatrixFree<GOP_SCHUR_APPLY, GOP_SCHUR_RHS, GOP_SCHUR_BACK_TRAFO, Decomposition, X>
schurOperatorMatrixFree(gop_schur_apply, gop_schur_rhs, gop_schur_back_trafo, decomp);
NeumannNeumann<GOP_NN, Decomposition, X> pre_nn(gop_nn, decomp);
CoarseGridCorrection<I, R, Decomposition, X> pre_coarse(interpolation, restriction, decomp, c_gfs, c_cc, c_lop);
AdditiveTwoLevel<decltype(pre_nn), decltype(pre_coarse), X> pre_twoLevel(pre_nn, pre_coarse);
SchurComplement schurComplement(schurOperatorMatrixFree, pre_twoLevel, initree.get<bool>("only_preconditioning", true));
Dune::PDELab::solveMatrixFree(gop, x_r, schurComplement);
// Maybe calculate errors for test results...
auto exact_solution = func_0000;
using DiscreteGridFunction_ = Dune::PDELab::DiscreteGridFunction<GFS, X>;
DiscreteGridFunction_ discreteGridFunction_(gfs, x_r);
using DifferenceSquaredAdapter_ = Dune::PDELab::DifferenceSquaredAdapter<decltype(exact_solution), DiscreteGridFunction_>;
DifferenceSquaredAdapter_ dsa_(exact_solution, discreteGridFunction_);
Dune::FieldVector<RangeType, 1> l2error(0.0);
{
// L2 error squared of difference between numerical
// solution and the interpolation of exact solution
// for treepath ()
typename decltype(dsa_)::Traits::RangeType err(0.0);
Dune::PDELab::integrateGridFunction(dsa_, err, 10);
l2error += err;
if (gv.comm().rank() == 0){
std::cout << "L2 Error for treepath : " << err << std::endl;
}
}
if (isnan(l2error[0]) or abs(l2error[0])>1e-7)
return 1;
else
return 0;
}