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
