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

[mpfa] do not allow mixed BC

The treatment of mixed BC is questionable for cc methods. Not allowing
for it enables us to not have to store n interaction volume seeds at the boundary.
parent 651643f9
...@@ -278,10 +278,8 @@ public: ...@@ -278,10 +278,8 @@ public:
{ return boundaryInfo_[eqIdx].isDirichlet; } { return boundaryInfo_[eqIdx].isDirichlet; }
/*! /*!
* \brief Returns true if an equation is used to specify a * \brief Returns true if all equations are used to specify a
* Dirichlet condition. * Dirichlet condition.
*
* \param eqIdx The index of the equation
*/ */
bool hasOnlyDirichlet() const bool hasOnlyDirichlet() const
{ {
...@@ -312,6 +310,18 @@ public: ...@@ -312,6 +310,18 @@ public:
bool isNeumann(unsigned eqIdx) const bool isNeumann(unsigned eqIdx) const
{ return boundaryInfo_[eqIdx].isNeumann; } { return boundaryInfo_[eqIdx].isNeumann; }
/*!
* \brief Returns true if all equations are used to specify a
* Neumann condition.
*/
bool hasOnlyNeumann() const
{
return std::all_of(boundaryInfo_.begin(),
boundaryInfo_.end(),
[](const BoundaryInfo& b){ return b.isNeumann; }
);
}
/*! /*!
* \brief Returns true if some equation is used to specify a * \brief Returns true if some equation is used to specify a
* Neumann condition. * Neumann condition.
......
...@@ -119,7 +119,7 @@ public: ...@@ -119,7 +119,7 @@ public:
Stencil stencil; Stencil stencil;
if(problem.model().globalFvGeometry().scvfTouchesBoundary(scvf)) if(problem.model().globalFvGeometry().scvfTouchesBoundary(scvf))
{ {
const auto& ivSeed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf, /*dummy*/0); const auto& ivSeed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
const auto& localScvSeeds = ivSeed.scvSeeds(); const auto& localScvSeeds = ivSeed.scvSeeds();
stencil.reserve(localScvSeeds.size()); stencil.reserve(localScvSeeds.size());
......
...@@ -78,13 +78,13 @@ public: ...@@ -78,13 +78,13 @@ public:
// update the flux var caches for this scvf // update the flux var caches for this scvf
if (problem.model().globalFvGeometry().scvfTouchesBoundary(scvf)) if (problem.model().globalFvGeometry().scvfTouchesBoundary(scvf))
{ {
const auto& boundarySeed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
BoundaryInteractionVolume iv(boundarySeed, problem, fvGeometry, elemVolVars);
iv.solveLocalSystem(permFunction);
// we assume phaseIdx = eqIdx here for purely advective problems // we assume phaseIdx = eqIdx here for purely advective problems
for (unsigned int eqIdx = 0; eqIdx < numEq; ++eqIdx) for (unsigned int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{ {
const auto& boundarySeed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf, eqIdx);
BoundaryInteractionVolume iv(boundarySeed, problem, fvGeometry, elemVolVars);
iv.solveLocalSystem(permFunction);
// lambda function defining the upwind factor of the advective flux // lambda function defining the upwind factor of the advective flux
auto advectionUpwindFunction = [eqIdx](const VolumeVariables& volVars) { return volVars.density(eqIdx)/volVars.viscosity(eqIdx); }; auto advectionUpwindFunction = [eqIdx](const VolumeVariables& volVars) { return volVars.density(eqIdx)/volVars.viscosity(eqIdx); };
iv.assembleNeumannFluxes(advectionUpwindFunction, eqIdx); iv.assembleNeumannFluxes(advectionUpwindFunction, eqIdx);
......
...@@ -107,8 +107,8 @@ public: ...@@ -107,8 +107,8 @@ public:
{ return globalInteractionVolumeSeeds_.seed(scvf); } { return globalInteractionVolumeSeeds_.seed(scvf); }
//! Get the boundary interaction volume seed corresponding to an scvf //! Get the boundary interaction volume seed corresponding to an scvf
const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const SubControlVolumeFace& scvf, const LocalIndexType eqIdx) const const BoundaryInteractionVolumeSeed& boundaryInteractionVolumeSeed(const SubControlVolumeFace& scvf) const
{ return globalInteractionVolumeSeeds_.boundarySeed(scvf, eqIdx); } { return globalInteractionVolumeSeeds_.boundarySeed(scvf); }
//! Returns whether or not a scvf touches the boundary (has to be called before getting an interaction volume) //! Returns whether or not a scvf touches the boundary (has to be called before getting an interaction volume)
bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const
......
...@@ -74,18 +74,13 @@ public: ...@@ -74,18 +74,13 @@ public:
const InteractionVolumeSeed& seed(const SubControlVolumeFace& scvf) const const InteractionVolumeSeed& seed(const SubControlVolumeFace& scvf) const
{ return seeds_[scvfIndexMap_[scvf.index()]]; } { return seeds_[scvfIndexMap_[scvf.index()]]; }
const BoundaryInteractionVolumeSeed& boundarySeed(const SubControlVolumeFace& scvf, const LocalIndexType eqIdx) const const BoundaryInteractionVolumeSeed& boundarySeed(const SubControlVolumeFace& scvf) const
{ { return boundarySeeds_[scvfIndexMap_[scvf.index()]]; }
assert(boundarySeeds_[scvfIndexMap_[scvf.index()]].size() > eqIdx);
return boundarySeeds_[scvfIndexMap_[scvf.index()]][eqIdx];
}
private: private:
void initializeBoundarySeeds_() void initializeBoundarySeeds_()
{ {
auto numBoundaryScvf = problem_().model().globalFvGeometry().numBoundaryScvf(); boundarySeeds_.reserve(problem_().model().globalFvGeometry().numBoundaryScvf());
boundarySeeds_.reserve(numBoundaryScvf);
IndexType boundarySeedIndex = 0; IndexType boundarySeedIndex = 0;
for (const auto& element : elements(gridView_)) for (const auto& element : elements(gridView_))
{ {
...@@ -97,19 +92,16 @@ private: ...@@ -97,19 +92,16 @@ private:
if (scvfIndexMap_[scvf.index()] != -1 || !scvf.boundary()) if (scvfIndexMap_[scvf.index()] != -1 || !scvf.boundary())
continue; continue;
// container to store the interaction volume seeds // the boundary interaction volume seed
std::vector<BoundaryInteractionVolumeSeed> seedVector; auto seed = Helper::makeBoundaryInteractionVolumeSeed(problem_(), element, fvGeometry, scvf);
seedVector.reserve(numEq);
for (LocalIndexType eqIdx = 0; eqIdx < numEq; ++eqIdx)
seedVector.emplace_back(Helper::makeBoundaryInteractionVolumeSeed(problem_(), element, fvGeometry, scvf, eqIdx));
// update the index map entries for the global scv faces in the interaction volume // update the index map entries for the global scv faces in the interaction volume
for (const auto& localScvfSeed : seedVector[0].scvfSeeds()) for (const auto& localScvfSeed : seed.scvfSeeds())
for (const auto scvfIdxGlobal : localScvfSeed.globalScvfIndices()) for (const auto scvfIdxGlobal : localScvfSeed.globalScvfIndices())
scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex; scvfIndexMap_[scvfIdxGlobal] = boundarySeedIndex;
// store interaction volume and increment counter // store interaction volume and increment counter
boundarySeeds_.emplace_back(std::move(seedVector)); boundarySeeds_.emplace_back(std::move(seed));
boundarySeedIndex++; boundarySeedIndex++;
} }
} }
...@@ -120,13 +112,13 @@ private: ...@@ -120,13 +112,13 @@ private:
void initializeInteriorSeeds_() void initializeInteriorSeeds_()
{ {
auto numScvf = problem_().model().globalFvGeometry().numScvf(); const auto& globalFvGeometry = problem_().model().globalFvGeometry();
seeds_.reserve(numScvf); seeds_.reserve(globalFvGeometry.numScvf() - globalFvGeometry.numBoundaryScvf());
IndexType seedIndex = 0; IndexType seedIndex = 0;
for (const auto& element : elements(gridView_)) for (const auto& element : elements(gridView_))
{ {
auto fvGeometry = localView(problem_().model().globalFvGeometry()); auto fvGeometry = localView(globalFvGeometry);
fvGeometry.bind(element); fvGeometry.bind(element);
for (const auto& scvf : scvfs(fvGeometry)) for (const auto& scvf : scvfs(fvGeometry))
{ {
...@@ -158,7 +150,7 @@ private: ...@@ -158,7 +150,7 @@ private:
GridView gridView_; GridView gridView_;
std::vector<IndexType> scvfIndexMap_; std::vector<IndexType> scvfIndexMap_;
std::vector<InteractionVolumeSeed> seeds_; std::vector<InteractionVolumeSeed> seeds_;
std::vector<std::vector<BoundaryInteractionVolumeSeed>> boundarySeeds_; std::vector<BoundaryInteractionVolumeSeed> boundarySeeds_;
}; };
} // end namespace } // end namespace
......
...@@ -95,25 +95,6 @@ public: ...@@ -95,25 +95,6 @@ public:
volumeVariables_[scvf.outsideScvIdx()].update(dirichletPriVars, problem, element, insideScv); volumeVariables_[scvf.outsideScvIdx()].update(dirichletPriVars, problem, element, insideScv);
} }
else if (bcTypes.hasDirichlet())
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto dirichletPriVars = problem.dirichlet(element, scvf);
PrimaryVariables priVars(0.0);
for (unsigned int eqIdx = 0; eqIdx < GET_PROP_VALUE(TypeTag, NumEq); ++eqIdx)
{
auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
if (bcTypes.isDirichlet(eqIdx))
priVars[pvIdx] = dirichletPriVars[pvIdx];
else
priVars[pvIdx] = volumeVariables_[scvf.insideScvIdx()].priVar(pvIdx);
}
volumeVariables_[scvf.outsideScvIdx()].update(priVars, problem, element, insideScv);
}
else else
{ {
volumeVariables_[scvf.outsideScvIdx()] = volumeVariables_[scvf.insideScvIdx()]; volumeVariables_[scvf.outsideScvIdx()] = volumeVariables_[scvf.insideScvIdx()];
......
...@@ -245,20 +245,22 @@ public: ...@@ -245,20 +245,22 @@ public:
// Returns the MpfaFaceType of an scv face // Returns the MpfaFaceType of an scv face
static MpfaFaceTypes getMpfaFaceType(const Problem& problem, static MpfaFaceTypes getMpfaFaceType(const Problem& problem,
const Element& element, const Element& element,
const SubControlVolumeFace& scvf, const SubControlVolumeFace& scvf)
const LocalIndexType eqIdx)
{ {
if (!scvf.boundary()) if (!scvf.boundary())
return MpfaFaceTypes::interior; return MpfaFaceTypes::interior;
auto bcTypes = problem.boundaryTypes(element, scvf); auto bcTypes = problem.boundaryTypes(element, scvf);
if (bcTypes.isNeumann(eqIdx)) if (bcTypes.hasOnlyNeumann())
return MpfaFaceTypes::neumann; return MpfaFaceTypes::neumann;
if (bcTypes.isDirichlet(eqIdx)) if (bcTypes.hasOnlyDirichlet())
return MpfaFaceTypes::dirichlet; return MpfaFaceTypes::dirichlet;
if (bcTypes.isOutflow(eqIdx)) // throw for outflow or mixed boundary conditions
if (bcTypes.hasOutflow())
DUNE_THROW(Dune::NotImplemented, "outflow BC for mpfa schemes"); DUNE_THROW(Dune::NotImplemented, "outflow BC for mpfa schemes");
if (bcTypes.hasDirichlet() && bcTypes.hasNeumann())
DUNE_THROW(Dune::InvalidStateException, "Mixed BC are not allowed for cellcentered schemes");
DUNE_THROW(Dune::InvalidStateException, "unknown boundary condition type"); DUNE_THROW(Dune::InvalidStateException, "unknown boundary condition type");
} }
......
...@@ -71,20 +71,19 @@ public: ...@@ -71,20 +71,19 @@ public:
std::vector<ScvSeed> scvSeeds; std::vector<ScvSeed> scvSeeds;
std::vector<ScvfSeed> scvfSeeds; std::vector<ScvfSeed> scvfSeeds;
fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, /*dummy*/0); fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf);
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false); return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), false);
} }
static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem, static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
const Element& element, const Element& element,
const FVElementGeometry& fvGeometry, const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf, const SubControlVolumeFace& scvf)
const LocalIndexType eqIdx)
{ {
std::vector<ScvSeed> scvSeeds; std::vector<ScvSeed> scvSeeds;
std::vector<ScvfSeed> scvfSeeds; std::vector<ScvfSeed> scvfSeeds;
fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, eqIdx); fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf);
return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true); return InteractionVolumeSeed(std::move(scvSeeds), std::move(scvfSeeds), true);
} }
...@@ -94,8 +93,7 @@ private: ...@@ -94,8 +93,7 @@ private:
const Problem& problem, const Problem& problem,
const Element& element, const Element& element,
const FVElementGeometry& fvGeometry, const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf, const SubControlVolumeFace& scvf)
const LocalIndexType eqIdx)
{ {
// Check whether or not we are touching the boundary here // Check whether or not we are touching the boundary here
bool onBoundary = problem.model().globalFvGeometry().scvfTouchesBoundary(scvf); bool onBoundary = problem.model().globalFvGeometry().scvfTouchesBoundary(scvf);
...@@ -107,7 +105,7 @@ private: ...@@ -107,7 +105,7 @@ private:
auto scvIdx0 = scvf.insideScvIdx(); auto scvIdx0 = scvf.insideScvIdx();
// rotate counter clockwise and create the entities // rotate counter clockwise and create the entities
performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, eqIdx); performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0);
if (onBoundary) if (onBoundary)
{ {
...@@ -115,7 +113,7 @@ private: ...@@ -115,7 +113,7 @@ private:
LocalIndexType storeIdx = scvfSeeds.size(); LocalIndexType storeIdx = scvfSeeds.size();
// clockwise rotation until hitting the boundary again // clockwise rotation until hitting the boundary again
performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, eqIdx, /*clockwise*/true); performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, /*clockwise*/true);
// Finish by creating the first scv // Finish by creating the first scv
scvSeeds.emplace(scvSeeds.begin(), ScvSeed(GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}), scvSeeds.emplace(scvSeeds.begin(), ScvSeed(GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
...@@ -136,7 +134,6 @@ private: ...@@ -136,7 +134,6 @@ private:
std::vector<ScvSeed>& scvSeeds, std::vector<ScvSeed>& scvSeeds,
std::vector<ScvfSeed>& scvfSeeds, std::vector<ScvfSeed>& scvfSeeds,
const GlobalIndexType scvIdx0, const GlobalIndexType scvIdx0,
const LocalIndexType eqIdx,
const bool clockWise = false) const bool clockWise = false)
{ {
// extract the actual local indices from the containers // extract the actual local indices from the containers
...@@ -162,7 +159,7 @@ private: ...@@ -162,7 +159,7 @@ private:
// the current element inside of the scv face // the current element inside of the scv face
auto insideElement = problem.model().globalFvGeometry().element(insideGlobalScvIdx); auto insideElement = problem.model().globalFvGeometry().element(insideGlobalScvIdx);
auto faceType = Implementation::getMpfaFaceType(problem, insideElement, curScvf, eqIdx); auto faceType = Implementation::getMpfaFaceType(problem, insideElement, curScvf);
// if the face touches the boundary, create a boundary scvf entity // if the face touches the boundary, create a boundary scvf entity
if (curScvf.boundary()) if (curScvf.boundary())
...@@ -287,8 +284,7 @@ public: ...@@ -287,8 +284,7 @@ public:
static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem, static InteractionVolumeSeed makeBoundaryInteractionVolumeSeed(const Problem& problem,
const Element& element, const Element& element,
const FVElementGeometry& fvGeometry, const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf, const SubControlVolumeFace& scvf)
const LocalIndexType eqIdx)
{ {
std::vector<ScvSeed> scvSeeds; std::vector<ScvSeed> scvSeeds;
std::vector<ScvfSeed> scvfSeeds; std::vector<ScvfSeed> scvfSeeds;
...@@ -301,7 +297,7 @@ public: ...@@ -301,7 +297,7 @@ public:
auto vIdxGlobal = scvf.vertexIndex(); auto vIdxGlobal = scvf.vertexIndex();
// create the scv entity seeds // create the scv entity seeds
fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, eqIdx, vIdxGlobal); fillEntitySeeds_(scvSeeds, scvfSeeds, problem, element, fvGeometry, scvf, vIdxGlobal);
// shrink containers to necessary size // shrink containers to necessary size
scvSeeds.shrink_to_fit(); scvSeeds.shrink_to_fit();
...@@ -317,7 +313,6 @@ private: ...@@ -317,7 +313,6 @@ private:
const Element& element, const Element& element,
const FVElementGeometry& fvGeometry, const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf, const SubControlVolumeFace& scvf,
const LocalIndexType eqIdx,
const GlobalIndexType vIdxGlobal) const GlobalIndexType vIdxGlobal)
{ {
// Get the three scv faces in the scv // Get the three scv faces in the scv
...@@ -347,7 +342,7 @@ private: ...@@ -347,7 +342,7 @@ private:
actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size()); actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size());
// create the scvf seed // create the scvf seed
auto faceType = Implementation::getMpfaFaceType(problem, element, *scvfVector[coordDir], eqIdx); auto faceType = Implementation::getMpfaFaceType(problem, element, *scvfVector[coordDir]);
scvfSeeds.emplace_back( actualScvf, scvfSeeds.emplace_back( actualScvf,
LocalIndexSet({actualLocalScvIdx}), LocalIndexSet({actualLocalScvIdx}),
GlobalIndexSet({scvfVector[coordDir]->index()}), GlobalIndexSet({scvfVector[coordDir]->index()}),
...@@ -388,13 +383,13 @@ private: ...@@ -388,13 +383,13 @@ private:
const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
// create scvf seed // create scvf seed
auto faceType = Implementation::getMpfaFaceType(problem, element, actualScvf, eqIdx); auto faceType = Implementation::getMpfaFaceType(problem, element, actualScvf);
LocalIndexSet scvIndicesLocal({actualLocalScvIdx, static_cast<LocalIndexType>(scvSeeds.size())}); LocalIndexSet scvIndicesLocal({actualLocalScvIdx, static_cast<LocalIndexType>(scvSeeds.size())});
GlobalIndexSet scvfIndicesGlobal({globalScvfIndex, outsideScvf.index()}); GlobalIndexSet scvfIndicesGlobal({globalScvfIndex, outsideScvf.index()});
scvfSeeds.emplace_back(actualScvf, std::move(scvIndicesLocal), std::move(scvfIndicesGlobal), faceType); scvfSeeds.emplace_back(actualScvf, std::move(scvIndicesLocal), std::move(scvfIndicesGlobal), faceType);
// make outside scv by recursion // make outside scv by recursion
fillEntitySeeds_(scvSeeds, scvfSeeds, problem, outsideElement, outsideFvGeometry, outsideScvf, eqIdx, vIdxGlobal); fillEntitySeeds_(scvSeeds, scvfSeeds, problem, outsideElement, outsideFvGeometry, outsideScvf, vIdxGlobal);
} }
// we have to find out if it is necessary to create a new scvf // we have to find out if it is necessary to create a new scvf
else else
...@@ -434,7 +429,7 @@ private: ...@@ -434,7 +429,7 @@ private:
const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx]; const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
// some data on the face // some data on the face
auto faceType = Implementation::getMpfaFaceType(problem, element, *scvfVector[coordDir], eqIdx); auto faceType = Implementation::getMpfaFaceType(problem, element, *scvfVector[coordDir]);
LocalIndexSet scvIndicesLocal( {actualLocalScvIdx, outsideLocalScvIdx} ); LocalIndexSet scvIndicesLocal( {actualLocalScvIdx, outsideLocalScvIdx} );
GlobalIndexSet scvfIndicesGlobal( {globalScvfIndex, outsideScvf.index()} ); GlobalIndexSet scvfIndicesGlobal( {globalScvfIndex, outsideScvf.index()} );
scvfSeeds.emplace_back(*scvfVector[coordDir], std::move(scvIndicesLocal), std::move(scvfIndicesGlobal), faceType); scvfSeeds.emplace_back(*scvfVector[coordDir], std::move(scvIndicesLocal), std::move(scvfIndicesGlobal), faceType);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment