Commit 04e7ec8f authored by Timo Koch's avatar Timo Koch Committed by Dennis Gläser
Browse files

[next] Implement fvproblem. Move geometry stuff to fvGridGeometry.

parent 844b8639
This diff is collapsed.
......@@ -38,6 +38,7 @@ namespace Properties
{
// Property forward declarations
NEW_PROP_TAG(ElementVolumeVariables);
NEW_PROP_TAG(FVGridGeometry);
NEW_PROP_TAG(FVElementGeometry);
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(ImplicitIsBox);
......@@ -346,26 +347,25 @@ private:
template<class TypeTag>
class BoundingBoxTreePointSourceHelper
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PointSource) PointSource;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using PointSource = typename GET_PROP_TYPE(TypeTag, PointSource);
static const int dim = GridView::dimension;
static const int dimworld = GridView::dimensionworld;
typedef Dumux::BoundingBoxTree<GridView> BoundingBoxTree;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
public:
//! calculate a DOF index to point source map from given vector of point sources
static void computePointSourceMap(const Problem& problem,
const BoundingBoxTree& boundingBoxTree,
static void computePointSourceMap(const FVGridGeometry& fvGridGeometry,
std::vector<PointSource>& sources,
std::map<std::pair<unsigned int, unsigned int>, std::vector<PointSource> >& pointSourceMap)
{
const auto& boundingBoxTree = fvGridGeometry.boundingBoxTree();
for (auto&& source : sources)
{
// compute in which elements the point source falls
......@@ -378,9 +378,8 @@ public:
if(isBox)
{
// check in which subcontrolvolume(s) we are
// TODO mapper/problem in bboxtree would allow to make this much better
const auto element = boundingBoxTree.entity(eIdx);
auto fvGeometry = localView(problem.model().fvGridGeometry());
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
const auto globalPos = source.position();
......
......@@ -25,7 +25,10 @@
#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_GRID_GEOMETRY_HH
#define DUMUX_DISCRETIZATION_CCTPFA_FV_GRID_GEOMETRY_HH
#include <dune/common/version.hh>
#include <dumux/common/elementmap.hh>
#include <dumux/common/boundingboxtree.hh>
#include <dumux/implicit/cellcentered/tpfa/properties.hh>
#include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
#include <dumux/discretization/cellcentered/connectivitymap.hh>
......@@ -58,10 +61,12 @@ class CCTpfaFVGridGeometry<TypeTag, true>
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Element = typename GridView::template Codim<0>::Entity;
using ConnectivityMap = CCSimpleConnectivityMap<TypeTag>;
using BoundingBoxTree = Dumux::BoundingBoxTree<GridView>;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using CoordScalar = typename GridView::ctype;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
//! The local class needs access to the scv, scvfs and the fv element geometry
//! as they are globally cached
......@@ -72,9 +77,38 @@ public:
CCTpfaFVGridGeometry(const GridView& gridView)
: gridView_(gridView)
, elementMap_(gridView)
#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
, elementMapper_(gridView, Dune::mcmgElementLayout())
, vertexMapper_(gridView, Dune::mcmgVertexLayout())
{}
#else
, elementMapper_(gridView)
, vertexMapper_(gridView)
#endif
, bBoxMin_(std::numeric_limits<double>::max())
, bBoxMax_(-std::numeric_limits<double>::max())
{
// calculate the bounding box of the local partition of the grid view
for (const auto& vertex : vertices(gridView))
{
for (int i=0; i<dimWorld; i++)
{
using std::min;
using std::max;
bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().corner(0)[i]);
bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().corner(0)[i]);
}
}
// communicate to get the bounding box of the whole domain
if (gridView.comm().size() > 1)
{
for (int i = 0; i < dimWorld; ++i)
{
bBoxMin_[i] = gridView.comm().min(bBoxMin_[i]);
bBoxMax_[i] = gridView.comm().max(bBoxMax_[i]);
}
}
}
//! The total number of sub control volumes
std::size_t numScv() const
......@@ -303,6 +337,42 @@ public:
const ConnectivityMap &connectivityMap() const
{ return connectivityMap_; }
/*!
* \brief Returns the bounding box tree of the grid
*/
BoundingBoxTree& boundingBoxTree()
{
if(!boundingBoxTree_)
boundingBoxTree_ = std::make_unique<BoundingBoxTree>(gridView_);
return *boundingBoxTree_;
}
/*!
* \brief Returns the bounding box tree of the grid
*/
const BoundingBoxTree& boundingBoxTree() const
{
if(!boundingBoxTree_)
boundingBoxTree_ = std::make_unique<BoundingBoxTree>(gridView_);
return *boundingBoxTree_;
}
/*!
* \brief The coordinate of the corner of the GridView's bounding
* box with the smallest values.
*/
const GlobalPosition &bBoxMin() const
{ return bBoxMin_; }
/*!
* \brief The coordinate of the corner of the GridView's bounding
* box with the largest values.
*/
const GlobalPosition &bBoxMax() const
{ return bBoxMax_; }
private:
// find the scvf that has insideScvIdx in its outsideScvIdx list and outsideScvIdx as its insideScvIdx
IndexType findFlippedScvfIndex_(IndexType insideScvIdx, IndexType outsideScvIdx)
......@@ -332,8 +402,16 @@ private:
std::vector<SubControlVolumeFace> scvfs_;
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
IndexType numBoundaryScvf_;
// needed for embedded surface and network grids (dim < dimWorld)
std::vector<std::vector<IndexType>> flipScvfIndices_;
// The bounding box tree of the grid for effecient element intersections
mutable std::unique_ptr<BoundingBoxTree> boundingBoxTree_;
// The bounding box of the global grid
GlobalPosition bBoxMin_;
GlobalPosition bBoxMax_;
};
// specialization in case the FVElementGeometries are not stored
......@@ -350,10 +428,14 @@ class CCTpfaFVGridGeometry<TypeTag, false>
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Element = typename GridView::template Codim<0>::Entity;
using ConnectivityMap = CCSimpleConnectivityMap<TypeTag>;
using BoundingBoxTree = Dumux::BoundingBoxTree<GridView>;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
using CoordScalar = typename GridView::ctype;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
// TODO is this necessary?
friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
......@@ -362,9 +444,38 @@ public:
CCTpfaFVGridGeometry(const GridView& gridView)
: gridView_(gridView)
, elementMap_(gridView)
#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
, elementMapper_(gridView, Dune::mcmgElementLayout())
, vertexMapper_(gridView, Dune::mcmgVertexLayout())
{}
#else
, elementMapper_(gridView)
, vertexMapper_(gridView)
#endif
, bBoxMin_(std::numeric_limits<double>::max())
, bBoxMax_(-std::numeric_limits<double>::max())
{
// calculate the bounding box of the local partition of the grid view
for (const auto& vertex : vertices(gridView))
{
for (int i=0; i<dimWorld; i++)
{
using std::min;
using std::max;
bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().corner(0)[i]);
bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().corner(0)[i]);
}
}
// communicate to get the bounding box of the whole domain
if (gridView.comm().size() > 1)
{
for (int i = 0; i < dimWorld; ++i)
{
bBoxMin_[i] = gridView.comm().min(bBoxMin_[i]);
bBoxMax_[i] = gridView.comm().max(bBoxMax_[i]);
}
}
}
//! The total number of sub control volumes
std::size_t numScv() const
......@@ -538,6 +649,42 @@ public:
const ConnectivityMap &connectivityMap() const
{ return connectivityMap_; }
/*!
* \brief Returns the bounding box tree of the grid
*/
BoundingBoxTree& boundingBoxTree()
{
if(!boundingBoxTree_)
boundingBoxTree_ = std::make_unique<BoundingBoxTree>(gridView_);
return *boundingBoxTree_;
}
/*!
* \brief Returns the bounding box tree of the grid
*/
const BoundingBoxTree& boundingBoxTree() const
{
if(!boundingBoxTree_)
boundingBoxTree_ = std::make_unique<BoundingBoxTree>(gridView_);
return *boundingBoxTree_;
}
/*!
* \brief The coordinate of the corner of the GridView's bounding
* box with the smallest values.
*/
const GlobalPosition &bBoxMin() const
{ return bBoxMin_; }
/*!
* \brief The coordinate of the corner of the GridView's bounding
* box with the largest values.
*/
const GlobalPosition &bBoxMax() const
{ return bBoxMax_; }
private:
const GridView gridView_;
......@@ -556,6 +703,13 @@ private:
// vectors that store the global data
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
std::vector<std::vector<std::vector<IndexType>>> neighborVolVarIndices_;
// The bounding box tree of the grid for effecient element intersections
mutable std::unique_ptr<BoundingBoxTree> boundingBoxTree_;
// The bounding box of the global grid
GlobalPosition bBoxMin_;
GlobalPosition bBoxMax_;
};
} // end namespace
......
......@@ -701,53 +701,6 @@ protected:
LocalResidual &localResidual_()
{ return localJacobian_.localResidual(); }
/*!
* \brief Applies the initial solution for all degrees of freedom of the grid.
*
*/
void applyInitialSolution_()
{
// first set the whole domain to zero
uCur_ = Scalar(0.0);
// set the initial values by forwarding to a specialized method
applyInitialSolutionImpl_(std::integral_constant<bool, isBox>());
// add up the primary variables which cross process borders
if (isBox && gridView_().comm().size() > 1)
{
VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
sumPVHandle(uCur_, vertexMapper());
gridView_().communicate(sumPVHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
}
}
/*!
* \brief Applies the initial solution for the box method
*/
void applyInitialSolutionImpl_(std::true_type)
{
for (const auto& vertex : vertices(problem_().gridView()))
{
const auto dofIdxGlobal = dofMapper().index(vertex);
uCur_[dofIdxGlobal] = problem_().initial(vertex);
}
}
/*!
* \brief Applies the initial solution for cell-centered methods
*/
void applyInitialSolutionImpl_(std::false_type)
{
for (const auto& element : elements(problem_().gridView()))
{
const auto dofIdxGlobal = dofMapper().index(element);
uCur_[dofIdxGlobal] = problem_().initial(element);
}
}
/*!
* \brief Find all indices of boundary vertices (box) / elements (cell centered).
*/
......
......@@ -71,11 +71,6 @@ public:
FVSpatialParamsOneP(const GridView &gridView)
{
}
~FVSpatialParamsOneP()
{
}
/*!
* \brief Averages the intrinsic permeability (Scalar).
* \param K1 intrinsic permeability of the first element
......
......@@ -62,8 +62,8 @@ class ImplicitSpatialParams: public ImplicitSpatialParamsOneP<TypeTag>
using GlobalPosition = Dune::FieldVector<CoordScalar,dimWorld>;
public:
ImplicitSpatialParams(const Problem &problem, const GridView &gridView)
: ImplicitSpatialParamsOneP<TypeTag>(problem, gridView)
ImplicitSpatialParams(const Problem& problem)
: ImplicitSpatialParamsOneP<TypeTag>(problem)
{}
/*!
......@@ -78,7 +78,7 @@ public:
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{
return asImp_().materialLawParamsAtPos(scv.center());
return this->asImp_().materialLawParamsAtPos(scv.center());
}
/*!
......@@ -93,13 +93,6 @@ public:
"The spatial parameters do not provide "
"a materialLawParamsAtPos() method.");
}
private:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
};
} // namespace Dumux
......
......@@ -68,7 +68,7 @@ class ImplicitSpatialParamsOneP
using GlobalPosition = Dune::FieldVector<CoordScalar,dimWorld>;
public:
ImplicitSpatialParamsOneP(const Problem& problem, const GridView &gridView)
ImplicitSpatialParamsOneP(const Problem& problem)
: problemPtr_(&problem)
{}
......@@ -306,7 +306,8 @@ public:
return forchCoeff;
}
const Problem& problem()
//! The problem we are associated with
const Problem& problem() const
{
return *problemPtr_;
}
......
// -*- 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 Base class for all porous media problems
*/
#ifndef DUMUX_POROUS_MEDIUM_FLOW_PROBLEM_HH
#define DUMUX_POROUS_MEDIUM_FLOW_PROBLEM_HH
#include <dumux/implicit/properties.hh>
#include <dumux/common/fvproblem.hh>
namespace Dumux
{
namespace Properties
{
//! Property forward declarations
NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters object
NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
}
/*!
* \ingroup ImplicitBaseProblems
* \brief Base class for all fully implicit porous media problems
* TODO: derive from base problem property?
*/
template<class TypeTag>
class PorousMediumFlowProblem : public FVProblem<TypeTag>
{
using ParentType = FVProblem<TypeTag>;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
using CoordScalar = typename GridView::ctype;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
/*!
* \brief The constructor
*
* \param fvGridGeometry The finite volume grid geometry
*/
PorousMediumFlowProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
, spatialParams_(std::make_shared<SpatialParams>(this->asImp_()))
{
// TODO: spatial params init?
if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
gravity_[dimWorld-1] = -9.81;
}
/*!
* \name Physical parameters for porous media problems
*/
// \{
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at a given global position.
*
* This is not specific to the discretization. By default it just
* calls temperature().
*
* \param globalPos The position in global coordinates where the temperature should be specified.
*/
Scalar temperatureAtPos(const GlobalPosition &globalPos) const
{ return this->asImp_().temperature(); }
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
*
* This is not specific to the discretization. By default it just
* throws an exception so it must be overloaded by the problem if
* no energy equation is used.
*/
Scalar temperature() const
{
DUNE_THROW(Dune::NotImplemented, "temperature() method not implemented by the user problem");
}
/*!
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
*
* This is discretization independent interface. By default it
* just calls gravity().
*/
const GlobalPosition &gravityAtPos(const GlobalPosition &pos) const
{ return this->asImp_().gravity(); }
/*!
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
*
* This method is used for problems where the gravitational
* acceleration does not depend on the spatial position. The
* default behaviour is that if the <tt>ProblemEnableGravity</tt>
* property is true, \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$ holds,
* else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$.
*/
const GlobalPosition &gravity() const
{ return gravity_; }
/*!
* \brief Returns the spatial parameters object.
*/
SpatialParams &spatialParams()
{ return *spatialParams_; }
/*!
* \brief Returns the spatial parameters object.
*/
const SpatialParams &spatialParams() const
{ return *spatialParams_; }
// \}
protected:
//! The gravity acceleration vector
GlobalPosition gravity_;
// material properties of the porous medium
std::shared_ptr<SpatialParams> spatialParams_;
};
} // end namespace Dumux
#endif
......@@ -23,6 +23,7 @@
#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
......@@ -34,11 +35,12 @@
#include <dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh>
#include <dumux/porousmediumflow/1p/implicit/propertydefaults.hh>
#include "spatialparams.hh"
namespace Dumux
{
// forward declarations
template<class TypeTag> class OnePTestProblem;
template<class TypeTag> class OnePTestSpatialParams;
namespace Properties
{
......@@ -73,7 +75,7 @@ SET_PROP(IncompressibleTestProblem, Fluid)
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = FluidSystems::LiquidPhase<Scalar, SimpleH2O<Scalar> >;
using type = FluidSystems::LiquidPhase<Scalar, TabulatedComponent<Scalar, H2O<Scalar>>>;
};
// Enable caching
......@@ -84,76 +86,48 @@ SET_BOOL_PROP(IncompressibleTestProblem, EnableFVGridGeometryCache, false);
} // end namespace Properties
template<class TypeTag>
class OnePTestSpatialParams
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
public:
using PermeabilityType = Scalar;
PermeabilityType permeability(const Element &element,
const SubControlVolume &scv,
const ElementSolutionVector &elemSol) const