2pproblem.hh

// -*- 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 twophase box model.
 */
#ifndef DUMUX_EXERCISE_FLUIDSYSTEM_A_PROBLEM_HH
#define DUMUX_EXERCISE_FLUIDSYSTEM_A_PROBLEM_HH

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

// The water component
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>

// The interface to create plots during simulation
#include <dumux/io/gnuplotinterface.hh>

namespace Dumux {

/*!
 * \ingroup TwoPBoxModel
 * \brief  Tutorial problem for a fully coupled twophase box model.
 */
template <class TypeTag>
class ExerciseFluidsystemProblemTwoP : 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 FluidState = GetPropType<TypeTag, Properties::FluidState>;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;

    enum {
        waterPressureIdx = Indices::pressureIdx,
        naplSaturationIdx = Indices::saturationIdx,
        contiWEqIdx = Indices::conti0EqIdx + FluidSystem::comp0Idx, // water transport equation index
        contiNEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx // napl transport equation index
    };

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

        // initialize the tables for the water properties
        Components::TabulatedComponent<Components::H2O<Scalar>>::init(/*tempMin=*/273.15,
                                                                      /*tempMax=*/623.15,
                                                                      /*numTemp=*/100,
                                                                      /*pMin=*/0.0,
                                                                      /*pMax=*/20e6,
                                                                      /*numP=*/200);

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

        // plot density over pressure of the phase consisting of your component
        if(getParam<bool>("Output.PlotDensity"))
            plotDensity_();
    }

    /*!
     * \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 the initial values as Dirichlet values
        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.
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
     */
    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];

        // influx of oil (30 g/m/s) over a segment of the top boundary
        if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40)
        {
            values[contiWEqIdx] = 0;
            values[contiNEqIdx] = -3e-2;
        }
        else
        {
            // no-flow on the remaining Neumann-boundaries.
            values[contiWEqIdx] = 0;
            values[contiNEqIdx] = 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);

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

        return values;
    }
    // \}

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

private:
    void plotDensity_()
    {
        FluidState fluidState;
        fluidState.setTemperature(temperature());
        int numberOfPoints = 100;
        Scalar xMin = 1e4;
        Scalar xMax = 1e7;
        Scalar spacing = std::pow((xMax/xMin), 1.0/(numberOfPoints-1));
        std::vector<double> x(numberOfPoints);
        std::vector<double> y(numberOfPoints);
        for (int i=0; i<numberOfPoints; ++i)
            x[i] = xMin*std::pow(spacing,i);
        for (int i=0; i<numberOfPoints; ++i)
        {
            fluidState.setPressure(FluidSystem::phase1Idx, x[i]);
            y[i] = FluidSystem::density(fluidState, FluidSystem::phase1Idx);
        }
        gnuplot_.resetPlot();
        gnuplot_.setXRange(xMin, xMax);
        gnuplot_.setOption("set logscale x 10");
        gnuplot_.setOption("set key left top");
        gnuplot_.setYRange(1440, 1480);
        gnuplot_.setXlabel("pressure [Pa]");
        gnuplot_.setYlabel("density [kg/m^3]");
        gnuplot_.addDataSetToPlot(x, y, "YourComponentPhase_density.dat", "w l t 'Phase consisting of your component'");
        gnuplot_.plot("YourComponentPhase_density");
    }

    Scalar eps_; //! small epsilon value
    Scalar depthBOR_; //! depth at the bottom of the reservoir
    Dumux::GnuplotInterface<double> gnuplot_; //! collects data for plotting
};

} // end namespace Dumux

#endif