1p/sequential/diffusion/cellcentered/pressure.hh

Go to the documentation of this file.

// -*- 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 3 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/>.   *
 *****************************************************************************/
#ifndef DUMUX_FVPRESSURE1P_HH
#define DUMUX_FVPRESSURE1P_HH
 
 
// dumux environment
#include <dumux/porousmediumflow/sequential/cellcentered/pressure.hh>
#include <dumux/porousmediumflow/1p/sequential/properties.hh>
 
namespace Dumux {
 
template<class TypeTag> class FVPressure1P: public FVPressure<TypeTag>
{
    using ParentType = FVPressure<TypeTag>;
 
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
 
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
    using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
    using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
    using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
    using CellData = GetPropType<TypeTag, Properties::CellData>;
    using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
 
    enum
    {
        dim = GridView::dimension, dimWorld = GridView::dimensionworld
    };
 
    enum
    {
        pressEqIdx = Indices::pressureEqIdx // only one equation!
    };
 
    enum
    {
        rhs = ParentType::rhs, matrix = ParentType::matrix
    };
 
    using Element = typename GridView::Traits::template Codim<0>::Entity;
    using Intersection = typename GridView::Intersection;
 
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
    using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
    using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
 
 
public:
    // Function which calculates the source entry
    void getSource(Dune::FieldVector<Scalar, 2>&, const Element&, const CellData&, const bool);
    // Function which calculates the storage entry
    void getStorage(Dune::FieldVector<Scalar, 2>& entry, const Element& element,
                    const CellData& cellData, const bool first)
    {
        entry = 0;
    }
    // Function which calculates the flux entry
    void getFlux(Dune::FieldVector<Scalar, 2>&, const Intersection&, const CellData&, const bool);
    // Function which calculates the boundary flux entry
    void getFluxOnBoundary(Dune::FieldVector<Scalar, 2>&,
    const Intersection&, const CellData&, const bool);
 
    void initialize(bool solveTwice = true)
    {
        ParentType::initialize();
        this->assemble(true);
        this->solve();
        if (solveTwice)
        {
            this->assemble(false);
            this-> solve();
        }
        storePressureSolution();
        return;
    }
 
    void update()
    {
        ParentType::update();
        storePressureSolution();
    }
 
    void storePressureSolution()
    {
        int size = problem_.gridView().size(0);
        for (int i = 0; i < size; i++)
        {
            CellData& cellData = problem_.variables().cellData(i);
            storePressureSolution(i, cellData);
            cellData.fluxData().resetVelocity();
        }
    }
 
    void storePressureSolution(int eIdxGlobal, CellData& cellData)
    {
            Scalar press = this->pressure()[eIdxGlobal];
 
            cellData.setPressure(press);
    }
 
    template<class MultiWriter>
    void addOutputVtkFields(MultiWriter &writer)
    {
        ScalarSolutionType *pressure = writer.allocateManagedBuffer (
                problem_.gridView().size(0));
 
        *pressure = this->pressure();
 
        writer.attachCellData(*pressure, "pressure");
 
        return;
    }
 
    FVPressure1P(Problem& problem) :
        ParentType(problem), problem_(problem),
        gravity_(
        problem.gravity())
    {
        auto element = *problem_.gridView().template begin<0> ();
        Scalar temperature = problem_.temperature(element);
        Scalar referencePress = problem_.referencePressure(element);
 
        density_ = FluidSystem::density(temperature, referencePress);
        viscosity_ = FluidSystem::viscosity(temperature, referencePress);
    }
 
private:
    Problem& problem_;
    const GravityVector& gravity_; 
    Scalar density_;
    Scalar viscosity_;
};
 
template<class TypeTag>
void FVPressure1P<TypeTag>::getSource(Dune::FieldVector<Scalar, 2>& entry, const Element& element
        , const CellData& cellData, const bool first)
{
    // cell volume, assume linear map here
    Scalar volume = element.geometry().volume();
 
    // get sources from problem
    PrimaryVariables sourcePhase(0.0);
    problem_.source(sourcePhase, element);
        sourcePhase /= density_;
 
    entry[rhs] = volume * sourcePhase;
 
    return;
}
 
template<class TypeTag>
void FVPressure1P<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entry, const Intersection& intersection
        , const CellData& cellData, const bool first)
{
    auto elementI = intersection.inside();
    auto elementJ = intersection.outside();
 
    // get global coordinates of cell centers
    const GlobalPosition& globalPosI = elementI.geometry().center();
    const GlobalPosition& globalPosJ = elementJ.geometry().center();
 
    //get face normal
    const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
 
    // get face area
    Scalar faceArea = intersection.geometry().volume();
 
    // distance vector between barycenters
    GlobalPosition distVec = globalPosJ - globalPosI;
 
    // compute distance between cell centers
    Scalar dist = distVec.two_norm();
 
    // compute vectorized permeabilities
    DimMatrix meanPermeability(0);
 
    problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
            problem_.spatialParams().intrinsicPermeability(elementJ));
 
    Dune::FieldVector<Scalar, dim> permeability(0);
    meanPermeability.mv(unitOuterNormal, permeability);
 
    permeability/=viscosity_;
 
    //calculate current matrix entry
    entry[matrix] = ((permeability * unitOuterNormal) / dist) * faceArea;
 
    //calculate right hand side
    entry[rhs] = density_ * (permeability * gravity_) * faceArea;
 
    return;
}
 
template<class TypeTag>
void FVPressure1P<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& entry,
const Intersection& intersection, const CellData& cellData, const bool first)
{
    auto element = intersection.inside();
 
    // get global coordinates of cell centers
    const GlobalPosition& globalPosI = element.geometry().center();
 
    // center of face in global coordinates
    const GlobalPosition& globalPosJ = intersection.geometry().center();
 
    //get face normal
    const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
 
    // get face area
    Scalar faceArea = intersection.geometry().volume();
 
    // distance vector between barycenters
    GlobalPosition distVec = globalPosJ - globalPosI;
 
    // compute distance between cell centers
    Scalar dist = distVec.two_norm();
 
    BoundaryTypes bcType;
    problem_.boundaryTypes(bcType, intersection);
    PrimaryVariables boundValues(0.0);
 
    if (bcType.isDirichlet(pressEqIdx))
    {
        problem_.dirichlet(boundValues, intersection);
 
        // permeability vector at boundary
        // compute vectorized permeabilities
        DimMatrix meanPermeability(0);
 
        problem_.spatialParams().meanK(meanPermeability,
                problem_.spatialParams().intrinsicPermeability(element));
 
        Dune::FieldVector<Scalar, dim> permeability(0);
        meanPermeability.mv(unitOuterNormal, permeability);
 
        permeability/= viscosity_;
 
        // get Dirichlet pressure boundary condition
        Scalar pressBound = boundValues;
 
        // calculate current matrix entry
        entry[matrix] = ((permeability * unitOuterNormal) / dist) * faceArea;
        entry[rhs] = entry[matrix] * pressBound;
 
        // calculate right hand side
        entry[rhs] -= density_ * (permeability * gravity_) * faceArea;
 
    }
    // set Neumann boundary condition
    else if (bcType.isNeumann(pressEqIdx))
    {
        problem_.neumann(boundValues, intersection);
        Scalar J = boundValues /= density_;
 
        entry[rhs] = -(J * faceArea);
    }
    else
    {
        DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
    }
 
    return;
}
 
} // end namespace Dumux
#endif