Hi Dennis,
Thanks a lot, it works now. Best regards, Nikolai ________________________________ From: Dumux <[email protected]> on behalf of Dennis Gläser <[email protected]> Sent: Tuesday, August 21, 2018 2:03:46 PM To: [email protected] Subject: Re: [DuMuX] Check global conservation error Hi Nikolai, you have to make sure that you only compute this on Dirichlet faces. On Neumann faces the fluxes are the ones you set in your neumann() or neumannAtPos() function and cannot be computed by the flux assembly. So if you have only no-flow Neumann boundaries just make sure to only compute boundary fluxes on Dirichlet faces and your calculations will be correct. If you also consider Neumann fluxes then integrate these fluxes by calling neumann() and multiplying the values with scvf.area() and the extrusion factor. I hope this helps! Dennis On 21.08.2018 13:59, Nikolai Andrianov wrote: Hi Dennis, Thanks a lot, I got through elemFluxVarsCache.bind runtime error! However, a new problem has popped up đ When I get to computing (or recovering) the advective flux with fluxVars.advectiveFlux(nPhaseIdx, upwindTermN), the execution brings me to caluclating Dracy's flux in CCTpfaDarcysLaw::flux(...) since - !advFluxIsCached_[phaseIdx] is true in PorousMediumFluxVariables::advectiveFlux(..) and - (!problem.couplingManager().isOnInteriorBoundary(element, scvf)) is true in CCTpfaFacetCouplingDarcysLawImpl::flux(..) So I am ending up in CCTpfaDarcysLaw::flux(...) and the outsideVolVars = elemVolVars[scvf.outsideScvIdx()] cannot be computed because scvf is a boundary face, resulting in the error message "Assertion `it != volVarIndices_.end() && "Could not find the current volume variables for volVarIdx!"' failed." Looking forward for your feedback. Thanks, Nikolai ________________________________ From: Dumux <[email protected]><mailto:[email protected]> on behalf of Dennis Gläser <[email protected]><mailto:[email protected]> Sent: Tuesday, August 21, 2018 11:39:34 AM To: [email protected]<mailto:[email protected]> Subject: Re: [DuMuX] Check global conservation error 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]<mailto:[email protected]> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux _______________________________________________ Dumux mailing list [email protected]<mailto:[email protected]> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________ Dumux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
