// -*- 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 The two-phase porousmediumflow problem for exercise runtime parameters
 */

#ifndef DUMUX_EXRUNTIMEPARAMS_INJECTION_PROBLEM_2P_HH
#define DUMUX_EXRUNTIMEPARAMS_INJECTION_PROBLEM_2P_HH

#include <dune/grid/yaspgrid.hh>

#include <dumux/discretization/cctpfa.hh>
#include <dumux/porousmediumflow/2p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>

#include "injection2pspatialparams.hh"

namespace Dumux {

// forward declare problem
template <class TypeTag>
class InjectionProblem2P;

namespace Properties {
// define the TypeTag for this problem with a cell-centered two-point flux approximation spatial discretization.
// Create new type tags
namespace TTag {
struct Injection2p { using InheritsFrom = std::tuple<TwoP>; };
struct Injection2pCC { using InheritsFrom = std::tuple<Injection2p, CCTpfaModel>; };
} // end namespace TTag

// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::Injection2p> { using type = Dune::YaspGrid<2>; };

// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::Injection2p> { using type = InjectionProblem2P<TypeTag>; };

// Set the spatial parameters
SET_TYPE_PROP(Injection2p, SpatialParams,
              InjectionSpatialParams<GetPropType<TypeTag, Properties::FVGridGeometry>,
                                     GetPropType<TypeTag, Properties::Scalar>>);

// Set fluid configuration
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::Injection2p> { using type = FluidSystems::H2ON2<GetPropType<TypeTag, Properties::Scalar>, FluidSystems::H2ON2DefaultPolicy</*fastButSimplifiedRelations=*/ true>>; };
} // end namespace Properties

/*!
 * \ingroup TwoPModel
 * \ingroup ImplicitTestProblems
 * \brief Gas injection problem where a gas (here  nitrogen) is injected into a fully
 *        water saturated medium. During buoyancy driven upward migration the gas
 *        passes a high temperature area.
 *
 * The domain is sized 60 m times 40 m.
 *
 * For the mass conservation equation neumann boundary conditions are used on
 * the top, on the bottom and on the right of the domain, while dirichlet conditions
 * apply on the left boundary.
 *
 * Gas is injected at the right boundary from 7 m to 15 m at a rate of
 * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
 * boundaries.
 *
 * At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
 *
 * This problem uses the \ref TwoPModel model.
 */
template<class TypeTag>
class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;

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

public:
    InjectionProblem2P(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
    : ParentType(fvGridGeometry)
    {
        // initialize the tables of the fluid system
        FluidSystem::init(/*tempMin=*/273.15,
                /*tempMax=*/423.15,
                /*numTemp=*/50,
                /*pMin=*/0.0,
                /*pMax=*/30e6,
                /*numP=*/300);

        // name of the problem and output file
        // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
        // of type TYPE given in the group GROUPNAME from the input file
        name_ = getParam<std::string>("Problem.Name");
        // depth of the aquifer, units: m
        aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
        // the duration of the injection, units: second
        injectionDuration_ = getParamFromGroup<Scalar>("Problem","InjectionDuration");
        //TODO: Task 2: Set a variable "TotalAreaSpecificInflow" to read in a value from the parameter tree via the input file
        totalAreaSpecificInflow_ = getParam<Scalar>("Problem.TotalAreaSpecificInflow");
        //TODO: Task 3: Set a default value for the above parameter.
//        totalAreaSpecificInflow_ = getParam<Scalar>("Problem.TotalAreaSpecificInflow", -1e-4);
        //TODO: Task 4: Provide output describing where the parameter value comes from using parameter bool functions.
//         if (hasParamInGroup("Problem","TotalAreaSpecificInflow"))
//             std::cout << "Parameter value is read from the input file." << std::endl;
//         else
//             std::cout << "Using the default parameter value." << std::endl;
    }

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

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

    /*!
     * \brief Returns the temperature \f$ K \f$
     */
    Scalar temperature() const
    {
        return 273.15 + 30; // [K]
    }

    // \}

    /*!
     * \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;
        // set the left of the domain (with the global position in "0 = x" direction as a Dirichlet boundary
        if (globalPos[0] < eps_)
            bcTypes.setAllDirichlet();
        // set all other as Neumann boundaries
        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
    {
        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.
     */
    PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
    {
        // initialize values to zero, i.e. no-flow Neumann boundary conditions
        PrimaryVariables values(0.0);

        // if we are inside the injection zone set inflow Neumann boundary conditions
        // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
        // than using <= or >= as it is robust with regard to imprecision introduced by rounding errors.
        if (time_ < injectionDuration_
            && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->fvGridGeometry().bBoxMax()[0])
        {
            // inject nitrogen. negative values mean injection
            // units kg/(s*m^2)
            //TODO: Task 2: incorporate "totalAreaSpecificInflow_"  into the injection boundary condition
            values[Indices::conti0EqIdx + FluidSystem::N2Idx] = totalAreaSpecificInflow_/FluidSystem::molarMass(FluidSystem::N2Idx);
            values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.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
     */
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
    {
        PrimaryVariables values(0.0);

        // get the water density at atmospheric conditions
        const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);

        // assume an intially hydrostatic liquid pressure profile
        // note: we subtract rho_w*g*h because g is defined negative
        const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);

        values[Indices::pressureIdx] = pw;
        values[Indices::saturationIdx] = 0.0;

        return values;
    }

    // \}

    //! set the time for the time dependent boundary conditions (called from main)
    void setTime(Scalar time)
    { time_ = time; }

private:
    static constexpr Scalar eps_ = 1e-6;
    std::string name_; //! Problem name
    Scalar aquiferDepth_; //! Depth of the aquifer in m
    Scalar injectionDuration_; //! Duration of the injection in seconds
    //TODO: Task 2: Set a variable "totalAreaSpecificInflow_" to read in a value from the parameter tree via the input file
    Scalar totalAreaSpecificInflow_; //! Rate of the Injection in kg/(s m^2)
    Scalar time_;
};

} //end namespace Dumux

#endif