diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
index 47bb3ddedb8364a1eb41fc426fce6257b923a0da..b4b6b14e98cd5f8c76caac992ca38e8939b1a7bb 100644
--- a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
@@ -137,7 +137,7 @@ public:
void bindElement(const Element& element)
{
elementPtr_ = &element;
- scvIndices_ = std::vector<IndexType>({fvGridGeometry().problem_().elementMapper().index(*elementPtr_)});
+ scvIndices_ = std::vector<IndexType>({fvGridGeometry().elementMapper().index(*elementPtr_)});
}
//! The global finite volume geometry we are a restriction of
@@ -299,21 +299,20 @@ public:
}
}
- //! Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
+ //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
//! on additional DOFs not included in the discretization schemes' occupation pattern
- const auto& problem = fvGridGeometry().problem_();
- const auto globalI = problem.elementMapper().index(element);
- const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
- if (!additionalDofDependencies.empty())
- {
- neighborScvs_.reserve(neighborScvs_.size() + additionalDofDependencies.size());
- neighborScvIndices_.reserve(neighborScvIndices_.size() + additionalDofDependencies.size());
- for (auto globalJ : additionalDofDependencies)
- {
- neighborScvs_.emplace_back(fvGridGeometry().element(globalJ).geometry(), globalJ);
- neighborScvIndices_.emplace_back(globalJ);
- }
- }
+ // const auto globalI = fvGridGeometry().elementMapper().index(element);
+ // const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
+ // if (!additionalDofDependencies.empty())
+ // {
+ // neighborScvs_.reserve(neighborScvs_.size() + additionalDofDependencies.size());
+ // neighborScvIndices_.reserve(neighborScvIndices_.size() + additionalDofDependencies.size());
+ // for (auto globalJ : additionalDofDependencies)
+ // {
+ // neighborScvs_.emplace_back(fvGridGeometry().element(globalJ).geometry(), globalJ);
+ // neighborScvIndices_.emplace_back(globalJ);
+ // }
+ // }
}
//! Binding of an element preparing the geometries only inside the element
@@ -361,7 +360,7 @@ private:
//! create scvs and scvfs of the bound element
void makeElementGeometries(const Element& element)
{
- auto eIdx = fvGridGeometry().problem_().elementMapper().index(element);
+ const auto eIdx = fvGridGeometry().elementMapper().index(element);
scvs_.emplace_back(element.geometry(), eIdx);
scvIndices_.emplace_back(eIdx);
@@ -378,8 +377,8 @@ private:
int scvfCounter = 0;
for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
{
- // check if intersection is on interior boundary
- const auto isInteriorBoundary = fvGridGeometry().problem_().isInteriorBoundary(element, intersection);
+ // TODO check if intersection is on interior boundary
+ const auto isInteriorBoundary = false;
if (dim < dimWorld)
if (handledScvf[intersection.indexInInside()])
@@ -408,7 +407,7 @@ private:
void makeNeighborGeometries(const Element& element)
{
// create the neighbor scv
- auto eIdx = fvGridGeometry().problem_().elementMapper().index(element);
+ const auto eIdx = fvGridGeometry().elementMapper().index(element);
neighborScvs_.emplace_back(element.geometry(), eIdx);
neighborScvIndices_.push_back(eIdx);
@@ -425,8 +424,8 @@ private:
int scvfCounter = 0;
for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
{
- // check if intersection is on interior boundary
- const auto isInteriorBoundary = fvGridGeometry().problem_().isInteriorBoundary(element, intersection);
+ // TODO check if intersection is on interior boundary
+ const auto isInteriorBoundary = false;
if (dim < dimWorld)
if (handledScvf[intersection.indexInInside()])
@@ -458,7 +457,7 @@ private:
{
for (unsigned outsideScvIdx = 0; outsideScvIdx < neighborVolVarIndices[scvfCounter].size(); ++outsideScvIdx)
{
- if (neighborVolVarIndices[scvfCounter][outsideScvIdx] == fvGridGeometry().problem_().elementMapper().index(*elementPtr_))
+ if (neighborVolVarIndices[scvfCounter][outsideScvIdx] == fvGridGeometry().elementMapper().index(*elementPtr_))
{
std::vector<IndexType> scvIndices({eIdx});
scvIndices.insert(scvIndices.end(), neighborVolVarIndices[scvfCounter].begin(), neighborVolVarIndices[scvfCounter].end());
diff --git a/dumux/implicit/cellcentered/assembler.hh b/dumux/implicit/cellcentered/assembler.hh
index ae80971092189c22109c76ef1d49559b62a0baa9..c00c49716ac61cad79389c79aecb7a72297ab689 100644
--- a/dumux/implicit/cellcentered/assembler.hh
+++ b/dumux/implicit/cellcentered/assembler.hh
@@ -32,15 +32,21 @@
namespace Dumux {
+template<class TypeTag, DiscretizationMethods DM>
+class ImplicitAssemblerImplementation;
+
+template<class TypeTag>
+class CCImplicitAssembler;
+
//! TPFA and MPFA use the same assembler class
template<class TypeTag>
-class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCTpfa> : CCImplicitAssembler<TypeTag>
+class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCTpfa> : public CCImplicitAssembler<TypeTag>
{
using CCImplicitAssembler<TypeTag>::CCImplicitAssembler;
};
template<class TypeTag>
-class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCMpfa> : CCImplicitAssembler<TypeTag>
+class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCMpfa> : public CCImplicitAssembler<TypeTag>
{
using CCImplicitAssembler<TypeTag>::CCImplicitAssembler;
};
@@ -54,21 +60,23 @@ template<class TypeTag>
class CCImplicitAssembler
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using AssemblyMap = typename GET_PROP_TYPE(TypeTag, AssemblyMap);
using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
- using GridFvGeometry = typename GET_PROP_TYPE(TypeTag, GridFvGeometry);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using LocalAssembler = typename GET_PROP_TYPE(TypeTag, LocalAssembler);
+ // using LocalAssembler = typename GET_PROP_TYPE(TypeTag, LocalAssembler);
+ using LocalAssembler = CCImplicitLocalAssembler<TypeTag, DifferentiationMethods::numeric>;
public:
using ResidualType = SolutionVector;
//! The constructor
CCImplicitAssembler(std::shared_ptr<const Problem> problem,
- std::shared_ptr<const GridFvGeometry> gridFvGeometry,
+ std::shared_ptr<const FVGridGeometry> gridFvGeometry,
std::shared_ptr<const SolutionVector> prevSol,
std::shared_ptr<const SolutionVector> curSol,
std::shared_ptr<GridVariables> gridVariables)
@@ -79,7 +87,7 @@ public:
, gridVariables_(gridVariables)
{
// build assembly map
- assemblyMap_.init(globalFvGeometry);
+ assemblyMap_.init(gridFvGeometry);
}
/*!
@@ -100,18 +108,18 @@ public:
// if we get here, everything worked well
succeeded = true;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
}
// throw exception if a problem ocurred
catch (NumericalProblem &e)
{
- std::cout << "rank " << problem_().gridView().comm().rank()
+ std::cout << "rank " << gridView().comm().rank()
<< " caught an exception while assembling:" << e.what()
<< "\n";
succeeded = false;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
}
if (!succeeded)
DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
@@ -134,8 +142,8 @@ public:
// if we get here, everything worked well
succeeded = true;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
}
// throw exception if a problem ocurred
catch (NumericalProblem &e)
@@ -144,8 +152,8 @@ public:
<< " caught an exception while assembling:" << e.what()
<< "\n";
succeeded = false;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
}
if (!succeeded)
DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
@@ -229,7 +237,7 @@ public:
// calculate the square norm of the residual
Scalar result2 = residual.two_norm2();
if (gridView().comm().size() > 1)
- result2 = gridView_().comm().sum(result2);
+ result2 = gridView().comm().sum(result2);
using std::sqrt;
return sqrt(result2);
@@ -238,7 +246,7 @@ public:
const Problem& problem() const
{ return *problem; }
- const GridFvGeometry& gridFvGeometry() const
+ const FVGridGeometry& gridFvGeometry() const
{ return *gridFvGeometry_; }
const GridView& gridView() const
@@ -302,7 +310,7 @@ private:
std::shared_ptr<const Problem> problem_;
// the finite volume geometry of the grid
- std::shared_ptr<const GridFvGeometry> gridFvGeometry_;
+ std::shared_ptr<const FVGridGeometry> gridFvGeometry_;
// previous and current solution to the problem
std::shared_ptr<const SolutionVector> prevSol_;
diff --git a/dumux/implicit/cellcentered/assemblymap.hh b/dumux/implicit/cellcentered/assemblymap.hh
index 78ad0e17cf0c7acdb8f056fd4fc184baf5b35079..597e3fa0c308e68de0f7f9c267abfaf74101f60b 100644
--- a/dumux/implicit/cellcentered/assemblymap.hh
+++ b/dumux/implicit/cellcentered/assemblymap.hh
@@ -47,7 +47,7 @@ class CCSimpleAssemblyMap
{
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
-
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
struct DataJ
diff --git a/dumux/implicit/cellcentered/localassembler.hh b/dumux/implicit/cellcentered/localassembler.hh
index a3461b54bed62a95c761848662a145ab1bc73bc7..04c856284cb613f0685fc0eb6b614c09e2e760fe 100644
--- a/dumux/implicit/cellcentered/localassembler.hh
+++ b/dumux/implicit/cellcentered/localassembler.hh
@@ -30,6 +30,11 @@
namespace Dumux {
+enum DifferentiationMethods
+{
+ numeric
+};
+
/*!
* \ingroup ImplicitModel
* \brief An assembler for the local contributions (per element) to the global
@@ -47,6 +52,10 @@ class CCImplicitLocalAssembler<TypeTag, DifferentiationMethods::numeric>
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalAssembler);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
public:
@@ -57,6 +66,7 @@ public:
template<class Assembler>
static void assemble(Assembler& assembler, SolutionVector& res, const Element& element)
{
+ const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
res[globalI] = Implementation::assemble_(assembler, element);
}
@@ -76,7 +86,7 @@ public:
*
* \param priVar The value of the primary variable
*/
- static Scalar numericEpsilon(const Scalar priVar) const
+ static Scalar numericEpsilon(const Scalar priVar)
{
// define the base epsilon as the geometric mean of 1 and the
// resolution of the scalar type. E.g. for standard 64 bit
@@ -115,7 +125,7 @@ private:
auto& A = assembler.matrix();
// prepare the local views
- auto fvGeometry = localView(globalFvGeometry());
+ auto fvGeometry = localView(assembler.fvGridGeometry());
fvGeometry.bind(element);
auto curElemVolVars = localView(gridVariables.curGlobalVolVars());
@@ -251,7 +261,7 @@ private:
neighborDeriv = origFlux;
}
- if (numericDifferenceMethod_ <= 0)
+ if (numericDifferenceMethod <= 0)
{
// we are not using forward differences, i.e. we
// need to calculate f(x - \epsilon)
@@ -266,7 +276,8 @@ private:
// calculate the residual with the deflected primary variables and subtract it
if (!isGhost)
- partialDeriv -= localResidual.eval(element,
+ partialDeriv -= localResidual.eval(problem,
+ element,
fvGeometry,
prevElemVolVars,
curElemVolVars,
@@ -275,12 +286,13 @@ private:
// calculate the fluxes into element with the deflected primary variables
for (std::size_t k = 0; k < numNeighbors; ++k)
- for (auto scvfIdx : assemblyMap_[globalI][k].scvfsJ)
- neighborDeriv[k] -= this->localResidual().evalFlux_(neighborElements[k],
- fvGeometry,
- curElemVolVars,
- fvGeometry.scvf(scvfIdx),
- elemFluxVarsCache);
+ for (auto scvfIdx : assemblyMap[globalI][k].scvfsJ)
+ neighborDeriv[k] -= localResidual.evalFlux(problem,
+ neighborElements[k],
+ fvGeometry,
+ curElemVolVars,
+ fvGeometry.scvf(scvfIdx),
+ elemFluxVarsCache);
}
else
{
@@ -312,7 +324,7 @@ private:
// off-diagonal entries
j = 0;
- for (const auto& dataJ : assemblyMap_[globalI])
+ for (const auto& dataJ : assemblyMap[globalI])
A[dataJ.globalJ][globalI][eqIdx][pvIdx] += neighborDeriv[j++][pvIdx];
}
}
@@ -325,101 +337,103 @@ private:
// //
//////////////////////////////////////////////////////////////////////////////////////////////
- const auto& additionalDofDepedencies = problem.getAdditionalDofDependencies(globalI);
- if (!additionalDofDepedencies.empty() && !isGhost)
- {
- // compute the source in the undeflected state
- auto source = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
- source *= -scv.volume()*curVolVarsI.extrusionFactor();
-
- // deflect solution at given dofs and recalculate the source
- for (auto globalJ : additionalDofDependencies)
- {
- const auto& scvJ = fvGeometry.scv(globalJ);
- auto& curVolVarsJ = curElemVolVars[scv];
- const auto& elementJ = gridFvGeometry.element(globalJ);
-
- // save a copy of the original privars and volvars
- // to restore original solution after deflection
- const auto origPriVars = curSol[globalJ];
- const auto origVolVarsJ = curVolVarsJ;
-
- // derivatives with repect to the additional DOF we depend on
- for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
- {
- // derivatives of element dof with respect to itself
- NumEqVector partialDeriv(0.0);
- const auto eps = Implementation::numericEpsilon(curVolVarsJ.priVar(pvIdx));
- Scalar delta = 0;
-
- if (numericDifferenceMethod_ >= 0)
- {
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- curSol[globalJ][pvIdx] += eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
-
- // calculate the source with the deflected primary variables
- auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
- deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
- partialDeriv = std::move(deflSource);
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) source f(x)
- partialDeriv = source;
- }
-
- if (numericDifferenceMethod_ <= 0)
- {
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- curSol[globalJ][pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
-
- // calculate the source with the deflected primary variables and subtract
- auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
- deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
- partialDeriv -= std::move(deflSource);
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) source f(x)
- partialDeriv -= source;
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- partialDeriv /= delta;
-
- // restore the original state of the dofs privars and the volume variables
- curSol[globalJ] = origPriVars;
- curVolVarsJ = origVolVarsJ;
-
- // add the current partial derivatives to the global jacobian matrix
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
- }
- }
- }
+ // const auto& additionalDofDepedencies = problem.getAdditionalDofDependencies(globalI);
+ // if (!additionalDofDepedencies.empty() && !isGhost)
+ // {
+ // // compute the source in the undeflected state
+ // auto source = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
+ // source *= -scv.volume()*curVolVarsI.extrusionFactor();
+
+ // // deflect solution at given dofs and recalculate the source
+ // for (auto globalJ : additionalDofDependencies)
+ // {
+ // const auto& scvJ = fvGeometry.scv(globalJ);
+ // auto& curVolVarsJ = curElemVolVars[scv];
+ // const auto& elementJ = gridFvGeometry.element(globalJ);
+
+ // // save a copy of the original privars and volvars
+ // // to restore original solution after deflection
+ // const auto origPriVars = curSol[globalJ];
+ // const auto origVolVarsJ = curVolVarsJ;
+
+ // // derivatives with repect to the additional DOF we depend on
+ // for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ // {
+ // // derivatives of element dof with respect to itself
+ // NumEqVector partialDeriv(0.0);
+ // const auto eps = Implementation::numericEpsilon(curVolVarsJ.priVar(pvIdx));
+ // Scalar delta = 0;
+
+ // if (numericDifferenceMethod >= 0)
+ // {
+ // // we are not using backward differences, i.e. we need to
+ // // calculate f(x + \epsilon)
+
+ // // deflect primary variables
+ // curSol[globalJ][pvIdx] += eps;
+ // delta += eps;
+
+ // // update the volume variables and the flux var cache
+ // curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
+
+ // // calculate the source with the deflected primary variables
+ // auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
+ // deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
+ // partialDeriv = std::move(deflSource);
+ // }
+ // else
+ // {
+ // // we are using backward differences, i.e. we don't need
+ // // to calculate f(x + \epsilon) and we can recycle the
+ // // (already calculated) source f(x)
+ // partialDeriv = source;
+ // }
+
+ // if (numericDifferenceMethod <= 0)
+ // {
+ // // we are not using forward differences, i.e. we
+ // // need to calculate f(x - \epsilon)
+
+ // // deflect the primary variables
+ // curSol[globalJ][pvIdx] -= delta + eps;
+ // delta += eps;
+
+ // // update the volume variables and the flux var cache
+ // curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
+
+ // // calculate the source with the deflected primary variables and subtract
+ // auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
+ // deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
+ // partialDeriv -= std::move(deflSource);
+ // }
+ // else
+ // {
+ // // we are using forward differences, i.e. we don't need to
+ // // calculate f(x - \epsilon) and we can recycle the
+ // // (already calculated) source f(x)
+ // partialDeriv -= source;
+ // }
+
+ // // divide difference in residuals by the magnitude of the
+ // // deflections between the two function evaluation
+ // partialDeriv /= delta;
+
+ // // restore the original state of the dofs privars and the volume variables
+ // curSol[globalJ] = origPriVars;
+ // curVolVarsJ = origVolVarsJ;
+
+ // // add the current partial derivatives to the global jacobian matrix
+ // for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ // A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ // }
+ // }
+ // }
// return the original residual
return residual;
}
};
-}
\ No newline at end of file
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh
index 56868fa1d154866060f7ed2e0bcd913663497927..b6d2abb8dcf42c076b6bfe3590113304790344d0 100644
--- a/dumux/implicit/cellcentered/localresidual.hh
+++ b/dumux/implicit/cellcentered/localresidual.hh
@@ -44,262 +44,62 @@ template<class TypeTag>
class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
{
using ParentType = ImplicitLocalResidual<TypeTag>;
- friend class ImplicitLocalResidual<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-
- using Element = typename GridView::template Codim<0>::Entity;
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using ElementResidualVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
public:
- // copying the local residual class is not a good idea
- CCLocalResidual(const CCLocalResidual &) = delete;
-
- CCLocalResidual() : ParentType() {}
-
-protected:
- void evalFluxes_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const ElementBoundaryTypes& bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ void evalFlux(ElementResidualVector& residual,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementBoundaryTypes& elemBcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
{
- // calculate the mass flux over the scv faces
- for (auto&& scvf : scvfs(fvGeometry))
- {
- this->residual_[0] += this->asImp_().computeFlux_(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
- }
- }
+ const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto localScvIdx = scv.indexInElement();
- PrimaryVariables computeFlux_(const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf,
- const ElementFluxVariablesCache& elemFluxVarsCache)
- {
+ // inner faces
if (!scvf.boundary())
- return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
- else
- return PrimaryVariables(0.0);
-
- }
-
- void evalBoundary_(const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const ElementBoundaryTypes& elemBcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
- {
- for (auto&& scvf : scvfs(fvGeometry))
{
- if (scvf.boundary())
- {
- auto bcTypes = this->problem().boundaryTypes(element, scvf);
- this->residual_[0] += evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache);
- }
+ residual[localScvIdx] += this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
}
- // additionally treat mixed D/N conditions in a strong sense
- if (elemBcTypes.hasDirichlet())
+ // boundary faces
+ else
{
- for (auto&& scvf : scvfs(fvGeometry))
- {
- if (scvf.boundary())
- this->asImp_().evalDirichlet_(element, fvGeometry, elemVolVars, scvf);
- }
- }
- }
-
- /*!
- * \brief Add all fluxes resulting from Neumann, outflow and pure Dirichlet
- * boundary conditions to the local residual.
- */
- PrimaryVariables evalBoundaryFluxes_(const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf,
- const BoundaryTypes& bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
- {
- // evaluate the Neumann conditions at the boundary face
- PrimaryVariables flux(0);
- if (bcTypes.hasNeumann())
- flux += this->asImp_().evalNeumannSegment_(element, fvGeometry, elemVolVars, scvf, bcTypes);
-
- // TODO: evaluate the outflow conditions at the boundary face
- //if (bcTypes.hasOutflow())
- // flux += this->asImp_().evalOutflowSegment_(&intersection, bcTypes);
-
- // evaluate the pure Dirichlet conditions at the boundary face
- if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
- flux += this->asImp_().evalDirichletSegment_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache);
-
- return flux;
- }
-
-
- /*!
- * \brief Evaluate Dirichlet conditions on faces that have mixed
- * Dirichlet/Neumann boundary conditions.
- */
- void evalDirichlet_(const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf)
- {
- BoundaryTypes bcTypes = this->problem().boundaryTypes(element, scvf);
-
- if (bcTypes.hasDirichlet() && bcTypes.hasNeumann())
- this->asImp_().evalDirichletSegmentMixed_(element, fvGeometry, elemVolVars, scvf, bcTypes);
- }
+ const auto& bcTypes = elemBcTypes[localScvIdx];
- /*!
- * \brief Add Neumann boundary conditions for a single scv face
- */
- PrimaryVariables evalNeumannSegment_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf,
- const BoundaryTypes &bcTypes)
- {
- // temporary vector to store the neumann boundary fluxes
- PrimaryVariables flux(0);
-
- auto neumannFluxes = this->problem().neumann(element, fvGeometry, elemVolVars, scvf);
-
- // multiply neumann fluxes with the area and the extrusion factor
- auto&& scv = fvGeometry.scv(scvf.insideScvIdx());
- neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
-
- // add fluxes to the residual
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- if (bcTypes.isNeumann(eqIdx))
- flux[eqIdx] += neumannFluxes[eqIdx];
-
- return flux;
- }
-
- /*!
- * \brief Treat Dirichlet boundary conditions in a weak sense for a single
- * intersection that only has Dirichlet boundary conditions
- */
- PrimaryVariables evalDirichletSegment_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf,
- const BoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
- {
- // temporary vector to store the Dirichlet boundary fluxes
- PrimaryVariables flux(0);
+ // Dirichlet boundaries
+ if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
+ residual[localScvIdx] += this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
- auto dirichletFlux = this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
-
- // add fluxes to the residual
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- if (bcTypes.isDirichlet(eqIdx))
- flux[eqIdx] += dirichletFlux[eqIdx];
-
- return flux;
- }
-
- /*!
- * \brief Treat Dirichlet boundary conditions in a strong sense for a
- * single intersection that has mixed D/N boundary conditions
- */
- void evalDirichletSegmentMixed_(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf,
- const BoundaryTypes &bcTypes)
- {
- // temporary vector to store the Dirichlet values
- PrimaryVariables dirichletValues = this->problem().dirichlet(element, scvf);
-
- // get the primary variables
- const auto& priVars = elemVolVars[scvf.insideScvIdx()].priVars();
-
- // set Dirichlet conditions in a strong sense
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (bcTypes.isDirichlet(eqIdx))
+ // Neumann and Robin ("solution dependent Neumann") boundary conditions
+ else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
{
- auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- this->residual_[0][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- }
- }
- }
-
- /*!
- * \brief Add outflow boundary conditions for a single intersection
- */
- /*template <class IntersectionIterator>
- void evalOutflowSegment_(const IntersectionIterator &isIt,
- const BoundaryTypes &bcTypes)
- {
- if (this->element_().geometry().type().isCube() == false)
- DUNE_THROW(Dune::InvalidStateException,
- "for cell-centered models, outflow BCs only work for cubes.");
-
- // store pointer to the current FVElementGeometry
- const FVElementGeometry *oldFVGeometryPtr = this->fvElemGeomPtr_;
-
- // copy the current FVElementGeometry to a local variable
- // and set the pointer to this local variable
- FVElementGeometry fvGeometry = this->fvGeometry_();
- this->fvElemGeomPtr_ = &fvGeometry;
-
- // get the index of the boundary face
- unsigned bfIdx = isIt->indexInInside();
- unsigned oppositeIdx = bfIdx^1;
+ auto neumannFluxes = this->problem().neumann(element, fvGeometry, elemVolVars, scvf);
- // manipulate the corresponding subcontrolvolume face
- SCVFace& boundaryFace = fvGeometry.boundaryFace[bfIdx];
+ // multiply neumann fluxes with the area and the extrusion factor
+ neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
- // set the second flux approximation index for the boundary face
- for (int nIdx = 0; nIdx < fvGeometry.numNeighbors-1; nIdx++)
- {
- // check whether the two faces are opposite of each other
- if (fvGeometry.subContVolFace[nIdx].fIdx == oppositeIdx)
- {
- boundaryFace.j = nIdx+1;
- break;
+ residual[localScvIdx] += neumannFluxes;
}
- }
- boundaryFace.fapIndices[1] = boundaryFace.j;
- boundaryFace.grad[0] *= -0.5;
- boundaryFace.grad[1] *= -0.5;
-
- // temporary vector to store the outflow boundary fluxes
- PrimaryVariables values;
- Valgrind::SetUndefined(values);
- this->asImp_().computeFlux(values, bfIdx, true);
- values *= this->curVolVars_(0).extrusionFactor();
-
- // add fluxes to the residual
- Valgrind::CheckDefined(values);
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (bcTypes.isOutflow(eqIdx))
- this->residual_[0][eqIdx] += values[eqIdx];
+ else
+ DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
}
- // restore the pointer to the FVElementGeometry
- this->fvElemGeomPtr_ = oldFVGeometryPtr;
- }*/
+ }
};
-}
+} // end namespace Dumux
-#endif // DUMUX_CC_LOCAL_RESIDUAL_HH
+#endif
diff --git a/dumux/implicit/cellcentered/propertydefaults.hh b/dumux/implicit/cellcentered/propertydefaults.hh
index 99fbb8925c9fadc81dfe84ac82ad9ae3929d5e99..dbc4672ce365ceac52f4b3ad5fd5d24868bee5ed 100644
--- a/dumux/implicit/cellcentered/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/propertydefaults.hh
@@ -57,7 +57,7 @@ SET_TYPE_PROP(CCModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper)
SET_TYPE_PROP(CCModel, LocalJacobian, CCLocalJacobian<TypeTag>);
//! Assembler for the global jacobian matrix
-SET_TYPE_PROP(CCModel, JacobianAssembler, CCAssembler<TypeTag>);
+// SET_TYPE_PROP(CCModel, JacobianAssembler, CCAssembler<TypeTag>);
//! The sub control volume
SET_PROP(CCModel, SubControlVolume)
diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh
index c86f3aff1580fed9ae31ed060c3647366f0f82d4..88a7add46147a75f3ce07514b2c52650f15b4865 100644
--- a/dumux/implicit/localresidual.hh
+++ b/dumux/implicit/localresidual.hh
@@ -42,18 +42,15 @@ namespace Dumux
template<class TypeTag>
class ImplicitLocalResidual
{
- friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Model = typename GET_PROP_TYPE(TypeTag, Model);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Element = typename GridView::template Codim<0>::Entity;
+ using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ElementResidualVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
@@ -62,29 +59,11 @@ class ImplicitLocalResidual
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
public:
- // copying the local residual class is not a good idea
- ImplicitLocalResidual(const ImplicitLocalResidual &) = delete;
-
- // the default constructor
- ImplicitLocalResidual() = default;
-
- /*!
- * \brief Initialize the local residual.
- *
- * This assumes that all objects of the simulation have been fully
- * allocated but not necessarily initialized completely.
- *
- * \param problem The representation of the physical problem to be
- * solved.
- */
- void init(Problem &problem)
- { problemPtr_ = &problem; }
-
/*!
* \name User interface
* \note The following methods are usually expensive to evaluate
- * They are useful for outputting residual information.
+ * They are useful for outputting / postprocessing residual information.
*/
// \{
@@ -95,26 +74,26 @@ public:
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
*/
- void eval(const Element &element)
- {
- // make sure FVElementGeometry and volume variables are bound to the element
- auto fvGeometry = localView(problem().model().fvGridGeometry());
- fvGeometry.bind(element);
+ // ElementResidualVector eval(const Problem& problem, const Element &element)
+ // {
+ // // make sure FVElementGeometry and volume variables are bound to the element
+ // auto fvGeometry = localView(problem.model().globalFvGeometry());
+ // fvGeometry.bind(element);
- auto curElemVolVars = localView(problem().model().curGlobalVolVars());
- curElemVolVars.bind(element, fvGeometry, problem().model().curSol());
+ // auto curElemVolVars = localView(problem.model().curGlobalVolVars());
+ // curElemVolVars.bind(element, fvGeometry, problem.model().curSol());
- auto prevElemVolVars = localView(problem().model().prevGlobalVolVars());
- prevElemVolVars.bindElement(element, fvGeometry, problem().model().prevSol());
+ // auto prevElemVolVars = localView(problem.model().prevGlobalVolVars());
+ // prevElemVolVars.bindElement(element, fvGeometry, problem.model().prevSol());
- auto elemFluxVarsCache = localView(problem().model().globalFluxVarsCache());
- elemFluxVarsCache.bindElement(element, fvGeometry, curElemVolVars);
+ // auto elemFluxVarsCache = localView(problem.model().globalFluxVarsCache());
+ // elemFluxVarsCache.bindElement(element, fvGeometry, curElemVolVars);
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem(), element, fvGeometry);
+ // ElementBoundaryTypes bcTypes;
+ // bcTypes.update(problem, element, fvGeometry);
- asImp_().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, bcTypes, elemFluxVarsCache);
- }
+ // return asImp_().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, bcTypes, elemFluxVarsCache);
+ // }
/*!
* \brief Compute the storage term for the current solution.
@@ -125,52 +104,33 @@ public:
* \param element The DUNE Codim<0> entity for which the storage
* term ought to be calculated
*/
- void evalStorage(const Element &element)
- {
- // make sure FVElementGeometry and volume variables are bound to the element
- auto fvGeometry = localView(problem().model().fvGridGeometry());
- fvGeometry.bindElement(element);
-
- auto curElemVolVars = localView(problem().model().curGlobalVolVars());
- curElemVolVars.bindElement(element, fvGeometry, problem().model().curSol());
-
- auto prevElemVolVars = localView(problem().model().prevGlobalVolVars());
- prevElemVolVars.bindElement(element, fvGeometry, problem().model().prevSol());
-
- asImp_().evalStorage_(fvGeometry, prevElemVolVars, curElemVolVars);
- }
-
- // !
- // * \brief Compute the flux term for the current solution.
- // *
- // * \param element The DUNE Codim<0> entity for which the residual
- // * ought to be calculated
- // * \param curVolVars The volume averaged variables for all
- // * sub-contol volumes of the element
-
- // void evalFluxes(const Element &element)
+ // ElementResidualVector evalStorage(const Problem& problem, const Element &element)
// {
- // elemPtr_ = &element;
-
// // make sure FVElementGeometry and volume variables are bound to the element
- // problem().model().fvGeometries_().bind(element);
- // problem().model().curVolVars_().bind(element);
- // problem().model().prevVolVars_().bindElement(element);
- // problem().model().fluxVariablesCache_().bindElement(element);
+ // auto fvGeometry = localView(problem.model().globalFvGeometry());
+ // fvGeometry.bindElement(element);
- // ElementBoundaryTypes bcTypes;
- // bcTypes.update(problem(), element, fvGeometry_());
+ // auto curElemVolVars = localView(problem.model().curGlobalVolVars());
+ // curElemVolVars.bindElement(element, fvGeometry, problem.model().curSol());
- // residual_.resize(fvGeometry_().numScv);
- // residual_ = 0;
+ // ElementResidualVector storage(fvGeometry.numScv());
+ // storage.resize(fvGeometry.numScv(), 0.0);
+
+ // // calculate the amount of conservation each quantity inside
+ // // all sub control volumes
+ // for (auto&& scv : scvs(fvGeometry))
+ // {
+ // auto localScvIdx = scv.indexInElement();
+ // const auto& volVars = elemVolVars[scv];
+ // storage[localScvIdx] = asImp_().computeStorage(scv, volVars);
+ // storage[localScvIdx] *= scv.volume() * volVars.extrusionFactor();
+ // }
- // bcTypesPtr_ = &bcTypes;
- // asImp_().evalFluxes_();
+ // return storage;
// }
// \}
-
/*!
* \name Main interface
* \note Methods used by the assembler to compute derivatives and residual
@@ -180,7 +140,8 @@ public:
/*!
* \brief Compute the local residual, i.e. the deviation of the
* equations from zero.
- *
+ * \param problem The problem to solve
+
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
* \param fvGeometry The finite-volume geometry of the element
@@ -193,230 +154,183 @@ public:
* \param bcTypes The types of the boundary conditions for all
* vertices of the element
*/
- void eval(const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementBoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ ElementResidualVector eval(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const ElementBoundaryTypes &bcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
- // resize the vectors for all terms
- auto numScv = fvGeometry.numScv();
- residual_.resize(numScv);
- storageTerm_.resize(numScv);
+ // initialize the residual vector for all scvs in this element
+ ElementResidualVector residual(fvGeometry.numScv(), 0.0);
- residual_ = 0.0;
- storageTerm_ = 0.0;
+ // evaluate the volume terms (storage + source terms)
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ //! foward to the local residual specialized for the discretization methods
+ asImp_().evalStorage(residual, problem, element, fvGeometry, curElemVolVars, prevElemVolVars, scv);
+ asImp_().evalSource(residual, problem, element, fvGeometry, curElemVolVars, scv);
+ }
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ //! foward to the local residual specialized for the discretization methods
+ asImp_().evalFlux(residual, problem, element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache, scvf);
+ asImp_().evalBoundary(residual, problem, element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache, scvf);
+ }
- asImp_().evalVolumeTerms_(element, fvGeometry, prevElemVolVars, curElemVolVars, bcTypes);
- asImp_().evalFluxes_(element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache);
- asImp_().evalBoundary_(element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache);
+ return residual;
}
+ // \}
+
+
+ /*!
+ * \name Model specific interface
+ * \note The following method are the model specific implementations of the actual residual
+ */
+ // \{
+
/*!
* \brief Calculate the source term of the equation
*
* \param scv The sub-control volume over which we integrate the source term
+ * \note has to be implemented by the model specific residual class
*
*/
- PrimaryVariables computeSource(const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolume &scv)
+ ResidualVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
{
- PrimaryVariables source(0);
-
- // add contributions from volume flux sources
- source += this->problem().source(element, fvGeometry, elemVolVars, scv);
-
- // add contribution from possible point sources
- source += this->problem().scvPointSources(element, fvGeometry, elemVolVars, scv);
-
- return source;
+ DUNE_THROW(Dune::NotImplemented, "This model does not implement a storage method!");
}
/*!
- * \brief Returns the local residual for all sub-control
- * volumes of the element.
- */
- const ElementSolutionVector &residual() const
- { return residual_; }
-
- /*!
- * \brief Returns the local residual for a given sub-control
- * volume of the element.
+ * \brief Calculate the source term of the equation
+ *
+ * \param scv The sub-control volume over which we integrate the source term
+ * \note This is the default implementation for all models as sources are computed
+ * in the user interface of the problem
*
- * \param localScvIdx The local index of the sub-control volume
*/
- const PrimaryVariables &residual(const int localScvIdx) const
- { return residual_[localScvIdx]; }
+ ResidualVector computeSource(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ ResidualVector source(0.0);
- /*!
- * \brief Returns the storage term for all sub-control volumes of the
- * element.
- */
- const ElementSolutionVector &storageTerm() const
- { return storageTerm_; }
+ // add contributions from volume flux sources
+ source += problem.source(element, fvGeometry, elemVolVars, scv);
- /*!
- * \brief Returns the storage term for a given sub-control volumes
- * of the element.
- */
- const PrimaryVariables &storageTerm(const int localScvIdx) const
- { return storageTerm_[localScvIdx]; }
+ // add contribution from possible point sources
+ source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv);
- /*!
- * \brief Return the problem we are solving. Only call this after init()!
- */
- const Problem& problem() const
- { return *problemPtr_; }
+ return source;
+ }
/*!
- * \brief Return the problem we are solving. Only call this after init()!
+ * \brief Calculate the source term of the equation
+ *
+ * \param scv The sub-control volume over which we integrate the source term
+ * \note has to be implemented by the model specific residual class
+ *
*/
- Problem& problem()
- { return *problemPtr_; }
-
- // \}
-
-protected:
- Implementation &asImp_()
- { return *static_cast<Implementation*>(this); }
-
- const Implementation &asImp_() const
- { return *static_cast<const Implementation*>(this); }
-
- PrimaryVariables evalFlux_(const Element &element,
+ ResidualVector computeFlux(const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
- const ElementFluxVariablesCache& elemFluxVarsCache)
- {
- return asImp_().computeFlux_(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
- }
-
- /*!
- * \brief Set the local residual to the storage terms of all
- * sub-control volumes of the current element.
- */
- void evalStorage_(const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
- storageTerm_.resize(fvGeometry.numScv());
- storageTerm_ = 0;
-
- // calculate the amount of conservation each quantity inside
- // all sub control volumes
- for (auto&& scv : scvs(fvGeometry))
- {
- auto localScvIdx = scv.indexInElement();
- const auto& volVars = curElemVolVars[scv];
- storageTerm_[localScvIdx] = asImp_().computeStorage(scv, volVars);
- storageTerm_[localScvIdx] *= scv.volume() * volVars.extrusionFactor();
- }
+ DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!");
}
- // PrimaryVariables evalSource_(const Element& element,
- // const FVElementGeometry& fvGeometry)
- // {
- // PrimaryVariables source(0);
- // const auto& fvGeometry = fvGeometry.fvElementGeometry();
- // for (auto&& scv : scvs(fvGeometry))
- // {
- // source += this->problem().source(element, scv);
-
- // // add contribution from possible point sources
- // source += this->problem().scvPointSources(element, scv);
- // }
-
- // return source;
- // }
+ // \}
/*!
- * \brief Add the change the source term for stationary problems
- * to the local residual of all sub-control volumes of the
- * current element.
+ * \name Discretization specific interface
+ * \note The following method are the discretization specific wrapper methods
*/
- template<class P = Problem>
- typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTerms_(const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementBoundaryTypes &bcTypes)
- {
- // evaluate the volume terms (storage + source terms)
- for (auto&& scv : scvs(fvGeometry))
- {
- auto localScvIdx = scv.indexInElement();
- auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor();
-
- // subtract the source term from the local rate
- PrimaryVariables source = asImp_().computeSource(element, fvGeometry, curElemVolVars, scv);
- source *= scv.volume()*curExtrusionFactor;
-
- residual_[localScvIdx] -= source;
- }
- }
+ // \{
- /*!
- * \brief Add the change in the storage terms and the source term
- * to the local residual of all sub-control volumes of the
- * current element.
- */
- template<class P = Problem>
- typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTerms_(const Element &element,
+ void evalStorage(ElementResidualVector& residual,
+ const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
const ElementVolumeVariables& curElemVolVars,
- const ElementBoundaryTypes &bcTypes)
+ const ElementVolumeVariables& prevElemVolVars,
+ const SubControlVolume& scv) const
{
- // evaluate the volume terms (storage + source terms)
- for (auto&& scv : scvs(fvGeometry))
- {
- auto localScvIdx = scv.indexInElement();
+ const auto& curVolVars = curElemVolVars[scv];
+ const auto& prevVolVars = prevElemVolVars[scv];
- const auto& curVolVars = curElemVolVars[scv];
- const auto& prevVolVars = prevElemVolVars[scv];
+ // mass balance within the element. this is the
+ // \f$\frac{m}{\partial t}\f$ term if using implicit
+ // euler as time discretization.
+ //
+ // We might need a more explicit way for
+ // doing the time discretization...
- // mass balance within the element. this is the
- // \f$\frac{m}{\partial t}\f$ term if using implicit
- // euler as time discretization.
- //
- // We might need a more explicit way for
- // doing the time discretization...
- PrimaryVariables prevStorage = asImp_().computeStorage(scv, prevVolVars);
- PrimaryVariables curStorage = asImp_().computeStorage(scv, curVolVars);
+ //! Compute storage with the model specific storage residual
+ ResidualVector prevStorage = asImp_().computeStorage(problem, scv, prevVolVars);
+ ResidualVector storage = asImp_().computeStorage(problem, scv, curVolVars);
- prevStorage *= prevVolVars.extrusionFactor();
- curStorage *= curVolVars.extrusionFactor();
+ prevStorage *= prevVolVars.extrusionFactor();
+ storage *= curVolVars.extrusionFactor();
- storageTerm_[localScvIdx] = std::move(curStorage);
- storageTerm_[localScvIdx] -= std::move(prevStorage);
- storageTerm_[localScvIdx] *= scv.volume();
- storageTerm_[localScvIdx] /= problem().timeManager().timeStepSize();
+ storage -= prevStorage;
+ storage *= scv.volume();
+ storage /= problem.timeManager().timeStepSize();
- // add the storage term to the residual
- residual_[localScvIdx] += storageTerm_[localScvIdx];
+ residual[scv.indexInElement()] += storage;
+ }
- // subtract the source term from the local rate
- PrimaryVariables source = asImp_().computeSource(element, fvGeometry, curElemVolVars, scv);
- source *= scv.volume()*curVolVars.extrusionFactor();
+ void evalSource(ElementResidualVector& residual,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& curElemVolVars,
+ const SubControlVolume& scv) const
+ {
+ //! Compute source with the model specific storage residual
+ const auto& curVolVars = curElemVolVars[scv];
+ ResidualVector source = asImp_().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
+ source *= scv.volume()*curVolVars.extrusionFactor();
- residual_[localScvIdx] -= source;
- }
+ //! subtract source from local rate (sign convention in user interface)
+ residual[scv.indexInElement()] -= source;
}
-protected:
- ElementSolutionVector storageTerm_;
- ElementSolutionVector residual_;
+ void evalFlux(ElementResidualVector& residual,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementBoundaryTypes& elemBcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const {}
+
+ void evalBoundary(ElementResidualVector& residual,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementBoundaryTypes& elemBcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const {}
-private:
- Problem* problemPtr_;
+ // \}
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
};
-}
+} // end namespace Dumux
#endif
diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh
index 87b8cf07e1840ac4cbd7780a53651f180417eeb1..e7f872434023c84eedb94d99abb326bfd7a7c430 100644
--- a/dumux/implicit/properties.hh
+++ b/dumux/implicit/properties.hh
@@ -65,6 +65,7 @@ NEW_PROP_TAG(FVGridGeometry); //!< The type of the global finite volume geometry
NEW_PROP_TAG(EnableFVGridGeometryCache); //! specifies if geometric data is be saved (faster, but more memory consuming)
NEW_PROP_TAG(JacobianAssembler); //!< Assembles the global jacobian matrix
+NEW_PROP_TAG(LocalAssembler); //!< Assembles the global jacobian matrix
NEW_PROP_TAG(JacobianMatrix); //!< Type of the global jacobian matrix
NEW_PROP_TAG(BoundaryTypes); //!< Stores the boundary types of a single degree of freedom
NEW_PROP_TAG(ElementBoundaryTypes); //!< Stores the boundary types on an element
@@ -73,6 +74,9 @@ NEW_PROP_TAG(PrimaryVariables); //!< A vector of primary variables within a sub-
NEW_PROP_TAG(SolutionVector); //!< Vector containing all primary variable vector of the grid
NEW_PROP_TAG(ElementSolutionVector); //!< A vector of primary variables within a sub-control volume
+// TODO GridVariables
+NEW_PROP_TAG(GridVariables); //!< Doc me
+
NEW_PROP_TAG(EnableGlobalVolumeVariablesCache); //!< If disabled, the volume variables are not stored (reduces memory, but is slower)
NEW_PROP_TAG(VolumeVariables); //!< The secondary variables within a sub-control volume
NEW_PROP_TAG(ElementVolumeVariables); //!< The type for a local (element/stencil) container for the volume variables
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index a611e62a4d321898fcef76fc432322deed9c136e..a617f44e76fd5688204386c96d523c9705c228b8 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -167,7 +167,9 @@ public:
/*!
* \brief Constructor without convergence writer
*/
- NewtonController(const Communicator& comm) : comm_(comm), endIterMsgStream_(std::ostringstream::out)
+ NewtonController(const Communicator& comm)
+ : comm_(comm)
+ , endIterMsgStream_(std::ostringstream::out)
{
enablePartialReassemble_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnablePartialReassemble);
enableJacobianRecycling_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableJacobianRecycling);
@@ -196,12 +198,8 @@ public:
numSteps_ = 0;
}
- /*!
- * \brief Constructor with a given convergence writer
- */
- NewtonController(const Communicator& comm, ConvergenceWriter&& writer)
- : NewtonController(comm), convergenceWriter_(std::move(writer))
- {}
+ Communicator& communicator() const
+ { return comm_; }
/*!
* \brief Set the maximum acceptable difference of any primary variable
@@ -372,8 +370,8 @@ public:
shift_ = max(shift_, shiftAtDof);
}
- if (comm().size() > 1)
- shift_ = comm().max(shift_);
+ if (communicator().size() > 1)
+ shift_ = communicator().max(shift_);
}
/*!
@@ -406,8 +404,8 @@ public:
if (numSteps_ == 0)
{
Scalar norm2 = b.two_norm2();
- if (comm().size() > 1)
- norm2 = comm().sum(norm2);
+ if (communicator().size() > 1)
+ norm2 = communicator().sum(norm2);
using std::sqrt;
initialResidual_ = sqrt(norm2);
@@ -432,8 +430,8 @@ public:
// make sure all processes converged
int convergedRemote = converged;
- if (comm().size() > 1)
- convergedRemote = comm().min(converged);
+ if (communicator().size() > 1)
+ convergedRemote = communicator().min(converged);
if (!converged) {
DUNE_THROW(NumericalProblem,
@@ -447,8 +445,8 @@ public:
catch (Dune::MatrixBlockError e) {
// make sure all processes converged
int converged = 0;
- if (comm().size() > 1)
- converged = comm().min(converged);
+ if (communicator().size() > 1)
+ converged = communicator().min(converged);
NumericalProblem p;
std::string msg;
@@ -460,8 +458,8 @@ public:
catch (const Dune::Exception &e) {
// make sure all processes converged
int converged = 0;
- if (comm().size() > 1)
- converged = comm().min(converged);
+ if (communicator().size() > 1)
+ converged = communicator().min(converged);
NumericalProblem p;
p.message(e.what());
@@ -495,7 +493,11 @@ public:
if (enableShiftCriterion_)
newtonUpdateShift(uLastIter, deltaU);
- writeConvergence_(uLastIter, deltaU);
+ if (writeConvergence_)
+ {
+ convergenceWriter_->advanceIteration();
+ convergenceWriter_->write(uLastIter, deltaU);
+ }
if (useLineSearch_)
{
@@ -544,7 +546,7 @@ public:
// classes directly as well?
// -> The assembler has all the data w.r.t. the problem etc...
// -> also, how do we call this?
- assembler.updateVariables()
+ assembler.updateVariables();
++numSteps_;
@@ -628,7 +630,7 @@ public:
* \brief Returns true if the Newton method ought to be chatty.
*/
bool verbose() const
- { return verbose_ && comm().rank() == 0; }
+ { return verbose_ && communicator().rank() == 0; }
protected:
@@ -641,16 +643,6 @@ protected:
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
- void writeConvergence_(const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- if (writeConvergence_)
- {
- convergenceWriter_->advanceIteration();
- convergenceWriter_->write(uLastIter, deltaU);
- }
- }
-
void lineSearchUpdate_(const JacobianAssembler& assembler,
SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
@@ -694,6 +686,7 @@ protected:
* \param priVars1 The first vector of primary variables
* \param priVars2 The second vector of primary variables
*/
+ template<class PrimaryVariables>
Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1,
const PrimaryVariables &priVars2)
{
diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh
index 0b142216c5f5613043ca94ad76eca74190f204c1..41e3fe14aff766121a8791b3273af33d316d8554 100644
--- a/dumux/nonlinear/newtonmethod.hh
+++ b/dumux/nonlinear/newtonmethod.hh
@@ -46,6 +46,9 @@ NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(NewtonController);
NEW_PROP_TAG(SolutionVector);
NEW_PROP_TAG(JacobianAssembler);
+NEW_PROP_TAG(JacobianMatrix);
+NEW_PROP_TAG(NewtonConvergenceWriter);
+NEW_PROP_TAG(WriteConvergence);
}
/*!
@@ -62,6 +65,7 @@ class NewtonMethod
using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
using ConvergenceWriter = typename GET_PROP_TYPE(TypeTag, NewtonConvergenceWriter);
using JacobianAssembler = typename GET_PROP_TYPE(TypeTag, JacobianAssembler);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
using LinearSolver = typename GET_PROP_TYPE(TypeTag, LinearSolver);
public:
@@ -69,18 +73,8 @@ public:
std::shared_ptr<LinearSolver> linearSolver)
: jacobianAssembler_(jacobianAssembler)
, linearSolver_(linearSolver)
- , matrix_ = std::make_shared<JacobianMatrix>()
- , residual_ = std::make_shared<SolutionVector>()
{
- // construct the newton controller with or without convergence writer
- if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence))
- {
- newtonController_ = std::make_unique<NewtonController>(assembler().gridView().comm(),
- NewtonConvergenceWriter(gridView, assembler().numDofs()));
- }
- else
- newtonController_ = std::make_unique<NewtonController>(assembler().gridView().comm());
-
+ newtonController_ = std::make_unique<NewtonController>(assembler().gridView().comm());
// set the linear system (matrix & residual) in the assembler
assembler().setLinearSystem(matrix(), residual());
}
@@ -95,9 +89,9 @@ public:
return execute_();
}
catch (const NumericalProblem &e) {
- if (ctl.verbose())
+ if (controller().verbose())
std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
- ctl.newtonFail();
+ controller().newtonFail();
return false;
}
}
@@ -144,7 +138,7 @@ private:
if (controller().newtonNumSteps() > 0)
uLastIter = uCurrentIter;
- if (ctl.verbose()) {
+ if (controller().verbose()) {
std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r";
std::cout.flush();
}
@@ -167,7 +161,7 @@ private:
// http://en.wikipedia.org/wiki/ANSI_escape_code
const char clearRemainingLine[] = { 0x1b, '[', 'K', 0 };
- if (ctl.verbose()) {
+ if (controller().verbose()) {
std::cout << "\rSolve: M deltax^k = r";
std::cout << clearRemainingLine;
std::cout.flush();
@@ -180,15 +174,15 @@ private:
deltaU = 0;
// ask the controller to solve the linearized system
controller().solveLinearSystem(linearSolver(),
- jacobianAsm.matrix(),
+ matrix(),
deltaU,
- jacobianAsm.residual());
+ residual());
solveTimer.stop();
///////////////
// update
///////////////
- if (ctl.verbose()) {
+ if (controller().verbose()) {
std::cout << "\rUpdate: x^(k+1) = x^k - deltax^k";
std::cout << clearRemainingLine;
std::cout.flush();
@@ -197,17 +191,17 @@ private:
updateTimer.start();
// update the current solution (i.e. uOld) with the delta
// (i.e. u). The result is stored in u
- ctl.newtonUpdate(assembler(), uCurrentIter, uLastIter, deltaU);
+ controller().newtonUpdate(assembler(), uCurrentIter, uLastIter, deltaU);
updateTimer.stop();
// tell the controller that we're done with this iteration
- ctl.newtonEndStep(assembler(), uCurrentIter, uLastIter);
+ controller().newtonEndStep(assembler(), uCurrentIter, uLastIter);
}
// tell the controller that we're done
- ctl.newtonEnd();
+ controller().newtonEnd();
- if (ctl.verbose()) {
+ if (controller().verbose()) {
Scalar elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
std::cout << "Assemble/solve/update time: "
<< assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
@@ -216,12 +210,12 @@ private:
<< "\n";
}
- if (!ctl.newtonConverged()) {
- ctl.newtonFail();
+ if (!controller().newtonConverged()) {
+ controller().newtonFail();
return false;
}
- ctl.newtonSucceed();
+ controller().newtonSucceed();
return true;
}
diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh
index c8f4361b3d4dba4d9de8804dc9fdc4676a6cd444..2365cefbddd2898d87695582d62c678f34eebba6 100644
--- a/dumux/porousmediumflow/immiscible/localresidual.hh
+++ b/dumux/porousmediumflow/immiscible/localresidual.hh
@@ -36,10 +36,9 @@ namespace Dumux
template<class TypeTag>
class ImmiscibleLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
-
using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
@@ -69,7 +68,8 @@ public:
* \note The volVars can be different to allow computing
* the implicit euler time derivative here
*/
- PrimaryVariables computeStorage(const SubControlVolume& scv,
+ PrimaryVariables computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
// partial time derivative of the phase mass
@@ -96,14 +96,15 @@ public:
* \brief Evaluate the mass flux over a face of a sub control volume
* \param scvf The sub control volume face to compute the flux on
*/
- PrimaryVariables computeFlux(const Element& element,
+ PrimaryVariables computeFlux(const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
- fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
PrimaryVariables flux;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@@ -124,15 +125,8 @@ public:
return flux;
}
-
-private:
- Implementation *asImp_()
- { return static_cast<Implementation *> (this); }
-
- const Implementation *asImp_() const
- { return static_cast<const Implementation *> (this); }
};
-}
+} // end namespace Dumux
#endif
diff --git a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
index 349b045b476f4feb39cfa0b428306367b0dc1b3c..ca43301d6f6050ba6ff2be7a7ad35618024d1c25 100644
--- a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
+++ b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
@@ -36,27 +36,6 @@
namespace Dumux
{
-template<class TypeTag>
-class MockProblem
-{
- using ElementMapper = typename GET_PROP_TYPE(TypeTag, DofMapper);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-public:
- MockProblem(const GridView& gridView) : mapper_(gridView) {}
-
- const ElementMapper& elementMapper() const
- { return mapper_; }
-
- template<class Element, class Intersection>
- bool isInteriorBoundary(const Element& e, const Intersection& i) const
- { return false; }
-
- std::vector<unsigned int> getAdditionalDofDependencies(unsigned int index) const
- { return std::vector<unsigned int>(); }
-private:
- ElementMapper mapper_;
-};
-
namespace Properties
{
NEW_TYPE_TAG(TestFVGeometry, INHERITS_FROM(CCTpfaModel));
diff --git a/test/porousmediumflow/1p/implicit/CMakeLists.txt b/test/porousmediumflow/1p/implicit/CMakeLists.txt
index 0997c152fe295579c19a28b42039a1834026ccb6..9713f6fd8116c8652db833919a15e9ea2630f87a 100644
--- a/test/porousmediumflow/1p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory(pointsources)
+add_subdirectory(incompressible)
add_input_file_links()
add_gstat_file_links()
diff --git a/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt b/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4f3ee5b2806725005411a3ce4bf0b6fceea476fb
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt
@@ -0,0 +1,5 @@
+dune_symlink_to_source_files(FILES "test_1pincompressible.input")
+
+# incompressible
+add_dumux_test(test_1pincompressible test_1pincompressible test_1pincompressible.cc
+ ${CMAKE_CURRENT_BINARY_DIR}/test_1pincompressible)
diff --git a/test/porousmediumflow/1p/implicit/incompressible/properties.hh b/test/porousmediumflow/1p/implicit/incompressible/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c9b2584d2b2a49e4986ac78e4b5be00cd1df8262
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/incompressible/properties.hh
@@ -0,0 +1,75 @@
+// -*- 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 The properties for the incompressible test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROPERTIES_HH
+#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROPERTIES_HH
+
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh>
+#include <dumux/porousmediumflow/1p/implicit/propertydefaults.hh>
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class MockProblem
+{
+ using ElementMapper = typename GET_PROP_TYPE(TypeTag, DofMapper);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+public:
+ MockProblem(const GridView& gridView) : mapper_(gridView) {}
+
+ const ElementMapper& elementMapper() const
+ { return mapper_; }
+
+ template<class Element, class Intersection>
+ bool isInteriorBoundary(const Element& e, const Intersection& i) const
+ { return false; }
+
+ std::vector<unsigned int> getAdditionalDofDependencies(unsigned int index) const
+ { return std::vector<unsigned int>(); }
+private:
+ ElementMapper mapper_;
+};
+
+namespace Properties
+{
+
+NEW_PROP_TAG(EnableFVGridGeometryCache);
+NEW_PROP_TAG(FVGridGeometry);
+
+NEW_TYPE_TAG(IncompressibleTestProblem, INHERITS_FROM(OneP, CCTpfaModel));
+
+// Set the grid type
+SET_TYPE_PROP(IncompressibleTestProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the finite volume grid geometry
+SET_TYPE_PROP(IncompressibleTestProblem, FVGridGeometry, CCTpfaFVGridGeometry<TypeTag, false>);
+
+// Set the problem type
+SET_TYPE_PROP(IncompressibleTestProblem, Problem, Dumux::MockProblem<TypeTag>);
+
+} // end namespace Properties
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/1p/implicit/incompressible/test_1pincompressible.cc b/test/porousmediumflow/1p/implicit/incompressible/test_1pincompressible.cc
new file mode 100644
index 0000000000000000000000000000000000000000..daaf557e3cff1c3fc2fcc255f91baa695ac63d41
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/incompressible/test_1pincompressible.cc
@@ -0,0 +1,126 @@
+// -*- 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 test for the one-phase CC model
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+
+#include <dumux/common/propertysystem.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+#include <dumux/common/parameterparser.hh>
+
+#include "properties.hh"
+//#include "problem.hh"
+
+template<class TypeTag>
+class MockProblem
+{
+ using ElementMapper = typename GET_PROP_TYPE(TypeTag, DofMapper);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+public:
+ MockProblem(const GridView& gridView) : mapper_(gridView) {}
+
+ const ElementMapper& elementMapper() const
+ { return mapper_; }
+
+ template<class Element, class Intersection>
+ bool isInteriorBoundary(const Element& e, const Intersection& i) const
+ { return false; }
+
+ std::vector<unsigned int> getAdditionalDofDependencies(unsigned int index) const
+ { return std::vector<unsigned int>(); }
+private:
+ ElementMapper mapper_;
+};
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ using TypeTag = TTAG(IncompressibleTestProblem);
+
+ // some aliases for better readability
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ // using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
+ using ParameterTree = typename GET_PROP(TypeTag, ParameterTree);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments and input file
+ ////////////////////////////////////////////////////////////
+
+ // parse command line arguments
+ ParameterParser::parseCommandLineArguments(argc, argv, ParameterTree::tree());
+
+ // parse the input file into the parameter tree
+ // check first if the user provided an input file through the command line, if not use the default
+ const auto parameterFileName = ParameterTree::tree().hasKey("ParameterFile") ? GET_RUNTIME_PARAM(TypeTag, std::string, ParameterFile) : "";
+ ParameterParser::parseInputFile(argc, argv, ParameterTree::tree(), parameterFileName);
+
+ //////////////////////////////////////////////////////////////////////
+ // try to create a grid (from the given grid file or the input file)
+ /////////////////////////////////////////////////////////////////////
+
+ try { GridCreator::makeGrid(); }
+ catch (...) {
+ std::cout << "\n\t -> Creation of the grid failed! <- \n\n";
+ throw;
+ }
+ GridCreator::loadBalance();
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // do some testing
+ for (const auto& element : elements(leafGridView))
+ {
+ auto fvGeometry = localView(*fvGridGeometry);
+ fvGeometry.bind(element);
+ for (const auto& scv : scvs(fvGeometry))
+ std::cout << scv.dofPosition() << " ";
+ std::cout << std::endl;
+ }
+
+ return 0;
+
+} // end main
diff --git a/test/porousmediumflow/1p/implicit/incompressible/test_1pincompressible.input b/test/porousmediumflow/1p/implicit/incompressible/test_1pincompressible.input
new file mode 100644
index 0000000000000000000000000000000000000000..bd21b5f555c67b6d7cf4624e7e8af6f29dc663aa
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/incompressible/test_1pincompressible.input
@@ -0,0 +1,19 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 1 # [s]
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 1 1
+Cells = 10 10
+
+[Problem]
+Name = 1ptestcc # name passed to the output routines
+
+[SpatialParams]
+LensLowerLeft = 0.2 0.2
+LensUpperRight = 0.8 0.8
+
+Permeability = 1e-10 # [m^2]
+PermeabilityLens = 1e-12 # [m^2]
+