Commit 9bc34068 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[multidomain][facet] Add coupling manager for conforming facet coupling (box/cc) with multidomain

parent 8c475f91
add_subdirectory(boundary)
add_subdirectory(embedded)
add_subdirectory(facet)
install(FILES
couplingjacobianpattern.hh
......
add_subdirectory(box)
add_subdirectory(cellcentered)
install(FILES
codimonegridadapter.hh
couplingmanager.hh
couplingmapper.hh
couplingmapperbase.hh
enrichmenthelper.hh
gmshreader.hh
gridcreator.hh
vertexmapper.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/multidomain/facet)
install(FILES
couplingmanager.hh
couplingmapper.hh
darcyslaw.hh
elementboundarytypes.hh
fvelementgeometry.hh
fvgridgeometry.hh
localresidual.hh
properties.hh
subcontrolvolumeface.hh
upwindscheme.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/multidomain/facet/box)
This diff is collapsed.
/*****************************************************************************
* 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
* \ingroup MixedDimension
* \ingroup MixedDimensionFacet
* \copydoc Dumux::FacetCouplingMapper
*/
#ifndef DUMUX_BOX_FACETCOUPLING_MAPPER_HH
#define DUMUX_BOX_FACETCOUPLING_MAPPER_HH
#include <dune/common/indices.hh>
#include <dune/geometry/referenceelements.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/multidomain/facet/couplingmapper.hh>
#include <dumux/multidomain/facet/couplingmapperbase.hh>
#include <dumux/multidomain/facet/codimonegridadapter.hh>
namespace Dumux {
/*!
* \ingroup MixedDimension
* \ingroup MixedDimensionFacet
* \brief Base class for the coupling mapper that sets up and stores
* the coupling maps between two domains of dimension d and (d-1).
* This specialization is for the bulk domain using the box scheme.
*
* \tparam BulkFVG The d-dimensional finite-volume grid geometry
* \tparam LowDimFVG The (d-1)-dimensional finite-volume grid geometry
* \tparam bulkId The domain id of the bulk problem
* \tparam lowDimId The domain id of the lower-dimensional problem
*/
template<class BulkFVG, class LowDimFVG, std::size_t bulkId, std::size_t lowDimId>
class FacetCouplingMapper<BulkFVG, LowDimFVG, bulkId, lowDimId, DiscretizationMethod::box>
: public virtual FacetCouplingMapperBase<BulkFVG, LowDimFVG, bulkId, lowDimId>
{
using ParentType = FacetCouplingMapperBase<BulkFVG, LowDimFVG, bulkId, lowDimId>;
using LowDimElement = typename LowDimFVG::GridView::template Codim<0>::Entity;
// dimensions of the two grids
using BulkGridView = typename BulkFVG::GridView;
using LowDimGridView = typename LowDimFVG::GridView;
static constexpr int bulkDim = BulkGridView::dimension;
static constexpr int lowDimDim = LowDimGridView::dimension;
public:
//! export domain ids
using ParentType::bulkGridId;
using ParentType::facetGridId;
/*!
* \brief Update coupling maps. This is the standard
* interface required by any mapper implementation.
*
* \param bulkFvGridGeometry The finite-volume grid geometry of the bulk grid
* \param lowDimFvGridGeometry The finite-volume grid geometry of the lower-dimensional grid
* \param gridCreator Class that contains the grid factories and embedments
*/
template< class GridCreator >
void update(const BulkFVG& bulkFvGridGeometry,
const LowDimFVG& lowDimFvGridGeometry,
const GridCreator& gridCreator)
{
// forward to update function with instantiated vertex adapter
using GridAdapter = CodimOneGridAdapter<GridCreator, bulkGridId, facetGridId>;
update(bulkFvGridGeometry, lowDimFvGridGeometry, gridCreator, GridAdapter{gridCreator});
}
/*!
* \brief Update coupling maps. This is the standard
* interface required by any mapper implementation.
*
* \param bulkFvGridGeometry The finite-volume grid geometry of the bulk grid
* \param lowDimFvGridGeometry The finite-volume grid geometry of the lower-dimensional grid
* \param gridCreator Class that contains the grid factories and embedments
* \param codimOneGridAdapter Allows direct access to data on the bulk grid for lowdim grid entities
*/
template< class GridCreator, class CodimOneGridAdapter >
void update(const BulkFVG& bulkFvGridGeometry,
const LowDimFVG& lowDimFvGridGeometry,
const GridCreator& gridCreator,
const CodimOneGridAdapter& codimOneGridAdapter)
{
// define the execution policy how to add map entries from an embedment
auto addEmbedmentPolicy = [&] (auto&& embedments,
const LowDimElement& lowDimElement,
const LowDimFVG& lowDimFvGridGeometry,
const BulkFVG& bulkFvGridGeometry)
{
using LowDimIndexType = typename LowDimGridView::IndexSet::IndexType;
using BulkIndexType = typename BulkGridView::IndexSet::IndexType;
using BulkReferenceElements = Dune::ReferenceElements<typename BulkGridView::ctype, bulkDim>;
const auto lowDimElemIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement);
auto& lowDimData = this->couplingMap_(facetGridId, bulkGridId)[lowDimElemIdx];
// determine corner indices (in bulk grid indices)
const auto& eg = lowDimElement.geometry();
const auto numElementCorners = eg.corners();
std::vector<BulkIndexType> elementCorners(numElementCorners);
for (int i = 0; i < numElementCorners; ++i)
elementCorners[i] = codimOneGridAdapter.bulkGridVertexIndex(lowDimElement.template subEntity<lowDimDim>(i));
// save unsorted set of corner indices and search scvfs in embedments
const auto unsortedElemCorners = elementCorners;
std::sort(elementCorners.begin(), elementCorners.end());
for (auto bulkElemIdx : embedments)
{
const auto bulkElement = bulkFvGridGeometry.element(bulkElemIdx);
const auto bulkRefElem = BulkReferenceElements::general(bulkElement.geometry().type());
// find the bulk element facet that lies on this low dim element (assumes conformity!)
bool found = false;
unsigned int coupledFacetIndex;
for (const auto& is : intersections(bulkFvGridGeometry.gridView(), bulkElement))
{
// determine if it lies on low dim element by comparing corner indices
const auto numCorners = is.geometry().corners();
std::vector<BulkIndexType> facetIndices(numCorners);
for (int i = 0; i < numCorners; ++i)
{
const auto vIdxLocal = bulkRefElem.subEntity(is.indexInInside(), 1, i, bulkDim);
facetIndices[i] = bulkFvGridGeometry.vertexMapper().vertexIndex(bulkElement.template subEntity<bulkDim>(vIdxLocal));
}
std::sort(facetIndices.begin(), facetIndices.end());
if ( std::equal(facetIndices.begin(), facetIndices.end(), elementCorners.begin(), elementCorners.end()) )
{
coupledFacetIndex = is.indexInInside();
found = true; break;
}
}
// ensure that we found the facet!
if (!found)
DUNE_THROW(Dune::InvalidStateException, "Could not find the bulk element coupling facet!");
// we should always find numElementCorners coupling scvfs
auto fvGeometry = localView(bulkFvGridGeometry);
fvGeometry.bindElement(bulkElement);
unsigned int foundCounter = 0;
std::vector<BulkIndexType> embeddedScvfIndices(numElementCorners);
for (const auto& scvf : scvfs(fvGeometry))
{
if (scvf.interiorBoundary() && scvf.facetIndexInElement() == coupledFacetIndex)
{
// we want to order the set of scvfs lying on the lower-dimensional element such that the i-th scvf
// coincides with the i-th low dim element corner. This will help later to identify which scvf fluxes
// enter which scv of the low dim element if the lower-dimensional domain uses the box scheme
const auto vIdxLocal = bulkRefElem.subEntity(coupledFacetIndex, 1, scvf.indexInElementFacet(), bulkDim);
const auto vIdxGlobal = bulkFvGridGeometry.vertexMapper().vertexIndex(bulkElement, vIdxLocal, bulkDim);
const auto it = std::find(unsortedElemCorners.begin(), unsortedElemCorners.end(), vIdxGlobal);
assert(it != unsortedElemCorners.end());
const auto lowDimElemLocalCornerIdx = std::distance(unsortedElemCorners.begin(), it);
embeddedScvfIndices[lowDimElemLocalCornerIdx] = scvf.index();
foundCounter++;
}
}
// ensure we found all scvfs
if (foundCounter != numElementCorners)
DUNE_THROW(Dune::InvalidStateException, "Could not find all coupling scvfs in the bulk element");
// add each dof in the low dim element to coupling stencil of the bulk element
auto& bulkData = this->couplingMap_(bulkGridId, facetGridId)[bulkElemIdx];
const auto lowDimElementDofs = LowDimFVG::discMethod == DiscretizationMethod::cctpfa
? std::vector<LowDimIndexType>({lowDimElemIdx})
: this->extractNodalDofs_(lowDimElement, lowDimFvGridGeometry);
for (auto lowDimDofIdx : lowDimElementDofs)
{
bulkData.couplingStencil.push_back(lowDimDofIdx);
auto& couplingScvfs = bulkData.dofToCouplingScvfMap[lowDimDofIdx];
couplingScvfs.insert(couplingScvfs.end(), embeddedScvfIndices.begin(), embeddedScvfIndices.end());
}
// add info on which scvfs coincide with which low dim element
bulkData.couplingElementStencil.push_back(lowDimElemIdx);
auto& elemToScvfMap = bulkData.elementToScvfMap[lowDimElemIdx];
elemToScvfMap.insert(elemToScvfMap.end(), embeddedScvfIndices.begin(), embeddedScvfIndices.end());
// add embedment and the bulk cell dofs to coupling stencil of low dim element
lowDimData.embedments.emplace_back(bulkElemIdx, std::move(embeddedScvfIndices));
const auto bulkElementDofs = this->extractNodalDofs_(bulkElement, bulkFvGridGeometry);
for (auto bulkDofIdx : bulkElementDofs)
lowDimData.couplingStencil.push_back(bulkDofIdx);
}
};
// let the parent do the update subject to the execution policy defined above
ParentType::update_(bulkFvGridGeometry, lowDimFvGridGeometry, gridCreator, addEmbedmentPolicy);
// coupling stencils might not be unique with the policy above
auto makeStencilUnique = [] (auto& data)
{
auto& cs = data.second.couplingStencil;
std::sort(cs.begin(), cs.end());
cs.erase( std::unique(cs.begin(), cs.end()), cs.end() );
};
auto& lowDimCouplingData = this->couplingMap_(facetGridId, bulkGridId);
std::for_each(lowDimCouplingData.begin(), lowDimCouplingData.end(), makeStencilUnique);
// bulk coupling stencil is only non-unique if box is used
if (LowDimFVG::discMethod == DiscretizationMethod::box)
{
auto& bulkCouplingData = this->couplingMap_(bulkGridId, facetGridId);
std::for_each(bulkCouplingData.begin(), bulkCouplingData.end(), makeStencilUnique);
}
}
};
} // end namespace Dumux
#endif // DUMUX_FACETCOUPLING_CCTPFA_COUPLING_MAPPER_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
* \copydoc Dumux::BoxFacetCouplingDarcysLaw
*/
#ifndef DUMUX_DISCRETIZATION_BOX_FACET_COUPLING_DARCYS_LAW_HH
#define DUMUX_DISCRETIZATION_BOX_FACET_COUPLING_DARCYS_LAW_HH
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/discretization/box/darcyslaw.hh>
namespace Dumux {
/*!
* \ingroup MixedDimension
* \ingroup MixedDimensionFacet
* \brief Darcy's law for the box scheme scheme in the context of coupled models
* where coupling occurs across the facets of the bulk domain elements
* with a lower-dimensional domain living on these facets.
*/
template<class Scalar, class FVGridGeometry>
class BoxFacetCouplingDarcysLaw
{
using DefaultBoxDarcysLaw = BoxDarcysLaw<Scalar, FVGridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using CoordScalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
public:
template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const int phaseIdx,
const ElementFluxVarsCache& elemFluxVarCache)
{
// if this scvf is not on an interior boundary, use the standard law
if (!scvf.interiorBoundary())
return DefaultBoxDarcysLaw::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarCache);
// get some references for convenience
const auto& fluxVarCache = elemFluxVarCache[scvf];
const auto& shapeValues = fluxVarCache.shapeValues();
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
// evaluate user-defined interior boundary types
const auto bcTypes = problem.interiorBoundaryTypes(element, scvf);
static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
// on interior Neumann boundaries, evaluate the flux using the facet permeability
if (bcTypes.hasOnlyNeumann())
{
// interpolate pressure/density to scvf integration point
Scalar p = 0.0;
Scalar rho = 0.0;
for (const auto& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
p += volVars.pressure(phaseIdx)*shapeValues[scv.indexInElement()][0];
rho += volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0];
}
// compute tpfa flux from integration point to facet centerline
const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
using std::sqrt;
// If this is a surface grid, use the square root of the facet extrusion factor
// as an approximate average distance from scvf ip to facet center
using std::sqrt;
const auto a = facetVolVars.extrusionFactor();
auto gradP = scvf.unitOuterNormal();
gradP *= dim == dimWorld ? 0.5*a : 0.5*sqrt(a);
gradP /= gradP.two_norm2();
gradP *= (facetVolVars.pressure(phaseIdx) - p);
if (enableGravity)
gradP.axpy(-rho, problem.gravityAtPos(scvf.center()));
// apply facet permeability and return the flux
return -1.0*scvf.area()
*insideVolVars.extrusionFactor()
*vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), gradP);
}
// on interior Dirichlet boundaries use the facet pressure and evaluate flux
else if (bcTypes.hasOnlyDirichlet())
{
// create vector with nodal pressures
std::vector<Scalar> pressures(element.subEntities(dim));
for (const auto& scv : scvs(fvGeometry))
pressures[scv.localDofIndex()] = elemVolVars[scv].pressure(phaseIdx);
// substitute with facet pressures for those scvs touching this facet
for (const auto& scvfJ : scvfs(fvGeometry))
if (scvfJ.interiorBoundary() && scvfJ.facetIndexInElement() == scvf.facetIndexInElement())
pressures[ fvGeometry.scv(scvfJ.insideScvIdx()).localDofIndex() ]
= problem.couplingManager().getLowDimVolVars(element, scvfJ).pressure(phaseIdx);
// evaluate gradP - rho*g at integration point
Scalar rho(0.0);
Dune::FieldVector<Scalar, dimWorld> gradP(0.0);
for (const auto& scv : scvs(fvGeometry))
{
rho += elemVolVars[scv].density(phaseIdx)*shapeValues[scv.indexInElement()][0];
gradP.axpy(pressures[scv.localDofIndex()], fluxVarCache.gradN(scv.indexInElement()));
}
if (enableGravity)
gradP.axpy(-rho, problem.gravityAtPos(scvf.center()));
// apply matrix permeability and return the flux
return -1.0*scvf.area()
*insideVolVars.extrusionFactor()
*vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), gradP);
}
// mixed boundary types are not supported
else
DUNE_THROW(Dune::NotImplemented, "Mixed boundary types are not supported");
}
// compute transmissibilities ti for analytical jacobians
template<class Problem, class ElementVolumeVariables, class FluxVarCache>
static std::vector<Scalar> calculateTransmissibilities(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const FluxVarCache& fluxVarCache)
{
DUNE_THROW(Dune::NotImplemented, "transmissibilty computation for BoxFacetCouplingDarcysLaw");
}
};
} // end 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 Boundary types gathered on an element for the box scheme
* with coupling occurring across the element facets.
*/
#ifndef DUMUX_FACET_COUPLING_BOX_ELEMENT_BOUNDARY_TYPES_HH
#define DUMUX_FACET_COUPLING_BOX_ELEMENT_BOUNDARY_TYPES_HH
#include <cassert>
#include <vector>
#include <dumux/discretization/box/elementboundarytypes.hh>
namespace Dumux {
/*!
* \ingroup MixedDimension
* \ingroup MixedDimensionFacet
* \brief This class stores an array of BoundaryTypes objects
* on an element for the box scheme with coupling occuring
* across the element facets.
*/
template<class BTypes>
class BoxFacetCouplingElementBoundaryTypes : public BoxElementBoundaryTypes<BTypes>
{
public:
/*!
* \brief Update the boundary types for all vertices of an element.
*
* \param problem The problem object which needs to be simulated
* \param element The DUNE Codim<0> entity for which the boundary
* types should be collected
* \param fvGeometry The element's finite volume geometry
*
* \note We need this overload so that we can call different
* problem interfaces for domain boundary and interior boundaries.
*/
template<class Problem, class Element, class FVElementGeometry>
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry)
{
using FVGridGeometry = typename FVElementGeometry::FVGridGeometry;
using GridView = typename FVGridGeometry::GridView;
this->vertexBCTypes_.resize( element.subEntities(GridView::dimension) );
this->hasDirichlet_ = false;
this->hasNeumann_ = false;
this->hasOutflow_ = false;
for (const auto& scv : scvs(fvGeometry))
{
int scvIdxLocal = scv.localDofIndex();
this->vertexBCTypes_[scvIdxLocal].reset();
// lambda to update the element boundary info
auto updateElemBCInfo = [&] (const BTypes& bcTypes)
{
this->hasDirichlet_ = this->hasDirichlet_ || this->vertexBCTypes_[scvIdxLocal].hasDirichlet();
this->hasNeumann_ = this->hasNeumann_ || this->vertexBCTypes_[scvIdxLocal].hasNeumann();
this->hasOutflow_ = this->hasOutflow_ || this->vertexBCTypes_[scvIdxLocal].hasOutflow();
};
// We have to have Neumann-type BCs on nodes that touch interior boundaries
if (fvGeometry.fvGridGeometry().dofOnInteriorBoundary(scv.dofIndex()))
{
this->vertexBCTypes_[scvIdxLocal].setAllNeumann();
updateElemBCInfo(this->vertexBCTypes_[scvIdxLocal]);
}
// otherwise, let the problem decide
else if (fvGeometry.fvGridGeometry().dofOnBoundary(scv.dofIndex()))
{
this->vertexBCTypes_[scvIdxLocal] = problem.boundaryTypes(element, scv);
updateElemBCInfo(this->vertexBCTypes_[scvIdxLocal]);
}
}
}
};
} // 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
* \ingroup MixedDimension
* \ingroup MixedDimensionFacet
* \ingroup BoxDiscretization
* \brief \copydoc Dumux::BoxFacetCouplingFVElementGeometry
*/
#ifndef DUMUX_FACETCOUPLING_BOX_FV_ELEMENT_GEOMETRY_HH
#define DUMUX_FACETCOUPLING_BOX_FV_ELEMENT_GEOMETRY_HH
#include <algorithm>
#include <dune/geometry/referenceelements.hh>
#include <dumux/discretization/scvandscvfiterators.hh>
#include <dumux/discretization/box/boxgeometryhelper.hh>
namespace Dumux {
/*!
* \ingroup MixedDimension
* \ingroup MixedDimensionFacet
* \ingroup BoxDiscretization
* \brief Base class for the element-local finite volume geometry for box models
* in the context of models considering coupling of different domains across the
* bulk grid facets. This builds up the sub control volumes and sub control volume