BlackoilModel.hpp

/*
  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
  Copyright 2014, 2015 Statoil ASA.
  Copyright 2015 NTNU
  Copyright 2015, 2016, 2017 IRIS AS

  This file is part of the Open Porous Media project (OPM).

  OPM 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.

  OPM 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 OPM.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED
#define OPM_BLACKOILMODEL_HEADER_INCLUDED

#include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>

#include <opm/simulators/flow/BlackoilModelConvergenceMonitor.hpp>
#include <opm/simulators/flow/BlackoilModelNldd.hpp>
#include <opm/simulators/flow/BlackoilModelProperties.hpp>
#include <opm/simulators/flow/FlowProblemBlackoilProperties.hpp>
#include <opm/simulators/flow/RSTConv.hpp>

#include <opm/simulators/linalg/ISTLSolver.hpp>

#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>

#include <opm/simulators/utils/ComponentName.hpp>

#include <opm/simulators/wells/BlackoilWellModel.hpp>

#include <memory>
#include <tuple>
#include <vector>

#include <fmt/format.h>

namespace Opm {

template <class TypeTag>
class BlackoilModel
{
public:
    // ---------  Types and enums  ---------
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Grid = GetPropType<TypeTag, Properties::Grid>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
    using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using ModelParameters = BlackoilModelParameters<Scalar>;

    static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();

    static constexpr int numEq = Indices::numEq;
    static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
    static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx;
    static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
    static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
    static constexpr int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
    static constexpr int contiFoamEqIdx = Indices::contiFoamEqIdx;
    static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
    static constexpr int contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
    static constexpr int contiOxygenEqIdx = Indices::contiOxygenEqIdx;
    static constexpr int contiUreaEqIdx = Indices::contiUreaEqIdx;
    static constexpr int contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
    static constexpr int contiCalciteEqIdx = Indices::contiCalciteEqIdx;
    static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
    static constexpr int zFractionIdx = Indices::zFractionIdx;
    static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
    static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
    static constexpr int temperatureIdx = Indices::temperatureIdx;
    static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
    static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
    static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
    static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
    static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
    static constexpr int biofilmVolumeFractionIdx = Indices::biofilmVolumeFractionIdx;
    static constexpr int calciteVolumeFractionIdx = Indices::calciteVolumeFractionIdx;

    using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
    using MatrixBlockType = typename SparseMatrixAdapter::MatrixBlock;
    using Mat = typename SparseMatrixAdapter::IstlMatrix;
    using BVector = Dune::BlockVector<VectorBlockType>;

    using ComponentName = ::Opm::ComponentName<FluidSystem,Indices>;

    // Helper structs
    struct CnvPvSplitData {
        std::pair<std::vector<double>, std::vector<int>> cnvPvSplit;
        std::vector<unsigned> ixCells;
    };

    struct MaxSolutionUpdateData {
        Scalar dPMax = 0.0;
        Scalar dSMax = 0.0;
        Scalar dRsMax = 0.0;
        Scalar dRvMax = 0.0;
    };

    // Output debug flags for which tolerances used
    enum class DebugFlags {
        STRICT = 0,
        RELAXED = 1,
        TUNINGDP = 2
    };

    // ---------  Public methods  ---------

    BlackoilModel(Simulator& simulator,
                  const ModelParameters& param,
                  BlackoilWellModel<TypeTag>& well_model,
                  const bool terminal_output);

    bool isParallel() const
    { return grid_.comm().size() > 1; }

    const EclipseState& eclState() const
    { return simulator_.vanguard().eclState(); }

    SimulatorReportSingle prepareStep(const SimulatorTimerInterface& timer);

    void initialLinearization(SimulatorReportSingle& report,
                              const int minIter,
                              const int maxIter,
                              const SimulatorTimerInterface& timer);

    template <class NonlinearSolverType>
    SimulatorReportSingle nonlinearIteration(const SimulatorTimerInterface& timer,
                                             NonlinearSolverType& nonlinear_solver);

    template <class NonlinearSolverType>
    SimulatorReportSingle nonlinearIterationNewton(const SimulatorTimerInterface& timer,
                                                   NonlinearSolverType& nonlinear_solver);

    SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface& timer);

    // compute the "relative" change of the solution between time steps
    Scalar relativeChange() const;

    int linearIterationsLastSolve() const
    { return simulator_.model().newtonMethod().linearSolver().iterations (); }

    // Obtain reference to linear solver setup time
    double& linearSolveSetupTime()
    { return linear_solve_setup_time_; }

    void solveJacobianSystem(BVector& x);

    void updateSolution(const BVector& dx);

    void prepareStoringSolutionUpdate();
    void storeSolutionUpdate(const BVector& dx);
    MaxSolutionUpdateData getMaxSolutionUpdate(const std::vector<unsigned>& ixCells);

    bool terminalOutputEnabled() const
    { return terminal_output_; }

    std::tuple<Scalar,Scalar>
    convergenceReduction(Parallel::Communication comm,
                         const Scalar pvSumLocal,
                         const Scalar numAquiferPvSumLocal,
                         std::vector<Scalar>& R_sum,
                         std::vector<Scalar>& maxCoeff,
                         std::vector<Scalar>& B_avg);

    std::pair<Scalar,Scalar>
    localConvergenceData(std::vector<Scalar>& R_sum,
                         std::vector<Scalar>& maxCoeff,
                         std::vector<Scalar>& B_avg,
                         std::vector<int>& maxCoeffCell);

    CnvPvSplitData characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt);

    void convergencePerCell(const std::vector<Scalar>& B_avg,
                            const double dt,
                            const double tol_cnv,
                            const double tol_cnv_energy);

    void updateTUNING(const Tuning& tuning);
    void updateTUNINGDP(const TuningDp& tuning_dp);

    ConvergenceReport
    getReservoirConvergence(const double reportTime,
                            const double dt,
                            const int maxIter,
                            std::vector<Scalar>& B_avg,
                            std::vector<Scalar>& residual_norms);

    ConvergenceReport
    getConvergence(const SimulatorTimerInterface& timer,
                   const int maxIter,
                   std::vector<Scalar>& residual_norms);

    int numPhases() const
    { return Indices::numPhases; }

    template<class T>
    std::vector<std::vector<Scalar> >
    computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
    { return this->computeFluidInPlace(fipnum); }

    std::vector<std::vector<Scalar> >
    computeFluidInPlace(const std::vector<int>& /*fipnum*/) const;

    const Simulator& simulator() const
    { return simulator_; }

    Simulator& simulator()
    { return simulator_; }

    const SimulatorReportSingle& failureReport() const
    { return failureReport_; }

    const SimulatorReport& localAccumulatedReports() const;

    const std::vector<SimulatorReport>& domainAccumulatedReports() const;

    void writeNonlinearIterationsPerCell(const std::filesystem::path& odir) const;

    const std::vector<StepReport>& stepReports() const
    { return convergence_reports_; }

    void popLastConvergenceReport()
    {
        convergence_reports_.back().report.pop_back();
        residual_norms_history_.pop_back();
    }

    void writePartitions(const std::filesystem::path& odir) const;

    BlackoilWellModel<TypeTag>&
    wellModel()
    { return well_model_; }

    const BlackoilWellModel<TypeTag>&
    wellModel() const
    { return well_model_; }

    void beginReportStep()
    { simulator_.problem().beginEpisode(); }

    void endReportStep()
    { simulator_.problem().endEpisode(); }

    template<class FluidState, class Residual>
    void getMaxCoeff(const unsigned cell_idx,
                     const IntensiveQuantities& intQuants,
                     const FluidState& fs,
                     const Residual& modelResid,
                     const Scalar pvValue,
                     std::vector<Scalar>& B_avg,
                     std::vector<Scalar>& R_sum,
                     std::vector<Scalar>& maxCoeff,
                     std::vector<int>& maxCoeffCell);

    const ModelParameters& param() const
    { return param_; }

    const ComponentName& compNames() const
    { return compNames_; }

    bool hasNlddSolver() const
    { return nlddSolver_ != nullptr; }

protected:
    // ---------  Data members  ---------
    Simulator& simulator_;
    const Grid& grid_;
    static constexpr bool has_solvent_ = getPropValue<TypeTag, Properties::EnableSolvent>();
    static constexpr bool has_extbo_ = getPropValue<TypeTag, Properties::EnableExtbo>();
    static constexpr bool has_polymer_ = getPropValue<TypeTag, Properties::EnablePolymer>();
    static constexpr bool has_polymermw_ = getPropValue<TypeTag, Properties::EnablePolymerMW>();
    static constexpr bool has_energy_ = getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal;
    static constexpr bool has_foam_ = getPropValue<TypeTag, Properties::EnableFoam>();
    static constexpr bool has_brine_ = getPropValue<TypeTag, Properties::EnableBrine>();
    static constexpr bool has_bioeffects_ = getPropValue<TypeTag, Properties::EnableBioeffects>();
    static constexpr bool has_micp_ = Indices::enableMICP;

    ModelParameters                 param_;
    SimulatorReportSingle failureReport_;

    // Well Model
    BlackoilWellModel<TypeTag>& well_model_;

    bool terminal_output_;
    long int global_nc_;

    std::vector<std::vector<Scalar>> residual_norms_history_;
    Scalar current_relaxation_;
    BVector dx_old_;

    SolutionVector solUpd_;

    std::vector<StepReport> convergence_reports_;
    ComponentName compNames_{};

    std::unique_ptr<BlackoilModelNldd<TypeTag>> nlddSolver_; 
    BlackoilModelConvergenceMonitor<Scalar> conv_monitor_;

private:
    Scalar dpMaxRel() const { return param_.dp_max_rel_; }
    Scalar dsMax() const { return param_.ds_max_; }
    Scalar drMaxRel() const { return param_.dr_max_rel_; }
    Scalar maxResidualAllowed() const { return param_.max_residual_allowed_; }
    double linear_solve_setup_time_;
    std::vector<bool> wasSwitched_;
};

} // namespace Opm

#include <opm/simulators/flow/BlackoilModel_impl.hpp>

#endif // OPM_BLACKOILMODEL_HEADER_INCLUDED