MimeticIPAnisoRelpermEvaluator.hpp

//===========================================================================
//
// File: MimeticIPAnisoRelpermEvaluator.hpp
//
// Created: Mon Oct 19 10:22:22 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            B�rd Skaflestad     <bard.skaflestad@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_MIMETICIPANISORELPERMEVALUATOR_HEADER
#define OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER


#include <algorithm>
#include <vector>

#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 MimeticIPAnisoRelpermEvaluator
    {
    public:
        static constexpr int dim = GridInterface::Dimension;
        typedef typename GridInterface::CellIterator CellIter;
        typedef typename CellIter::Scalar Scalar;


        MimeticIPAnisoRelpermEvaluator()
            : max_nf_(-1),
              prock_(0)
        {}


        explicit MimeticIPAnisoRelpermEvaluator(const int max_nf)
            : max_nf_       (max_nf         ),
              fa_           (max_nf * max_nf),
              t1_           (max_nf * dim   ),
              t2_           (max_nf * dim   ),
              second_term_  (               ),
              n_            (               ),
              Kg_           (               ),
              prock_        (        0      )
        {}


        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; });

            second_term_.allocate(sz2.begin(), sz2.end());

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

            n_.allocate(sz2.begin(), sz2.end());

            std::ranges::fill(sz2, vt(dim));
            Kg_.allocate(sz2.begin(), sz2.end());
        }


        void buildStaticContrib(const CellIter& c,
                                const RockInterface& r,
                                const typename CellIter::Vector& grav,
                                const int nf)
        {
            // Binv = (N*lambda*K*N'   +   t*diag(A)*(I - Q*Q')*diag(A))/vol
            //         ^                     ^^^^^^^^^^^^^^^^^^^^^^^^^^
            //         precompute: n_        precompute: second_term_
            // t = 6/dim * trace(lambda*K)

            typedef typename CellIter::FaceIterator FI;
            typedef typename CellIter::Vector       CV;
            typedef typename FI      ::Vector       FV;

            // Now we need to remember the rocks, since we will need
            // the permeability for dynamic assembly.
            prock_ = &r;

            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 T2  (nf, dim, &t2_      [0]);
            SharedFortranMatrix fa  (nf, nf , &fa_      [0]);
            SharedFortranMatrix second_term(nf, nf, &second_term_[ci][0]);
            SharedFortranMatrix n(nf, dim, &n_[ci][0]);

            // Clear matrices of any residual data.
            zero(second_term);  zero(n);  zero(T2);  zero(fa);

            // Setup: second_term <- I, n <- N, T2 <- C
            const CV cc = c->centroid();
            int i = 0;
            for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
                second_term(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) {
                    n (i,j) = fn[j];
                    T2(i,j) = fc[j];
                }
            }
            assert (i == nf);

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

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

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

            // Gravity term: Kg_ = K * grav
            vecMulAdd_N(Scalar(1.0), r.permeability(ci), &grav[0],
                        Scalar(0.0), &Kg_[ci][0]);
        }


        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, dim * dim> lambda_t;
            std::array<Scalar, dim * dim> pmob_data;

            SharedFortranMatrix pmob(dim, dim, &pmob_data[0]);
            SharedFortranMatrix Kg  (dim, 1  , &Kg_[ci][0]);

            std::array<Scalar, FluidInterface::NumberOfPhases> rho;
            fl.phaseDensities(ci, rho);

            std::ranges::fill(dyn_Kg_, Scalar(0.0));
            std::ranges::fill(lambda_t, 0.0);

            for (int phase = 0; phase < FluidInterface::NumberOfPhases; ++phase) {
                fl.phaseMobility(phase, ci, s[ci], pmob);

                // dyn_Kg_ += (\rho_phase \lambda_phase) Kg
                vecMulAdd_N(rho[phase], pmob, Kg.data(), Scalar(1.0), dyn_Kg_.data());

                // \lambda_t += \lambda_phase
                std::ranges::transform(lambda_t, pmob_data, lambda_t.begin(), std::plus<Scalar>());
            }

            // lambdaK_ = (\sum_i \lambda_i) K
            SharedFortranMatrix lambdaT(dim, dim, lambda_t.data());
            SharedFortranMatrix lambdaK(dim, dim, lambdaK_.data());
            prod(lambdaT, prock_->permeability(ci), lambdaK);
        }


        template<template<typename> class SP>
        void getInverseMatrix(const CellIter&                        c,
                              FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
        {
            // Binv = (N*lambda*K*N'   +   t*diag(A)*(I - Q*Q')*diag(A))/vol
            //         ^                     ^^^^^^^^^^^^^^^^^^^^^^^^^^
            //         precomputed: n_       precomputed: second_term_
            // t = 6/dim * trace(lambda*K)
            int ci = c->index();
            int nf = Binv.numRows();
            ImmutableFortranMatrix n(nf, dim, &n_[ci][0]);
            ImmutableFortranMatrix t2(nf, nf, &second_term_[ci][0]);
            Binv = t2;
            ImmutableFortranMatrix lambdaK(dim, dim, lambdaK_.data());
            SharedFortranMatrix T2(nf, dim, &t2_[0]);

            // T2 <- N*lambda*K
            matMulAdd_NN(Scalar(1.0), n, lambdaK, Scalar(0.0), T2);

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


        template<class Vector>
        void gravityFlux(const CellIter& c,
                         Vector&         gflux) const
        {
            const int ci = c->index();
            const int nf = n_.rowSize(ci) / dim;

            ImmutableFortranMatrix N(nf, dim, &n_[ci][0]);

            // gflux = N (\sum_i \rho_i \lambda_i) Kg
            vecMulAdd_N(Scalar(1.0), N, &dyn_Kg_[0],
                        Scalar(0.0), &gflux[0]);
        }

    private:
        int                           max_nf_      ;
        mutable std::vector<Scalar>   fa_, t1_, t2_;
        Opm::SparseTable<Scalar>           second_term_ ;
        Opm::SparseTable<Scalar>           n_           ;
        Opm::SparseTable<Scalar>           Kg_          ;
        std::array<Scalar, dim>     dyn_Kg_      ;
        std::array<double, dim*dim> lambdaK_     ;
        const RockInterface*          prock_       ;
    };
} // namespace Opm


#endif // OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER