Operators.hpp

Go to the documentation of this file.

#ifndef OPM_MIXED_OPERATORS_HEADER_INCLUDED
#define OPM_MIXED_OPERATORS_HEADER_INCLUDED
 
 
namespace Opm
{
 
template<class M, class V, class C>
class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,V,V>
{
public:
    MixedGhostLastMatrixAdapter (const M& A, const C& comm) : A_( A ), comm_(comm) {}
 
    // y = A * x
    void apply( const V& x, V& y ) const override
    {
        A_.mv(x,y);
        ghostLast_project(y);
    }
 
    // y += \alpha * A * x
    void applyscaleadd (double alpha, const V& x, V& y) const override
    {
        A_.usmv(alpha,x,y);
        ghostLast_project(y);
    }
 
    // accessor to matix object
    const M& getmat() const override { return A_; }
 
    // solver category
    Dune::SolverCategory::Category category() const override
    {
        return Dune::SolverCategory::overlapping;
    }
 
private:
    // extract block size
    static constexpr auto block_size = V::block_type::dimension;
 
    void ghostLast_project( V& y ) const
    {
        double *yy = &y[0][0];
        int n = block_size*A_.nrows();
        int N = block_size*y.N();
        for (int i=n;i<N;i++) yy[i] = 0;
    }
 
    const M& A_;
    const C& comm_;
};
 
 
template<class M, class V, class C>
class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,V,V>
{
public:
    using field_type = typename V::field_type;
    using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
 
    // extract block size
    static constexpr auto block_size = V::block_type::dimension;
 
    WellModelMixedGhostLastMatrixAdapter (const M& A,
                                     const LinearOperatorExtra<V, V>& wellOper,
                                     const C& comm
                                    )
        : A_( A ), wellOper_( wellOper ), comm_ ( comm )
    {}
 
    // y = A * x
    void apply( const V& x, V& y ) const override
    {
        A_.mv(x,y);
        wellOper_.apply(x, y);
        //comm_.project(y);
        ghostLast_project(y);
    }
 
    // y += \alpha * A * x
    void applyscaleadd (double alpha, const V& x, V& y) const override
    {
        A_.usmv(alpha,x,y);
        wellOper_.applyscaleadd(alpha, x, y);
        //comm_.project(y);
        ghostLast_project(y);
    }
 
    // accessor to matix object
    const M& getmat() const override { return A_; }
 
    void addWellPressureEquations(PressureMatrix& jacobian,
                                  const V& weights,
                                  const bool use_well_weights) const
    {
        OPM_TIMEBLOCK(addWellPressureEquations);
        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
    }
 
    void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
    {
        OPM_TIMEBLOCK(addWellPressureEquationsStruct);
        wellOper_.addWellPressureEquationsStruct(jacobian);
    }
 
    int getNumberOfExtraEquations() const
    {
        return wellOper_.getNumberOfExtraEquations();
    }
 
    // solver category
    Dune::SolverCategory::Category category() const override
    {
        return Dune::SolverCategory::overlapping;
    }
 
 
protected:
 
    const M& A_ ;
    const LinearOperatorExtra<V, V>& wellOper_;
    const C& comm_ ;
 
private:
 
    void ghostLast_project( V& y ) const
    {
        double *yy = &y[0][0];
        int n = block_size*A_.nrows();
        int N = block_size*y.N();
        for (int i=n;i<N;i++) yy[i] = 0;
    }
 
};
 
} // namespace Opm
#endif // OPM_MIXED_OPERATORS_HEADER_INCLUDED