/*
  Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
  Copyright 2016 Statoil 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_BLACKOILSEQUENTIALMODEL_HEADER_INCLUDED
#define OPM_BLACKOILSEQUENTIALMODEL_HEADER_INCLUDED


#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>

namespace Opm {

    struct BlackoilSequentialModelParameters : public BlackoilModelParameters
    {
        bool iterate_to_fully_implicit;
        explicit BlackoilSequentialModelParameters( const ParameterGroup& param )
            : BlackoilModelParameters(param),
              iterate_to_fully_implicit(param.getDefault("iterate_to_fully_implicit", false))
        {
        }
    };


    /// A sequential splitting model implementation for three-phase black oil.
    template<class Grid, class WellModel,
             template <class G, class W> class PressureModelT,
             template <class G, class W> class TransportModelT>
    class BlackoilSequentialModel
    {
    public:
        typedef BlackoilState ReservoirState;
        typedef WellStateFullyImplicitBlackoil WellState;
        typedef BlackoilSequentialModelParameters ModelParameters;
        typedef DefaultBlackoilSolutionState SolutionState;

        typedef PressureModelT<Grid, WellModel> PressureModel;
        typedef TransportModelT<Grid, WellModel> TransportModel;
        typedef NonlinearSolver<PressureModel> PressureSolver;
        typedef NonlinearSolver<TransportModel> TransportSolver;

        typedef typename TransportModel::SimulatorData SimulatorData;
        typedef typename TransportModel::FIPDataType   FIPDataType;

        /// Construct the model. It will retain references to the
        /// arguments of this functions, and they are expected to
        /// remain in scope for the lifetime of the solver.
        /// \param[in] param            parameters
        /// \param[in] grid             grid data structure
        /// \param[in] fluid            fluid properties
        /// \param[in] geo              rock properties
        /// \param[in] rock_comp_props  if non-null, rock compressibility properties
        /// \param[in] wells            well structure
        /// \param[in] vfp_properties   Vertical flow performance tables
        /// \param[in] linsolver        linear solver
        /// \param[in] eclState         eclipse state
        /// \param[in] has_disgas       turn on dissolved gas
        /// \param[in] has_vapoil       turn on vaporized oil feature
        /// \param[in] terminal_output  request output to cout/cerr
        BlackoilSequentialModel(const ModelParameters&          param,
                                const Grid&                     grid ,
                                const BlackoilPropsAdFromDeck& fluid,
                                const DerivedGeology&           geo  ,
                                const RockCompressibility*      rock_comp_props,
                                const WellModel                 well_model,
                                const NewtonIterationBlackoilInterface& linsolver,
                                std::shared_ptr< const EclipseState > eclState,
                                std::shared_ptr< const Schedule> schedule,
                                std::shared_ptr< const SummaryConfig> summary_config,
                                const bool has_disgas,
                                const bool has_vapoil,
                                const bool terminal_output)
        : pressure_model_(new PressureModel(param, grid, fluid, geo, rock_comp_props, well_model,
                                            linsolver, eclState, schedule, summary_config, has_disgas, has_vapoil, terminal_output)),
          transport_model_(new TransportModel(param, grid, fluid, geo, rock_comp_props, well_model,
                                              linsolver, eclState, schedule, summary_config, has_disgas, has_vapoil, terminal_output)),
          // TODO: fix solver parameters for pressure and transport solver.
          pressure_solver_(typename PressureSolver::SolverParameters(), std::move(pressure_model_)),
          transport_solver_(typename TransportSolver::SolverParameters(), std::move(transport_model_)),
          initial_reservoir_state_(0, 0, 0), // will be overwritten
          iterate_to_fully_implicit_(param.iterate_to_fully_implicit)
        {
            typename PressureSolver::SolverParameters pp;
            pp.min_iter_ = 0;
            pressure_solver_.setParameters(pp);
            typename TransportSolver::SolverParameters tp;
            tp.min_iter_ = 0;
            transport_solver_.setParameters(tp);
        }





        /// Called once before each time step.
        /// \param[in] timer             simulation timer
        /// \param[in] reservoir_state   reservoir state variables
        /// \param[in] well_state        well state variables
        void prepareStep(const SimulatorTimerInterface& /*timer*/,
                         const ReservoirState& reservoir_state,
                         const WellState& well_state)
        {
            initial_reservoir_state_ = reservoir_state;
            initial_well_state_ = well_state;
        }





        /// Called once per nonlinear iteration.
        /// This model will first solve the pressure model to convergence, then the
        /// transport model.
        /// \param[in] iteration              should be 0 for the first call of a new timestep
        /// \param[in] timer             simulation timer
        /// \param[in] nonlinear_solver       nonlinear solver used (for oscillation/relaxation control)
        /// \param[in, out] reservoir_state   reservoir state variables
        /// \param[in, out] well_state        well state variables
        template <class NonlinearSolverType>
        SimulatorReport nonlinearIteration(const int iteration,
                                           const SimulatorTimerInterface& timer,
                                           NonlinearSolverType& /* nonlinear_solver */,
                                           ReservoirState& reservoir_state,
                                           WellState& well_state)
        {
            if (!iterate_to_fully_implicit_) {
                // Do a single pressure solve, followed by a single transport solve.
                if (terminalOutputEnabled()) {
                    OpmLog::info("Using sequential model.");
                }

                // Pressure solve.
                if (terminalOutputEnabled()) {
                    OpmLog::info("Solving the pressure equation.");
                }
                ReservoirState initial_state = reservoir_state;
                const SimulatorReport pressure_report = pressure_solver_.step(timer, reservoir_state, well_state);
                const int pressure_liniter = pressure_report.total_linear_iterations;
                if (pressure_liniter == -1) {
                    OPM_THROW(std::runtime_error, "Pressure solver failed to converge.");
                }

                // Transport solve.
                if (terminalOutputEnabled()) {
                    OpmLog::info("Solving the transport equations.");
                }
                const SimulatorReport transport_report = transport_solver_.step(timer, initial_state, well_state, reservoir_state, well_state);
                const int transport_liniter = transport_report.total_linear_iterations;
                if (transport_liniter == -1) {
                    OPM_THROW(std::runtime_error, "Transport solver failed to converge.");
                }

                // Report and return.
                SimulatorReport report;
                report.converged = true;
                report.total_linear_iterations = pressure_liniter + transport_liniter;
                return report;
            } else {
                // Iterate to fully implicit solution.
                // This call is just for a single iteration (one pressure and one transport solve),
                // we return a 'false' converged status if more are needed
                if (terminalOutputEnabled()) {
                    OpmLog::info("Using sequential model in iterative mode, outer iteration " + std::to_string(iteration));
                }

                // Pressure solve.
                if (terminalOutputEnabled()) {
                    OpmLog::info("Solving the pressure equation.");
                }

                const SimulatorReport pressure_report = pressure_solver_.step(timer, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
                const int pressure_liniter = pressure_report.total_linear_iterations;
                if (pressure_liniter == -1) {
                    OPM_THROW(std::runtime_error, "Pressure solver failed to converge.");
                }

                // Transport solve.
                if (terminalOutputEnabled()) {
                    OpmLog::info("Solving the transport equations.");
                }
                const SimulatorReport transport_report = transport_solver_.step(timer, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
                const int transport_liniter = transport_report.total_linear_iterations;
                if (transport_liniter == -1) {
                    OPM_THROW(std::runtime_error, "Transport solver failed to converge.");
                }

                // Revisit pressure equation to check if it is still converged.
                bool done = false;
                {
                    auto rstate = reservoir_state;
                    auto wstate = well_state;
                    pressure_solver_.model().prepareStep(timer, initial_reservoir_state_, initial_well_state_);
                    SimulatorReport rep = pressure_solver_.model().nonlinearIteration(0, timer, pressure_solver_, rstate, wstate);
                    if (rep.converged && rep.total_newton_iterations == 0) {
                        done = true;
                    }
                }

                SimulatorReport report;
                report.converged = done;
                report.total_linear_iterations = pressure_liniter + transport_liniter;
                return report;
            }
        }





        /// Called once after each time step.
        /// In this class, this function does nothing.
        /// \param[in] timer                  simulation timer
        /// \param[in, out] reservoir_state   reservoir state variables
        /// \param[in, out] well_state        well state variables
        void afterStep(const SimulatorTimerInterface& /* timer */,
                       ReservoirState& /* reservoir_state */,
                       WellState& /* well_state */)
        {
        }





        /// \brief Set threshold pressures that prevent or reduce flow.
        /// This prevents flow across faces if the potential
        /// difference is less than the threshold. If the potential
        /// difference is greater, the threshold value is subtracted
        /// before calculating flow. This is treated symmetrically, so
        /// flow is prevented or reduced in both directions equally.
        /// \param[in]  threshold_pressures_by_face   array of size equal to the number of faces
        ///                                   of the grid passed in the constructor.
        void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face)
        {
            pressure_solver_.model().setThresholdPressures(threshold_pressures_by_face);
            transport_solver_.model().setThresholdPressures(threshold_pressures_by_face);
        }





        /// Return true if output to cout is wanted.
        bool terminalOutputEnabled() const
        {
            return pressure_solver_.model().terminalOutputEnabled();
        }





        /// Return the relative change in variables relevant to this model.
        double relativeChange(const SimulationDataContainer& previous,
                              const SimulationDataContainer& current ) const
        {
            // TODO: this is a quick stopgap implementation, and should be evaluated more carefully.
            return std::max(pressure_solver_.model().relativeChange(previous, current),
                            transport_solver_.model().relativeChange(previous, current));
        }

        /// Return the well model
        const WellModel& wellModel() const
        {
            return pressure_solver_.model().wellModel();
        }


        /// Compute fluid in place.
        /// \param[in]    ReservoirState
        /// \param[in]    FIPNUM for active cells not global cells.
        /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
        std::vector<std::vector<double> >
        computeFluidInPlace(const ReservoirState& x,
                            const std::vector<int>& fipnum) const
        {
            return transport_solver_.computeFluidInPlace(x, fipnum);
        }


        /// Return reservoir simulation data (for output functionality)
        const SimulatorData& getSimulatorData(const SimulationDataContainer& localState) const {
            return transport_solver_.model().getSimulatorData(localState);
        }

        /// Return fluid-in-place data (for output functionality)
        FIPDataType getFIPData() const {
            return transport_solver_.model().getFIPData();
        }

        /// return the statistics if the nonlinearIteration() method failed.
        ///
        /// NOTE: for the flow_legacy simulator family this method is a stub, i.e. the
        /// failure report object will *not* contain any meaningful data.
        const SimulatorReport& failureReport() const
        { return failureReport_; }

    protected:
        SimulatorReport failureReport_;

        std::unique_ptr<PressureModel> pressure_model_;
        std::unique_ptr<TransportModel> transport_model_;
        PressureSolver pressure_solver_;
        TransportSolver transport_solver_;

        ReservoirState initial_reservoir_state_;
        WellState initial_well_state_;

        bool iterate_to_fully_implicit_;
    };

} // namespace Opm



#endif // OPM_BLACKOILSEQUENTIALMODEL_HEADER_INCLUDED