OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

changed: ewoms/models/discretefracture -> opm/models/discretefacture

Browse Source
This commit is contained in:
Arne Morten Kvarving 2019-09-16 11:39:30 +02:00
parent 474ae4ded8
commit 16f4bdcf02
8 changed files with 1121 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

138
opm/models/discretefracture/discretefractureextensivequantities.hh Normal file
View File

@ -0,0 +1,138 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::DiscreteFractureExtensiveQuantities
*/
#ifndef EWOMS_DISCRETE_FRACTURE_EXTENSIVE_QUANTITIES_HH
#define EWOMS_DISCRETE_FRACTURE_EXTENSIVE_QUANTITIES_HH
#include <ewoms/models/immiscible/immiscibleextensivequantities.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
* \ingroup ExtensiveQuantities
*
* \brief This class expresses all intensive quantities of the discrete fracture model.
*/
template <class TypeTag>
class DiscreteFractureExtensiveQuantities : public ImmiscibleExtensiveQuantities<TypeTag>
{
typedef ImmiscibleExtensiveQuantities<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum { dimWorld = GridView::dimensionworld };
enum { numPhases = FluidSystem::numPhases };
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
public:
/*!
* \copydoc MultiPhaseBaseExtensiveQuantities::update()
*/
void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{
ParentType::update(elemCtx, scvfIdx, timeIdx);
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.interiorFace(scvfIdx);
unsigned insideScvIdx = scvf.interiorIndex();
unsigned outsideScvIdx = scvf.exteriorIndex();
unsigned globalI = elemCtx.globalSpaceIndex(insideScvIdx, timeIdx);
unsigned globalJ = elemCtx.globalSpaceIndex(outsideScvIdx, timeIdx);
const auto& fractureMapper = elemCtx.problem().fractureMapper();
if (!fractureMapper.isFractureEdge(globalI, globalJ))
// do nothing if no fracture goes though the current edge
return;
// average the intrinsic permeability of the fracture
elemCtx.problem().fractureFaceIntrinsicPermeability(fractureIntrinsicPermeability_,
elemCtx, scvfIdx, timeIdx);
auto distDirection = elemCtx.pos(outsideScvIdx, timeIdx);
distDirection -= elemCtx.pos(insideScvIdx, timeIdx);
distDirection /= distDirection.two_norm();
const auto& problem = elemCtx.problem();
fractureWidth_ = problem.fractureWidth(elemCtx, insideScvIdx,
outsideScvIdx, timeIdx);
assert(fractureWidth_ < scvf.area());
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto& pGrad = extQuants.potentialGrad(phaseIdx);
unsigned upstreamIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx);
// multiply with the fracture mobility of the upstream vertex
fractureIntrinsicPermeability_.mv(pGrad,
fractureFilterVelocity_[phaseIdx]);
fractureFilterVelocity_[phaseIdx] *= -up.fractureMobility(phaseIdx);
// divide the volume flux by two. This is required because
// a fracture is always shared by two sub-control-volume
// faces.
fractureVolumeFlux_[phaseIdx] = 0;
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
fractureVolumeFlux_[phaseIdx] +=
(fractureFilterVelocity_[phaseIdx][dimIdx] * distDirection[dimIdx])
* (fractureWidth_ / 2.0) / scvf.area();
}
}
public:
const DimMatrix& fractureIntrinsicPermeability() const
{ return fractureIntrinsicPermeability_; }
Scalar fractureVolumeFlux(unsigned phaseIdx) const
{ return fractureVolumeFlux_[phaseIdx]; }
Scalar fractureWidth() const
{ return fractureWidth_; }
const DimVector& fractureFilterVelocity(unsigned phaseIdx) const
{ return fractureFilterVelocity_[phaseIdx]; }
private:
DimMatrix fractureIntrinsicPermeability_;
DimVector fractureFilterVelocity_[numPhases];
Scalar fractureVolumeFlux_[numPhases];
Scalar fractureWidth_;
};
} // namespace Opm
#endif

244
opm/models/discretefracture/discretefractureintensivequantities.hh Normal file
View File

@ -0,0 +1,244 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::DiscreteFractureIntensiveQuantities
*/
#ifndef EWOMS_DISCRETE_FRACTURE_INTENSIVE_QUANTITIES_HH
#define EWOMS_DISCRETE_FRACTURE_INTENSIVE_QUANTITIES_HH
#include "discretefractureproperties.hh"
#include <ewoms/models/immiscible/immiscibleintensivequantities.hh>
#include <opm/material/common/Valgrind.hpp>
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
* \ingroup IntensiveQuantities
*
* \brief Contains the quantities which are are constant within a
* finite volume in the discret fracture immiscible multi-phase
* model.
*/
template <class TypeTag>
class DiscreteFractureIntensiveQuantities : public ImmiscibleIntensiveQuantities<TypeTag>
{
typedef ImmiscibleIntensiveQuantities<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum { numPhases = FluidSystem::numPhases };
enum { dimWorld = GridView::dimensionworld };
static_assert(dimWorld == 2, "The fracture module currently is only "
"implemented for the 2D case!");
static_assert(numPhases == 2, "The fracture module currently is only "
"implemented for two fluid phases!");
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
enum { wettingPhaseIdx = MaterialLaw::wettingPhaseIdx };
enum { nonWettingPhaseIdx = MaterialLaw::nonWettingPhaseIdx };
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
typedef Opm::ImmiscibleFluidState<Scalar, FluidSystem,
/*storeEnthalpy=*/enableEnergy> FluidState;
public:
DiscreteFractureIntensiveQuantities()
{ }
DiscreteFractureIntensiveQuantities(const DiscreteFractureIntensiveQuantities& other) = default;
DiscreteFractureIntensiveQuantities& operator=(const DiscreteFractureIntensiveQuantities& other) = default;
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned vertexIdx, unsigned timeIdx)
{
ParentType::update(elemCtx, vertexIdx, timeIdx);
const auto& problem = elemCtx.problem();
const auto& fractureMapper = problem.fractureMapper();
unsigned globalVertexIdx = elemCtx.globalSpaceIndex(vertexIdx, timeIdx);
Opm::Valgrind::SetUndefined(fractureFluidState_);
Opm::Valgrind::SetUndefined(fractureVolume_);
Opm::Valgrind::SetUndefined(fracturePorosity_);
Opm::Valgrind::SetUndefined(fractureIntrinsicPermeability_);
Opm::Valgrind::SetUndefined(fractureRelativePermeabilities_);
// do nothing if there is no fracture within the current degree of freedom
if (!fractureMapper.isFractureVertex(globalVertexIdx)) {
fractureVolume_ = 0;
return;
}
// Make sure that the wetting saturation in the matrix fluid
// state does not get larger than 1
Scalar SwMatrix =
std::min<Scalar>(1.0, this->fluidState_.saturation(wettingPhaseIdx));
this->fluidState_.setSaturation(wettingPhaseIdx, SwMatrix);
this->fluidState_.setSaturation(nonWettingPhaseIdx, 1 - SwMatrix);
// retrieve the facture width and intrinsic permeability from
// the problem
fracturePorosity_ =
problem.fracturePorosity(elemCtx, vertexIdx, timeIdx);
fractureIntrinsicPermeability_ =
problem.fractureIntrinsicPermeability(elemCtx, vertexIdx, timeIdx);
// compute the fracture volume for the current sub-control
// volume. note, that we don't take overlaps of fractures into
// account for this.
fractureVolume_ = 0;
const auto& vertexPos = elemCtx.pos(vertexIdx, timeIdx);
for (unsigned vertex2Idx = 0; vertex2Idx < elemCtx.numDof(/*timeIdx=*/0); ++ vertex2Idx) {
unsigned globalVertex2Idx = elemCtx.globalSpaceIndex(vertex2Idx, timeIdx);
if (vertexIdx == vertex2Idx ||
!fractureMapper.isFractureEdge(globalVertexIdx, globalVertex2Idx))
continue;
Scalar fractureWidth =
problem.fractureWidth(elemCtx, vertexIdx, vertex2Idx, timeIdx);
auto distVec = elemCtx.pos(vertex2Idx, timeIdx);
distVec -= vertexPos;
Scalar edgeLength = distVec.two_norm();
// the fracture is always adjacent to two sub-control
// volumes of the control volume, so when calculating the
// volume of the fracture which gets attributed to one
// SCV, the fracture width needs to divided by 2. Also,
// only half of the edge is located in the current control
// volume, so its length also needs to divided by 2.
fractureVolume_ += (fractureWidth / 2) * (edgeLength / 2);
}
//////////
// set the fluid state for the fracture.
//////////
// start with the same fluid state as in the matrix. This
// implies equal saturations, pressures, temperatures,
// enthalpies, etc.
fractureFluidState_.assign(this->fluidState_);
// ask the problem for the material law parameters of the
// fracture.
const auto& fractureMatParams =
problem.fractureMaterialLawParams(elemCtx, vertexIdx, timeIdx);
// calculate the fracture saturations which would be required
// to be consistent with the pressures
Scalar saturations[numPhases];
MaterialLaw::saturations(saturations, fractureMatParams,
fractureFluidState_);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fractureFluidState_.setSaturation(phaseIdx, saturations[phaseIdx]);
// Make sure that the wetting saturation in the fracture does
// not get negative
Scalar SwFracture =
std::max<Scalar>(0.0, fractureFluidState_.saturation(wettingPhaseIdx));
fractureFluidState_.setSaturation(wettingPhaseIdx, SwFracture);
fractureFluidState_.setSaturation(nonWettingPhaseIdx, 1 - SwFracture);
// calculate the relative permeabilities of the fracture
MaterialLaw::relativePermeabilities(fractureRelativePermeabilities_,
fractureMatParams,
fractureFluidState_);
// make sure that valgrind complains if the fluid state is not
// fully defined.
fractureFluidState_.checkDefined();
}
public:
/*!
* \brief Returns the effective mobility of a given phase within
* the control volume.
*
* \param phaseIdx The phase index
*/
Scalar fractureRelativePermeability(unsigned phaseIdx) const
{ return fractureRelativePermeabilities_[phaseIdx]; }
/*!
* \brief Returns the effective mobility of a given phase within
* the control volume.
*
* \param phaseIdx The phase index
*/
Scalar fractureMobility(unsigned phaseIdx) const
{
return fractureRelativePermeabilities_[phaseIdx]
/ fractureFluidState_.viscosity(phaseIdx);
}
/*!
* \brief Returns the average porosity within the fracture.
*/
Scalar fracturePorosity() const
{ return fracturePorosity_; }
/*!
* \brief Returns the average intrinsic permeability within the
* fracture.
*/
const DimMatrix& fractureIntrinsicPermeability() const
{ return fractureIntrinsicPermeability_; }
/*!
* \brief Return the volume [m^2] occupied by fractures within the
* given sub-control volume.
*/
Scalar fractureVolume() const
{ return fractureVolume_; }
/*!
* \brief Returns a fluid state object which represents the
* thermodynamic state of the fluids within the fracture.
*/
const FluidState& fractureFluidState() const
{ return fractureFluidState_; }
protected:
FluidState fractureFluidState_;
Scalar fractureVolume_;
Scalar fracturePorosity_;
DimMatrix fractureIntrinsicPermeability_;
Scalar fractureRelativePermeabilities_[numPhases];
};
} // namespace Opm
#endif

159
opm/models/discretefracture/discretefracturelocalresidual.hh Normal file
View File

@ -0,0 +1,159 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::DiscreteFractureLocalResidual
*/
#ifndef EWOMS_DISCRETE_FRACTURE_LOCAL_RESIDUAL_BASE_HH
#define EWOMS_DISCRETE_FRACTURE_LOCAL_RESIDUAL_BASE_HH
#include <ewoms/models/immiscible/immisciblelocalresidual.hh>
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
*
* \brief Calculates the local residual of the discrete fracture
* immiscible multi-phase model.
*/
template <class TypeTag>
class DiscreteFractureLocalResidual : public ImmiscibleLocalResidual<TypeTag>
{
typedef ImmiscibleLocalResidual<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
public:
/*!
* \brief Adds the amount all conservation quantities (e.g. phase
* mass) within a single fluid phase
*
* \copydetails Doxygen::storageParam
* \copydetails Doxygen::dofCtxParams
* \copydetails Doxygen::phaseIdxParam
*/
void addPhaseStorage(EqVector& storage,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx,
unsigned phaseIdx) const
{
EqVector phaseStorage(0.0);
ParentType::addPhaseStorage(phaseStorage, elemCtx, dofIdx, timeIdx, phaseIdx);
const auto& problem = elemCtx.problem();
const auto& fractureMapper = problem.fractureMapper();
unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
if (!fractureMapper.isFractureVertex(globalIdx)) {
// don't do anything in addition to the immiscible model for degrees of
// freedom that do not feature fractures
storage += phaseStorage;
return;
}
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
const auto& scv = elemCtx.stencil(timeIdx).subControlVolume(dofIdx);
// reduce the matrix storage by the fracture volume
phaseStorage *= 1 - intQuants.fractureVolume()/scv.volume();
// add the storage term inside the fractures
const auto& fsFracture = intQuants.fractureFluidState();
phaseStorage[conti0EqIdx + phaseIdx] +=
intQuants.fracturePorosity()*
fsFracture.saturation(phaseIdx) *
fsFracture.density(phaseIdx) *
intQuants.fractureVolume()/scv.volume();
EnergyModule::addFracturePhaseStorage(phaseStorage, intQuants, scv,
phaseIdx);
// add the result to the overall storage term
storage += phaseStorage;
}
/*!
* \copydoc FvBaseLocalResidual::computeFlux
*/
void computeFlux(RateVector& flux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx) const
{
ParentType::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned i = extQuants.interiorIndex();
unsigned j = extQuants.exteriorIndex();
unsigned I = elemCtx.globalSpaceIndex(i, timeIdx);
unsigned J = elemCtx.globalSpaceIndex(j, timeIdx);
const auto& fractureMapper = elemCtx.problem().fractureMapper();
if (!fractureMapper.isFractureEdge(I, J))
// do nothing if the edge from i to j is not part of a
// fracture
return;
const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
Scalar scvfArea = scvf.area();
// advective mass fluxes of all phases
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!elemCtx.model().phaseIsConsidered(phaseIdx))
continue;
// reduce the matrix mass flux by the width of the scv
// face that is occupied by the fracture. As usual, the
// fracture is shared between two SCVs, so the its width
// needs to be divided by two.
flux[conti0EqIdx + phaseIdx] *=
1 - extQuants.fractureWidth() / (2 * scvfArea);
// intensive quantities of the upstream and the downstream DOFs
unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
flux[conti0EqIdx + phaseIdx] +=
extQuants.fractureVolumeFlux(phaseIdx) * up.fractureFluidState().density(phaseIdx);
}
EnergyModule::handleFractureFlux(flux, elemCtx, scvfIdx, timeIdx);
}
};
} // namespace Opm
#endif

156
opm/models/discretefracture/discretefracturemodel.hh Normal file
View File

@ -0,0 +1,156 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::DiscreteFractureModel
*/
#ifndef EWOMS_DISCRETE_FRACTURE_MODEL_HH
#define EWOMS_DISCRETE_FRACTURE_MODEL_HH
#include <opm/material/densead/Math.hpp>
#include "discretefractureproperties.hh"
#include "discretefractureprimaryvariables.hh"
#include "discretefractureintensivequantities.hh"
#include "discretefractureextensivequantities.hh"
#include "discretefracturelocalresidual.hh"
#include "discretefractureproblem.hh"
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <ewoms/io/vtkdiscretefracturemodule.hh>
#include <opm/material/common/Exceptions.hpp>
#include <string>
namespace Opm {
template <class TypeTag>
class DiscreteFractureModel;
}
BEGIN_PROPERTIES
//! The generic type tag for problems using the immiscible multi-phase model
NEW_TYPE_TAG(DiscreteFractureModel, INHERITS_FROM(ImmiscibleTwoPhaseModel, VtkDiscreteFracture));
//! The class for the model
SET_TYPE_PROP(DiscreteFractureModel, Model, Opm::DiscreteFractureModel<TypeTag>);
//! The class for the model
SET_TYPE_PROP(DiscreteFractureModel, BaseProblem, Opm::DiscreteFractureProblem<TypeTag>);
//! Use the immiscible multi-phase local jacobian operator for the immiscible multi-phase model
SET_TYPE_PROP(DiscreteFractureModel, LocalResidual, Opm::DiscreteFractureLocalResidual<TypeTag>);
// The type of the base base class for actual problems.
// TODO!?
//SET_TYPE_PROP(DiscreteFractureModel BaseProblem, DiscreteFractureBaseProblem<TypeTag>);
//! the PrimaryVariables property
SET_TYPE_PROP(DiscreteFractureModel, PrimaryVariables,
Opm::DiscreteFracturePrimaryVariables<TypeTag>);
//! the IntensiveQuantities property
SET_TYPE_PROP(DiscreteFractureModel, IntensiveQuantities,
Opm::DiscreteFractureIntensiveQuantities<TypeTag>);
//! the ExtensiveQuantities property
SET_TYPE_PROP(DiscreteFractureModel, ExtensiveQuantities,
Opm::DiscreteFractureExtensiveQuantities<TypeTag>);
//! For the discrete fracture model, we need to use two-point flux approximation or it
//! will converge very poorly
SET_BOOL_PROP(DiscreteFractureModel, UseTwoPointGradients, true);
// The intensive quantity cache cannot be used by the discrete fracture model, because
// the intensive quantities of a control degree of freedom are not identical to the
// intensive quantities of the other intensive quantities of the same of the same degree
// of freedom. This is because the fracture properties (volume, permeability, etc) are
// specific for each...
SET_BOOL_PROP(DiscreteFractureModel, EnableIntensiveQuantityCache, false);
END_PROPERTIES
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
* \brief A fully-implicit multi-phase flow model which assumes
* immiscibility of the phases and is able to include fractures
* in the domain.
*
* This model implements multi-phase flow of \f$M > 0\f$ immiscible
* fluids \f$\alpha\f$. It also can consider edges of the
* computational grid as fractures i.e. as a porous medium with
* different higher permeability than the rest of the domain.
*
* \todo So far, the discrete fracture model only works for 2D grids
* and without energy. Also only the Darcy velocity model is
* supported for the fractures.
*
* \sa ImmiscibleModel
*/
template <class TypeTag>
class DiscreteFractureModel : public ImmiscibleModel<TypeTag>
{
typedef ImmiscibleModel<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
public:
DiscreteFractureModel(Simulator& simulator)
: ParentType(simulator)
{
if (EWOMS_GET_PARAM(TypeTag, bool, EnableIntensiveQuantityCache)) {
throw std::runtime_error("The discrete fracture model does not work in conjunction "
"with intensive quantities caching");
}
}
/*!
* \brief Register all run-time parameters for the immiscible model.
*/
static void registerParameters()
{
ParentType::registerParameters();
// register runtime parameters of the VTK output modules
Opm::VtkDiscreteFractureModule<TypeTag>::registerParameters();
}
/*!
* \copydoc FvBaseDiscretization::name
*/
static std::string name()
{ return "discretefracture"; }
void registerOutputModules_()
{
ParentType::registerOutputModules_();
this->addOutputModule(new Opm::VtkDiscreteFractureModule<TypeTag>(this->simulator_));
}
};
} // namespace Opm
#endif

124
opm/models/discretefracture/discretefractureprimaryvariables.hh Normal file
View File

@ -0,0 +1,124 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::DiscreteFracturePrimaryVariables
*/
#ifndef EWOMS_DISCRETE_FRACTURE_PRIMARY_VARIABLES_HH
#define EWOMS_DISCRETE_FRACTURE_PRIMARY_VARIABLES_HH
#include "discretefractureproperties.hh"
#include <ewoms/models/immiscible/immiscibleprimaryvariables.hh>
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
*
* \brief Represents the primary variables used by the discrete fracture
* multi-phase model.
*/
template <class TypeTag>
class DiscreteFracturePrimaryVariables
: public ImmisciblePrimaryVariables<TypeTag>
{
typedef ImmisciblePrimaryVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
public:
/*!
* \brief Default constructor
*/
DiscreteFracturePrimaryVariables() : ParentType()
{}
/*!
* \brief Constructor with assignment from scalar
*
* \param value The scalar value to which all entries of the vector will be set.
*/
DiscreteFracturePrimaryVariables(Scalar value) : ParentType(value)
{}
/*!
* \brief Copy constructor
*
* \param value The primary variables that will be duplicated.
*/
DiscreteFracturePrimaryVariables(const DiscreteFracturePrimaryVariables& value) = default;
DiscreteFracturePrimaryVariables& operator=(const DiscreteFracturePrimaryVariables& value) = default;
/*!
* \brief Directly retrieve the primary variables from an
* arbitrary fluid state of the fractures.
*
* \param fractureFluidState The fluid state of the fractures
* which should be represented by the
* primary variables. The temperatures,
* pressures and compositions of all
* phases must be defined.
* \param matParams The parameters for the capillary-pressure law
* which apply for the fracture.
*/
template <class FluidState>
void assignNaiveFromFracture(const FluidState& fractureFluidState,
const MaterialLawParams& matParams)
{
FluidState matrixFluidState;
fractureToMatrixFluidState_(matrixFluidState, fractureFluidState,
matParams);
ParentType::assignNaive(matrixFluidState);
}
private:
template <class FluidState>
void fractureToMatrixFluidState_(FluidState& matrixFluidState,
const FluidState& fractureFluidState,
const MaterialLawParams& matParams) const
{
// start with the same fluid state as in the fracture
matrixFluidState.assign(fractureFluidState);
// the condition for the equilibrium is that the pressures are
// the same in the fracture and in the matrix. This means that
// we have to find saturations for the matrix which result in
// the same pressures as in the fracture. this can be done by
// inverting the capillary pressure-saturation curve.
Scalar saturations[numPhases];
MaterialLaw::saturations(saturations, matParams, matrixFluidState);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
matrixFluidState.setSaturation(phaseIdx, saturations[phaseIdx]);
}
};
} // namespace Opm
#endif

149
opm/models/discretefracture/discretefractureproblem.hh Normal file
View File

@ -0,0 +1,149 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::DiscreteFractureProblem
*/
#ifndef EWOMS_DISCRETE_FRACTURE_PROBLEM_HH
#define EWOMS_DISCRETE_FRACTURE_PROBLEM_HH
#include "discretefractureproperties.hh"
#include <ewoms/models/common/multiphasebaseproblem.hh>
#include <opm/material/common/Means.hpp>
#include <opm/material/common/Unused.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
BEGIN_PROPERTIES
NEW_PROP_TAG(ThermalConductionLawParams);
NEW_PROP_TAG(EnableGravity);
NEW_PROP_TAG(FluxModule);
END_PROPERTIES
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
* \brief The base class for the problems of ECFV discretizations which deal
* with a multi-phase flow through a porous medium.
*/
template<class TypeTag>
class DiscreteFractureProblem
: public MultiPhaseBaseProblem<TypeTag>
{
typedef Opm::MultiPhaseBaseProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
enum { dimWorld = GridView::dimensionworld };
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
public:
/*!
* \copydoc Problem::FvBaseProblem(Simulator& )
*/
DiscreteFractureProblem(Simulator& simulator)
: ParentType(simulator)
{}
/*!
* \brief Returns the intrinsic permeability of a face due to a fracture.
*
* This method is specific to the finite volume discretizations. If left unspecified,
* it calls the intrinsicPermeability() methods for the face's interior and exterior
* finite volume and averages them harmonically. Note that if this function is
* defined, the intrinsicPermeability() method does not need to be defined by the
* problem (if a finite-volume discretization is used).
*/
template <class Context>
void fractureFaceIntrinsicPermeability(DimMatrix& result,
const Context& context,
unsigned localFaceIdx,
unsigned timeIdx) const
{
const auto& scvf = context.stencil(timeIdx).interiorFace(localFaceIdx);
unsigned interiorElemIdx = scvf.interiorIndex();
unsigned exteriorElemIdx = scvf.exteriorIndex();
const DimMatrix& K1 = asImp_().fractureIntrinsicPermeability(context, interiorElemIdx, timeIdx);
const DimMatrix& K2 = asImp_().fractureIntrinsicPermeability(context, exteriorElemIdx, timeIdx);
// entry-wise harmonic mean. this is almost certainly wrong if
// you have off-main diagonal entries in your permeabilities!
for (unsigned i = 0; i < dimWorld; ++i)
for (unsigned j = 0; j < dimWorld; ++j)
result[i][j] = Opm::harmonicMean(K1[i][j], K2[i][j]);
}
/*!
* \brief Returns the intrinsic permeability tensor \f$[m^2]\f$ at a given position due to a fracture
*
* \param context Reference to the object which represents the
* current execution context.
* \param spaceIdx The local index of spatial entity defined by the context
* \param timeIdx The index used by the time discretization.
*/
template <class Context>
const DimMatrix& fractureIntrinsicPermeability(const Context& context OPM_UNUSED,
unsigned spaceIdx OPM_UNUSED,
unsigned timeIdx OPM_UNUSED) const
{
throw std::logic_error("Not implemented: Problem::fractureIntrinsicPermeability()");
}
/*!
* \brief Returns the porosity [] inside fractures for a given control volume.
*
* \param context Reference to the object which represents the
* current execution context.
* \param spaceIdx The local index of spatial entity defined by the context
* \param timeIdx The index used by the time discretization.
*/
template <class Context>
Scalar fracturePorosity(const Context& context OPM_UNUSED,
unsigned spaceIdx OPM_UNUSED,
unsigned timeIdx OPM_UNUSED) const
{
throw std::logic_error("Not implemented: Problem::fracturePorosity()");
}
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation& asImp_()
{ return *static_cast<Implementation *>(this); }
//! \copydoc asImp_()
const Implementation& asImp_() const
{ return *static_cast<const Implementation *>(this); }
};
} // namespace Opm
#endif

43
opm/models/discretefracture/discretefractureproperties.hh Normal file
View File

@ -0,0 +1,43 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \ingroup DiscreteFractureModel
*
* \brief Defines the properties required for the immiscible
* multi-phase model which considers discrete fractures.
*/
#ifndef EWOMS_DISCRETE_FRACTIRE_PROPERTIES_HH
#define EWOMS_DISCRETE_FRACTIRE_PROPERTIES_HH
#include <ewoms/models/immiscible/immiscibleproperties.hh>
#include <ewoms/io/vtkdiscretefracturemodule.hh>
BEGIN_PROPERTIES
NEW_PROP_TAG(UseTwoPointGradients);
END_PROPERTIES
#endif

108
opm/models/discretefracture/fracturemapper.hh Normal file
View File

@ -0,0 +1,108 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \copydoc Opm::FractureMapper
*/
#ifndef EWOMS_FRACTURE_MAPPER_HH
#define EWOMS_FRACTURE_MAPPER_HH
#include <opm/models/utils/propertysystem.hh>
#include <algorithm>
#include <set>
namespace Opm {
/*!
* \ingroup DiscreteFractureModel
* \brief Stores the topology of fractures.
*/
template <class TypeTag>
class FractureMapper
{
struct FractureEdge
{
FractureEdge(unsigned edgeVertex1Idx, unsigned edgeVertex2Idx)
: i_(std::min(edgeVertex1Idx, edgeVertex2Idx)),
j_(std::max(edgeVertex1Idx, edgeVertex2Idx))
{}
bool operator<(const FractureEdge& e) const
{ return i_ < e.i_ || (i_ == e.i_ && j_ < e.j_); }
bool operator==(const FractureEdge& e) const
{ return i_ == e.i_ && j_ == e.j_; }
unsigned i_;
unsigned j_;
};
public:
/*!
* \brief Constructor
*/
FractureMapper()
{}
/*!
* \brief Marks an edge as having a fracture.
*
* \param vertexIdx1 The index of the edge's first vertex.
* \param vertexIdx2 The index of the edge's second vertex.
*/
void addFractureEdge(unsigned vertexIdx1, unsigned vertexIdx2)
{
fractureEdges_.insert(FractureEdge(vertexIdx1, vertexIdx2));
fractureVertices_.insert(vertexIdx1);
fractureVertices_.insert(vertexIdx2);
}
/*!
* \brief Returns true iff a fracture cuts through a given vertex.
*
* \param vertexIdx The index of the vertex.
*/
bool isFractureVertex(unsigned vertexIdx) const
{ return fractureVertices_.count(vertexIdx) > 0; }
/*!
* \brief Returns true iff a fracture is associated with a given edge.
*
* \param vertex1Idx The index of the first vertex of the edge.
* \param vertex2Idx The index of the second vertex of the edge.
*/
bool isFractureEdge(unsigned vertex1Idx, unsigned vertex2Idx) const
{
FractureEdge tmp(vertex1Idx, vertex2Idx);
return fractureEdges_.count(tmp) > 0;
}
private:
std::set<FractureEdge> fractureEdges_;
std::set<unsigned> fractureVertices_;
};
} // namespace Opm
#endif // EWOMS_FRACTURE_MAPPER_HH