MimeticIPEvaluator.hpp

//===========================================================================
//
// File: MimeticIPEvaluator.hpp
//
// Created: Thu Jun 25 13:09:23 2009
//
// Author(s): B�rd Skaflestad     <bard.skaflestad@sintef.no>
//            Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================

/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010 Statoil ASA.

  This file is part of The Open Reservoir Simulator Project (OpenRS).

  OpenRS is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OpenRS is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPENRS_MIMETICIPEVALUATOR_HEADER
#define OPENRS_MIMETICIPEVALUATOR_HEADER

#include <algorithm>
#include <vector>

#include <array>

#include <opm/common/ErrorMacros.hpp>
#include <opm/grid/utility/SparseTable.hpp>

#include <opm/porsol/common/fortran.hpp>
#include <opm/porsol/common/blas_lapack.hpp>
#include <opm/porsol/common/Matrix.hpp>

namespace Opm {
    template <class GridInterface, class RockInterface>
    class MimeticIPEvaluator
    {
    public:
        static constexpr int dim = GridInterface::Dimension;
        typedef typename GridInterface::CellIterator CellIter;
        typedef typename CellIter::Scalar Scalar;


        MimeticIPEvaluator()
            : max_nf_(-1),
              fa_    (  ),
              t1_    (  ),
              t2_    (  ),
              Binv_  (  )
        {}


        explicit MimeticIPEvaluator(const int max_nf)
            : max_nf_(max_nf         ),
              fa_    (max_nf * max_nf),
              t1_    (max_nf * dim   ),
              t2_    (max_nf * dim   ),
              Binv_  (               ),
              gflux_ (               )
        {}


        void init(const int max_nf)
        {
            max_nf_ = max_nf;
            std::vector<double>(max_nf * max_nf).swap(fa_);
            std::vector<double>(max_nf * dim   ).swap(t1_);
            std::vector<double>(max_nf * dim   ).swap(t2_);
        }


        template<class Vector>
        void reserveMatrices(const Vector& sz)
        {
            typedef typename Vector::value_type vt;

            Vector sz2(sz.size());

            std::ranges::transform(sz, sz2.begin(),
                                   [](const vt& input) { return input*input; });

            Binv_ .allocate(sz2.begin(), sz2.end());
            gflux_.allocate(sz .begin(), sz .end());
        }


        void buildStaticContrib(const CellIter& c,
                                const RockInterface& r,
                                const typename CellIter::Vector& grav,
                                const int nf)
        {
            typedef typename CellIter::FaceIterator FI;
            typedef typename CellIter::Vector       CV;
            typedef typename FI      ::Vector       FV;

            const int ci = c->index();

            static_assert (FV::dimension == int(dim), "");
            assert (int(t1_.size()) >= nf * dim);
            assert (int(t2_.size()) >= nf * dim);
            assert (int(fa_.size()) >= nf * nf);

            SharedFortranMatrix T1  (nf, dim, &t1_      [0]);
            SharedFortranMatrix T2  (nf, dim, &t2_      [0]);
            SharedFortranMatrix fa  (nf, nf , &fa_      [0]);
            SharedFortranMatrix Binv(nf, nf , &Binv_[ci][0]);

            // Clear matrices of any residual data.
            zero(Binv);  zero(T1);  zero(T2);  zero(fa);

            typename RockInterface::PermTensor K  = r.permeability(ci);
            const    CV          Kg = prod(K, grav);

            // Setup:    Binv  <- I, T1 <- N, T2 <- C
            // Complete: gflux <- N*K*g
            const CV cc = c->centroid();
            int i = 0;
            for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
                Binv(i,i) = Scalar(1.0);
                fa(i,i)   = f->area();

                FV fc = f->centroid();  fc -= cc;  fc *= fa(i,i);
                FV fn = f->normal  ();             fn *= fa(i,i);

                gflux_[ci][i] = fn * Kg;

                for (int j = 0; j < dim; ++j) {
                    T1(i,j) = fn[j];
                    T2(i,j) = fc[j];
                }
            }
            assert (i == nf);

            // T2 <- orth(T2)
            if (orthogonalizeColumns(T2) != 0) {
                assert (false);
            }

            // Binv <- Binv - T2*T2' == I - Q*Q'
            symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);

            // Binv <- diag(A) * Binv * diag(A)
            symmetricUpdate(fa, Binv);

            // T2 <- N*K
            matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);

            // Binv <- (T2*N' + Binv) / vol(c)
            //      == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
            //
            // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
            //
            Scalar t = Scalar(6.0) * trace(K) / dim;
            matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
                         t           / c->volume(), Binv  );
        }


        template<class FluidInterface, class Sat>
        void computeDynamicParams(const CellIter&         c,
                                  const FluidInterface&   fl,
                                  const std::vector<Sat>& s)
        {
            const int ci = c->index();

            std::array<Scalar, FluidInterface::NumberOfPhases> mob ;
            std::array<Scalar, FluidInterface::NumberOfPhases> rho ;
            fl.phaseMobilities(ci, s[ci], mob);
            fl.phaseDensities (ci, rho);

            totmob_   = std::accumulate   (mob.begin(), mob.end(), Scalar(0.0));
            mob_dens_ = std::inner_product(rho.begin(), rho.end(), mob.begin(),
                                           Scalar(0.0));
        }


        template<template<typename> class SP>
        void getInverseMatrix(const CellIter&                        c,
                              FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
        {
            getInverseMatrix(c, totmob_, Binv);
        }

        template<template<typename> class SP>
        void getInverseMatrix(const CellIter&                        c,
                              const Scalar                           totmob,
                              FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
        {
            const int ci = c->index();
            std::transform(Binv_[ci].begin(), Binv_[ci].end(), Binv.data(),
                           [totmob](const Scalar& input) { return input*totmob; });
        }

        template<class PermTensor, template<typename> class SP>
        void evaluate(const CellIter&                        c,
                      const PermTensor&                      K,
                      FullMatrix<Scalar,SP,FortranOrdering>& Binv)
        {
            typedef typename CellIter::FaceIterator FI;
            typedef typename CellIter::Vector       CV;
            typedef typename FI      ::Vector       FV;

            const int nf = Binv.numRows();

            static_assert(FV::dimension == int(dim), "");
            assert(Binv.numRows()  <= max_nf_);
            assert(Binv.numRows()  == Binv.numCols());
            assert(int(t1_.size()) >= nf * dim);
            assert(int(t2_.size()) >= nf * dim);
            assert(int(fa_.size()) >= nf * nf);

            SharedFortranMatrix T1(nf, dim, &t1_[0]);
            SharedFortranMatrix T2(nf, dim, &t2_[0]);
            SharedFortranMatrix fa(nf, nf , &fa_[0]);

            // Clear matrices of any residual data.
            zero(Binv);  zero(T1);  zero(T2);  zero(fa);

            // Setup: Binv <- I, T1 <- N, T2 <- C
            const CV cc = c->centroid();
            int i = 0;
            for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
                Binv(i,i) = Scalar(1.0);
                fa(i,i)   = f->area();

                FV fc = f->centroid();  fc -= cc;  fc *= fa(i,i);
                FV fn = f->normal  ();             fn *= fa(i,i);

                for (int j = 0; j < dim; ++j) {
                    T1(i,j) = fn[j];
                    T2(i,j) = fc[j];
                }
            }
            assert(i == nf);

            // T2 <- orth(T2)
            if (orthogonalizeColumns(T2) != 0) {
                assert (false);
            }

            // Binv <- Binv - T2*T2' == I - Q*Q'
            symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);

            // Binv <- diag(A) * Binv * diag(A)
            symmetricUpdate(fa, Binv);

            // T2 <- N*K
            matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);

            // Binv <- (T2*N' + Binv) / vol(c)
            //      == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
            //
            // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
            //
            Scalar t = Scalar(6.0) * trace(K) / dim;
            matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
                         t           / c->volume(), Binv  );
        }


        template<class Vector>
        void gravityTerm(const CellIter& c,
                         const typename CellIter::Vector& grav,
                         const Scalar    omega,
                         Vector&         gterm) const
        {
            typedef typename CellIter::FaceIterator FI;
            typedef typename CellIter::Vector Point;

            assert (gterm.size() <= max_nf_);

            const Point cc = c->centroid();
            int i = 0;
            for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
                Point fc = f->centroid();
                fc -= cc;
                gterm[i] = omega * (fc * grav);
            }
        }

        template<class Vector>
        void gravityTerm(const CellIter& c,
                         const typename CellIter::Vector&    grav,
                         Vector&         gterm) const
        {
            gravityTerm(c, grav, mob_dens_ / totmob_, gterm);
        }

        template<class FluidInterface, class Sat, class Vector>
        void gravityTerm(const CellIter&         c,
                         const FluidInterface&   fl,
                         const std::vector<Sat>& s,
                         const typename CellIter::Vector& grav,
                         Vector&                 gterm) const
        {
            const int ci = c->index();

            std::array<Scalar, FluidInterface::NumberOfPhases> mob;
            std::array<Scalar, FluidInterface::NumberOfPhases> rho;
            fl.phaseMobilities(ci, s[ci], mob);
            fl.phaseDensities (ci, rho);

            Scalar totmob = std::accumulate   (mob.begin(), mob.end(), Scalar(0.0));
            Scalar omega  = std::inner_product(rho.begin(), rho.end(), mob.begin(),
                                               Scalar(0.0)) / totmob;

            gravityTerm(c, grav, omega, gterm);
        }


        template<class Vector>
        void gravityFlux(const CellIter& c,
                         Vector&         gflux) const
        {
            Scalar mob_dens = mob_dens_;
            std::transform(gflux_[c->index()].begin(), gflux_[c->index()].end(),
                           gflux.begin(),
                           [mob_dens](const Scalar& input) { return input*mob_dens; });
        }

    private:
        int                 max_nf_      ;
        Scalar              totmob_      ;
        Scalar              mob_dens_    ;
        std::vector<Scalar> fa_, t1_, t2_;
        Opm::SparseTable<Scalar> Binv_        ;
        Opm::SparseTable<Scalar> gflux_       ;
    };
} // namespace Opm

#endif // OPENRS_MIMETICIPEVALUATOR_HEADER