Commit 486d0918 authored by Eike Mueller's avatar Eike Mueller
Browse files

Pass grid operator directly to (partially-) matrixfree preconditions

This avoids wrapping the operator outside the class and also makes the constructors more consistent.
parent 32e0fe1e
......@@ -32,29 +32,36 @@ namespace Dune {
* IterativeInverseBlockDiagonalLocalOperatorWrapper. The operator \f$A\f$
* is also applied in a matrix-free way.
*/
template <class OP, class IDOP>
template <class GO, class GOInvDiag>
class SeqMatrixfreeBlockJac
: public Dune::Preconditioner<typename OP::domain_type,
typename OP::range_type> {
: public Dune::Preconditioner<typename GO::Traits::Domain,
typename GO::Traits::Range> {
public:
/** \brief preconditioner category */
enum { category=SolverCategory::sequential };
/** \brief Domain type */
typedef typename OP::domain_type domain_type;
typedef typename GO::Traits::Domain domain_type;
/** \brief Range type */
typedef typename OP::range_type range_type;
typedef typename GO::Traits::Domain range_type;
/** \brief matrixfree operator */
typedef Dune::PDELab::OnTheFlyOperator<domain_type,range_type,GO> OP;
/** \brief matrixfree operator for inverse diagonal */
typedef Dune::PDELab::OnTheFlyOperator<domain_type,range_type,GOInvDiag> OPInvDiag;
/** \brief create a new instance
*
* \param[in] opA_ Operator \f$A\f$
* \param[in] opinvD_ Inverse diagonal operator \f$D^{-1}\f$
* \param[in] go_ Grid operator \f$A\f$
* \param[in] goinvdiag_ Inverse diagonal gridoperator \f$D^{-1}\f$
* \param[in] niter_ Number of block-Jacobi iterations
* \param[in] weight_ Scalar weight \f$w\f$
*/
SeqMatrixfreeBlockJac(const OP& opA_,
const IDOP& opinvD_,
SeqMatrixfreeBlockJac(const GO& go_,
const GOInvDiag& goinvdiag_,
const size_t niter_,
const double weight_) :
niter(niter_), weight(weight_), opA(opA_), opinvD(opinvD_) {}
niter(niter_), weight(weight_),
op(go_), opinvdiag(goinvdiag_) {}
/** \brief Prepare preconditioner (do nothing)
*
......@@ -81,10 +88,10 @@ namespace Dune {
if (i>0) r = d; // Reassign d
if ( (i==0) and (solution_is_zero) ) {
r *= weight;
opinvD.apply(r,v);
opinvdiag.apply(r,v);
} else {
opA.applyscaleadd(-1.0,v,r); // r <- d - A.v
opinvD.applyscaleadd(weight,r,v); // v <- v + w.D^{-1}.r
op.applyscaleadd(-1.0,v,r); // r <- d - A.v
opinvdiag.applyscaleadd(weight,r,v); // v <- v + w.D^{-1}.r
}
}
}
......@@ -112,9 +119,9 @@ namespace Dune {
/** \brief Jacobi weight \f$w\f$ */
const double weight;
/** \brief Operator \f$A\f$ */
const OP& opA;
const OP op;
/** \brief Inverse diagonal part \f$D^{-1}\f$ of operator \f$A\f$ */
const IDOP& opinvD;
const OPInvDiag opinvdiag;
};
/** \class SeqPartiallyMatrixfreeBlockJac
......@@ -133,32 +140,34 @@ namespace Dune {
* assembled, factorised and stored explicitly. The operator \f$A\f$ is
* still applied in a matrix-free way.
*/
template <class OP, class M>
template <class GO, class M>
class SeqPartiallyMatrixfreeBlockJac
: public Dune::Preconditioner<typename OP::domain_type,
typename OP::range_type> {
: public Dune::Preconditioner<typename GO::Traits::Domain,
typename GO::Traits::Range> {
public:
/** \brief preconditioner category */
enum { category=SolverCategory::sequential };
/** \brief Domain type */
typedef typename OP::domain_type domain_type;
typedef typename GO::Traits::Domain domain_type;
/** \brief Range type */
typedef typename OP::range_type range_type;
typedef typename GO::Traits::Range range_type;
/** \brief matrixfree operator */
typedef Dune::PDELab::OnTheFlyOperator<domain_type,range_type,GO> OP;
/** \brief create a new instance
*
* \param[in] opA_ Operator \f$A\f$
* \param[in] go_ Grid operator \f$A\f$
* \param[in] D_ Diagonal matrix \f$D\f$
* \param[in] niter_ Number of iterations
* \param[in] weight_ Jacobi weight \f$w\f$
* \param[in] symmetric_blockdiag_ Are the block-diagonal entries
* symmetric?
*/
SeqPartiallyMatrixfreeBlockJac(const OP& opA_,
SeqPartiallyMatrixfreeBlockJac(const GO& go_,
const M& D_,
const size_t niter_,
const double weight_,
const bool symmetric_blockdiag_) :
opA(opA_), D(D_), niter(niter_), weight(weight_),
op(go_), D(D_), niter(niter_), weight(weight_),
symmetric_blockdiag(symmetric_blockdiag_) {
// check that matrix is square
assert (D.N() == D.M());
......@@ -201,7 +210,7 @@ namespace Dune {
PDELab::Backend::native(v));
} else {
// r <- d - A.v
opA.applyscaleadd(-1.0,v,r);
op.applyscaleadd(-1.0,v,r);
// v <- v + w*D^{-1}.r = v + w*D^{-1}.(d-A.v)
invD->usmv(weight,
PDELab::Backend::native(r),
......@@ -229,7 +238,7 @@ namespace Dune {
private:
/** \brief Operator \f$A\f$ */
const OP& opA;
const OP op;
/** \brief Matrix for diagonal part \f$D\f$ of operator \f$A\f$ */
const M& D;
/** \brief Number of iterations */
......
......@@ -128,18 +128,9 @@ namespace Dune {
typename GO::Traits::Range,
GO> Op;
Op op(cc,go);
// Construct sequential operators for sequential preconditioner
typedef OnTheFlyOperator<typename GO::Traits::Domain,
typename GO::Traits::Range,
GO> SEQ_Op;
SEQ_Op seq_op(go);
typedef OnTheFlyOperator<typename GOInvDiag::Traits::Domain,
typename GOInvDiag::Traits::Range,
GOInvDiag> SEQ_OpInvDiag;
SEQ_OpInvDiag seq_opinvdiag(goinvdiag);
// Construct sequential preconditioner
typedef Dune::SeqMatrixfreeBlockJac<SEQ_Op,SEQ_OpInvDiag> SEQ_Prec;
SEQ_Prec seq_prec(seq_op,seq_opinvdiag,steps,1.0);
typedef Dune::SeqMatrixfreeBlockJac<GO,GOInvDiag> SEQ_Prec;
SEQ_Prec seq_prec(go,goinvdiag,steps,1.0);
// Construct parallel preconditioner
typedef MatrixFreeOVLPPreconditionerWrapper<CC,GFS,SEQ_Prec> Prec;
Prec prec(gfs,seq_prec,cc);
......@@ -383,10 +374,6 @@ namespace Dune {
GO> Op;
Op pop(cc,go);
// Construct sequential operator for sequential preconditioner
typedef OnTheFlyOperator<typename GO::Traits::Domain,
typename GO::Traits::Range,
GO> SEQ_Op;
SEQ_Op op(go);
typedef OVLPScalarProduct<GFS,V> PSP;
PSP psp(*this);
Dune::Timer watch;
......@@ -397,8 +384,8 @@ namespace Dune {
typename GODiag::Traits::Jacobian djac(godiag);
godiag.jacobian(u,djac);
// sequential preconditioner type
typedef Preconditioner<SEQ_Op,typename GODiag::Traits::Jacobian> SeqPrec;
SeqPrec seq_prec(op,djac,steps,1.0,symmetric_blockdiag);
typedef Preconditioner<GO,typename GODiag::Traits::Jacobian> SeqPrec;
SeqPrec seq_prec(go,djac,steps,1.0,symmetric_blockdiag);
double prec_setup_time = watch.elapsed();
if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "=== preconditioner setup (max) " << prec_setup_time << " s" << std::endl;
typedef MatrixFreeOVLPPreconditionerWrapper<CC,GFS,SeqPrec> WPREC;
......
......@@ -439,18 +439,10 @@ namespace Dune {
void apply(const GO& go, V& z, W& r,
typename V::ElementType reduction) {
// Construct DG preconditioner for preconditioner
typedef OnTheFlyOperator<typename GOPrec::Traits::Domain,
typename GOPrec::Traits::Range,
GOPrec> SeqPrecOp;
typedef OnTheFlyOperator<typename GOPrecInvDiag::Traits::Domain,
typename GOPrecInvDiag::Traits::Range,
GOPrecInvDiag> SeqPrecOpInvDiag;
typedef DGPrec<SeqPrecOp,SeqPrecOpInvDiag> DGPrecType;
typedef DGPrec<GOPrec,GOPrecInvDiag> DGPrecType;
Dune::Timer watch;
watch.reset();
SeqPrecOp seq_precop(goprec);
SeqPrecOpInvDiag seq_precopinvdiag(goprecinvdiag);
DGPrecType dgprec(seq_precop,seq_precopinvdiag,1,rho);
DGPrecType dgprec(goprec,goprecinvdiag,1,rho);
double dgsmoother_setup_time = watch.elapsed();
if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "=== DG smoother setup time " << dgsmoother_setup_time << " s" << std::endl;
Base::apply(go,dgprec,z,r,reduction);
......@@ -529,15 +521,11 @@ namespace Dune {
void apply(const GO& go, V& z, W& r,
typename V::ElementType reduction) {
// Construct DG preconditioner
typedef OnTheFlyOperator<typename GOPrec::Traits::Domain,
typename GOPrec::Traits::Range,
GOPrec> SeqPrecOp;
typedef typename GOPrecDiag::Traits::Jacobian Matrix;
typedef DGPrec<SeqPrecOp,Matrix> DGPrecType;
typedef DGPrec<GO,Matrix> DGPrecType;
Dune::Timer watch;
watch.reset();
SeqPrecOp seq_precop(goprec);
Matrix djac(goprecdiag);
typedef typename GOPrecDiag::Traits::Domain U;
U u(gfs,0.0);
......@@ -548,7 +536,7 @@ namespace Dune {
const unsigned int idx = col.first[col.first.size()-1];
PDELab::Backend::native(djac)[idx][idx] = 0.0;
}
DGPrecType dgprec(seq_precop,djac,1,rho,symmetric_blockdiag);
DGPrecType dgprec(goprec,djac,1,rho,symmetric_blockdiag);
double dgsmoother_setup_time = watch.elapsed();
if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "=== DG smoother setup time " << dgsmoother_setup_time << " s" << std::endl;
Base::apply(go,dgprec,z,r,reduction);
......
......@@ -869,12 +869,8 @@ void test_Jacobi(TestHarness& test_harness,
const double omega=0.9;
// *** Matrix-free application ***
typedef Dune::PDELab::OnTheFlyOperator<U,V,GO> SeqOP;
SeqOP seq_op(go);
typedef Dune::PDELab::OnTheFlyOperator<U,V,IDGO> SeqIDOP;
SeqIDOP seq_idop(idgo);
typedef typename Dune::SeqMatrixfreeBlockJac<SeqOP,SeqIDOP> SeqPrecMF;
SeqPrecMF seq_prec_mf(seq_op,seq_idop,niter,omega);
typedef typename Dune::SeqMatrixfreeBlockJac<GO,IDGO> SeqPrecMF;
SeqPrecMF seq_prec_mf(go,idgo,niter,omega);
typedef typename Dune::PDELab::MatrixFreeOVLPPreconditionerWrapper<CC, GFS, SeqPrecMF> PrecMF;
PrecMF prec_mf(gfs,seq_prec_mf,cc);
......@@ -890,8 +886,8 @@ void test_Jacobi(TestHarness& test_harness,
typename DGO::Traits::Jacobian djac(dgo);
typedef typename Dune::PDELab::Backend::Native<typename DGO::Traits::Jacobian> DMat;
dgo.jacobian(u,djac);
typedef typename Dune::SeqPartiallyMatrixfreeBlockJac<SeqOP,DMat> SeqPrecPMF;
SeqPrecPMF seq_prec_pmf(seq_op,native(djac),niter,omega,true);
typedef typename Dune::SeqPartiallyMatrixfreeBlockJac<GO,DMat> SeqPrecPMF;
SeqPrecPMF seq_prec_pmf(go,native(djac),niter,omega,true);
typedef typename Dune::PDELab::MatrixFreeOVLPPreconditionerWrapper<CC, GFS, SeqPrecPMF> PrecPMF;
PrecPMF prec_pmf(gfs,seq_prec_pmf,cc);
// Apply preconditioner
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment