diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index ec65166b1..f8b9e97e1 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -131,6 +131,15 @@ public: enableBrine, enableSaltPrecipitation, Indices::numPhases>; + using ScalarFluidState = BlackOilFluidState; using Problem = GetPropType; BlackOilIntensiveQuantities() diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index e5fc3242b..4537c9d48 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -61,6 +61,7 @@ class BlackOilLocalResidualTPFA : public GetPropType; using GridView = GetPropType; using Problem = GetPropType; + using FluidState = typename IntensiveQuantities::FluidState; enum { conti0EqIdx = Indices::conti0EqIdx }; enum { numEq = getPropValue() }; @@ -414,6 +415,114 @@ public: } + template + static void computeBoundaryFlux(RateVector& bdyFlux, + const Problem& problem, + const BoundaryConditionData& bdyInfo, + const IntensiveQuantities& insideIntQuants, + unsigned globalSpaceIdx) + { + if (bdyInfo.type == BCType::RATE) { + computeBoundaryFluxRate(bdyFlux, bdyInfo); + } else if (bdyInfo.type == BCType::FREE) { + computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx); + } else { + throw std::logic_error("Unknown boundary condition type " + std::to_string(static_cast(bdyInfo.type)) + " in computeBoundaryFlux()." ); + } + } + + template + static void computeBoundaryFluxRate(RateVector& bdyFlux, + const BoundaryConditionData& bdyInfo) + { + bdyFlux.setMassRate(bdyInfo.massRate, bdyInfo.pvtRegionIdx); + } + + template + static void computeBoundaryFluxFree(const Problem& problem, + RateVector& bdyFlux, + const BoundaryConditionData& bdyInfo, + const IntensiveQuantities& insideIntQuants, + unsigned globalSpaceIdx) + { + std::array upIdx; + std::array dnIdx; + RateVector volumeFlux; + RateVector pressureDifference; + ExtensiveQuantities::calculateBoundaryGradients_(problem, + globalSpaceIdx, + insideIntQuants, + bdyInfo.boundaryFaceIndex, + bdyInfo.faceArea, + bdyInfo.faceZCoord, + bdyInfo.exFluidState, + upIdx, + dnIdx, + volumeFlux, + pressureDifference); + + //////// + // advective fluxes of all components in all phases + //////// + bdyFlux = 0.0; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + const auto& pBoundary = bdyInfo.exFluidState.pressure(phaseIdx); + const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx); + const unsigned pvtRegionIdx = insideIntQuants.pvtRegionIndex(); + + RateVector tmp; + + // mass conservation + if (pBoundary < pInside) { + // outflux + const auto& invB = getInvB_(insideIntQuants.fluidState(), phaseIdx, pvtRegionIdx); + Evaluation surfaceVolumeFlux = invB * volumeFlux[phaseIdx]; + evalPhaseFluxes_(tmp, + phaseIdx, + insideIntQuants.pvtRegionIndex(), + surfaceVolumeFlux, + insideIntQuants.fluidState()); + } else if (pBoundary > pInside) { + // influx + using ScalarFluidState = decltype(bdyInfo.exFluidState); + const auto& invB = getInvB_(bdyInfo.exFluidState, phaseIdx, pvtRegionIdx); + Evaluation surfaceVolumeFlux = invB * volumeFlux[phaseIdx]; + evalPhaseFluxes_(tmp, + phaseIdx, + insideIntQuants.pvtRegionIndex(), + surfaceVolumeFlux, + bdyInfo.exFluidState); + } + + for (unsigned i = 0; i < tmp.size(); ++i) { + bdyFlux[i] += tmp[i]; + } + + static_assert(!enableEnergy, "Relevant treatment of boundary conditions must be implemented before enabling."); + // Add energy flux treatment per phase here. + } + + static_assert(!enableSolvent, "Relevant treatment of boundary conditions must be implemented before enabling."); + static_assert(!enablePolymer, "Relevant treatment of boundary conditions must be implemented before enabling."); + static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling."); + + // make sure that the right mass conservation quantities are used + adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex()); + + // heat conduction + static_assert(!enableEnergy, "Relevant treatment of boundary conditions must be implemented before enabling."); + +#ifndef NDEBUG + for (unsigned i = 0; i < numEq; ++i) { + Valgrind::CheckDefined(bdyFlux[i]); + } + Valgrind::CheckDefined(bdyFlux); +#endif + } + static void computeSource(RateVector& source, const Problem& problem, unsigned globalSpaceIdex, diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 809d37dee..250a6e41e 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -91,6 +91,7 @@ class TpfaLinearizer enum { numEq = getPropValue() }; enum { historySize = getPropValue() }; + enum { dimWorld = GridView::dimensionworld }; using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock; using VectorBlock = Dune::FieldVector; @@ -343,6 +344,28 @@ private: } } neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end()); + for (unsigned bfIndex = 0; bfIndex < stencil.numBoundaryFaces(); ++bfIndex) { + const auto& bf = stencil.boundaryFace(bfIndex); + const int dir_id = bf.dirId(); + const auto [free, massrateAD] = problem_().boundaryCondition(myIdx, dir_id); + // Strip the unnecessary (and zero anyway) derivatives off massrate. + VectorBlock massrate(0.0); + for (size_t ii = 0; ii < massrate.size(); ++ii) { + massrate[ii] = massrateAD[ii].value(); + } + const bool nonzero_massrate = massrate != VectorBlock(0.0); + if (free || nonzero_massrate) { + const auto& exFluidState = problem_().initialFluidState(myIdx); + BoundaryConditionData bcdata{free ? BCType::FREE : BCType::RATE, + massrate, + exFluidState.pvtRegionIndex(), + bfIndex, + bf.area(), + bf.integrationPos()[dimWorld - 1], + exFluidState}; + boundaryInfo_.push_back({myIdx, bcdata}); + } + } } } @@ -458,6 +481,23 @@ private: residual_[globI] += res; jacobian_->addToBlock(globI, globI, bMat); } + } // end of loop for cell globI. + + // Boundary terms. Only looping over cells with nontrivial bcs. + for (const auto& bdyInfo : boundaryInfo_) { + VectorBlock res(0.0); + MatrixBlock bMat(0.0); + ADVectorBlock adres(0.0); + const unsigned globI = bdyInfo.cell; + const IntensiveQuantities* insideIntQuants = model_().cachedIntensiveQuantities(globI, /*timeIdx*/ 0); + if (insideIntQuants == nullptr) { + throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI)); + } + LocalResidual::computeBoundaryFlux(adres, problem_(), bdyInfo.bcdata, *insideIntQuants, globI); + adres *= bdyInfo.bcdata.faceArea; + setResAndJacobi(res, bMat, adres); + residual_[globI] += res; + jacobian_->addToBlock(globI, globI, bMat); } } @@ -479,6 +519,24 @@ private: FaceDir::DirEnum faceDirection; }; SparseTable neighborInfo_; + + using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState; + struct BoundaryConditionData + { + BCType type; + VectorBlock massRate; + unsigned pvtRegionIdx; + unsigned boundaryFaceIndex; + double faceArea; + double faceZCoord; + ScalarFluidState exFluidState; + }; + struct BoundaryInfo + { + unsigned int cell; + BoundaryConditionData bcdata; + }; + std::vector boundaryInfo_; }; } // namespace Opm diff --git a/opm/models/discretization/ecfv/ecfvstencil.hh b/opm/models/discretization/ecfv/ecfvstencil.hh index 343692436..425485f4d 100644 --- a/opm/models/discretization/ecfv/ecfvstencil.hh +++ b/opm/models/discretization/ecfv/ecfvstencil.hh @@ -212,7 +212,10 @@ public: { return area_; } /*! - * \brief Returns the direction of the face + * \brief Returns the direction id of the face w.r.t the cell. + * + * For corner point grids, this is 0-5 for I-, I+, J-, J+, K- and K+ faces, + * and -1 for NNC faces. */ int dirId() const {