// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   Copyright (C) 2012 by                                                   *
 *   Institute for Modelling Hydraulic and Environmental Systems             *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
 *   This program 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 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 *
 * \brief Definition of an injection problem for the generic leakage scnearios
 * of the ECO2 project
 *
 * Up to now five different scenarios have been developed:
 * a chimney scenario
 * a fault scenario
 * a leaky well scenario
 * a catastrophic blowout scneario
 *
 * This file contains the boundary conditions of the chimney scenario
 * the boundary conditions for the three remaining scenarios are given in comments
 */
#ifndef DUMUX_ECO2_PROBLEM_HH
#define DUMUX_ECO2_PROBLEM_HH

#include <dune/grid/io/file/dgfparser/dgfs.hh>

#include "eco2model.hh"
#include <dumux/implicit/co2/co2volumevariables.hh>
#include <dumux/material/fluidsystems/brineco2fluidsystem.hh>
#include <dumux/implicit/common/implicitporousmediaproblem.hh>
#include <dumux/implicit/box/intersectiontovertexbc.hh>

#include "eco2spatialparams.hh"
#include "eco2co2tables.hh"

namespace Dumux
{

template <class TypeTag>
class ECO2Problem;

namespace Properties
{
NEW_TYPE_TAG(ECO2Problem, INHERITS_FROM(BoxTwoPTwoC, ECO2SpatialParams));

// Set the grid type
SET_PROP(ECO2Problem, Grid)
{
    typedef Dune::ALUSimplexGrid<3,3> type;
};

// Set the problem property
SET_PROP(ECO2Problem, Problem)
{
    typedef Dumux::ECO2Problem<TTAG(ECO2Problem)> type;
};

// Set fluid configuration
SET_PROP(ECO2Problem, FluidSystem)
{
    typedef Dumux::BrineCO2FluidSystem<TypeTag> type;
};

// Set the CO2 table to be used; in this case not the the default table
SET_TYPE_PROP(ECO2Problem, CO2Table, Dumux::ECO2CO2Tables::CO2Tables);

// Set the salinity mass fraction of the brine in the reservoir
SET_SCALAR_PROP(ECO2Problem, ProblemSalinity, 1e-1);

//! the CO2 Model and VolumeVariables properties
SET_TYPE_PROP(ECO2Problem, Model, CO2Model<TypeTag>);
SET_TYPE_PROP(ECO2Problem, VolumeVariables, CO2VolumeVariables<TypeTag>);

// Enable gravity
SET_BOOL_PROP(ECO2Problem, ProblemEnableGravity, true);
SET_BOOL_PROP(ECO2Problem, ImplicitEnableJacobianRecycling, false);
SET_BOOL_PROP(ECO2Problem, VtkAddVelocity, true);
}


/*!
 * \ingroup TwoPTwoCModel
 * \ingroup ECO2Problem
 * \brief Definition of an injection problem for the generic leakage scnearios of the ECO2 project
 *
 * To run the simulation execute the following line in shell:
 * <tt>./eco2_scenarios -parameterFile ./eco2_scenarios.input</tt>
 */
template <class TypeTag = TTAG(ECO2Problem)>
class ECO2Problem : public ImplicitPorousMediaProblem<TypeTag>
{
    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;

    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;

    enum {
        // Grid and world dimension
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };

    // copy some indices for convenience
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    enum {
        lPhaseIdx = Indices::wPhaseIdx,
        gPhaseIdx = Indices::nPhaseIdx,


        BrineIdx = FluidSystem::BrineIdx,
        CO2Idx = FluidSystem::CO2Idx,

        conti0EqIdx = Indices::conti0EqIdx,
        contiCO2EqIdx = conti0EqIdx + CO2Idx
    };

    typedef Dune::FieldVector<Scalar, dim> DimVector;

    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;

    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::Intersection Intersection;

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(CO2Table)) CO2Table;
    typedef Dumux::CO2<Scalar, CO2Table> CO2;

public:
    /*!
     * \brief The constructor
     *
     * \param timeManager The time manager
     * \param gridView The grid view
     */
    ECO2Problem(TimeManager &timeManager,
                     const GridView &gridView)
        : ParentType(timeManager, GridCreator::grid().leafView())
    {
        try
        {
            nTemperature_       = GET_RUNTIME_PARAM(TypeTag, int, FluidSystem.NTemperature);
            nPressure_          = GET_RUNTIME_PARAM(TypeTag, int, FluidSystem.NPressure);
            pressureLow_        = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.PressureLow);
            pressureHigh_       = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.PressureHigh);
            temperatureLow_     = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.TemperatureLow);
            temperatureHigh_    = GET_RUNTIME_PARAM(TypeTag, Scalar, FluidSystem.TemperatureHigh);
            depthBOR_           = GET_RUNTIME_PARAM(TypeTag, Scalar, InitialConditions.DepthBOR);
            depthSeafloor_		= GET_RUNTIME_PARAM(TypeTag, Scalar, InitialConditions.DepthSeafloor);
            name_               = GET_RUNTIME_PARAM(TypeTag, std::string, SimulationControl.Name);
            chimneyScenario_	= GET_RUNTIME_PARAM(TypeTag, bool, Problem.ChimneyScenario);
            catastrophicBlowoutScenario_ = GET_RUNTIME_PARAM(TypeTag, bool, Problem.CatastrophicBlowoutScenario);
            faultScenario_	= GET_RUNTIME_PARAM(TypeTag, bool, Problem.FaultScenario);
            leakyWellScenario_	= GET_RUNTIME_PARAM(TypeTag, bool, Problem.LeakyWellScenario);
        }
        catch (Dumux::ParameterException &e) {
            std::cerr << e << ". Abort!\n";
            exit(1) ;
        }
        catch (...) {
            std::cerr << "Unknown exception thrown!\n";
            exit(1);
        }

        eps_ = 1e-6;

        // initialize the tables of the fluid system
        //FluidSystem::init();
        FluidSystem::init(/*Tmin=*/temperatureLow_,
                          /*Tmax=*/temperatureHigh_,
                          /*nT=*/nTemperature_,
                          /*pmin=*/pressureLow_,
                          /*pmax=*/pressureHigh_,
                          /*np=*/nPressure_);
        this->timeManager().startNextEpisode(3.1536e7);
    }

    /*!
     * \brief Returns true if the current solution should be written to
     *        disk (i.e. as a VTK file)
     *
     * The default behaviour is to write out every the solution for
     * very time step. This file is intented to be overwritten by the
     * implementation.
     */
    bool shouldWriteOutput() const
    {
    	// output files are written according to episode lengths
        return
            this->timeManager().timeStepIndex() == 0 ||
            this->timeManager().episodeWillBeOver() ||
            this->timeManager().willBeFinished();
    }

    /*!
     * \brief Returns true if a restart file should be written to
     *        disk.
     *
     * The default behaviour is to write one restart file every 5 time
     * steps. This file is intented to be overwritten by the
     * implementation.
     */
    bool shouldWriteRestartFile() const
    {
    	// restart files are only written for certain episodes (250s, 1, 10, 20, 30, 40, 50, 100, 150, 200 years)
        return 
            this->timeManager().timeStepIndex() == 0 ||
            (this->timeManager().episodeIndex() == 1 && this->timeManager().episodeIsOver()) ||
            (this->timeManager().episodeIndex() == 10 && this->timeManager().episodeIsOver()) ||
            (this->timeManager().episodeIndex() == 11 && this->timeManager().episodeIsOver()) ||
            (this->timeManager().episodeIndex() == 12 && this->timeManager().episodeIsOver()) ||
            (this->timeManager().episodeIndex() == 13 && this->timeManager().episodeIsOver()) ||
            (this->timeManager().episodeIndex() == 14 && this->timeManager().episodeIsOver()) ||
            (this->timeManager().episodeIndex() == 15 && this->timeManager().episodeIsOver()) ||
            (this->timeManager().episodeIndex() == 16 && this->timeManager().episodeIsOver()) ||
            this->timeManager().willBeFinished();
    }

    /*!
     * \brief Returns true if a data file should be written to
     *        disk.
     *
     * The data files contain cell-wise pressure and saturation information
     * for the geomechanical evaluation performed by TNO.
     */
    bool shouldWriteDataFile() const
    {
    	// restart files are only written for certain episodes (250s, 1, 2, 5, 10, 30, 40, 50, 100, 150, 200 years)
        return
            this->timeManager().timeStepIndex() == 0 ||
            (this->timeManager().episodeIndex() == 1 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 2 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 5 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 10 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 11 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 12 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 13 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 14 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 15 && this->timeManager().episodeWillBeOver()) ||
            (this->timeManager().episodeIndex() == 16 && this->timeManager().episodeWillBeOver()) ||
            this->timeManager().willBeFinished();
    }

    /*!
     * \brief Called directly after the time integration.
     */
    void postTimeStep()
    {
        // Calculate storage terms
        PrimaryVariables storageL, storageG;
        Dune::FieldVector<Scalar, 2> flux(0.0);
        this->model().globalPhaseStorage(storageL, lPhaseIdx);
        this->model().globalPhaseStorage(storageG, gPhaseIdx);

        //calculating flux across the layer of the coordinate 2 (Z axis)
        //intercepted at coordinate value = 999m (just below the top of the cap rock)
        this->model().calculateFluxAcrossLayer(flux, 2, 999);

        // Write mass balance information for rank 0
        if (this->gridView().comm().rank() == 0) {
                    std::cout<<"Storage: liquid=[" << storageL << "]"
                             << " gas=[" << storageG << "]"<< "flux=[" <<flux << "]\n";
        }

    }

    void episodeEnd()
    {
        // Start new episode if episode is over
    	// for first 10 year episode length is 1 year
        if(this->timeManager().time() < 3.1536e8 ){
            this->timeManager().startNextEpisode(3.1536e7);
            std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
        }
        // until 50 years episode length is 10 years
        else if(this->timeManager().time() < 1.5768e9){
        	this->timeManager().startNextEpisode(3.1536e8);
        	std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
        }
        // until end of simulation (200 years) episode length is 50 years
        else{
            this->timeManager().startNextEpisode(1.5768e9);
            std::cout<<"Episode index is set to: "<<this->timeManager().episodeIndex()<<std::endl;
        }
    }


    /*!
     * \name Problem parameters
     */
    // \{

    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     */
    const std::string name() const
    { return name_; }

    /*!
     * \brief Returns the temperature within the domain.
     *
     * This problem assumes a temperature of 10 degrees Celsius.
     */
    Scalar temperatureAtPos(const GlobalPosition &pos) const
    {
    	return 273.15 + 4 + (depthBOR_ - depthSeafloor_ - pos[2]) * 0.032;
    };

    void sourceAtPos(PrimaryVariables &values,
                const GlobalPosition &globalPos) const
    {
        values = 0;
    }

    // \}

    /*!
     * \name Boundary conditions
     */
    // \{

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
     * \param values The boundary types for the conservation equations
     * \param vertex The vertex for which the boundary type is set
     */
    void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
    {
        const GlobalPosition globalPos = vertex.geometry().center();

        //set dirichlet on the top boundary for evaluation of flux across layer
        if (onLeftBoundary(globalPos) || onRightBoundary(globalPos) || onBackBoundary(globalPos) || onUpperBoundary(globalPos))
            values.setAllDirichlet();

        else if(onLowerBoundary(globalPos))
            values.setAllNeumann();

        else
            values.setAllNeumann();
    }

    /*!
     * \brief Evaluate the boundary conditions for a dirichlet
     *        boundary segment.
     *
     * \param values The dirichlet values for the primary variables
     * \param vertex The vertex for which the boundary type is set
     *
     * For this method, the \a values parameter stores primary variables.
     */
    void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
    {
        const GlobalPosition globalPos = vertex.geometry().center();

        initial_(values, globalPos);
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
     * \param values The neumann values for the conservation equations
     * \param element The finite element
     * \param fvElemGeom The finite-volume geometry in the box scheme
     * \param is The intersection between element and boundary
     * \param scvIdx The local vertex index
     * \param boundaryFaceIdx The index of the boundary face
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
     */
    void neumann(PrimaryVariables &values,
                 const Element &element,
                 const FVElementGeometry &fvElemGeom,
                 const Intersection &is,
                 int scvIdx,
                 int boundaryFaceIdx) const
    {
        const GlobalPosition &globalPos = element.geometry().corner(scvIdx);

        values = 0;

        // injection well boundary condition for the chimney scenario
        if (this->atWell(globalPos) && this->timeManager().time() <= 1.26144e9 )
        {
        	if(chimneyScenario_)
        		 values[contiCO2EqIdx] = -1.692912e-2; // kg/(s*m^2)

        	else if(catastrophicBlowoutScenario_)
        		values[contiCO2EqIdx] = -1.692912e-2; // kg/(s*m^2)

        	else if(faultScenario_)
        		values[contiCO2EqIdx] = -2.28557e-2; // kg/(s*m^2)

        	else if(leakyWellScenario_)
        		values[contiCO2EqIdx] = -4.14319e-2; // kg/(s*m^2)

        	else
        		std::cout<<"no scenario specified in input file!"<<std::endl;
        }
    }

    // \}

    /*!
     * \name Volume terms
     */
    // \{

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
     * \param values The initial values for the primary variables
     * \param element The finite element
     * \param fvElemGeom The finite-volume geometry in the box scheme
     * \param scvIdx The local vertex index
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
    void initial(PrimaryVariables &values,
                 const Element &element,
                 const FVElementGeometry &fvElemGeom,
                 int scvIdx) const
    {
        const GlobalPosition &globalPos = element.geometry().corner(scvIdx);

        initial_(values, globalPos);
    }

    /*!
     * \brief Return the initial phase state inside a control volume.
     *
     * \param vert The vertex
     * \param globalIdx The index of the global vertex
     * \param globalPos The global position
     */
    int initialPhasePresence(const Vertex &vert,
                             int &globalIdx,
                             const GlobalPosition &globalPos) const
    { return Indices::wPhaseOnly; }


    bool inLeakageFeature(const GlobalPosition &globalPos) const
	{
    	if(chimneyScenario_)
    		return (globalPos[0] > 3000 && globalPos[0] < 3500 && globalPos[1] < 250 && globalPos[2] > 100);

    	else if(catastrophicBlowoutScenario_)
    		return ((globalPos[0] > 2500 && globalPos[0] <= 3000 && globalPos[1] < 1.e-6 && globalPos[2] > 100 && globalPos[2] < 550) ||
    				(globalPos[0] > 3000 && globalPos[0] < 3500 && globalPos[1] < 250 && globalPos[2] > 100));

    	else if(faultScenario_)
    		return (globalPos[0] > 2999 && globalPos[0] < 3011) && (globalPos[2] > 100);

    	else if(leakyWellScenario_)
    		return ((globalPos[0] >2999 && globalPos[0] < 3002) && globalPos[1] < 2 && globalPos[2] > 100);

    	else{
    		std::cout<<"no scenario specified in input file!"<<std::endl;
    		return false;
    	}
    }

    bool atWell(const GlobalPosition &globalPos) const
    {
    	if(chimneyScenario_)
    		return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2485 && globalPos[0] < 2501);

    	else if(catastrophicBlowoutScenario_)
    		return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2485 && globalPos[0] < 2501);

    	else if(faultScenario_)
    		return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2499 && globalPos[0] < 2508);

    	else if(leakyWellScenario_)
    		return (globalPos[2] < 100 && globalPos[1] < eps_ && globalPos[0] > 2495 && globalPos[0] < 2503);

    	else{
    		std::cout<<"no scenario specified in input file!"<<std::endl;
    		return false;
    	}
    }


    bool onLeftBoundary(const GlobalPosition &globalPos) const
    { return globalPos[0] < this->bboxMin()[0] + eps_; }

    bool onRightBoundary(const GlobalPosition &globalPos) const
    { return globalPos[0] > this->bboxMax()[0] - eps_; }

    bool onLowerBoundary(const GlobalPosition &globalPos) const
    { return globalPos[2] < this->bboxMin()[2] + eps_; }

    bool onUpperBoundary(const GlobalPosition &globalPos) const
    { return globalPos[2] > this->bboxMax()[2] - eps_; }

    bool onFrontBoundary(const GlobalPosition &globalPos) const
    {return globalPos[1] < this->bboxMin()[1] + eps_; }

    bool onBackBoundary(const GlobalPosition &globalPos) const
    {return globalPos[1] > this->bboxMax()[1] - eps_; }
    // \}

private:
    // the internal method for the initial condition
    void initial_(PrimaryVariables &values,
                  const GlobalPosition &globalPos) const
    {
    	Scalar temperature = this->temperatureAtPos(globalPos);
        Scalar pressure = (1.013e5 + (depthBOR_ - globalPos[2]) * 1100 * 9.81);
        Scalar densityW = FluidSystem::Brine::liquidDensity(temperature, pressure);

        values[Indices::pressureIdx] = (1.013e5 + (depthBOR_ - globalPos[2]) * densityW * 9.81);
        values[Indices::switchIdx] = 0.0;
    }

    Scalar depthBOR_;
    Scalar depthSeafloor_;
    Scalar eps_;

    bool chimneyScenario_;
    bool catastrophicBlowoutScenario_;
    bool faultScenario_;
    bool leakyWellScenario_;

    int nTemperature_;
    int nPressure_;

    std::string name_ ;

    Scalar pressureLow_, pressureHigh_;
    Scalar temperatureLow_, temperatureHigh_;




};
} //end namespace

#endif