porousmediumflow/2pve/spatialparams.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:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
 
#ifndef DUMUX_TWOPVE_SPATIAL_PARAMS_HH
#define DUMUX_TWOPVE_SPATIAL_PARAMS_HH
 
#include <cmath>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/grid/common/rangegenerators.hh>
 
#include <dumux/porousmediumflow/fvspatialparamsmp.hh>
#include <dumux/porousmediumflow/2pve/columnmapping.hh>
 
namespace Dumux {

template<class GridGeometry, class FineGridGeometry, class Scalar, class SpatialParamsFine, class Implementation>
class TwoPVESpatialParams
: public FVPorousMediumFlowSpatialParamsMP<GridGeometry, Scalar, Implementation>
{
    using ParentType = FVPorousMediumFlowSpatialParamsMP<GridGeometry, Scalar, Implementation>;
    using GridView = typename GridGeometry::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using SubControlVolume = typename GridGeometry::LocalView::SubControlVolume;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
    using ColumnMapping = TwoPVEColumnMapping<GridGeometry, FineGridGeometry, Scalar>;
 
public:
    using PermeabilityType = Scalar;

    TwoPVESpatialParams(std::shared_ptr<const GridGeometry> gridGeometry,
                        const ColumnMapping& columnMapping,
                        std::shared_ptr<const SpatialParamsFine> spatialParamsFine,
                        const Scalar fineCellHeight)
    : ParentType(gridGeometry)
    , spatialParamsFine_(spatialParamsFine)
    , permeabilityIntegral_(gridGeometry->elementMapper().size(), 0.0)
    , porosityIntegral_(gridGeometry->elementMapper().size(), 0.0)
    , columnHeight_(gridGeometry->elementMapper().size(), 0.0)
    {
        using FinePermeability = std::decay_t<decltype(spatialParamsFine_->permeabilityAtElement(std::declval<const typename FineGridGeometry::GridView::template Codim<0>::Entity&>()))>;
        static_assert(std::is_convertible_v<FinePermeability, Scalar>, "The VE upscaling requires a scalar fine-level permeability");
 
        for (const auto& element : elements(gridGeometry->gridView()))
        {
            const auto columnIdx = gridGeometry->elementMapper().index(element);
            const auto& column = columnMapping.column(columnIdx);
            const Scalar columnPorosity = spatialParamsFine_->porosityAtElement(column.front());
            for (const auto& fineElement : column)
            {
                using std::abs;
                const Scalar porosityFine = spatialParamsFine_->porosityAtElement(fineElement);
                // the gas plume distance follows from the water balance of the unweighted saturation profile
                if (abs(porosityFine - columnPorosity) > 1e-10*abs(columnPorosity))
                    DUNE_THROW(Dune::InvalidStateException, "The VE model requires a uniform porosity within each column, column " << columnIdx << " has porosities " << columnPorosity << " and " << porosityFine);
 
                permeabilityIntegral_[columnIdx] += spatialParamsFine_->permeabilityAtElement(fineElement)*fineCellHeight;
                porosityIntegral_[columnIdx] += porosityFine*fineCellHeight;
            }
 
            const Scalar columnHeight = column.size()*fineCellHeight;
            columnHeight_[columnIdx] = columnHeight;
        }
    }

    const SpatialParamsFine& spatialParamsFine() const
    { return *spatialParamsFine_; }

    template<class ElementSolution>
    PermeabilityType permeability(const Element& element,
                                  const SubControlVolume& scv,
                                  const ElementSolution& elemSol) const
    { return permeabilityIntegralAtElement(element)/columnHeightAtElement(element); }

    template<class ElementSolution>
    Scalar porosity(const Element& element,
                    const SubControlVolume& scv,
                    const ElementSolution& elemSol) const
    { return porosityIntegralAtElement(element)/columnHeightAtElement(element); }

    template<class ElementSolution>
    Scalar extrusionFactor(const Element& element,
                           const SubControlVolume&,
                           const ElementSolution&) const
    { return columnHeightAtElement(element); }

    template<class FluidSystem>
    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
    { return FluidSystem::phase0Idx; }

    Scalar permeabilityIntegralAtElement(const Element& element) const
    { return permeabilityIntegral_[this->gridGeometry().elementMapper().index(element)]; }

    Scalar porosityIntegralAtElement(const Element& element) const
    { return porosityIntegral_[this->gridGeometry().elementMapper().index(element)]; }

    Scalar columnHeightAtElement(const Element& element) const
    { return columnHeight_[this->gridGeometry().elementMapper().index(element)]; }

    GlobalPosition gravity(const GlobalPosition&) const
    { return GlobalPosition(0.0); }
 
private:
    std::shared_ptr<const SpatialParamsFine> spatialParamsFine_;
    std::vector<PermeabilityType> permeabilityIntegral_;
    std::vector<Scalar> porosityIntegral_;
    std::vector<Scalar> columnHeight_;
};
 
} // end namespace Dumux
 
#endif