diff --git a/opm/models/discretefracture/discretefractureextensivequantities.hh b/opm/models/discretefracture/discretefractureextensivequantities.hh new file mode 100644 index 000000000..678b12a3e --- /dev/null +++ b/opm/models/discretefracture/discretefractureextensivequantities.hh @@ -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 . + + 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 + +#include +#include + +namespace Opm { + +/*! + * \ingroup DiscreteFractureModel + * \ingroup ExtensiveQuantities + * + * \brief This class expresses all intensive quantities of the discrete fracture model. + */ +template +class DiscreteFractureExtensiveQuantities : public ImmiscibleExtensiveQuantities +{ + typedef ImmiscibleExtensiveQuantities 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 DimMatrix; + typedef Dune::FieldVector 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(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 diff --git a/opm/models/discretefracture/discretefractureintensivequantities.hh b/opm/models/discretefracture/discretefractureintensivequantities.hh new file mode 100644 index 000000000..2c89d1cec --- /dev/null +++ b/opm/models/discretefracture/discretefractureintensivequantities.hh @@ -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 . + + 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 + +#include + +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 DiscreteFractureIntensiveQuantities : public ImmiscibleIntensiveQuantities +{ + typedef ImmiscibleIntensiveQuantities 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 DimMatrix; + typedef Opm::ImmiscibleFluidState 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(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(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 diff --git a/opm/models/discretefracture/discretefracturelocalresidual.hh b/opm/models/discretefracture/discretefracturelocalresidual.hh new file mode 100644 index 000000000..b056484cb --- /dev/null +++ b/opm/models/discretefracture/discretefracturelocalresidual.hh @@ -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 . + + 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 + +namespace Opm { + +/*! + * \ingroup DiscreteFractureModel + * + * \brief Calculates the local residual of the discrete fracture + * immiscible multi-phase model. + */ +template +class DiscreteFractureLocalResidual : public ImmiscibleLocalResidual +{ + typedef ImmiscibleLocalResidual 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 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(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 diff --git a/opm/models/discretefracture/discretefracturemodel.hh b/opm/models/discretefracture/discretefracturemodel.hh new file mode 100644 index 000000000..92b724558 --- /dev/null +++ b/opm/models/discretefracture/discretefracturemodel.hh @@ -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 . + + 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 + +#include "discretefractureproperties.hh" +#include "discretefractureprimaryvariables.hh" +#include "discretefractureintensivequantities.hh" +#include "discretefractureextensivequantities.hh" +#include "discretefracturelocalresidual.hh" +#include "discretefractureproblem.hh" + +#include +#include + +#include + +#include + +namespace Opm { +template +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); + +//! The class for the model +SET_TYPE_PROP(DiscreteFractureModel, BaseProblem, Opm::DiscreteFractureProblem); + +//! Use the immiscible multi-phase local jacobian operator for the immiscible multi-phase model +SET_TYPE_PROP(DiscreteFractureModel, LocalResidual, Opm::DiscreteFractureLocalResidual); + +// The type of the base base class for actual problems. +// TODO!? +//SET_TYPE_PROP(DiscreteFractureModel BaseProblem, DiscreteFractureBaseProblem); + +//! the PrimaryVariables property +SET_TYPE_PROP(DiscreteFractureModel, PrimaryVariables, + Opm::DiscreteFracturePrimaryVariables); + +//! the IntensiveQuantities property +SET_TYPE_PROP(DiscreteFractureModel, IntensiveQuantities, + Opm::DiscreteFractureIntensiveQuantities); + +//! the ExtensiveQuantities property +SET_TYPE_PROP(DiscreteFractureModel, ExtensiveQuantities, + Opm::DiscreteFractureExtensiveQuantities); + +//! 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 DiscreteFractureModel : public ImmiscibleModel +{ + typedef ImmiscibleModel 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::registerParameters(); + } + + /*! + * \copydoc FvBaseDiscretization::name + */ + static std::string name() + { return "discretefracture"; } + + void registerOutputModules_() + { + ParentType::registerOutputModules_(); + + this->addOutputModule(new Opm::VtkDiscreteFractureModule(this->simulator_)); + } +}; +} // namespace Opm + +#endif diff --git a/opm/models/discretefracture/discretefractureprimaryvariables.hh b/opm/models/discretefracture/discretefractureprimaryvariables.hh new file mode 100644 index 000000000..454fc300b --- /dev/null +++ b/opm/models/discretefracture/discretefractureprimaryvariables.hh @@ -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 . + + 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 + +namespace Opm { +/*! + * \ingroup DiscreteFractureModel + * + * \brief Represents the primary variables used by the discrete fracture + * multi-phase model. + */ +template +class DiscreteFracturePrimaryVariables + : public ImmisciblePrimaryVariables +{ + typedef ImmisciblePrimaryVariables 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 + void assignNaiveFromFracture(const FluidState& fractureFluidState, + const MaterialLawParams& matParams) + { + FluidState matrixFluidState; + fractureToMatrixFluidState_(matrixFluidState, fractureFluidState, + matParams); + + ParentType::assignNaive(matrixFluidState); + } + +private: + template + 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 diff --git a/opm/models/discretefracture/discretefractureproblem.hh b/opm/models/discretefracture/discretefractureproblem.hh new file mode 100644 index 000000000..c7c9aa5dd --- /dev/null +++ b/opm/models/discretefracture/discretefractureproblem.hh @@ -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 . + + 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 + +#include +#include +#include + +#include +#include + +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 DiscreteFractureProblem + : public MultiPhaseBaseProblem +{ + typedef Opm::MultiPhaseBaseProblem 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 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 + 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 + 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 + 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(this); } + //! \copydoc asImp_() + const Implementation& asImp_() const + { return *static_cast(this); } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/discretefracture/discretefractureproperties.hh b/opm/models/discretefracture/discretefractureproperties.hh new file mode 100644 index 000000000..a2b96cb13 --- /dev/null +++ b/opm/models/discretefracture/discretefractureproperties.hh @@ -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 . + + 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 + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(UseTwoPointGradients); + +END_PROPERTIES + +#endif diff --git a/opm/models/discretefracture/fracturemapper.hh b/opm/models/discretefracture/fracturemapper.hh new file mode 100644 index 000000000..c18fb1cb1 --- /dev/null +++ b/opm/models/discretefracture/fracturemapper.hh @@ -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 . + + 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 + +#include +#include + +namespace Opm { + +/*! + * \ingroup DiscreteFractureModel + * \brief Stores the topology of fractures. + */ +template +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 fractureEdges_; + std::set fractureVertices_; +}; + +} // namespace Opm + +#endif // EWOMS_FRACTURE_MAPPER_HH