Commit 539bcad8 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

port mpfa to new param interface and fix some issues

parent c407e15b
......@@ -94,6 +94,7 @@ class CCMpfaElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
......@@ -131,7 +132,10 @@ public:
// update the volume variables of the element at hand
auto eIdx = problem.elementMapper().index(element);
auto&& scvI = fvGeometry.scv(eIdx);
volumeVariables_[localIdx].update(sol[eIdx], problem, element, scvI);
volumeVariables_[localIdx].update(problem.model().elementSolution(element, sol),
problem,
element,
scvI);
volVarIndices_[localIdx] = scvI.index();
++localIdx;
......@@ -140,7 +144,10 @@ public:
{
const auto& elementJ = globalFvGeometry.element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
volumeVariables_[localIdx].update(sol[globalJ], problem, elementJ, scvJ);
volumeVariables_[localIdx].update(problem.model().elementSolution(elementJ, sol),
problem,
elementJ,
scvJ);
volVarIndices_[localIdx] = scvJ.index();
++localIdx;
}
......@@ -166,8 +173,10 @@ public:
{
// boundary volume variables
VolumeVariables dirichletVolVars;
const auto dirichletPriVars = problem.dirichlet(element, scvf);
dirichletVolVars.update(dirichletPriVars, problem, element, scvI);
dirichletVolVars.update(ElementSolution({problem.dirichlet(element, scvf)}),
problem,
element,
scvI);
volumeVariables_.emplace_back(std::move(dirichletVolVars));
volVarIndices_.push_back(scvf.outsideScvIdx());
......@@ -206,8 +215,10 @@ public:
// boundary volume variables
VolumeVariables dirichletVolVars;
auto&& ivScv = fvGeometry.scv(insideScvIdx);
const auto dirichletPriVars = problem.dirichlet(insideElement, ivScvf);
dirichletVolVars.update(dirichletPriVars, problem, insideElement, ivScv);
dirichletVolVars.update(ElementSolution({problem.dirichlet(insideElement, ivScvf)}),
problem,
insideElement,
ivScv);
volumeVariables_.emplace_back(std::move(dirichletVolVars));
volVarIndices_.push_back(ivScvf.outsideScvIdx());
......@@ -234,7 +245,10 @@ public:
// update the volume variables of the element
auto&& scv = fvGeometry.scv(eIdx);
volumeVariables_[0].update(sol[eIdx], globalVolVars().problem_(), element, scv);
volumeVariables_[0].update(globalVolVars().problem_().model().elementSolution(element, sol),
globalVolVars().problem_(),
element,
scv);
volVarIndices_[0] = scv.index();
}
......
......@@ -52,6 +52,7 @@ class CCMpfaAdvectionCacheFiller
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using Element = typename GridView::template Codim<0>::Entity;
......@@ -73,7 +74,7 @@ public:
auto getK = [&problem](const Element& element,
const VolumeVariables& volVars,
const SubControlVolume& scv)
{ return problem.spatialParams().intrinsicPermeability(scv, volVars); };
{ return volVars.permeability(); };
// update flux var caches
if (problem.model().globalFvGeometry().scvfTouchesBoundary(scvf))
......@@ -332,11 +333,10 @@ public:
//! set the neumann boundary conditions in case we do not use tpfa on the boundary
if (!useTpfaBoundary)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
iv.assembleNeumannFluxes(energyEqIdx);
cache.updateHeatNeumannFlux(scvf, iv, energyEqIdx);
}
{
iv.assembleNeumannFluxes(energyEqIdx);
cache.updateHeatNeumannFlux(scvf, iv);
}
for (const auto scvfIdxJ : iv.globalScvfs())
{
......@@ -349,11 +349,10 @@ public:
//! set the neumann boundary conditions in case we do not use tpfa on the boundary
if (!useTpfaBoundary)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
iv.assembleNeumannFluxes(energyEqIdx);
cache.updateHeatNeumannFlux(scvf, iv, energyEqIdx);
}
{
iv.assembleNeumannFluxes(energyEqIdx);
cache.updateHeatNeumannFlux(scvf, iv);
}
}
}
}
......@@ -371,14 +370,6 @@ public:
cache.updateHeatConduction(scvf, iv);
cache.setUpdateStatus(true);
//! set the neumann boundary conditions in case we do not use tpfa on the boundary
if (!useTpfaBoundary)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
iv.assembleNeumannFluxes(energyEqIdx);
cache.updateHeatNeumannFlux(scvf, iv, energyEqIdx);
}
for (const auto scvfIdxJ : iv.globalScvfs())
{
if (scvfIdxJ != scvfIdx)
......@@ -387,14 +378,6 @@ public:
auto& cacheJ = fluxVarsCacheContainer[scvfIdxJ];
cacheJ.updateHeatConduction(scvfJ, iv);
cacheJ.setUpdateStatus(true);
//! set the neumann boundary conditions in case we do not use tpfa on the boundary
if (!useTpfaBoundary)
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
iv.assembleNeumannFluxes(energyEqIdx);
cache.updateHeatNeumannFlux(scvf, iv, energyEqIdx);
}
}
}
}
......
......@@ -239,7 +239,7 @@ public:
continue;
// store info on which vertices are on the domain boundary
if (boundary)
if (boundary && !boundaryVertices_[vIdxGlobal])
{
boundaryVertices_[vIdxGlobal] = true;
numBoundaryVertices_++;
......@@ -573,7 +573,7 @@ public:
continue;
// store info on which vertices are on the domain boundary
if (boundary)
if (boundary && !boundaryVertices_[vIdxGlobal])
{
boundaryVertices_[vIdxGlobal] = true;
numBoundaryVertices_++;
......
......@@ -64,8 +64,10 @@ public:
// reserve memory
const auto numScvf = this->problem().model().globalFvGeometry().numScvf();
const auto numBoundaryScvf = this->problem().model().globalFvGeometry().numBoundaryScvf();
const int numInteriorScvf = (numScvf-numBoundaryScvf)/2;
seeds.reserve( std::size_t((numScvf-numBoundaryScvf)/2) );
if (numInteriorScvf > 0)
seeds.reserve( numInteriorScvf );
boundarySeeds.reserve(numBoundaryScvf);
scvfIndexMap.resize(numScvf);
......
......@@ -67,7 +67,7 @@ public:
// reserve memory
const auto numScvf = this->problem().model().globalFvGeometry().numScvf();
const auto numBoundaryVertices = this->problem().model().globalFvGeometry().numBoundaryVertices();
const auto numInteriorVertices = this->gridView().size(dim) - numBoundaryVertices;
const int numInteriorVertices = this->gridView().size(dim) - numBoundaryVertices;
if (numInteriorVertices > 0)
seeds.reserve(numInteriorVertices);
......
......@@ -432,12 +432,13 @@ private:
auto posLocalScvIdx = localScvf.insideLocalScvIndex();
auto&& posLocalScv = localScv_(posLocalScvIdx);
auto&& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
auto&& posVolVars = elemVolVars_()[posGlobalScv];
auto element = localElement_(posLocalScvIdx);
auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv);
auto tensor = getTensor(element, posVolVars, posGlobalScv);
// the omega factors of the "positive" sub volume
auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
posWijk *= problem_().boxExtrusionFactor(element, posGlobalScv);
posWijk *= posVolVars.extrusionFactor();
// Check the local directions of the positive sub volume
for (int localDir = 0; localDir < dim; localDir++)
......@@ -484,8 +485,9 @@ private:
{
auto&& negLocalScv = localScv_(negLocalScvIdx);
auto&& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
auto&& negVolVars = elemVolVars_()[negGlobalScv];
auto negElement = localElement_(negLocalScvIdx);
auto negTensor = getTensor(negElement, elemVolVars_()[negGlobalScv], negGlobalScv);
auto negTensor = getTensor(negElement, negVolVars, negGlobalScv);
// the omega factors of the "negative" sub volume
DimVector negWijk;
......@@ -503,7 +505,7 @@ private:
negWijk = calculateOmegas_(negLocalScv, localScvf.unitOuterNormal(), localScvf.area(), negTensor);
// scale by extrusion factpr
negWijk *= problem_().boxExtrusionFactor(negElement, negGlobalScv);
negWijk *= negVolVars.extrusionFactor();
// Check local directions of negative sub volume
for (int localDir = 0; localDir < dim; localDir++)
......@@ -564,12 +566,13 @@ private:
auto posLocalScvIdx = localScvf.insideLocalScvIndex();
auto&& posLocalScv = localScv_(posLocalScvIdx);
auto&& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
auto&& posVolVars = elemVolVars_()[posGlobalScv];
auto element = localElement_(posLocalScvIdx);
auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv);
auto tensor = getTensor(element, posVolVars, posGlobalScv);
// the omega factors of the "positive" sub volume
auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
posWijk *= problem_().boxExtrusionFactor(element, posGlobalScv);
posWijk *= posVolVars.extrusionFactor();
for (int localDir = 0; localDir < dim; localDir++)
{
......
......@@ -311,6 +311,15 @@ public:
heatConductionTij_ = interactionVolume.getTransmissibilities(localFaceData);
}
// update cached objects for heat conduction
template<typename InteractionVolume>
void updateHeatNeumannFlux(const SubControlVolumeFace &scvf,
const InteractionVolume& interactionVolume)
{
const auto& localFaceData = interactionVolume.getLocalFaceData(scvf);
heatNeumannFlux_ = interactionVolume.getNeumannFlux(localFaceData);
}
//! Returns the volume variables indices necessary for heat conduction flux
//! computation. This includes all participating boundary volume variables
//! and it can be different for the phases & components.
......
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