diff --git a/dumux/assembly/staggeredfvassembler.hh b/dumux/assembly/staggeredfvassembler.hh
index 28020a6beddcf309029e1ac025a05fb4902538ca..d3dc3203c8e16ea360b8d7a99a2a19090d1853cf 100644
--- a/dumux/assembly/staggeredfvassembler.hh
+++ b/dumux/assembly/staggeredfvassembler.hh
@@ -109,9 +109,7 @@ public:
if(!jacobian_)
{
- jacobian_ = std::make_shared<JacobianMatrix>();
- jacobian_->setBuildMode(JacobianMatrix::random);
- setJacobianPattern();
+ setLinearSystem();
}
if(!residual_)
@@ -129,12 +127,17 @@ public:
{
// let the local assembler add the element contributions
for (const auto element : elements(gridView()))
- LocalAssembler::assemble(*this, *jacobian_, *residual_, element, curSol);
+ LocalAssembler::assembleJacobianAndResidual(*this, *jacobian_, *residual_, element, curSol);
// if we get here, everything worked well
succeeded = true;
if (gridView().comm().size() > 1)
succeeded = gridView().comm().min(succeeded);
+
+ // printmatrix(std::cout, (*jacobian_)[cellCenterIdx][cellCenterIdx], "A11", "");
+ // printmatrix(std::cout, (*jacobian_)[cellCenterIdx][faceIdx], "A12", "");
+ // printmatrix(std::cout, (*jacobian_)[faceIdx][cellCenterIdx], "A21", "");
+ // printmatrix(std::cout, (*jacobian_)[faceIdx][faceIdx], "A22", "");
}
// throw exception if a problem ocurred
catch (NumericalProblem &e)
@@ -220,7 +223,9 @@ public:
//! computes the residual norm
Scalar residualNorm(const SolutionVector& curSol) const
{
- ResidualType residual(numDofs());
+ ResidualType residual;
+ residual[cellCenterIdx].resize(numCellCenterDofs());
+ residual[faceIdx].resize(numFaceDofs());
assembleResidual(residual, curSol);
// calculate the square norm of the residual
@@ -240,11 +245,23 @@ public:
void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
std::shared_ptr<SolutionVector> r)
{
- DUNE_THROW(Dune::NotImplemented, "Setting linear system with existing matrix not implemented yet.");
- // jacobian_ = A;
- // residual_ = r;
+ jacobian_ = A;
+ residual_ = r;
//
- // // check and/or set the BCRS matrix's build mode
+ // check and/or set the BCRS matrix's build mode
+ // convenience references
+ CCToCCMatrixBlock& A11 = (*jacobian_)[cellCenterIdx][cellCenterIdx];
+ CCToFaceMatrixBlock& A12 = (*jacobian_)[cellCenterIdx][faceIdx];
+ FaceToFaceMatrixBlock& A21 = (*jacobian_)[faceIdx][cellCenterIdx];
+ FaceToCCMatrixBlock& A22 = (*jacobian_)[faceIdx][faceIdx];
+
+ A11.setBuildMode(CCToCCMatrixBlock::random);
+ A12.setBuildMode(CCToFaceMatrixBlock::random);
+ A21.setBuildMode(FaceToFaceMatrixBlock::random);
+ A22.setBuildMode(FaceToCCMatrixBlock::random);
+
+ setJacobianPattern();
+ setResidualSize();
// if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown)
// jacobian_->setBuildMode(JacobianMatrix::random);
// else if (jacobian_->buildMode() != JacobianMatrix::BuildMode::random)
diff --git a/dumux/assembly/staggeredlocalassembler.hh b/dumux/assembly/staggeredlocalassembler.hh
index 272d4f2d7669d6e7c6f139e4a62ca37d73538e40..42aef23ca2606a796f06479fa3e54603b201318a 100644
--- a/dumux/assembly/staggeredlocalassembler.hh
+++ b/dumux/assembly/staggeredlocalassembler.hh
@@ -30,6 +30,8 @@
#include <dumux/implicit/properties.hh>
#include <dumux/assembly/diffmethod.hh>
+#include <dumux/implicit/staggered/primaryvariables.hh>
+
namespace Dumux {
/*!
@@ -51,6 +53,9 @@ class StaggeredLocalAssembler<TypeTag,
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
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);
@@ -60,8 +65,22 @@ class StaggeredLocalAssembler<TypeTag,
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using NumCellCenterEqVector = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using NumFaceEqVector = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using PriVarIndices = typename Dumux::PriVarIndices<TypeTag>;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ // using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ // typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ // typename DofTypeIndices::FaceIdx faceIdx;
+
static constexpr bool enableGlobalFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache);
public:
@@ -71,71 +90,13 @@ public:
* to the global matrix. The element residual is written into the right hand side.
*/
template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
+ static void assembleJacobianAndResidual(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
const Element& element, const SolutionVector& curSol)
{
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
- }
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
- }
-
- /*!
- * \brief Computes the epsilon used for numeric differentiation
- * for a given value of a primary variable.
- *
- * \param priVar The value of the primary variable
- */
- 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
- // floating point values, the resolution is about 10^-16 and
- // the base epsilon is thus approximately 10^-8.
- /*
- static const Scalar baseEps
- = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
- */
- static const Scalar baseEps = 1e-10;
- assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
- // the epsilon value used for the numeric differentiation is
- // now scaled by the absolute value of the primary variable...
- return baseEps*(std::abs(priVar) + 1.0);
- }
-
-private:
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
- {
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
// get some references for convenience
const auto& problem = assembler.problem();
@@ -152,473 +113,73 @@ private:
auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
- const bool isStationary = localResidual.isStationary();
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- if (!isStationary)
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
-
- // for compatibility with box models
- ElementBoundaryTypes elemBcTypes;
-
- // the actual element's current residual
- NumEqVector residual(0.0);
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
-
- return residual;
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- const auto& connectivityMap = fvGridGeometry.connectivityMap();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+ auto& curGlobalFaceVars = gridVariables.curGridFaceVars();
+ auto& prevGlobalFaceVars = gridVariables.prevGridFaceVars();
const bool isStationary = localResidual.isStationary();
auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
if (!isStationary)
prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
+ // for compatibility with box models
ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- // the actual element's current residual
- NumEqVector residual(0.0);
- if (!isGhost)
- {
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
-
-
- // TODO Do we really need this??????????
- // this->model_().updatePVWeights(fvGeometry);
+ const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
+ res[cellCenterIdx][cellCenterGlobalI] = localResidual.evalCellCenter(problem,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ curGlobalFaceVars,
+ prevGlobalFaceVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
- //////////////////////////////////////////////////////////////////////////////////////////////////
- // //
- // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
- // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
- // actual element. In the actual element we evaluate the derivative of the entire residual. //
- // //
- //////////////////////////////////////////////////////////////////////////////////////////////////
- static const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
- static const int numericDifferenceMethod = getParamFromGroup<int>(group, "Implicit.NumericDifferenceMethod");
+ // std::cout << "at elem: " << cellCenterGlobalI << std::endl;
+ // std::cout << "orig res: " << res[cellCenterIdx][cellCenterGlobalI] << std::endl;
- // get stencil informations
- const auto numNeighbors = connectivityMap[globalI].size();
+ // treat the local residua of the face dofs:
+ // create a cache to reuse some results for the calculation of the derivatives
- // container to store the neighboring elements
- std::vector<Element> neighborElements;
- neighborElements.reserve(numNeighbors);
+ FaceSolutionVector faceResidualCache;
+ faceResidualCache.resize(fvGeometry.numScvf());
+ faceResidualCache = 0.0;
- // get the elements in which we need to evaluate the fluxes
- // and calculate these in the undeflected state
- Dune::BlockVector<NumEqVector> origFlux(numNeighbors);
- origFlux = 0.0;
- unsigned int j = 0;
- for (const auto& dataJ : connectivityMap[globalI])
+ for(auto&& scvf : scvfs(fvGeometry))
{
- neighborElements.emplace_back(fvGridGeometry.element(dataJ.globalJ));
- for (const auto scvfIdx : dataJ.scvfsJ)
- {
- origFlux[j] += localResidual.evalFlux(problem,
- neighborElements.back(),
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- fvGeometry.scvf(scvfIdx));
- }
- // increment neighbor counter
- ++j;
+ // res[faceIdx][scvf.dofIndex()] = 0.0;
+ faceResidualCache[scvf.localFaceIdx()] = localResidual.evalFace(problem,
+ element,
+ fvGeometry,
+ scvf,
+ prevElemVolVars,
+ curElemVolVars,
+ curGlobalFaceVars,
+ prevGlobalFaceVars,
+ elemBcTypes,
+ elemFluxVarsCache)[0];
+
+ res[faceIdx][scvf.dofIndex()] += faceResidualCache[scvf.localFaceIdx()] ;
}
- // reference to the element's scv (needed later) and corresponding vol vars
- const auto& scv = fvGeometry.scv(globalI);
- auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
-
- // save a copy of the original privars and vol vars in order
- // to restore the original solution after deflection
- const auto origPriVars = curSol[globalI];
- const auto origVolVars = curVolVars;
-
- // element solution container to be deflected
- ElementSolutionVector elemSol({origPriVars});
-
- // derivatives in the neighbors with repect to the current elements
- Dune::BlockVector<NumEqVector> neighborDeriv(numNeighbors);
- for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
- {
- // reset derivatives of element dof with respect to itself
- // as well as neighbor derivatives
- NumEqVector partialDeriv(0.0);
- neighborDeriv = 0.0;
-
- if (isGhost)
- partialDeriv[pvIdx] = 1.0;
-
- Scalar eps = numericEpsilon(curVolVars.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
- elemSol[0][pvIdx] += eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
- if (enableGlobalFluxVarsCache)
- gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
- else
- elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
-
- // calculate the residual with the deflected primary variables
- if (!isGhost)
- {
- if (isStationary)
- {
- partialDeriv = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- partialDeriv = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
-
- // calculate the fluxes in the neighbors with the deflected primary variables
- for (std::size_t k = 0; k < numNeighbors; ++k)
- for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
- {
- neighborDeriv[k] += localResidual.evalFlux(problem,
- neighborElements[k],
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- fvGeometry.scvf(scvfIdx));
- }
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv = residual;
- neighborDeriv = origFlux;
- }
-
- if (numericDifferenceMethod <= 0)
- {
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- elemSol[0][pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
- if (enableGlobalFluxVarsCache)
- gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
- else
- elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
-
- // calculate the residual with the deflected primary variables and subtract it
- if (!isGhost)
- {
- if (isStationary)
- {
- partialDeriv -= localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- partialDeriv -= localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
-
- // calculate the fluxes into element with the deflected primary variables
- for (std::size_t k = 0; k < numNeighbors; ++k)
- for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
- {
- neighborDeriv[k] += localResidual.evalFlux(problem,
- neighborElements[k],
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- fvGeometry.scvf(scvfIdx));
- }
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv -= residual;
- neighborDeriv -= origFlux;
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- if (!isGhost)
- partialDeriv /= delta;
- neighborDeriv /= delta;
-
- // restore the original state of the scv's volume variables
- curVolVars = origVolVars;
-
- // restore the current element solution
- elemSol[0][pvIdx] = origPriVars[pvIdx];
-
- // add the current partial derivatives to the global jacobian matrix
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- // the diagonal entries
- A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
-
- // off-diagonal entries
- j = 0;
- for (const auto& dataJ : connectivityMap[globalI])
- A[dataJ.globalJ][globalI][eqIdx][pvIdx] += neighborDeriv[j++][eqIdx];
- }
- }
-
- //////////////////////////////////////////////////////////////////////////////////////////////
- // //
- // Calculate derivatives of the dofs in the element with respect to user-defined additional //
- // dof dependencies. We do so by evaluating the change in the source term of the current //
- // element with respect to the primary variables at the given additional dofs. //
- // //
- //////////////////////////////////////////////////////////////////////////////////////////////
-
- // 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 = fvGridGeometry.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 = 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;
- }
-private:
- template<class T = TypeTag>
- static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return elemVolVars[scv]; }
+ // calculate derivatives of all dofs in stencil with respect to the dofs in the element
+ evalPartialDerivatives_(assembler,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ prevGlobalFaceVars,
+ curGlobalFaceVars,
+ elemFluxVarsCache,
+ elemBcTypes,
+ jac,
+ res[cellCenterIdx][cellCenterGlobalI],
+ faceResidualCache);
- template<class T = TypeTag>
- static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return gridVolVars.volVars(scv); }
-};
-//! Explicit assembler with numeric differentiation
-template<class TypeTag>
-class StaggeredLocalAssembler<TypeTag,
- DiffMethod::numeric,
- /*implicit=*/false>
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- 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);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-
-public:
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix. The element residual is written into the right hand side.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
}
/*!
@@ -629,6 +190,7 @@ public:
static void assemble(Assembler& assembler, JacobianMatrix& jac,
const Element& element, const SolutionVector& curSol)
{
+ std::cout << "calling wrong \n";
assemble_(assembler, jac, element, curSol);
}
@@ -639,8 +201,39 @@ public:
static void assemble(Assembler& assembler, SolutionVector& res,
const Element& element, const SolutionVector& curSol)
{
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
+ std::cout << "calling wrong \n";
+ // using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ // typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ // typename DofTypeIndices::FaceIdx faceIdx;
+
+ // const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
+ // res[cellCenterIdx][cellCenterGlobalI] = localResidual.evalCellCenter(problem,
+ // element,
+ // fvGeometry,
+ // prevElemVolVars,
+ // curElemVolVars,
+ // curGlobalFaceVars,
+ // prevGlobalFaceVars,
+ // elemBcTypes,
+ // elemFluxVarsCache)[0];
+
+
+
+ // treat the local residua of the face dofs:
+ // create a cache to reuse some results for the calculation of the derivatives
+ // FaceSolutionVector faceResidualCache;
+ // faceResidualCache.resize(assembler.fvGridGeometry().numScvf());
+ // faceResidualCache = 0.0;
+ //
+ // auto fvGeometry = localView(assembler.fvGridGeometry());
+ // fvGeometry.bind(element);
+ // auto faceResiduals = assembleFace_(assembler, element, curSol);
+ //
+ // for(auto&& scvf : scvfs(fvGeometry))
+ // {
+ // res[faceIdx][scvf.dofIndex()] += faceResiduals[scvf.localFaceIdx()];
+ // faceResidualCache[scvf.localFaceIdx()] = faceResiduals[scvf.localFaceIdx()];
+ // }
}
/*!
@@ -666,730 +259,362 @@ public:
return baseEps*(std::abs(priVar) + 1.0);
}
-private:
-
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
- {
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
-
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // compatibility with box method
- ElementBoundaryTypes elemBcTypes;
-
- // the actual element's previous time step residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
-
- return residual;
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- *
- * \return The element residual at the current solution.
- */
+protected:
template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
+ static void evalPartialDerivatives_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevGlobalFaceVars,
+ GlobalFaceVars& curGlobalFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const NumCellCenterEqVector& ccResidual,
+ const FaceSolutionVector& faceResidualCache)
+{
+ // compute the derivatives of the cell center residuals with respect to cell center dofs
+ dCCdCC_(assembler, element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
+ // compute the derivatives of the cell center residuals with respect to face dofs
+ dCCdFace_(assembler, element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
+ // compute the derivatives of the face residuals with respect to cell center dofs
+ dFacedCC_(assembler, element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
+ // compute the derivatives of the face residuals with respect to face dofs
+ dFacedFace_(assembler, element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
+}
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
+/*!
+* \brief Computes the derivatives of the cell center residuals with respect to cell center dofs
+*/
+template<class Assembler>
+static void dCCdCC_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevGlobalFaceVars,
+ const GlobalFaceVars& curGlobalFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const NumCellCenterEqVector& ccResidual)
+{
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
- // the actual element's previous time step residual
- NumEqVector residual(0.0), storageResidual(0.0);
- if (!isGhost)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
- }
+ // build derivatives with for cell center dofs w.r.t. cell center dofs
+ const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
- //////////////////////////////////////////////////////////////////////////////////////////////////
- // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
- // neighboring elements all derivatives are zero. For the assembled element only the storage //
- // derivatives are non-zero. //
- //////////////////////////////////////////////////////////////////////////////////////////////////
+ const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
- static const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
- static const int numericDifferenceMethod = getParamFromGroup<int>(group, "Implicit.NumericDifferenceMethod");
+ for(const auto& globalJ : connectivityMap(cellCenterIdx, cellCenterIdx, cellCenterGlobalI))
+ {
+ // get the volVars of the element with respect to which we are going to build the derivative
+ auto&& scvJ = fvGeometry.scv(globalJ);
+ const auto elementJ = fvGeometry.fvGridGeometry().element(globalJ);
+ auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
+ // auto& curVolVars = getCurVolVars(curElemVolVars, scvJ);
+ VolumeVariables origVolVars(curVolVars);
- // reference to the element's scv (needed later) and corresponding vol vars
- const auto& scv = fvGeometry.scv(globalI);
- auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+ for(auto pvIdx : PriVarIndices(cellCenterIdx))
+ {
+ PrimaryVariables priVars(CellCenterPrimaryVariables(localResidual.prevSol()[cellCenterIdx][globalJ]),
+ FacePrimaryVariables(0.0));
- // save a copy of the original privars and vol vars in order
- // to restore the original solution after deflection
- const auto origPriVars = curSol[globalI];
- const auto origVolVars = curVolVars;
+ const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, cellCenterIdx);
+ priVars[pvIdx] += eps;
+ ElementSolutionVector elemSol{std::move(priVars)};
+ curVolVars.update(elemSol, problem, elementJ, scvJ);
- // element solution container to be deflected
- ElementSolutionVector elemSol({origPriVars});
+ auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry, prevElemVolVars, curElemVolVars,
+ prevGlobalFaceVars, curGlobalFaceVars,
+ elemBcTypes, elemFluxVarsCache);
- // derivatives in the neighbors with repect to the current elements
- for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
- {
- // reset derivatives of element dof with respect to itself
- // as well as neighbor derivatives
- NumEqVector partialDeriv(0.0);
-
- if (isGhost)
- partialDeriv[pvIdx] = 1.0;
-
- Scalar eps = numericEpsilon(curVolVars.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
- elemSol[0][pvIdx] += eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
-
- // calculate the residual with the deflected primary variables
- if (!isGhost)
- {
- partialDeriv = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv = storageResidual;
- }
-
- if (numericDifferenceMethod <= 0)
- {
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- elemSol[0][pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
-
- // calculate the residual with the deflected primary variables and subtract it
- if (!isGhost)
- {
- partialDeriv -= localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv -= storageResidual;
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- if (!isGhost)
- partialDeriv /= delta;
-
- // restore the original state of the scv's volume variables
- curVolVars = origVolVars;
-
- // restore the current element solution
- elemSol[0][pvIdx] = origPriVars[pvIdx];
-
- // add the current partial derivatives to the global jacobian matrix
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- // the diagonal entries
- A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
- }
- }
+ auto partialDeriv = (deflectedResidual - ccResidual);
+ partialDeriv /= eps;
- // return the original residual
- return residual;
- }
-private:
- template<class T = TypeTag>
- static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return elemVolVars[scv]; }
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[cellCenterIdx][cellCenterIdx], cellCenterGlobalI, globalJ, pvIdx, partialDeriv);
- template<class T = TypeTag>
- static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return gridVolVars.volVars(scv); }
-};
+ // restore the original volVars
+ curVolVars = origVolVars;
+ }
+ }
+}
-//! implicit assembler using analytic differentiation
-template<class TypeTag>
-class StaggeredLocalAssembler<TypeTag,
- DiffMethod::analytic,
- /*implicit=*/true>
+/*!
+* \brief Computes the derivatives of the cell center residuals with respect to face dofs
+*/
+template<class Assembler>
+static void dCCdFace_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevGlobalFaceVars,
+ GlobalFaceVars& curGlobalFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const NumCellCenterEqVector& ccResidual)
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using IndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ // build derivatives with for cell center dofs w.r.t. face dofs
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ // auto& gridVariables = assembler.gridVariables();
-public:
+ // build derivatives with for cell center dofs w.r.t. cell center dofs
+ const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix. The element residual is written into the right hand side.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
- }
+ const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
- }
+ for(const auto& globalJ : connectivityMap(cellCenterIdx, faceIdx, cellCenterGlobalI))
+ {
+ // get the faceVars of the face with respect to which we are going to build the derivative
+ auto origFaceVars = curGlobalFaceVars.faceVars(globalJ);
+ auto& curFaceVars = curGlobalFaceVars.faceVars(globalJ);
-private:
+ for(auto pvIdx : PriVarIndices(faceIdx))
+ {
+ PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(localResidual.prevSol()[faceIdx][globalJ]));
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
- {
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
-
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
+ // std::cout << "orig velo : " << curGlobalFaceVars.faceVars(globalJ).velocity() << std::endl;
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
+ const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, faceIdx);
+ priVars[pvIdx] += eps;
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+ // std::cout << "deflecting " << globalJ << std::endl;
+ // std::cout << "eps. " << eps << std::endl;
+ curFaceVars.update(priVars[faceIdx]);
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ // std::cout << "deflected velo: " << curGlobalFaceVars.faceVars(globalJ).velocity() << std::endl;
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
+ auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry,
+ prevElemVolVars, curElemVolVars,
+ prevGlobalFaceVars, curGlobalFaceVars,
+ elemBcTypes, elemFluxVarsCache);
- NumEqVector residual(0.0);
- if (!localResidual.isStationary())
- {
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+ auto partialDeriv = (deflectedResidual - ccResidual);
+ partialDeriv /= eps;
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
+ // std::cout << "ccResidual dccdface " << ccResidual << std::endl;
+ // std::cout << "deflectedResidual dccdface " << deflectedResidual << std::endl;
+ // std::cout << "partialDeriv dccdface " << partialDeriv << std::endl;
- return residual;
- }
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], cellCenterGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
+ // restore the original faceVars
+ curFaceVars = origFaceVars;
+ }
+ }
+}
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
+/*!
+* \brief Computes the derivatives of the face residuals with respect to cell center dofs
+*/
+template<class Assembler>
+static void dFacedCC_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevGlobalFaceVars,
+ const GlobalFaceVars& curGlobalFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const FaceSolutionVector& cachedResidual)
+{
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ auto& gridVariables = assembler.gridVariables();
+ const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
+
+ // set the actual dof index
+ const auto faceGlobalI = scvf.dofIndex();
+
+ // build derivatives with for face dofs w.r.t. cell center dofs
+ for(const auto& globalJ : connectivityMap(faceIdx, cellCenterIdx, scvf.index()))
+ {
+ // get the volVars of the element with respect to which we are going to build the derivative
+ auto&& scvJ = fvGeometry.scv(globalJ);
+ const auto elementJ = fvGeometry.fvGridGeometry().element(globalJ);
+ auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
+ VolumeVariables origVolVars(curVolVars);
+
+ for(auto pvIdx : PriVarIndices(cellCenterIdx))
+ {
+ PrimaryVariables priVars(CellCenterPrimaryVariables(localResidual.prevSol()[cellCenterIdx][globalJ]),
+ FacePrimaryVariables(0.0));
+
+ const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, cellCenterIdx);
+ priVars[pvIdx] += eps;
+ ElementSolutionVector elemSol{std::move(priVars)};
+ curVolVars.update(elemSol, problem, elementJ, scvJ);
+
+ auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
+ prevElemVolVars, curElemVolVars,
+ prevGlobalFaceVars, curGlobalFaceVars,
+ elemBcTypes, elemFluxVarsCache);
+
+ auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]);
+ partialDeriv /= eps;
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[faceIdx][cellCenterIdx], faceGlobalI, globalJ, pvIdx, partialDeriv);
+
+ // restore the original volVars
+ curVolVars = origVolVars;
+ }
+ }
+ }
+}
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
+/*!
+* \brief Computes the derivatives of the face residuals with respect to cell center dofs
+*/
+template<class Assembler>
+static void dFacedFace_(Assembler& assembler,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevGlobalFaceVars,
+ GlobalFaceVars& curGlobalFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ const FaceSolutionVector& cachedResidual)
+{
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::FaceIdx faceIdx;
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+ const auto& problem = assembler.problem();
+ auto& localResidual = assembler.localResidual();
+ const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
- const bool isStationary = localResidual.isStationary();
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- if (!isStationary)
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ // set the actual dof index
+ const auto faceGlobalI = scvf.dofIndex();
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
+ // build derivatives with for face dofs w.r.t. cell center dofs
+ for(const auto& globalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
+ {
+ // get the faceVars of the face with respect to which we are going to build the derivative
+ auto origFaceVars = curGlobalFaceVars.faceVars(globalJ);
+ auto& curFaceVars = curGlobalFaceVars.faceVars(globalJ);
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
+ for(auto pvIdx : PriVarIndices(faceIdx))
+ {
+ PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(localResidual.prevSol()[faceIdx][globalJ]));
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+ const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, faceIdx);
+ priVars[pvIdx] += eps;
+ curFaceVars.update(priVars[faceIdx]);
- // the actual element's current residual (will be returned by this function)
- NumEqVector residual(0.0);
- if (!isGhost)
- {
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
+ auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
+ prevElemVolVars, curElemVolVars,
+ prevGlobalFaceVars, curGlobalFaceVars,
+ elemBcTypes, elemFluxVarsCache);
- // get reference to the element's current vol vars
- const auto& scv = fvGeometry.scv(globalI);
- const auto& volVars = curElemVolVars[scv];
+ auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]);
+ partialDeriv /= eps;
- // if the problem is instationary, add derivative of storage term
- if (!isStationary)
- localResidual.addStorageDerivatives(A[globalI][globalI],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
-
- // add source term derivatives
- localResidual.addSourceDerivatives(A[globalI][globalI],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
-
- // add flux derivatives for each scvf
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (!scvf.boundary())
- {
- localResidual.addFluxDerivatives(A[globalI],
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
- }
- else
- {
- const auto& bcTypes = problem.boundaryTypes(element, scvf);
-
- // add Dirichlet boundary flux derivatives
- if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
- {
- localResidual.addCCDirichletFluxDerivatives(A[globalI],
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
- }
- // add Robin ("solution dependent Neumann") boundary flux derivatives
- else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
- {
- localResidual.addRobinFluxDerivatives(A[globalI],
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
- }
- else
- DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
- }
- }
+ // update the global jacobian matrix with the current partial derivatives
+ updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
- // TODO Do we really need this??????????
- // this->model_().updatePVWeights(fvGeometry);
+ // restore the original faceVars
+ curFaceVars = origFaceVars;
+ }
+ }
+ }
+}
- // TODO: Additional dof dependencies???
- // return element residual
- return residual;
- }
-};
-//! explicit assembler using analytic differentiation
-template<class TypeTag>
-class StaggeredLocalAssembler<TypeTag,
- DiffMethod::analytic,
- /*implicit=*/false>
+static Scalar numericEpsilon(const Scalar priVar, const int idx1, const int idx2)
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using IndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ // define the base epsilon as the geometric mean of 1 and the
+ // resolution of the scalar type. E.g. for standard 64 bit
+ // floating point values, the resolution is about 10^-16 and
+ // the base epsilon is thus approximately 10^-8.
+ /*
+ static const Scalar baseEps
+ = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
+ */
+ using BaseEpsilon = typename GET_PROP(TypeTag, BaseEpsilon);
+ const std::array<std::array<Scalar, 2>, 2> baseEps_ = BaseEpsilon::getEps();
+
+
+ static const Scalar baseEps = baseEps_[idx1][idx2];
+ assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
+ // the epsilon value used for the numeric differentiation is
+ // now scaled by the absolute value of the primary variable...
+ return baseEps*(std::abs(priVar) + 1.0);
+}
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-public:
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix. The element residual is written into the right hand side.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
- }
-
-private:
-
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
+/*!
+ * \brief Updates the current global Jacobian matrix with the
+ * partial derivatives of all equations in regard to the
+ * primary variable 'pvIdx' at dof 'col'. Specialization for cc methods.
+ */
+template<class SubMatrix, class CCOrFacePrimaryVariables>
+static void updateGlobalJacobian_(SubMatrix& matrix,
+ const int globalI,
+ const int globalJ,
+ const int pvIdx,
+ const CCOrFacePrimaryVariables &partialDeriv)
+{
+ for (int eqIdx = 0; eqIdx < partialDeriv.size(); eqIdx++)
{
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
-
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's previous time step residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
-
- return residual;
+ // A[i][col][eqIdx][pvIdx] is the rate of change of
+ // the residual of equation 'eqIdx' at dof 'i'
+ // depending on the primary variable 'pvIdx' at dof
+ // 'col'.
+
+ assert(pvIdx >= 0);
+ assert(eqIdx < matrix[globalI][globalJ].size());
+ assert(pvIdx < matrix[globalI][globalJ][eqIdx].size());
+ matrix[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ Valgrind::CheckDefined(matrix[globalI][globalJ][eqIdx][pvIdx]);
}
+}
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- // the actual element's previous time step residual
- NumEqVector residual(0.0), storageResidual(0.0);
- if (!isGhost)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
- }
- // get reference to the element's current vol vars
- const auto& scv = fvGeometry.scv(globalI);
- const auto& volVars = curElemVolVars[scv];
-
- // add derivative of storage term
- localResidual.addStorageDerivatives(A[globalI][globalI],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
+private:
+ template<class T = TypeTag>
+ static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return elemVolVars[scv]; }
- // return the original residual
- return residual;
- }
+ template<class T = TypeTag>
+ static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
+ getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return gridVolVars.volVars(scv); }
};
} // end namespace Dumux
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index 4635bb46e161dbf51ca74f36755801e45dbd0228..86a2e12b126d17e51a17b5bf7f26fc4ce602aa7c 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -117,7 +117,7 @@ public:
const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
- const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
+ const Scalar upWindWeight = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
flux = (upWindWeight * upstreamVolVars.density() +
(1.0 - upWindWeight) * downstreamVolVars.density()) * velocity;
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index 1a4a112db5d87b3e5cea03d275711c5fc3105b58..61caa336cbad9d1224a4142b52618a3cc001ab73 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -96,6 +96,10 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false> : public Dumux::Staggere
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
+ using CellCenterResidual = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FaceResidual = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using FaceResidualVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+
enum {
// grid and world dimension
dim = GridView::dimension,
@@ -119,29 +123,32 @@ public:
using ParentType::ParentType;
- CellCenterPrimaryVariables computeFluxForCellCenter(const Element &element,
+ CellCenterPrimaryVariables computeFluxForCellCenter(const Problem& problem,
+ const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars,
const SubControlVolumeFace &scvf,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
- CellCenterPrimaryVariables flux = fluxVars.computeFluxForCellCenter(this->problem(), element, fvGeometry, elemVolVars,
+ CellCenterPrimaryVariables flux = fluxVars.computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars,
globalFaceVars, scvf, elemFluxVarsCache[scvf]);
- EnergyFluxVariables::energyFlux(flux, this->problem(), element, fvGeometry, elemVolVars, globalFaceVars, scvf, elemFluxVarsCache[scvf]);
+ EnergyFluxVariables::energyFlux(flux, problem, element, fvGeometry, elemVolVars, globalFaceVars, scvf, elemFluxVarsCache[scvf]);
return flux;
}
- CellCenterPrimaryVariables computeSourceForCellCenter(const Element &element,
+ CellCenterPrimaryVariables computeSourceForCellCenter(const Problem& problem,
+ const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars,
- const SubControlVolume &scv)
+ const SubControlVolume &scv) const
{
return CellCenterPrimaryVariables(0.0);
+ // TODO sources
}
@@ -184,16 +191,17 @@ public:
return storage;
}
- FacePrimaryVariables computeSourceForFace(const SubControlVolumeFace& scvf,
+ FacePrimaryVariables computeSourceForFace(const Problem& problem,
+ const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars)
+ const GlobalFaceVars& globalFaceVars) const
{
FacePrimaryVariables source(0.0);
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
- source += this->problem().gravity()[scvf.directionIndex()] * insideVolVars.density();
+ source += problem.gravity()[scvf.directionIndex()] * insideVolVars.density();
- source += this->problem().sourceAtPos(scvf.center())[faceIdx][scvf.directionIndex()];
+ source += problem.sourceAtPos(scvf.center())[faceIdx][scvf.directionIndex()];
return source;
}
@@ -205,18 +213,19 @@ public:
* \param elemVolVars All volume variables for the element
* \param globalFaceVars The face variables
*/
- FacePrimaryVariables computeFluxForFace(const Element& element,
+ FacePrimaryVariables computeFluxForFace(const Problem& problem,
+ const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FacePrimaryVariables flux(0.0);
FluxVariables fluxVars;
- flux += fluxVars.computeNormalMomentumFlux(this->problem(), element, scvf, fvGeometry, elemVolVars, globalFaceVars);
- flux += fluxVars.computeTangetialMomentumFlux(this->problem(), element, scvf, fvGeometry, elemVolVars, globalFaceVars);
- flux += computePressureTerm_(element, scvf, fvGeometry, elemVolVars, globalFaceVars);
+ flux += fluxVars.computeNormalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars);
+ flux += fluxVars.computeTangetialMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars);
+ flux += computePressureTerm_(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars);
return flux;
}
@@ -242,21 +251,23 @@ protected:
/*!
* \brief Evaluate boundary conditions for a cell center dof
*/
- void evalBoundaryForCellCenter_(const Element& element,
+ void evalBoundaryForCellCenter_(CellCenterResidual& residual,
+ const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& faceVars,
const ElementBoundaryTypes& elemBcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
{
- auto boundaryFlux = computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
+ auto boundaryFlux = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
// handle the actual boundary conditions:
- const auto bcTypes = this->problem().boundaryTypes(element, scvf);
+ const auto bcTypes = problem.boundaryTypes(element, scvf);
if(bcTypes.hasNeumann())
{
@@ -265,14 +276,14 @@ protected:
if(bcTypes.isNeumann(eqIdx))
{
const auto extrusionFactor = 1.0; //TODO: get correct extrusion factor
- boundaryFlux[eqIdx] = this->problem().neumannAtPos(scvf.center())[cellCenterIdx][eqIdx]
+ boundaryFlux[eqIdx] = problem.neumannAtPos(scvf.center())[cellCenterIdx][eqIdx]
* extrusionFactor * scvf.area();
}
}
- this->ccResidual_ += boundaryFlux;
+ residual += boundaryFlux;
- asImp_().setFixedCell_(fvGeometry.scv(scvf.insideScvIdx()), elemVolVars, bcTypes);
+ asImp_().setFixedCell_(residual, problem, fvGeometry.scv(scvf.insideScvIdx()), elemVolVars, bcTypes);
}
}
}
@@ -285,40 +296,44 @@ protected:
* \param elemVolVars The current or previous element volVars
* \param bcTypes The boundary types
*/
- void setFixedCell_(const SubControlVolume& insideScv,
+ void setFixedCell_(CellCenterResidual& residual,
+ const Problem& problem,
+ const SubControlVolume& insideScv,
const ElementVolumeVariables& elemVolVars,
- const BoundaryTypes& bcTypes)
+ const BoundaryTypes& bcTypes) const
{
// set a fixed pressure for cells adjacent to a wall
if(bcTypes.isDirichletCell(massBalanceIdx))
{
const auto& insideVolVars = elemVolVars[insideScv];
- this->ccResidual_[pressureIdx] = insideVolVars.pressure() - this->problem().dirichletAtPos(insideScv.center())[cellCenterIdx][pressureIdx];
+ residual[pressureIdx] = insideVolVars.pressure() - problem.dirichletAtPos(insideScv.center())[cellCenterIdx][pressureIdx];
}
}
/*!
* \brief Evaluate boundary conditions for a face dof
*/
- void evalBoundaryForFace_(const Element& element,
+ void evalBoundaryForFace_(FaceResidual& residual,
+ const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& faceVars,
const ElementBoundaryTypes& elemBcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
if (scvf.boundary())
{
// handle the actual boundary conditions:
- const auto bcTypes = this->problem().boundaryTypes(element, scvf);
+ const auto bcTypes = problem.boundaryTypes(element, scvf);
// set a fixed value for the velocity for Dirichlet boundary conditions
if(bcTypes.isDirichlet(momentumBalanceIdx))
{
const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
- const Scalar dirichletValue = this->problem().dirichlet(element, scvf)[faceIdx][scvf.directionIndex()];
- this->faceResiduals_[scvf.localFaceIdx()] = velocity - dirichletValue;
+ const Scalar dirichletValue = problem.dirichlet(element, scvf)[faceIdx][scvf.directionIndex()];
+ residual = velocity - dirichletValue;
}
// For symmetry boundary conditions, there is no flow accross the boundary and
@@ -327,14 +342,14 @@ protected:
{
const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
const Scalar fixedValue = 0.0;
- this->faceResiduals_[scvf.localFaceIdx()] = velocity - fixedValue;
+ residual = velocity - fixedValue;
}
// outflow condition for the momentum balance equation
if(bcTypes.isOutflow(momentumBalanceIdx))
{
if(bcTypes.isDirichlet(massBalanceIdx))
- this->faceResiduals_[scvf.localFaceIdx()] += computeFluxForFace(element, scvf, fvGeometry, elemVolVars, faceVars, elemFluxVarsCache);
+ residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, faceVars, elemFluxVarsCache);
else
DUNE_THROW(Dune::InvalidStateException, "Face at " << scvf.center() << " has an outflow BC for the momentum balance but no Dirichlet BC for the pressure!");
}
@@ -353,23 +368,24 @@ private:
* \param elemVolVars All volume variables for the element
* \param globalFaceVars The face variables
*/
- FacePrimaryVariables computePressureTerm_(const Element& element,
+ FacePrimaryVariables computePressureTerm_(const Problem& problem,
+ const Element& element,
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars)
+ const GlobalFaceVars& globalFaceVars) const
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
- const Scalar deltaP = normalizePressure ? this->problem().initialAtPos(scvf.center())[cellCenterIdx][pressureIdx] : 0.0;
+ const Scalar deltaP = normalizePressure ? problem.initialAtPos(scvf.center())[cellCenterIdx][pressureIdx] : 0.0;
Scalar result = (insideVolVars.pressure() - deltaP) * scvf.area() * -1.0 * sign(scvf.outerNormalScalar());
// treat outflow BCs
if(scvf.boundary())
{
- const Scalar pressure = this->problem().dirichlet(element, scvf)[cellCenterIdx][pressureIdx] - deltaP;
+ const Scalar pressure = problem.dirichlet(element, scvf)[cellCenterIdx][pressureIdx] - deltaP;
result += pressure * scvf.area() * sign(scvf.outerNormalScalar());
}
return result;
diff --git a/dumux/freeflow/staggered/vtkoutputfields.hh b/dumux/freeflow/staggered/vtkoutputfields.hh
index 6708e6a8bdd08eb26bf8c6bea193b3ef3c734cf5..5d142ea1acca34be679a944090837b12d25b0432 100644
--- a/dumux/freeflow/staggered/vtkoutputfields.hh
+++ b/dumux/freeflow/staggered/vtkoutputfields.hh
@@ -45,7 +45,8 @@ public:
// vtk.addSecondaryVariable("Sn", [](const VolumeVariables& v){ return v.saturation(Indices::nPhaseIdx); });
// vtk.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(Indices::wPhaseIdx); });
// vtk.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(Indices::nPhaseIdx); });
- vtk.addSecondaryVariable("p", [](const VolumeVariables& v){ return v.pressure(); });
+ // vtk.addSecondaryVariable("p", [](const VolumeVariables& v){ return v.pressure(); });
+ vtk.addVolumeVariable([](const VolumeVariables& v){ return v.pressure(); }, "p");
// vtk.addSecondaryVariable("pc", [](const VolumeVariables& v){ return v.capillaryPressure(); });
// vtk.addSecondaryVariable("rhoW", [](const VolumeVariables& v){ return v.density(Indices::wPhaseIdx); });
// vtk.addSecondaryVariable("rhoN", [](const VolumeVariables& v){ return v.density(Indices::nPhaseIdx); });
diff --git a/dumux/implicit/staggered/localresidual.hh b/dumux/implicit/staggered/localresidual.hh
index 4b8d1be3e9b9edf1f41dfd84eb4787ccf9cf4376..e2f3caf56bb3a04505e51c143fa75bf96ed513af 100644
--- a/dumux/implicit/staggered/localresidual.hh
+++ b/dumux/implicit/staggered/localresidual.hh
@@ -68,6 +68,10 @@ class StaggeredLocalResidual
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using CellCenterResidual = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using FaceResidual = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using FaceResidualVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
@@ -98,89 +102,6 @@ public:
*/
// \{
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \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(this->problem().model().fvGridGeometry());
- // fvGeometry.bind(element);
- //
- // 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 elemFluxVarsCache = localView(problem().model().globalFluxVarsCache());
- // elemFluxVarsCache.bindElement(element, fvGeometry, curElemVolVars);
- //
- // ElementBoundaryTypes bcTypes;
- // bcTypes.update(problem(), element, fvGeometry);
- //
- // auto& curGlobalFaceVars = problem().model().curGlobalFaceVars();
- // auto& prevGlobalFaceVars = problem().model().prevGlobalFaceVars();
- //
- //
- // asImp_().eval(element, fvGeometry,
- // prevElemVolVars, curElemVolVars,
- // prevGlobalFaceVars, curGlobalFaceVars,
- // bcTypes, elemFluxVarsCache);
- // }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param fvGeometry The finite-volume geometry of the element
- * \param prevVolVars The volume averaged variables for all
- * sub-control volumes of the element at the previous
- * time level
- * \param curVolVars The volume averaged variables for all
- * sub-control volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
- */
- void eval(const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevFaceVars,
- const GlobalFaceVars& curFaceVars,
- const ElementBoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
- {
- // resize and reset all face terms
- const auto numScvf = fvGeometry.numScvf();
-
- faceResiduals_.resize(numScvf, false /*copyOldValues*/);
- faceStorageTerms_.resize(numScvf, false /*copyOldValues*/);
- faceResiduals_ = 0.0;
- faceStorageTerms_ = 0.0;
-
- // evaluate the volume terms (storage + source terms)
- for (auto&& scv : scvs(fvGeometry))
- {
- // treat the cell center dof
- evalCellCenter(problem, element, fvGeometry, scv, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes, elemFluxVarsCache);
-
- // now, treat the dofs on the facets:
- for(auto&& scvf : scvfs(fvGeometry))
- {
- evalFace(problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes, elemFluxVarsCache);
- }
- }
- }
-
/*!
* \brief Compute the local residual, i.e. the deviation of the
* equations from zero.
@@ -197,24 +118,26 @@ public:
* \param bcTypes The types of the boundary conditions for all
* vertices of the element
*/
- void evalCellCenter(const Problem& problem,
+ auto evalCellCenter(const Problem& problem,
const Element &element,
const FVElementGeometry& fvGeometry,
- const SubControlVolume& scv,
const ElementVolumeVariables& prevElemVolVars,
const ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevFaceVars,
const GlobalFaceVars& curFaceVars,
const ElementBoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
// reset all terms
- ccResidual_ = 0.0;
- ccStorageTerm_ = 0.0;
+ CellCenterResidual residual;
+ residual = 0.0;
+ // ccStorageTerm_ = 0.0;
+
+ asImp_().evalVolumeTermForCellCenter_(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
+ asImp_().evalFluxesForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ asImp_().evalBoundaryForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
- asImp_().evalVolumeTermForCellCenter_(problem, element, fvGeometry, scv, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
- asImp_().evalFluxesForCellCenter_(problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
- asImp_().evalBoundaryForCellCenter_(problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ return residual;
}
/*!
@@ -233,7 +156,7 @@ public:
* \param bcTypes The types of the boundary conditions for all
* vertices of the element
*/
- void evalFace(const Problem& problem,
+ auto evalFace(const Problem& problem,
const Element &element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
@@ -243,32 +166,38 @@ public:
const GlobalFaceVars& curFaceVars,
const ElementBoundaryTypes &bcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache,
- const bool resizeResidual = false)
+ const bool resizeResidual = false) const
{
- if(resizeResidual)
- {
- const auto numScvf = fvGeometry.numScvf();
- faceResiduals_.resize(numScvf);
- faceStorageTerms_.resize(numScvf);
- }
+ FaceResidual residual;
- faceResiduals_[scvf.localFaceIdx()] = 0.0;
- faceStorageTerms_[scvf.localFaceIdx()] = 0.0;
+ asImp_().evalVolumeTermForFace_(residual, problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
+ asImp_().evalFluxesForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ asImp_().evalBoundaryForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
- asImp_().evalVolumeTermForFace_(problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
- asImp_().evalFluxesForFace_(problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
- asImp_().evalBoundaryForFace_(problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+ return residual;
}
+ /*!
+ * \brief Sets the solution from which to start the time integration. Has to be
+ * called prior to assembly for time-dependent problems.
+ */
+ void setPreviousSolution(const SolutionVector& u)
+ { prevSol_ = &u; }
- const auto& ccResidual() const
- { return ccResidual_; }
-
- const auto& faceResiduals() const
- { return faceResiduals_; }
+ /*!
+ * \brief Return the solution that has been set as the previous one.
+ */
+ const SolutionVector& prevSol() const
+ {
+ assert(prevSol_ && "no solution set for storage term evaluation");
+ return *prevSol_;
+ }
- const auto& faceResidual(const int fIdx) const
- { return faceResiduals_[fIdx]; }
+ /*!
+ * \brief If no solution has been set, we treat the problem as stationary.
+ */
+ bool isStationary() const
+ { return !prevSol_; }
protected:
@@ -276,35 +205,37 @@ protected:
/*!
* \brief Evaluate the flux terms for cell center dofs
*/
- void evalFluxesForCellCenter_(const Problem& problem,
+ void evalFluxesForCellCenter_(CellCenterResidual& residual,
+ const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& faceVars,
const ElementBoundaryTypes& bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
for (auto&& scvf : scvfs(fvGeometry))
{
if(!scvf.boundary())
- ccResidual_ += asImp_().computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
+ residual += asImp_().computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
}
}
/*!
* \brief Evaluate the flux terms for face dofs
*/
- void evalFluxesForFace_(const Problem& problem,
+ void evalFluxesForFace_(FaceResidual& residual,
+ const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars,
const ElementBoundaryTypes& bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
if(!scvf.boundary())
- faceResiduals_[scvf.localFaceIdx()] += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars, elemFluxVarsCache);
+ residual += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars, elemFluxVarsCache);
}
/*!
@@ -328,23 +259,25 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForCellCenter_(const Problem& problem,
+ evalVolumeTermForCellCenter_(CellCenterResidual& residual,
+ const Problem& problem,
const Element &element,
const FVElementGeometry& fvGeometry,
- const SubControlVolume& scv,
const ElementVolumeVariables& prevElemVolVars,
const ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevFaceVars,
const GlobalFaceVars& curFaceVars,
- const ElementBoundaryTypes &bcTypes)
+ const ElementBoundaryTypes &bcTypes) const
{
- const auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor();
-
- // subtract the source term from the local rate
- CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curFaceVars, scv);
- source *= scv.volume()*curExtrusionFactor;
+ for(auto&& scv : scvs(fvGeometry))
+ {
+ const auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor();
- ccResidual_ -= source;
+ // subtract the source term from the local rate
+ CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curFaceVars, scv);
+ source *= scv.volume()*curExtrusionFactor;
+ residual -= source;
+ }
}
/*!
@@ -352,7 +285,8 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForFace_(const Problem& problem,
+ evalVolumeTermForFace_(FaceResidual& residual,
+ const Problem& problem,
const Element &element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
@@ -360,14 +294,14 @@ protected:
const ElementVolumeVariables& curElemVolVars,
const GlobalFaceVars& prevFaceVars,
const GlobalFaceVars& curFaceVars,
- const ElementBoundaryTypes &bcTypes)
+ const ElementBoundaryTypes &bcTypes) const
{
// the source term:
auto faceSource = asImp_().computeSourceForFace(problem, scvf, curElemVolVars, curFaceVars);
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
const auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor();
faceSource *= 0.5*scv.volume()*curExtrusionFactor;
- faceResiduals_[scvf.localFaceIdx()] -= faceSource;
+ residual -= faceSource;
}
/*!
@@ -375,44 +309,49 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForCellCenter_(const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolume& scv,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevFaceVars,
- const GlobalFaceVars& curFaceVars,
- const ElementBoundaryTypes &bcTypes)
+ evalVolumeTermForCellCenter_(CellCenterResidual& residual,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevFaceVars,
+ const GlobalFaceVars& curFaceVars,
+ const ElementBoundaryTypes &bcTypes) const
{
- const auto& curVolVars = curElemVolVars[scv];
- const auto& prevVolVars = prevElemVolVars[scv];
+ for(auto&& scv : scvs(fvGeometry))
+ {
+ 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...
- auto prevCCStorage = asImp_().computeStorageForCellCenter(problem, scv, prevVolVars);
- auto curCCStorage = asImp_().computeStorageForCellCenter(problem, scv, curVolVars);
+ // 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...
+ auto prevCCStorage = asImp_().computeStorageForCellCenter(problem, scv, prevVolVars);
+ auto curCCStorage = asImp_().computeStorageForCellCenter(problem, scv, curVolVars);
- prevCCStorage *= prevVolVars.extrusionFactor();
- curCCStorage *= curVolVars.extrusionFactor();
+ prevCCStorage *= prevVolVars.extrusionFactor();
+ curCCStorage *= curVolVars.extrusionFactor();
- ccStorageTerm_ = std::move(curCCStorage);
- ccStorageTerm_ -= std::move(prevCCStorage);
- ccStorageTerm_ *= scv.volume();
- ccStorageTerm_ /= timeLoop_->timeStepSize();
+ CellCenterResidual storageTerm = 0.0;
- // add the storage term to the residual
- ccResidual_ += ccStorageTerm_;
+ storageTerm = std::move(curCCStorage);
+ storageTerm -= std::move(prevCCStorage);
+ storageTerm *= scv.volume();
+ storageTerm /= timeLoop_->timeStepSize();
- // subtract the source term from the local rate
- CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curFaceVars, scv);
- source *= scv.volume()*curVolVars.extrusionFactor();
+ // add the storage term to the residual
+ residual += storageTerm;
- ccResidual_ -= source;
+ // subtract the source term from the local rate
+ CellCenterPrimaryVariables source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curFaceVars, scv);
+ source *= scv.volume()*curVolVars.extrusionFactor();
+
+ residual -= source;
+ }
}
/*!
@@ -420,15 +359,16 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTermForFace_(const Problem& problem,
- const Element &element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const GlobalFaceVars& prevFaceVars,
- const GlobalFaceVars& curFaceVars,
- const ElementBoundaryTypes &bcTypes)
+ evalVolumeTermForFace_(FaceResidual& residual,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevFaceVars,
+ const GlobalFaceVars& curFaceVars,
+ const ElementBoundaryTypes &bcTypes) const
{
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
const auto& curVolVars = curElemVolVars[scv];
@@ -437,16 +377,16 @@ protected:
auto curFaceStorage = asImp_().computeStorageForFace(problem, scvf, curVolVars, curFaceVars);
// the storage term
- faceStorageTerms_[scvf.localFaceIdx()] = std::move(curFaceStorage);
- faceStorageTerms_[scvf.localFaceIdx()] -= std::move(prevFaceStorage);
- faceStorageTerms_[scvf.localFaceIdx()] *= (scv.volume()/2.0);
- faceStorageTerms_[scvf.localFaceIdx()] /= timeLoop_->timeStepSize();
- faceResiduals_[scvf.localFaceIdx()] += faceStorageTerms_[scvf.localFaceIdx()];
+ residual = std::move(curFaceStorage);
+ residual -= std::move(prevFaceStorage);
+ residual *= (scv.volume()/2.0);
+ residual /= timeLoop_->timeStepSize();
+ // residuals[scvf.localFaceIdx()] += faceStorageTerms_[scvf.localFaceIdx()];
// the source term:
auto faceSource = asImp_().computeSourceForFace(problem, scvf, curElemVolVars, curFaceVars);
faceSource *= 0.5*scv.volume()*curVolVars.extrusionFactor();
- faceResiduals_[scvf.localFaceIdx()] -= faceSource;
+ residual -= faceSource;
}
Implementation &asImp_()
@@ -455,36 +395,7 @@ protected:
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
- CellCenterPrimaryVariables ccResidual_;
- CellCenterPrimaryVariables ccStorageTerm_;
- FaceSolutionVector faceResiduals_;
- FaceSolutionVector faceStorageTerms_;
-
-
-
- /*!
- * \brief Sets the solution from which to start the time integration. Has to be
- * called prior to assembly for time-dependent problems.
- */
- void setPreviousSolution(const SolutionVector& u)
- { prevSol_ = &u; }
-
- /*!
- * \brief Return the solution that has been set as the previous one.
- */
- const SolutionVector& prevSol() const
- {
- assert(prevSol_ && "no solution set for storage term evaluation");
- return *prevSol_;
- }
-
- /*!
- * \brief If no solution has been set, we treat the problem as stationary.
- */
- bool isStationary() const
- { return !prevSol_; }
-protected:
TimeLoop& timeLoop()
{ return *timeLoop_; }
diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh
index d142371dce7c2e81bb6f75dd318a8d635adeb6a3..c10185d12e2594e4980f65b219794acceefcfcb3 100644
--- a/dumux/implicit/staggered/newtoncontroller.hh
+++ b/dumux/implicit/staggered/newtoncontroller.hh
@@ -65,6 +65,9 @@ class StaggeredNewtonController : public NewtonController<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Communicator = typename GridView::CollectiveCommunication;
+
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
@@ -75,8 +78,12 @@ class StaggeredNewtonController : public NewtonController<TypeTag>
};
public:
- StaggeredNewtonController(const Problem &problem)
- : ParentType(problem)
+ StaggeredNewtonController(const Communicator& comm)
+ : ParentType(comm)
+ {}
+
+ StaggeredNewtonController(const Communicator& comm, std::shared_ptr<TimeLoop<Scalar>> timeLoop)
+ : ParentType(comm, timeLoop)
{}
/*!
@@ -92,11 +99,16 @@ public:
* \param x The vector which solves the linear system
* \param b The right hand side of the linear system
*/
- template<typename T = TypeTag>
- typename std::enable_if<!LinearSolverAcceptsMultiTypeMatrix<T>::value, void>::type
- newtonSolveLinear(JacobianMatrix &A,
- SolutionVector &x,
- SolutionVector &b)
+ // template<typename T = TypeTag>
+ // typename std::enable_if<!LinearSolverAcceptsMultiTypeMatrix<T>::value, void>::type
+ // newtonSolveLinear(JacobianMatrix &A,
+ // SolutionVector &x,
+ // SolutionVector &b)
+ template<class LinearSolver, class JacobianMatrix, class SolutionVector>
+ void solveLinearSystem(LinearSolver& ls,
+ JacobianMatrix& A,
+ SolutionVector& x,
+ SolutionVector& b)
{
try
{
@@ -127,7 +139,7 @@ public:
// printmatrix(std::cout, M, "", "");
// solve
- const bool converged = this->linearSolver_.solve(M, y, bTmp);
+ const bool converged = ls.solve(M, y, bTmp);
// copy back the result y into x
VectorConverter<SolutionVector>::retrieveValues(x, y);
@@ -153,29 +165,29 @@ public:
* \param x The vector which solves the linear system
* \param b The right hand side of the linear system
*/
- template<typename T = TypeTag>
- typename std::enable_if<LinearSolverAcceptsMultiTypeMatrix<T>::value, void>::type
- newtonSolveLinear(JacobianMatrix &A,
- SolutionVector &x,
- SolutionVector &b)
- {
- try
- {
- if (this->numSteps_ == 0)
- this->initialResidual_ = b.two_norm();
-
- bool converged = this->linearSolver_.solve(A, x, b);
-
- if (!converged)
- DUNE_THROW(NumericalProblem, "Linear solver did not converge");
- }
- catch (const Dune::Exception &e)
- {
- Dumux::NumericalProblem p;
- p.message(e.what());
- throw p;
- }
- }
+ // template<typename T = TypeTag>
+ // typename std::enable_if<LinearSolverAcceptsMultiTypeMatrix<T>::value, void>::type
+ // newtonSolveLinear(JacobianMatrix &A,
+ // SolutionVector &x,
+ // SolutionVector &b)
+ // {
+ // try
+ // {
+ // if (this->numSteps_ == 0)
+ // this->initialResidual_ = b.two_norm();
+ //
+ // bool converged = this->linearSolver_.solve(A, x, b);
+ //
+ // if (!converged)
+ // DUNE_THROW(NumericalProblem, "Linear solver did not converge");
+ // }
+ // catch (const Dune::Exception &e)
+ // {
+ // Dumux::NumericalProblem p;
+ // p.message(e.what());
+ // throw p;
+ // }
+ // }
/*!
* \brief Update the current solution with a delta vector.
@@ -194,18 +206,20 @@ public:
* system of equations. This parameter also stores
* the updated solution.
*/
- void newtonUpdate(SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
+ template<class JacobianAssembler, class SolutionVector>
+ void newtonUpdate(JacobianAssembler& assembler,
+ SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
{
if (this->enableShiftCriterion_)
this->newtonUpdateShift(uLastIter, deltaU);
- this->writeConvergence_(uLastIter, deltaU);
+ // this->writeConvergence_(uLastIter, deltaU);
if (this->useLineSearch_)
{
- this->lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
+ this->lineSearchUpdate_(assembler, uCurrentIter, uLastIter, deltaU);
}
else {
for (unsigned int i = 0; i < uLastIter[cellCenterIdx].size(); ++i) {
@@ -219,8 +233,8 @@ public:
if (this->enableResidualCriterion_)
{
- SolutionVector tmp(uLastIter);
- this->reduction_ = this->method().model().globalResidual(tmp, uCurrentIter);
+ this->residualNorm_ = assembler.residualNorm(uCurrentIter);
+ this->reduction_ = this->residualNorm_;
this->reduction_ /= this->initialResidual_;
}
}
@@ -242,7 +256,7 @@ public:
auto uNewI = uLastIter[cellCenterIdx][i];
uNewI -= deltaU[cellCenterIdx][i];
- Scalar shiftAtDof = this->model_().relativeShiftAtDof(uLastIter[cellCenterIdx][i],
+ Scalar shiftAtDof = this->relativeShiftAtDof_(uLastIter[cellCenterIdx][i],
uNewI);
this->shift_ = std::max(this->shift_, shiftAtDof);
}
@@ -250,13 +264,13 @@ public:
auto uNewI = uLastIter[faceIdx][i];
uNewI -= deltaU[faceIdx][i];
- Scalar shiftAtDof = this->model_().relativeShiftAtDof(uLastIter[faceIdx][i],
+ Scalar shiftAtDof = this->relativeShiftAtDof_(uLastIter[faceIdx][i],
uNewI);
this->shift_ = std::max(this->shift_, shiftAtDof);
}
- if (this->gridView_().comm().size() > 1)
- this->shift_ = this->gridView_().comm().max(this->shift_);
+ if (this->communicator().size() > 1)
+ this->shift_ = this->communicator().max(this->shift_);
}
};
diff --git a/test/freeflow/staggered/doneatestproblem.hh b/test/freeflow/staggered/doneatestproblem.hh
index 5aaff6433a23ab73c5947fcaddbfcde3c16a72c6..019360ee5b13e16644c11345acc0122ea1980b2b 100644
--- a/test/freeflow/staggered/doneatestproblem.hh
+++ b/test/freeflow/staggered/doneatestproblem.hh
@@ -274,8 +274,8 @@ public:
InitialValues initialAtPos(const GlobalPosition& globalPos) const
{
InitialValues values;
- values[pressureIdx] = 123.0;
- values[velocityXIdx] = 44.0;
+ values[pressureIdx] = 0.0;
+ values[velocityXIdx] = 0.0;
values[velocityYIdx] = 0.0;
return values;
diff --git a/test/freeflow/staggered/test_donea.cc b/test/freeflow/staggered/test_donea.cc
index a0a1d7044925dc149010bf00077f883fa06c4704..14e4edcf567a277c0a2645d18d29eddc3fd3f931 100644
--- a/test/freeflow/staggered/test_donea.cc
+++ b/test/freeflow/staggered/test_donea.cc
@@ -124,7 +124,7 @@ int main(int argc, char** argv)
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ // using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
@@ -167,5 +167,73 @@ int main(int argc, char** argv)
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
- assembler->setLinearSystem();
+ // assembler->setLinearSystem();
+
+ // the linear solver
+ using LinearSolver = Dumux::UMFPackBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
+ using NewtonMethod = Dumux::NewtonMethod<TypeTag, NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
}