Commit 35f692be authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[cornerpoint] add capability to use cornerpoint grids

Introduce a grid creator that reads from a Petrel output / Eclipse input
file and generates a CpGrid that is offered by the OPM module
dune-cornerpoint. Enhance the fully-implicit cell-centered models to be
able to deal with cornerpoint grids, test the functionality. This patch
consists of several parts:

- patches/ and dune.module: the CMake buildsystems from Dune and OPM are
  not compatible. The patches change some CMake files in the required
  OPM modules opm-parser, opm-core and dune-cornerpoint in such a way
  that those modules can serve as a dependency for Dumux. In addition,
  some specific CMake options have to be passed to dunecontrol. See
  patches/README for details.

- dumux/io/: defines the CpGridCreator which behaves like any other grid
  creator that reads from a file. In addition, it offers the Eclipse
  input deck in form of a static function deck(). This deck can be used
  to extract further simulation parameters, most prominently porosity
  and permeability values.

- dumux/implicit/cornerpoint/: changes some of the cell-centered
  classes to be able to deal with a CpGrid. The changes in the
  ElementVolumeVariables and the FVElementGeometry are minor and
  required by the facts that a CpGrid doesn't offer codim-1/2 entities
  and that the geometries don't offer local-to-global mappings. The
  changes in the DarcyFluxVariables are made for enhancing the
  robustness of the two-point flux approximation by forcing the
  inter-element transmissibilities to be positive.

- test/implicit/2p: runs a fully-implicit two-phase immiscible model on
  a CpGrid. An element index with an injection rate can be given in the
  input file. The grid provided in grids/ consists of only four
  elements, but exhibits typical cornerpoint features like a
  nonconforming situation and general hexahedrons. A more realistic grid
  will be provided in dumux-lecture soon.

Reviewed by gruenich.



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15183 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 82f99e83
// -*- 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 This file contains the data which is required to calculate
* volume fluxes of fluid phases over a face of a finite volume by means
* of the Darcy approximation.
*
*/
#ifndef DUMUX_CP_DARCY_FLUX_VARIABLES_HH
#define DUMUX_CP_DARCY_FLUX_VARIABLES_HH
#include <dune/common/float_cmp.hh>
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/implicit/common/implicitproperties.hh>
namespace Dumux
{
namespace Properties
{
// forward declaration of properties
NEW_PROP_TAG(ImplicitMobilityUpwindWeight);
NEW_PROP_TAG(SpatialParams);
NEW_PROP_TAG(NumPhases);
NEW_PROP_TAG(ProblemEnableGravity);
}
/*!
* \ingroup ImplicitFluxVariables
* \brief Evaluates the normal component of the Darcy velocity
* on a (sub)control volume face.
*/
template <class TypeTag>
class CpDarcyFluxVariables
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
enum { dim = GridView::dimension} ;
enum { dimWorld = GridView::dimensionworld} ;
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
public:
/*!
* \brief The constructor
*
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param fIdx The local index of the SCV (sub-control-volume) face
* \param elemVolVars The volume variables of the current element
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
CpDarcyFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int fIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: fvGeometry_(fvGeometry), faceIdx_(fIdx), onBoundary_(onBoundary)
{
mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
calculateVolumeFlux_(problem, element, elemVolVars);
}
public:
/*!
* \brief Return the volumetric flux over a face of a given phase.
*
* This is the calculated velocity multiplied by the unit normal
* and the area of the face. face().normal has already the
* magnitude of the area.
*
* \param phaseIdx index of the phase
*/
Scalar volumeFlux(const unsigned int phaseIdx) const
{ return volumeFlux_[phaseIdx]; }
/*!
* \brief Return the velocity of a given phase.
*
* This is the full velocity vector on the
* face (without being multiplied with normal).
*
* \param phaseIdx index of the phase
*/
GlobalPosition velocity(const unsigned int phaseIdx) const
{ return velocity_[phaseIdx] ; }
/*!
* \brief Return the local index of the downstream control volume
* for a given phase.
*
* \param phaseIdx index of the phase
*/
const unsigned int downstreamIdx(const unsigned phaseIdx) const
{ return downstreamIdx_[phaseIdx]; }
/*!
* \brief Return the local index of the upstream control volume
* for a given phase.
*
* \param phaseIdx index of the phase
*/
const unsigned int upstreamIdx(const unsigned phaseIdx) const
{ return upstreamIdx_[phaseIdx]; }
/*!
* \brief Return the SCV (sub-control-volume) face. This may be either
* a face within the element or a face on the element boundary,
* depending on the value of onBoundary_.
*/
const SCVFace &face() const
{
if (onBoundary_)
return fvGeometry_.boundaryFace[faceIdx_];
else
return fvGeometry_.subContVolFace[faceIdx_];
}
protected:
/*!
* \brief Actual calculation of the volume flux.
*
* \param problem The problem
* \param element The finite element
* \param elemVolVars The volume variables of the current element
*/
void calculateVolumeFlux_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
// calculate the transmissibilities
const SpatialParams &spatialParams = problem.spatialParams();
const Element& elementI = *fvGeometry_.neighbors[face().i];
FVElementGeometry fvGeometryI;
fvGeometryI.subContVol[0].global = elementI.geometry().center();
auto ki = spatialParams.intrinsicPermeability(elementI, fvGeometryI, 0);
Dune::FieldVector<Scalar, dimWorld> kin;
ki.mv(face().normal, kin);
kin /= face().area;
auto di = face().ipGlobal;
di -= elementI.geometry().center();
auto ti = std::abs(di*kin*face().area/(2*di.two_norm2()));
auto tij = ti;
if (!onBoundary_)
{
const Element& elementJ = *fvGeometry_.neighbors[face().j];
FVElementGeometry fvGeometryJ;
fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
auto kj = spatialParams.intrinsicPermeability(elementJ, fvGeometryJ, 0);
Dune::FieldVector<Scalar, dimWorld> kjn;
kj.mv(face().normal, kjn);
kjn /= face().area;
auto dj = face().ipGlobal;
dj -= elementJ.geometry().center();
auto tj = std::abs(dj*kjn*face().area/(2*dj.two_norm2()));
tij = harmonicMean(ti, tj);
}
// loop over all phases
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
{
const auto& volVarsI = elemVolVars[face().i];
auto potentialI = volVarsI.pressure(phaseIdx);
const auto& volVarsJ = elemVolVars[face().j];
auto potentialJ = volVarsJ.pressure(phaseIdx);
if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
{
// calculate the phase density at the integration point. we
// only do this if the phase is present in both cells
Scalar SI = volVarsI.fluidState().saturation(phaseIdx);
Scalar SJ = volVarsJ.fluidState().saturation(phaseIdx);
Scalar rhoI = volVarsI.fluidState().density(phaseIdx);
Scalar rhoJ = volVarsJ.fluidState().density(phaseIdx);
Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(fI + fJ, 0.0, 1.0e-30))
// doesn't matter because no wetting phase is present in
// both cells!
fI = fJ = 0.5;
const Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
auto globalPosI = elementI.geometry().center();
potentialI -= density*(problem.gravityAtPos(globalPosI)*globalPosI);
if (onBoundary_)
{
potentialJ -= density*(problem.gravityAtPos(face().ipGlobal)*face().ipGlobal);
}
else
{
const Element& elementJ = *fvGeometry_.neighbors[face().j];
auto globalPosJ = elementJ.geometry().center();
potentialJ -= density*(problem.gravityAtPos(globalPosJ)*globalPosJ);
}
}
auto potentialDiff = potentialI - potentialJ;
// determine the upwind direction
if (potentialDiff > 0)
{
upstreamIdx_[phaseIdx] = face().i;
downstreamIdx_[phaseIdx] = face().j;
}
else
{
upstreamIdx_[phaseIdx] = face().j;
downstreamIdx_[phaseIdx] = face().i;
}
// obtain the upwind volume variables
const VolumeVariables& upVolVars = elemVolVars[ upstreamIdx(phaseIdx) ];
const VolumeVariables& downVolVars = elemVolVars[ downstreamIdx(phaseIdx) ];
// set the volume flux
volumeFlux_[phaseIdx] = tij*potentialDiff;
volumeFlux_[phaseIdx] *= mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
+ (1.0 - mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx);
velocity_[phaseIdx] = face().normal;
velocity_[phaseIdx] *= volumeFlux_[phaseIdx]/face().area;
} // over loop all phases
}
const FVElementGeometry &fvGeometry_; //!< Information about the geometry of discretization
const unsigned int faceIdx_; //!< The index of the sub control volume face
const bool onBoundary_; //!< Specifying whether we are currently on the boundary of the simulation domain
unsigned int upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
Scalar volumeFlux_[numPhases] ; //!< Velocity multiplied with normal (magnitude=area)
GlobalPosition velocity_[numPhases] ; //!< The velocity as determined by Darcy's law or by the Forchheimer relation
Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
};
} // end namespace
#endif // DUMUX_CP_DARCY_FLUX_VARIABLES_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 Volume variables gathered on an element
*/
#ifndef DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH
#define DUMUX_CP_ELEMENT_VOLUME_VARIABLES_HH
#include <dune/common/version.hh>
#include <dumux/implicit/cellcentered/ccproperties.hh>
namespace Dumux
{
/*!
* \ingroup CCModel
* \brief This class stores an array of VolumeVariables objects, one
* volume variables object for each of the element's vertices
*/
template<class TypeTag>
class CpElementVolumeVariables : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables) >
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
public:
/*!
* \brief The constructor.
*/
CpElementVolumeVariables()
{ }
/*!
* \brief Construct the volume variables for all of vertices of an element.
*
* \param problem The problem which needs to be simulated.
* \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
* \param fvGeometry The finite volume geometry of the element
* \param oldSol Tells whether the model's previous or current solution should be used.
*/
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
bool oldSol)
{
const SolutionVector &globalSol =
oldSol?
problem.model().prevSol():
problem.model().curSol();
int numNeighbors = fvGeometry.numNeighbors;
this->resize(numNeighbors);
for (int i = 0; i < numNeighbors; i++)
{
const Element& neighbor = *(fvGeometry.neighbors[i]);
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
const PrimaryVariables &solI
= globalSol[problem.elementMapper().index(neighbor)];
#else
const PrimaryVariables &solI
= globalSol[problem.elementMapper().map(neighbor)];
#endif
FVElementGeometry neighborFVGeom;
neighborFVGeom.updateInner(neighbor);
(*this)[i].update(solI,
problem,
neighbor,
neighborFVGeom,
/*scvIdx=*/0,
oldSol);
}
// only treat boundary if current solution is evaluated
if (!oldSol)
{
// check if element intersects with the boundary
ElementBoundaryTypes elemBCTypes;
elemBCTypes.update(problem, element);
if (elemBCTypes.hasDirichlet()
|| elemBCTypes.hasNeumann()
|| elemBCTypes.hasOutflow())
{
const int numFaces = 6;
this->resize(numNeighbors + numFaces);
// add volume variables for the boundary faces
IntersectionIterator isIt = problem.gridView().ibegin(element);
IntersectionIterator isEndIt = problem.gridView().iend(element);
for (; isIt != isEndIt; ++isIt) {
if (!isIt->boundary())
continue;
BoundaryTypes bcTypes;
problem.boundaryTypes(bcTypes, *isIt);
int fIdx = isIt->indexInInside();
int indexInVariables = numNeighbors + fIdx;
if (bcTypes.hasDirichlet())
{
PrimaryVariables dirichletValues;
problem.dirichlet(dirichletValues, *isIt);
(*this)[indexInVariables].update(dirichletValues,
problem,
element,
fvGeometry,
/*scvIdx=*/0,
oldSol);
}
else
{
(*this)[indexInVariables] = (*this)[0];
}
}
}
}
}
};
} // namespace Dumux
#endif
// -*- 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 Represents the finite volume geometry of a single element in
* the cell-centered fv scheme.
*/
#ifndef DUMUX_CP_FV_ELEMENTGEOMETRY_HH
#define DUMUX_CP_FV_ELEMENTGEOMETRY_HH
#include <dune/common/version.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/intersectioniterator.hh>
#include <dumux/common/propertysystem.hh>
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(Scalar);
}
/*!
* \ingroup CCModel
* \brief Represents the finite volume geometry of a single element in
* the cell-centered fv scheme.
*/
template<class TypeTag>
class CpFVElementGeometry
{
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum{dim = GridView::dimension};
enum{dimWorld = GridView::dimensionworld};
enum{maxNFAP = 2}; //! maximum number of flux approximation points (two-point flux)
enum{maxNE = 50}; //! maximum number of neighbors
enum{maxBF = 2*dim}; //! maximum number of boundary faces
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GridView::ctype CoordScalar;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
typedef typename Element::Geometry Geometry;
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
typedef typename GridView::IntersectionIterator IntersectionIterator;
public:
struct SubControlVolume //! FV intersected with element
{
LocalPosition local; //!< local position
GlobalPosition global; //!< global position
Scalar volume; //!< volume of scv
bool inner;
};
struct SubControlVolumeFace //! interior face of a sub control volume
{
int i,j; //!< scvf seperates corner i and j of elem
LocalPosition ipLocal; //!< integration point in local coords
GlobalPosition ipGlobal; //!< integration point in global coords
GlobalPosition normal; //!< normal on face pointing to CV j or outward of the domain with length equal to |scvf|
Scalar area; //!< area of face
Dune::FieldVector<GlobalPosition, maxNFAP> grad; //!< derivatives of shape functions at ip
Dune::FieldVector<Scalar, maxNFAP> shapeValue; //!< value of shape functions at ip
Dune::FieldVector<int, maxNFAP> fapIndices; //!< indices w.r.t.neighbors of the flux approximation points
unsigned numFap; //!< number of flux approximation points
unsigned fIdx; //!< the index (w.r.t. the element) of the face (codim 1 entity) that the scvf is part of
};
typedef SubControlVolumeFace BoundaryFace; //!< compatibility typedef
LocalPosition elementLocal; //!< local coordinate of element center
GlobalPosition elementGlobal; //!< global coordinate of element center
Scalar elementVolume; //!< element volume
SubControlVolume subContVol[1]; //!< data of the sub control volumes
SubControlVolumeFace subContVolFace[maxNE]; //!< data of the sub control volume faces
BoundaryFace boundaryFace[maxBF]; //!< data of the boundary faces
int numScv; //!< number of subcontrol volumes
int numScvf; //!< number of inner-domain subcontrolvolume faces
int numNeighbors; //!< number of neighboring elements including the element itself
std::vector<ElementPointer> neighbors; //!< stores pointers for the neighboring elements
void updateInner(const Element& element)
{
const Geometry geometry = element.geometry();
elementVolume = geometry.volume();
elementGlobal = geometry.center();
//elementLocal = geometry.local(elementGlobal);
numScv = 1;
numScvf = 0;
subContVol[0].local = elementLocal;
subContVol[0].global = elementGlobal;
subContVol[0].inner = true;
subContVol[0].volume = elementVolume;
// initialize neighbors list with self:
numNeighbors = 1;
neighbors.clear();
neighbors.reserve(maxNE);
ElementPointer elementPointer(element);
neighbors.push_back(elementPointer);
}
void update(const GridView& gridView, const Element& element)
{
updateInner(element);
const Geometry geometry = element.geometry();
bool onBoundary = false;
// fill neighbor information and control volume face data:
IntersectionIterator isEndIt = gridView.iend(element);
for (IntersectionIterator isIt = gridView.ibegin(element); isIt != isEndIt; ++isIt)
{
const auto isGeometry = isIt->geometry();
// neighbor information and inner cvf data:
if (isIt->neighbor())
{
numNeighbors++;
ElementPointer elementPointer(isIt->outside());
neighbors.push_back(elementPointer);
int scvfIdx = numNeighbors - 2;
SubControlVolumeFace& scvFace = subContVolFace[scvfIdx];
scvFace.i = 0;
scvFace.j = scvfIdx + 1;
scvFace.ipGlobal = isGeometry.center();
//scvFace.ipLocal = geometry.local(scvFace.ipGlobal);
scvFace.normal = isIt->centerUnitOuterNormal();
auto di = scvFace.ipGlobal;
di -= elementGlobal;
if (scvFace.normal*di < 0)