UpscalerBase_impl.hpp

//===========================================================================
//
// File: UpscalerBase_impl.hpp
//
// Created: Thu Apr 29 10:22:06 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================

/*
  Copyright 2010 SINTEF ICT, Applied Mathematics.
  Copyright 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 OPM_UPSCALERBASE_IMPL_HEADER
#define OPM_UPSCALERBASE_IMPL_HEADER

#include <opm/porsol/common/setupGridAndProps.hpp>
#include <opm/porsol/common/setupBoundaryConditions.hpp>
#include <opm/porsol/common/ReservoirPropertyTracerFluid.hpp>

#include <iostream>

namespace Opm
{

    template <class Traits>
    inline UpscalerBase<Traits>::UpscalerBase()
        : bctype_(Fixed),
          twodim_hack_(false),
          residual_tolerance_(1e-8),
          linsolver_maxit_(0),
          linsolver_prolongate_factor_(1.0),
          linsolver_verbosity_(0),
          linsolver_type_(3),
          linsolver_smooth_steps_(1),
          gravity_(0.0)
#if HAVE_MPI
         , grid_(MPI_COMM_SELF)
#endif
    {
    }




    template <class Traits>
    inline void UpscalerBase<Traits>::init(const Opm::ParameterGroup& param)
    {
        initImpl(param);
        initFinal(param);
    }




    template <class Traits>
    inline void UpscalerBase<Traits>::initImpl(const Opm::ParameterGroup& param)
    {
        // Request the boundary condition type parameter early since,
        // depending on the actual type, we may have to manufacture
        // and insert other parameters into the ParameterGroup prior
        // to constructing the grid and associated properties.
        //
        int bct = param.get<int>("boundary_condition_type");
        bctype_ = static_cast<BoundaryConditionType>(bct);

        // Import control parameters pertaining to reduced physical
        // dimensionality ("2d_hack = true" precludes periodic
        // boundary conditions in the Z direction), and linear solves.
        //
        twodim_hack_ = param.getDefault("2d_hack", twodim_hack_);
        residual_tolerance_ = param.getDefault("residual_tolerance", residual_tolerance_);
        linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_);
        linsolver_type_ = param.getDefault("linsolver_type", linsolver_type_);
        linsolver_maxit_ = param.getDefault("linsolver_max_iterations", linsolver_maxit_);
        linsolver_prolongate_factor_ = param.getDefault("linsolver_prolongate_factor", linsolver_prolongate_factor_);
        linsolver_smooth_steps_ = param.getDefault("linsolver_smooth_steps", linsolver_smooth_steps_);
        gravity_ = param.getDefault("gravity", gravity_);

        // Ensure sufficient grid support for requested boundary
        // condition type.
        //
        Opm::ParameterGroup temp_param = param;
        if (bctype_ == Linear || bctype_ == Periodic) {
            if (!temp_param.has("use_unique_boundary_ids")) {
                temp_param.insertParameter("use_unique_boundary_ids", "true");
            }
            if (!temp_param.has("clip_z")) {
                temp_param.insertParameter("clip_z", "true");
            }
        }
        if (bctype_ == Periodic) {
            if (!temp_param.has("periodic_extension")) {
                temp_param.insertParameter("periodic_extension", "true");
            }
        }

        setupGridAndProps(temp_param, grid_, res_prop_);
        ginterf_.init(grid_);
    }




    template <class Traits>
    inline void UpscalerBase<Traits>::initFinal(const Opm::ParameterGroup& param)
    {
        // Report any unused parameters.
        std::cout << "====================   Unused parameters:   ====================\n";
        param.displayUsage();
        std::cout << "================================================================\n";
    }




    template <class Traits>
    inline void UpscalerBase<Traits>::init(const Opm::Deck& deck,
                                           BoundaryConditionType bctype,
                                           double perm_threshold,
                                           double residual_tolerance,
                                           int linsolver_verbosity,
                                           int linsolver_type,
                                           bool twodim_hack,
                                           int linsolver_maxit,
                                           double linsolver_prolongate_factor,
                                           int linsolver_smooth_steps,
                                           const double gravity)
    {
        bctype_ = bctype;
        residual_tolerance_ = residual_tolerance;
        linsolver_verbosity_ = linsolver_verbosity;
        linsolver_type_ = linsolver_type;
            linsolver_maxit_ = linsolver_maxit;
         linsolver_prolongate_factor_ = linsolver_prolongate_factor;
        linsolver_smooth_steps_ = linsolver_smooth_steps;
        twodim_hack_ = twodim_hack;
        gravity_ = gravity;


        // Faking some parameters depending on bc type.
        bool periodic_ext = (bctype_ == Periodic);
        bool turn_normals = false;
        bool clip_z = (bctype_ == Linear || bctype_ == Periodic);
        bool unique_bids = (bctype_ == Linear || bctype_ == Periodic);
        std::string rock_list("no_list");
        setupGridAndPropsEclipse(deck,
                                 periodic_ext, turn_normals, clip_z, unique_bids,
                                 perm_threshold, rock_list,
                                 useJ<ResProp>(), 1.0, 0.0,
                                 grid_, res_prop_);
        ginterf_.init(grid_);
    }




    template <class Traits>
    inline const typename UpscalerBase<Traits>::GridType&
    UpscalerBase<Traits>::grid() const
    {
        return grid_;
    }




    template <class Traits>
    inline void
    UpscalerBase<Traits>::setBoundaryConditionType(BoundaryConditionType type)
    {
        if ((type == Periodic && bctype_ != Periodic)
            || (type != Periodic && bctype_ == Periodic)) {
            OPM_THROW(std::runtime_error, "Cannot switch to or from Periodic boundary condition, "
                  "periodic must be set in init() params.");
        } else {
            bctype_ = type;
            if (type == Periodic || type == Linear) {
                grid_.setUniqueBoundaryIds(true);
            } else {
                grid_.setUniqueBoundaryIds(false);
            }
        }
    }




    template <class Traits>
    inline void
    UpscalerBase<Traits>::setPermeability(const int cell_index, const permtensor_t& k)
    {
        res_prop_.permeabilityModifiable(cell_index) = k;
    }




    template <class Traits>
    inline typename UpscalerBase<Traits>::permtensor_t
    UpscalerBase<Traits>::upscaleSinglePhase()
    {
        ReservoirPropertyTracerFluid fluid;
        return upscaleEffectivePerm(fluid);
    }




    template <class Traits>
    template <class FluidInterface>
    inline typename UpscalerBase<Traits>::permtensor_t
    UpscalerBase<Traits>::upscaleEffectivePerm(const FluidInterface& fluid)
    {
        int num_cells = ginterf_.numberOfCells();
        // No source or sink.
        std::vector<double> src(num_cells, 0.0);
        // Just water.
        std::vector<double> sat(num_cells, 1.0);
        // Gravity.
        Dune::FieldVector<double, 3> gravity(0.0);
        gravity[2] = gravity_;

        permtensor_t upscaled_K(3, 3, (double*)0);
        for (int pdd = 0; pdd < Dimension; ++pdd) {
            setupUpscalingConditions(ginterf_, bctype_, pdd, 1.0, 1.0, twodim_hack_, bcond_);
            if (pdd == 0) {
                // Only on first iteration, since we do not change the
                // structure of the system, the way the flow solver is
                // implemented.
                flow_solver_.init(ginterf_, res_prop_, gravity, bcond_);
            }

            // Run pressure solver.
            bool same_matrix = (bctype_ != Fixed) && (pdd != 0);
            flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_,
                               linsolver_verbosity_,
                               linsolver_type_, same_matrix,
                               linsolver_maxit_, linsolver_prolongate_factor_,
                               linsolver_smooth_steps_);
            double max_mod = flow_solver_.postProcessFluxes();
            std::cout << "Max mod = " << max_mod << std::endl;

            // Compute upscaled K.
            double Q[Dimension] =  { 0 };
            switch (bctype_) {
            case Fixed:
                Q[pdd] = computeAverageVelocity(flow_solver_.getSolution(), pdd, pdd);
                break;
            case Linear:
            case Periodic:
                for (int i = 0; i < Dimension; ++i) {
                    Q[i] = computeAverageVelocity(flow_solver_.getSolution(), i, pdd);
                }
                break;
            default:
        OPM_THROW(std::runtime_error, "Unknown boundary type: " + std::to_string(bctype_));
            }
            double delta = computeDelta(pdd);
            for (int i = 0; i < Dimension; ++i) {
                upscaled_K(i, pdd) = Q[i] * delta;
            }
        }
        return upscaled_K;
    }




    template <class Traits>
    template <class FlowSol>
    double UpscalerBase<Traits>::computeAverageVelocity(const FlowSol& flow_solution,
                                                               const int flow_dir,
                                                               const int pdrop_dir) const
    {
        double side1_flux = 0.0;
        double side2_flux = 0.0;
        double side1_area = 0.0;
        double side2_area = 0.0;

        for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
            for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
                if (f->boundary()) {
                    int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
                    if ((canon_bid - 1)/2 == flow_dir) {
                        double flux = flow_solution.outflux(f);
                        double area = f->area();
                        double norm_comp = f->normal()[flow_dir];
                        // std::cout << "bid " << f->boundaryId() << "   area " << area << "   n " << norm_comp << std::endl;
                        if (canon_bid - 1 == 2*flow_dir) {
                            if (flow_dir == pdrop_dir && flux > 0.0) {
#ifdef VERBOSE
                              std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()<<" (canonical: "<<canon_bid
                                          << ") Magnitude: " << std::fabs(flux) << std::endl;
#endif
                                // OPM_THROW(std::runtime_error, "Detected outflow at entry face: " << face);
                            }
                            side1_flux += flux*norm_comp;
                            side1_area += area;
                        } else {
                            assert(canon_bid - 1 == 2*flow_dir + 1);
                            if (flow_dir == pdrop_dir && flux < 0.0) {
#ifdef VERBOSE
                                std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()
                                          << " Magnitude: " << std::fabs(flux) << std::endl;
#endif
                                // OPM_THROW(std::runtime_error, "Detected inflow at exit face: " << face);
                            }
                            side2_flux += flux*norm_comp;
                            side2_area += area;
                        }
                    }
                }
            }
        }
        // q is the average velocity.
        return 0.5*(side1_flux/side1_area + side2_flux/side2_area);
    }




    template <class Traits>
    inline double UpscalerBase<Traits>::computeDelta(const int flow_dir) const
    {
        double side1_pos = 0.0;
        double side2_pos = 0.0;
        double side1_area = 0.0;
        double side2_area = 0.0;
        for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
            for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
                if (f->boundary()) {
                    int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
                    if ((canon_bid - 1)/2 == flow_dir) {
                        double area = f->area();
                        double pos_comp = f->centroid()[flow_dir];
                        if (canon_bid - 1 == 2*flow_dir) {
                            side1_pos += area*pos_comp;
                            side1_area += area;
                        } else {
                            side2_pos += area*pos_comp;
                            side2_area += area;
                        }
                    }
                }
            }
        }
        // delta is the average length.
        return  side2_pos/side2_area - side1_pos/side1_area;
    }




    template <class Traits>
    double UpscalerBase<Traits>::upscalePorosity() const
    {
        double total_vol = 0.0;
        double total_pore_vol = 0.0;
        for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
            total_vol += c->volume();
            total_pore_vol += c->volume()*res_prop_.porosity(c->index());
        }
        return total_pore_vol/total_vol;
    }


    template <class Traits>
    double UpscalerBase<Traits>::upscaleNetPorosity() const
    {
        double total_net_vol = 0.0;
        double total_pore_vol = 0.0;
        for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
            total_net_vol += c->volume()*res_prop_.ntg(c->index());
            total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
        }
        if (total_net_vol>0.0) return total_pore_vol/total_net_vol;
        else return 0.0;
    }

    template <class Traits>
    double UpscalerBase<Traits>::upscaleNTG() const
    {
        double total_vol = 0.0;
        double total_net_vol = 0.0;
        for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
            total_vol += c->volume();
            total_net_vol += c->volume()*res_prop_.ntg(c->index());
        }
        return total_net_vol/total_vol;
    }

    template <class Traits>
    double UpscalerBase<Traits>::upscaleSWCR(const bool NTG) const
    {
        double total_swcr = 0.0;
        double total_pore_vol = 0.0;
        if (NTG) {
            for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
                total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.swcr(c->index());
                total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
            }
        }
        else {
            for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
                total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.swcr(c->index());
                total_pore_vol += c->volume()*res_prop_.porosity(c->index());
            }
        }
        return total_swcr/total_pore_vol;
    }

    template <class Traits>
    double UpscalerBase<Traits>::upscaleSOWCR(const bool NTG) const
    {
        double total_sowcr = 0.0;
        double total_pore_vol = 0.0;
        if (NTG) {
            for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
                total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.sowcr(c->index());
                total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
            }
        }
        else {
            for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
                total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.sowcr(c->index());
                total_pore_vol += c->volume()*res_prop_.porosity(c->index());
            }
        }
        return total_sowcr/total_pore_vol;
    }

} // namespace Opm



#endif // OPM_UPSCALERBASE_IMPL_HEADER