1 // -*- tab-width: 2; indent-tabs-mode: nil -*- 2 // vi: set et ts=2 sw=2 sts=2: 3 4 #ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_BACKENDS_HH 5 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_BACKENDS_HH 6 7 #include<dune/grid/common/rangegenerators.hh> 8 #include<dune/istl/preconditioner.hh> 9 #include<dune/istl/operators.hh> 10 #include<dune/istl/solvers.hh> 11 12 #include<dune/pdelab/backend/istl/matrixfree/gridoperatorpreconditioner.hh> 13 #include<dune/pdelab/gridoperator/fastdg.hh> 14 #include<dune/pdelab/backend/istl/matrixfree/blocksorpreconditioner.hh> 15 #include<dune/pdelab/backend/istl/matrixfree/checklopinterface.hh> 16 17 18 namespace Dune { 19 namespace PDELab { 20 21 namespace impl{ 22 23 // Can be used to check if a local operator is a 24 // BlockSORPreconditionerLocalOperator 25 template <typename> 26 struct isBlockSORPreconditionerLocalOperator : std::false_type {}; 27 28 template <typename JacobianLOP, typename BlockOffDiagonalLOP, typename GridFunctionSpace> 29 struct isBlockSORPreconditionerLocalOperator<BlockSORPreconditionerLocalOperator<JacobianLOP, 30 BlockOffDiagonalLOP, 31 GridFunctionSpace>> 32 : std::true_type {}; 33 34 // Can be used to check if a grid operator is a FastDGGridOperator 35 template <typename> 36 struct isFastDGGridOperator : std::false_type {}; 37 38 template<typename GFSU, typename GFSV, typename LOP, 39 typename MB, typename DF, typename RF, typename JF, 40 typename CU, typename CV> 41 struct isFastDGGridOperator<FastDGGridOperator<GFSU, GFSV, LOP, MB, DF, RF, JF, CU, CV >> 42 : std::true_type {}; 43 44 } 45 46 /**\brief Sequential matrix-free solver backend 47 * 48 * This can be used with a combination of {CGSolver, BiCGSTABSolver, 49 * MINRESSolver} as a solver and grid operator build from one of the 50 * following local operators: 51 * 52 * - AssembledBlockJacobiPreconditionerLocalOperator for partial matrix-free Jacobi 53 * - IterativeBlockJacobiPreconditionerLocalOperator for fully matrix-free Jacobi 54 * - BlockSORPreconditionerLocalOperator for matrix-free SOR 55 * 56 * Note: If you use BlockSORPreconditionerLocalOperator you need to use 57 * FastDGGridOperator! 58 59 * \tparam GO Grid operator implementing the matrix-free operator application 60 * \tparam PrecGO Grid operator implementing matrix-free preconditioning 61 */ 62 template<class GO, class PrecGO, 63 template<class> class Solver> 64 class ISTLBackend_SEQ_MatrixFree_Base 65 : public Dune::PDELab::SequentialNorm 66 , public Dune::PDELab::LinearResultStorage 67 { 68 using V = typename GO::Traits::Domain; 69 using W = typename GO::Traits::Range; 70 using Operator = Dune::PDELab::OnTheFlyOperator<V,W,GO>; 71 using SeqPrec = GridOperatorPreconditioner<PrecGO>; 72 73 public : ISTLBackend_SEQ_MatrixFree_Base(const GO & go,const PrecGO & precgo,unsigned maxiter=5000,int verbose=1)74 ISTLBackend_SEQ_MatrixFree_Base(const GO& go, const PrecGO& precgo, 75 unsigned maxiter=5000, int verbose=1) 76 : opa_(go), seqprec_(precgo) 77 , maxiter_(maxiter), verbose_(verbose) 78 { 79 // First lets brake that down: The case we want to avoid is having SOR 80 // with regular grid operator. We check for SOR and not fast grid 81 // operator and assert that this is not the case. 82 // 83 // For more information why this is not working have a look at the 84 // documentation of the class BlockSORPreconditionerLocalOperator 85 static_assert(not(impl::isBlockSORPreconditionerLocalOperator< 86 typename PrecGO::Traits::LocalAssembler::LocalOperator>::value 87 and not impl::isFastDGGridOperator<PrecGO>::value), 88 "If you use the BlockSORPreconditioner you need to use FastDGGridOperator!"); 89 90 91 // We need to assert that the local operator uses the new interface and not the old one 92 static_assert(!decltype(Dune::PDELab::impl::hasOldLOPInterface(go.localAssembler().localOperator()))::value, 93 "\n" 94 "In order to use the matrix-free solvers your grid operator must follow the new\n" 95 "local operator interface for nonlinear jacobian applys, where the current\n" 96 "solution and the linearization point can have different types. This is best\n" 97 "explained by an example.\n\n" 98 "old: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const X& z, const LFSV& lfsv, Y& y) const {}\n" 99 "new: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const {}\n\n" 100 "Note that in the new interface the type of z is Z and doesn't have to be X.\n" 101 "You need to adjust all the nonlinear jacobian apply methods in your local operator\n" 102 "in a similar way."); 103 static_assert(!decltype(Dune::PDELab::impl::hasOldLOPInterface(precgo.localAssembler().localOperator()))::value, 104 "\n" 105 "In order to use the matrix-free solvers your grid operator must follow the new\n" 106 "local operator interface for nonlinear jacobian applys, where the current\n" 107 "solution and the linearization point can have different types. This is best\n" 108 "explained by an example.\n\n" 109 "old: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const X& z, const LFSV& lfsv, Y& y) const {}\n" 110 "new: void jacobian_apply_volume (const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const {}\n\n" 111 "Note that in the new interface the type of z is Z and doesn't have to be X.\n" 112 "You need to adjust all the nonlinear jacobian apply methods in your local operator\n" 113 "in a similar way."); 114 } 115 apply(V & z,W & r,typename V::ElementType reduction)116 void apply(V& z, W& r, typename V::ElementType reduction) 117 { 118 Solver<V> solver(opa_,seqprec_,reduction,maxiter_,verbose_); 119 Dune::InverseOperatorResult stat; 120 solver.apply(z, r, stat); 121 res.converged = stat.converged; 122 res.iterations = stat.iterations; 123 res.elapsed = stat.elapsed; 124 res.reduction = stat.reduction; 125 res.conv_rate = stat.conv_rate; 126 } 127 128 //! Set position of jacobian. 129 //! Must be called before apply(). setLinearizationPoint(const V & u)130 void setLinearizationPoint(const V& u) 131 { 132 opa_.setLinearizationPoint(u); 133 seqprec_.setLinearizationPoint(u); 134 } 135 136 constexpr static bool isMatrixFree {true}; 137 138 private : 139 Operator opa_; 140 SeqPrec seqprec_; 141 unsigned maxiter_; 142 int verbose_; 143 }; // end class ISTLBackend_SEQ_MatrixFree_Base 144 145 146 // A matrix based SOR backend for comparing the matrix-free version 147 class ISTLBackend_SEQ_BCGS_SOR 148 : public ISTLBackend_SEQ_Base<Dune::SeqSOR, Dune::BiCGSTABSolver> 149 { 150 public: 151 /*! \brief make a linear solver object 152 153 \param[in] maxiter_ maximum number of iterations to do 154 \param[in] verbose_ print messages if true 155 */ ISTLBackend_SEQ_BCGS_SOR(unsigned maxiter_=5000,int verbose_=1)156 explicit ISTLBackend_SEQ_BCGS_SOR (unsigned maxiter_=5000, int verbose_=1) 157 : ISTLBackend_SEQ_Base<Dune::SeqSOR, Dune::BiCGSTABSolver>(maxiter_, verbose_) 158 {} 159 }; 160 161 } // namespace PDELab 162 } // namespace Dune 163 #endif 164