// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   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 Tutorial problem for a fully coupled two phase-two component box model.
 */
#ifndef DUMUX_EXERCISE_FLUIDSYSTEM_B_PROBLEM_HH
#define DUMUX_EXERCISE_FLUIDSYSTEM_B_PROBLEM_HH

// The base porous media box problem
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/common/properties.hh>

namespace Dumux {

/*!
 * \ingroup TwoPBoxModel
 *
 * \brief  Tutorial problem for a fully coupled two phase-two component box model.
 */
template <class TypeTag>
class ExerciseFluidsystemProblemTwoPTwoC : public PorousMediumFlowProblem<TypeTag>
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;

    // Grid dimension
    enum { dim = GridView::dimension,
           dimWorld = GridView::dimensionworld
    };
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;

    // Dumux specific types
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
    using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;

public:
    ExerciseFluidsystemProblemTwoPTwoC(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
    : ParentType(fvGridGeometry)
        , eps_(3e-6)
    {

        // initialize the fluid system
        FluidSystem::init();

        // set the depth of the bottom of the reservoir
        depthBOR_ = this->gridGeometry().bBoxMax()[dimWorld-1];
    }

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

    /*!
     * \brief Returns the temperature \f$ K \f$
     */
    Scalar temperature() const
    { return 283.15; }

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

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     *
     * \param globalPos The position for which the bc type should be evaluated
     */
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
         BoundaryTypes bcTypes;

        if (globalPos[0] < eps_ || globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
           bcTypes.setAllDirichlet();
        else
            bcTypes.setAllNeumann();

        return bcTypes;
    }

    /*!
     * \brief Evaluates the boundary conditions for a Dirichlet
     *        boundary segment
     *
     * \param globalPos The global position
     */
    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
    {
        // use initial as Dirichlet conditions
        return initialAtPos(globalPos);
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
     * \param globalPos The position of the integration point of the boundary segment.
     */
    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
    {
        // initialize values to zero, i.e. no-flow Neumann boundary conditions
        NumEqVector values(0.0);

        Scalar up = this->gridGeometry().bBoxMax()[dimWorld-1];
        // extraction of oil (30 g/m/s) on a segment of the upper boundary
        if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40)
        {
            // we solve for the mole balance, so we have to divide by the molar mass
            values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
            values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = -3e-2/FluidSystem::MyCompressibleComponent::molarMass();
        }
        else
        {
            // no-flow on the remaining Neumann-boundaries.
            values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
            values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = 0;
        }

        return values;
    }

    // \}

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

    /*!
     * \brief Evaluate the initial value for a control volume.
     * \param globalPos The position for which the initial condition should be evaluated
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const

    {
        PrimaryVariables values(0.0);

        // tell the primary variables the phase state, i.e. which phase/phases
        // is/are present, because this changes the meaning of the primary variable
        // value at the index Indices::switchIdx
        values.setState(Indices::firstPhaseOnly);

        // use hydrostatic pressure distribution with 2 bar at the top and zero saturation
        values[Indices::pressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar
        values[Indices::switchIdx] = 0.0;

        return values;
    }

    // \}

    /*!
     * \brief Returns the source term
     * \param globalPos The global position
     */
    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
    {
        // we do not define any sources
        NumEqVector values(0.0);
        return values;
    }


private:
    Scalar eps_; //! small epsilon value
    Scalar depthBOR_; //! depth at the bottom of the reservoir
};

} // end namespace Dumux

#endif