Hi Nikolai,
I think your implementation is almost correct. If you want to balance
masses, you have to add the density to your upwind terms though, i.e.:
auto upwindTermN = [](const auto& volVars) { return
volVars.mobility(nPhaseIdx)*volVars.density(nPhaseIdx); };
...
Also, you have to call bind() instead of bindElement() on the geometry
and the volume variables, i.e.:
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bind(element);
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bind(element, fvGeometry, curSol);
This has the following background:
If you want to do computations on the element without assembling fluxes,
you do not have to prepare all variables within the entire stencil. For
these purposes you call bindElement(). This only prepares
primary/secondary variables on the element. This is enough e.g. to
evaluate the mass inside the domain after a time step. But, if you want
to compute fluxes, you need to prepare the entire stencil. This stencil
includes neighboring elements as well for cell-centered schemes. In this
case you prepare you local views accordingly by calling bind() instead
of bindElement().
I hope this makes it work!
Best wishes,
Dennis
On 21.08.2018 10:48, Nikolai Andrianov wrote:
Dear DuMuX experts,
Please advise how can I implement a global conservation error check,
that is, evaluate the residual ( mass(t^{n+1}) - mass(t^n) ) / dt +
sum of fluxes over the domain boundaries.
So far I can evaluate the mass in the domain by looping through the
elements as in void getVolumes of
https://git.iws.uni-stuttgart.de/andrian/rate-sens-nofrac/blob/master/src/rate_problem.hh.
I tried to use the following to access the fluxes in rate_problem.hh
(inspired from porousmediumflow/velocityoutput.hh and
porousmediumflow/nonequilibrium/gridvariables.hh):
// the upwind term to be used for the volume flux evaluation
auto upwindTermN = [](const auto& volVars) { return
volVars.mobility(nPhaseIdx); };
auto upwindTermW = [](const auto& volVars) { return
volVars.mobility(wPhaseIdx); };
auto elemFluxVarsCache =
localView(gridVariables.gridFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
for (auto&& scvf : scvfs(fvGeometry)) {
if (scvf.boundary()) {
//auto bcTypes = problemBoundaryTypes(element, scvf);
//if (bcTypes.hasOnlyDirichlet())
{
FluxVariables fluxVars;
fluxVars.init(*this, element, fvGeometry,
elemVolVars, scvf, elemFluxVarsCache);
Scalar extrusionFactor =
problem_.extrusionFactor(element, fvGeometry.scv(scvf.insideScvIdx()),
elementSolution(element, sol_, fvGridGeometry_));
FO += fluxVars.advectiveFlux(nPhaseIdx,
upwindTermN);// / extrusionFactor;
FW += fluxVars.advectiveFlux(wPhaseIdx,
upwindTermW);// / extrusionFactor;
}
}
}
However, I get compilation error at elemFluxVarsCache.bind(element,
fvGeometry, elemVolVars) saying that Assertion `it != indices.end() &&
"Could not find the scv/scvf! Make sure to properly bind this class!"'
failed.
Your feedback is greatly appreciated.
Many thanks,
Nikolai
PS: please disregard my previous mail with the same subject..
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
