Re: [DuMux] Computing Boundary Fluxes
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(curSol, gridVars, this->gridGeometry());
// Do whatever you want with the result struct
}
On Mon, Oct 17, 2022 at 3:53 PM Timo Koch wrote:
> Dear Rafael,
>
> welcome to the Dumux mailing list!
>
> > On 14 Oct 2022, at 23:38, Rafael March 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
Re: [DuMux] Computing Boundary Fluxes
Dear Rafael,
welcome to the Dumux mailing list!
> On 14 Oct 2022, at 23:38, Rafael March 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
> DuMux@listserv.uni-stuttgart.de
> 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
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(, 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 +=
[DuMux] Computing Boundary Fluxes
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? Thank you, Rafael March. ___ DuMux mailing list DuMux@listserv.uni-stuttgart.de https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
