ScalarProducts.hpp

Go to the documentation of this file.

#ifndef OPM_SCALAR_PRODUCTS_HEADER_INCLUDED
#define OPM_SCALAR_PRODUCTS_HEADER_INCLUDED
 
#include <opm/simulators/linalg/mixed/dot.h>
 
namespace Dune
{
 
template<class Vector, class Comm>
class GhostLastScalarProduct : public ScalarProduct<Vector>
{
    public:
 
    GhostLastScalarProduct (std::shared_ptr<const Comm> com, SolverCategory::Category cat)
        : _communication(com), _category(cat)
    {
        count_ = getLocalCount();        // number or local cells
        int verify = verifyLocalCount(); // redundant check on numbef of local cells
        if (count_ != verify) OPM_THROW(std::runtime_error, "Inconsistent local node count!!\n");
    }
 
    GhostLastScalarProduct (const Comm& com, SolverCategory::Category cat)
        : GhostLastScalarProduct(stackobject_to_shared_ptr(com), cat)
    {}
 
    double dot (const Vector& vx, const Vector& vy) const override
    {
 
        // access underlying data
        double const *x = &vx[0][0];
        double const *y = &vy[0][0];
 
        // total array length
        int NN = block_size*count_;
 
        auto cc = _communication->communicator();
        return cc.sum(vec_dot(x,y,NN));
    }
 
    double norm (const Vector& vx) const override
    {
        return sqrt(dot(vx,vx));
    }
 
    virtual SolverCategory::Category category() const override
    {
        return _category;
    }
 
    private:
 
    static constexpr auto block_size = Vector::block_type::dimension;
 
    std::shared_ptr<const Comm> _communication;
    SolverCategory::Category _category;
    int count_;
 
    int getLocalCount() const
    {
        int count = 0;
        // Loop over index set
        auto indexSet = _communication->indexSet();
        for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) {
            if (idx->local().attribute()==1) count++; // count non-local indices
        }
        return count;
    }
 
    int verifyLocalCount() const
    {
        auto indexSet = _communication->indexSet();
 
        size_t is = 0;
        // Loop over index set
        for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) {
            //Only take "owner" indices
            if (idx->local().attribute()==1) {
                //get local index
                auto loc = idx->local().local();
                // if loc is higher than "old interior size", update it
                if (loc > is) {
                    is = loc;
                }
            }
        }
        return is + 1; //size is plus 1 since we start at 0
    }
 
};
 
 
 
template<class Vector>
class SeqOptmizedProduct : public Dune::SeqScalarProduct<Vector>
{
public:
 
    double dot(const Vector& vx, const Vector& vy) const override
    {
        // access underlying data
        double const *x = &vx[0][0];
        double const *y = &vy[0][0];
 
        // total array length
        int NN = block_size*vx.N();
 
        return vec_dot(x,y,NN);
    }
 
    double norm(const Vector& vx) const override {
        return std::sqrt(this->dot(vx, vx));
    }
 
    private:
 
    // extract block size
    static constexpr auto block_size = Vector::block_type::dimension;
 
 
};
 
}
 
#endif //OPM_SCALAR_PRODUCTS_HEADER_INCLUDED