Thanks Timo! I really appreciate your help. It would have taken me a while
to figure this out. For reference, here's how I adapted your code to my own
problem. Maybe this will help someone else in the future.
1* - *Included Results struct and computeOutput as public members of the
Problem class (they could actually be private since I don't call them from
the main file)
2* - *Added fields in the Results struct to accumulate Oil flux (second
position of the massFlux array)
3* - *Accumulated oil flux in the new struct fields (see LISTING 1 below)
4* - *Created the public method computeBoundaryFluxes in the Problem class
(see LISTING 2 below)
5* - *Called computeBoundaryFluxes from the time integration loop inside
the main function:
problem->computeBoundaryFluxes(x, *gridVariables);
I still don't get the same total flux in and out, but I suppose this is due
to compressibility. Next step is to change the fluid system to an
incompressible one.
Thanks again!
Rafael.
LISTING 1 - EXCERPT INSIDE computeOutput METHOD TO ACCUMULATE OIL FLUX
if (isRightBoundary(scvf.ipGlobal()))
{
const auto& insideVolVars =
elemVolVars[scvf.insideScvIdx()];
output.topArea +=
scvf.area()*insideVolVars.extrusionFactor();
if (problem.boundaryTypes(element,
scvf).hasNeumann())
{
const auto massFlux = problem.neumann(element,
fvGeometry, elemVolVars, elemFluxVarsCache,
scvf)*scvf.area()*insideVolVars.extrusionFactor();
output.topBoundaryWaterFlux += massFlux[0];
output.topBoundaryOilFlux += massFlux[1];
}
else
{
const auto massFlux =
localResidual.computeFlux(problem, element, fvGeometry, elemVolVars, scvf,
elemFluxVarsCache);
output.topBoundaryWaterFlux += massFlux[0];
output.topBoundaryOilFlux += massFlux[1];
}
}
else if (isLeftBoundary(scvf.ipGlobal()))
{
const auto& insideVolVars =
elemVolVars[scvf.insideScvIdx()];
output.bottomArea +=
scvf.area()*insideVolVars.extrusionFactor();
if (problem.boundaryTypes(element,
scvf).hasNeumann())
{
const auto massFlux = problem.neumann(element,
fvGeometry, elemVolVars, elemFluxVarsCache,
scvf)*scvf.area()*insideVolVars.extrusionFactor();
output.bottomBoundaryWaterFlux += massFlux[0];
output.bottomBoundaryOilFlux += massFlux[1];
}
else
{
const auto massFlux =
localResidual.computeFlux(problem, element, fvGeometry, elemVolVars, scvf,
elemFluxVarsCache);
output.bottomBoundaryWaterFlux += massFlux[0];
output.bottomBoundaryOilFlux += massFlux[1];
}
}
LISTING 2 - computeBoundaryFluxes METHOD
void computeBoundaryFluxes(const SolutionVector& curSol, const
GridVariables& gridVars)
{
Results res;
res = computeOutput<LocalResidual, SolutionVector, GridVariables,
GridGeometry>(curSol, gridVars, this->gridGeometry());
// Do whatever you want with the result struct
}
On Mon, Oct 17, 2022 at 3:53 PM Timo Koch <[email protected]> wrote:
> Dear Rafael,
>
> welcome to the Dumux mailing list!
>
> > On 14 Oct 2022, at 23:38, Rafael March <[email protected]> wrote:
> >
> > Hello,
> >
> > I hope this email finds you well.
> >
> > I'm trying to use DuMux to perform single and multiphase flow-based
> upscaling (Darcy scale).
> >
> > This involves running two-phase flow simulations until steady-state is
> reached, then using the steady-state saturation field to compute effective
> permeabilities for different fractional flow levels.
> >
> > I managed to implement cell-wise relative permeability and capillary
> pressures, as well as heterogeneous permeability and porosity fields.
> >
> > However, I'm still struggling to compute boundary fluxes for pressure
> dirichlet boundary faces. These boundary fluxes will determine the
> effective permeability for the whole model. Could you point me to an
> example or code excerpt that computes the total volumetric fluxes for the
> inlet and outlet boundaries?
>
> This will depend a bit on the discretisation method. It seems for your
> case cell-entered TPFA finite volumes might be suitable, then computing
> fluxes over the boundaries can be achieved the same way it is done in the
> local residual.
> Admittedly that can be a bit more complicated then it should be.
>
> I attached a snippet below that should roughly do what you want. On
> Dirichlet boundaries the LocalResidual is reused to assemble the flux.
> Neumann and computeFlux will return a vector with a value for each
> equation in case of two-phase flow, so you might need to adapt the code a
> bit to sum up to fluxes of the phase that you are interested in.
>
> Let me know if you have questions.
>
> Best wishes
> Timo
>
> >
> > Thank you,
> > Rafael March.
> >
> >
> > _______________________________________________
> > DuMux mailing list
> > [email protected]
> > https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
>
>
> struct Results
> {
> double topBoundaryWaterFlux = 0.0;
> double bottomBoundaryWaterFlux = 0.0;
> double topArea = 0.0;
> double bottomArea = 0.0;
> };
>
> template<class LocalResidual, class SolutionVector, class GridVariables,
> class GridGeometry>
> Results computeOutput(const SolutionVector& curSol, const GridVariables&
> gridVars, const GridGeometry& gridGeometry)
> {
> static constexpr int dimWorld = GridGeometry::GridView::dimensionworld;
>
> Results output;
> const auto& problem = gridVars.curGridVolVars().problem();
> for (const auto& element : elements(gridGeometry.gridView(),
> Dune::Partitions::interior))
> {
> auto fvGeometry = localView(gridGeometry);
> fvGeometry.bind(element);
>
> auto elemVolVars = localView(gridVars.curGridVolVars());
> elemVolVars.bind(element, fvGeometry, curSol);
>
> const auto isTopBoundary = [&](const auto& pos) -> bool { return
> pos[dimWorld-1] > gridGeometry.bBoxMax()[dimWorld-1] - 1e-10; };
> const auto isBottomBoundary = [&](const auto& pos) -> bool {
> return pos[dimWorld-1] < gridGeometry.bBoxMin()[dimWorld-1] + 1e-10; };
> LocalResidual localResidual(&problem, nullptr);
>
> if (!element.hasBoundaryIntersections())
> continue;
>
> auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
> elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
>
> for (const auto& scvf : scvfs(fvGeometry))
> {
> if (scvf.boundary())
> {
> if (isTopBoundary(scvf.ipGlobal()))
> {
> const auto& insideVolVars =
> elemVolVars[scvf.insideScvIdx()];
> output.topArea +=
> scvf.area()*insideVolVars.extrusionFactor();
> if (problem.boundaryTypes(element, scvf).hasNeumann())
> {
> const auto massFlux = problem.neumann(element,
> fvGeometry, elemVolVars, elemFluxVarsCache,
> scvf)*scvf.area()*insideVolVars.extrusionFactor();
> output.topBoundaryWaterFlux += massFlux[0];
> }
> else
> {
> const auto massFlux =
> localResidual.computeFlux(problem, element, fvGeometry, elemVolVars, scvf,
> elemFluxVarsCache);
> output.topBoundaryWaterFlux += massFlux[0];
> }
> }
>
> else if (isBottomBoundary(scvf.ipGlobal()))
> {
> const auto& insideVolVars =
> elemVolVars[scvf.insideScvIdx()];
> output.bottomArea +=
> scvf.area()*insideVolVars.extrusionFactor();
> if (problem.boundaryTypes(element, scvf).hasNeumann())
> {
> const auto massFlux = problem.neumann(element,
> fvGeometry, elemVolVars, elemFluxVarsCache,
> scvf)*scvf.area()*insideVolVars.extrusionFactor();
> output.bottomBoundaryWaterFlux += massFlux[0];
> }
> else
> {
> const auto massFlux =
> localResidual.computeFlux(problem, element, fvGeometry, elemVolVars, scvf,
> elemFluxVarsCache);
> output.bottomBoundaryWaterFlux += massFlux[0];
> }
> }
> }
> }
> }
>
> const auto& comm = gridGeometry.gridView().comm();
> output.topBoundaryWaterFlux = comm.sum(output.topBoundaryWaterFlux);
> output.bottomBoundaryWaterFlux =
> comm.sum(output.bottomBoundaryWaterFlux);
> output.topArea = comm.sum(output.topArea);
> output.bottomArea = comm.sum(output.bottomArea);
> output.totalVolume = comm.sum(output.totalVolume);
>
> return output;
> }
>
>
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
