move segment_perforations_ and perforation_segment_depth_diffs_ to MultisegmentWellSegments

This commit is contained in:
Arne Morten Kvarving 2022-12-19 13:52:21 +01:00
parent aa684a10b8
commit 2766427df0
8 changed files with 153 additions and 141 deletions

View File

@ -50,7 +50,8 @@ template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
init(const int num_cells,
const int numPerfs,
const std::vector<int>& cells)
const std::vector<int>& cells,
const std::vector<std::vector<int>>& perforations)
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise);
@ -98,7 +99,7 @@ init(const int num_cells,
for (auto row = duneC_.createbegin(),
end = duneC_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : well_.segmentPerforations()[row.index()]) {
for (const int& perf : perforations[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}
@ -108,7 +109,7 @@ init(const int num_cells,
for (auto row = duneB_.createbegin(),
end = duneB_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : well_.segmentPerforations()[row.index()]) {
for (const int& perf : perforations[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}

View File

@ -71,9 +71,11 @@ public:
//! \param num_cells Total number of cells
//! \param numPerfs Number of perforations
//! \param cells Cell indices for perforations
//! \param perforations Cell indices for segment perforations
void init(const int num_cells,
const int numPerfs,
const std::vector<int>& cells);
const std::vector<int>& cells,
const std::vector<std::vector<int>>& perforations);
//! \brief Set all coefficients to 0.
void clear();

View File

@ -57,7 +57,7 @@ MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
, baseif_(baseif)
, linSys_(*this)
, primary_variables_(baseif)
, segments_(this->numberOfSegments(), baseif_.numComponents())
, segments_(this->numberOfSegments(), baseif)
, cell_perforation_depth_diffs_(baseif_.numPerfs(), 0.0)
, cell_perforation_pressure_diffs_(baseif_.numPerfs(), 0.0)
{
@ -68,7 +68,8 @@ void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
initMatrixAndVectors(const int num_cells)
{
linSys_.init(num_cells, baseif_.numPerfs(), baseif_.cells());
linSys_.init(num_cells, baseif_.numPerfs(),
baseif_.cells(), segments_.perforations_);
primary_variables_.resize(this->numberOfSegments());
}

View File

@ -46,42 +46,9 @@ template<typename Scalar>
MultisegmentWellGeneric<Scalar>::
MultisegmentWellGeneric(WellInterfaceGeneric& baseif)
: baseif_(baseif)
, segment_perforations_(numberOfSegments())
, segment_inlets_(numberOfSegments())
, segment_depth_diffs_(numberOfSegments(), 0.0)
, perforation_segment_depth_diffs_(baseif_.numPerfs(), 0.0)
{
// since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
// for other facilities needed but not available from parser, we need to process them here
// initialize the segment_perforations_ and update perforation_segment_depth_diffs_
const WellConnections& completion_set = baseif_.wellEcl().getConnections();
// index of the perforation within wells struct
// there might be some perforations not active, which causes the number of the perforations in
// well_ecl_ and wells struct different
// the current implementation is a temporary solution for now, it should be corrected from the parser
// side
int i_perf_wells = 0;
baseif.perfDepth().resize(baseif_.numPerfs(), 0.);
for (size_t perf = 0; perf < completion_set.size(); ++perf) {
const Connection& connection = completion_set.get(perf);
if (connection.state() == Connection::State::OPEN) {
const int segment_index = segmentNumberToIndex(connection.segment());
if ( segment_index == -1) {
OPM_THROW(std::logic_error,
fmt::format("COMPSEGS: Well {} has connection in cell {}, {}, {} "
"without associated segment.", baseif_.wellEcl().name(),
connection.getI() + 1, connection.getJ() + 1,
connection.getK() + 1));
}
segment_perforations_[segment_index].push_back(i_perf_wells);
baseif.perfDepth()[i_perf_wells] = connection.depth();
const double segment_depth = segmentSet()[segment_index].depth();
perforation_segment_depth_diffs_[i_perf_wells] = baseif.perfDepth()[i_perf_wells] - segment_depth;
i_perf_wells++;
}
}
// initialize the segment_inlets_
for (int seg = 0; seg < numberOfSegments(); ++seg) {
const Segment& segment = segmentSet()[seg];
@ -108,7 +75,8 @@ MultisegmentWellGeneric(WellInterfaceGeneric& baseif)
template<typename Scalar>
void
MultisegmentWellGeneric<Scalar>::
scaleSegmentRatesWithWellRates(WellState& well_state) const
scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_perforations,
WellState& well_state) const
{
auto& ws = well_state.well(baseif_.indexOfWell());
auto& segments = ws.segments;
@ -137,7 +105,10 @@ scaleSegmentRatesWithWellRates(WellState& well_state) const
}
std::vector<double> rates;
WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, num_single_phase, 0, rates);
WellState::calculateSegmentRates(segment_inlets_,
segment_perforations,
perforation_rates,
num_single_phase, 0, rates);
for (int seg = 0; seg < numberOfSegments(); ++seg) {
segment_rates[baseif_.numPhases() * seg + phase] = rates[seg];
}
@ -187,14 +158,6 @@ segmentInlets() const
return segment_inlets_;
}
template<typename Scalar>
const std::vector<std::vector<int>>&
MultisegmentWellGeneric<Scalar>::
segmentPerforations() const
{
return segment_perforations_;
}
template<typename Scalar>
int
MultisegmentWellGeneric<Scalar>::

View File

@ -44,9 +44,6 @@ public:
//! \brief Returns the inlet segments for each segment.
const std::vector<std::vector<int>>& segmentInlets() const;
//! \brief Returns the perforation index for each segment.
const std::vector<std::vector<int>>& segmentPerforations() const;
// get the WellSegments from the well_ecl_
const WellSegments& segmentSet() const;
@ -61,7 +58,8 @@ protected:
MultisegmentWellGeneric(WellInterfaceGeneric& baseif);
// scale the segment rates and pressure based on well rates and bhp
void scaleSegmentRatesWithWellRates(WellState& well_state) const;
void scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_perforations,
WellState& well_state) const;
void scaleSegmentPressuresWithBhp(WellState& well_state) const;
// components of the pressure drop to be included
@ -78,26 +76,10 @@ protected:
const WellInterfaceGeneric& baseif_;
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// the completions that is related to each segment
// the completions's ids are their index in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
std::vector<std::vector<int>> segment_perforations_;
// the inlet segments for each segment. It is for convenience and efficiency reason
std::vector<std::vector<int>> segment_inlets_;
std::vector<double> segment_depth_diffs_;
// depth difference between the segment and the perforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std::vector<double> perforation_segment_depth_diffs_;
};
}

View File

@ -21,6 +21,8 @@
#include <config.h>
#include <opm/simulators/wells/MultisegmentWellSegments.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
@ -28,21 +30,60 @@
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <fmt/format.h>
namespace Opm
{
template<class FluidSystem, class Indices, class Scalar>
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
MultisegmentWellSegments(const int numSegments,
const int numComponents)
: densities_(numSegments, 0.0)
WellInterfaceGeneric& well)
: perforations_(numSegments)
, perforation_depth_diffs_(well.numPerfs(), 0.0)
, densities_(numSegments, 0.0)
, mass_rates_(numSegments, 0.0)
, viscosities_(numSegments, 0.0)
, upwinding_segments_(numSegments, 0)
, phase_densities_(numSegments, std::vector<EvalWell>(numComponents, 0.0)) // number of phase here?
, phase_fractions_(numSegments, std::vector<EvalWell>(numComponents, 0.0)) // number of phase here?
, phase_viscosities_(numSegments, std::vector<EvalWell>(numComponents, 0.0)) // number of phase here?
{}
, phase_densities_(numSegments, std::vector<EvalWell>(well.numComponents(), 0.0)) // number of phase here?
, phase_fractions_(numSegments, std::vector<EvalWell>(well.numComponents(), 0.0)) // number of phase here?
, phase_viscosities_(numSegments, std::vector<EvalWell>(well.numComponents(), 0.0)) // number of phase here?
, well_(well)
{
// since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
// for other facilities needed but not available from parser, we need to process them here
// initialize the segment_perforations_ and update perforation_segment_depth_diffs_
const WellConnections& completion_set = well_.wellEcl().getConnections();
// index of the perforation within wells struct
// there might be some perforations not active, which causes the number of the perforations in
// well_ecl_ and wells struct different
// the current implementation is a temporary solution for now, it should be corrected from the parser
// side
int i_perf_wells = 0;
well.perfDepth().resize(well_.numPerfs(), 0.);
const auto& segment_set = well_.wellEcl().getSegments();
for (size_t perf = 0; perf < completion_set.size(); ++perf) {
const Connection& connection = completion_set.get(perf);
if (connection.state() == Connection::State::OPEN) {
const int segment_index = segment_set.segmentNumberToIndex(connection.segment());
if (segment_index == -1) {
OPM_THROW(std::logic_error,
fmt::format("COMPSEGS: Well {} has connection in cell {}, {}, {} "
"without associated segment.", well_.wellEcl().name(),
connection.getI() + 1, connection.getJ() + 1,
connection.getK() + 1));
}
perforations_[segment_index].push_back(i_perf_wells);
well.perfDepth()[i_perf_wells] = connection.depth();
const double segment_depth = segment_set[segment_index].depth();
perforation_depth_diffs_[i_perf_wells] = well_.perfDepth()[i_perf_wells] - segment_depth;
i_perf_wells++;
}
}
}
#define INSTANCE(...) \
template class MultisegmentWellSegments<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

View File

@ -29,6 +29,8 @@
namespace Opm
{
class WellInterfaceGeneric;
template<typename FluidSystem, typename Indices, typename Scalar>
class MultisegmentWellSegments
{
@ -37,7 +39,23 @@ class MultisegmentWellSegments
public:
MultisegmentWellSegments(const int numSegments,
const int numComponents);
WellInterfaceGeneric& well);
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// the completions that is related to each segment
// the completions's ids are their index in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
std::vector<std::vector<int>> perforations_;
// depth difference between the segment and the perforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std::vector<double> perforation_depth_diffs_;
// the densities of segment fluids
// we should not have this member variable
@ -55,6 +73,9 @@ public:
std::vector<std::vector<EvalWell>> phase_densities_;
std::vector<std::vector<EvalWell>> phase_fractions_;
std::vector<std::vector<EvalWell>> phase_viscosities_;
private:
const WellInterfaceGeneric& well_;
};
}

View File

@ -161,7 +161,7 @@ namespace Opm
Base::updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
// scale segment rates based on the wellRates
// and segment pressure based on bhp
this->scaleSegmentRatesWithWellRates(well_state);
this->scaleSegmentRatesWithWellRates(this->segments_.perforations_, well_state);
this->scaleSegmentPressuresWithBhp(well_state);
}
@ -368,7 +368,7 @@ namespace Opm
segments_copy.scale_pressure(bhp);
const auto& segment_pressure = segments_copy.pressure;
for (int seg = 0; seg < nseg; ++seg) {
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
// flux for each perforation
@ -441,7 +441,8 @@ namespace Opm
ws.surface_rates[phase] = sign * ws.well_potentials[phase];
}
}
well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
well_copy.scaleSegmentRatesWithWellRates(this->segments_.perforations_,
well_state_copy);
well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
const double dt = ebosSimulator.timeStepSize();
@ -770,7 +771,7 @@ namespace Opm
DeferredLogger& deferred_logger) const
{
// pressure difference between the segment and the perforation
const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf];
const Value perf_seg_press_diff = this->gravity() * segment_density * this->segments_.perforation_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
@ -1245,76 +1246,76 @@ namespace Opm
dp += seg_dp[outlet_segment_index];
}
seg_dp[seg] = dp;
for (const int perf : this->segment_perforations_[seg]) {
std::vector<Scalar> mob(this->num_components_, 0.0);
for (const int perf : this->segments_.perforations_[seg]) {
std::vector<Scalar> mob(this->num_components_, 0.0);
// TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobilityScalar(ebos_simulator, perf, mob);
// TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobilityScalar(ebos_simulator, perf, mob);
const int cell_idx = this->well_cells_[perf];
const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = int_quantities.fluidState();
// pressure difference between the segment and the perforation
const double perf_seg_press_diff = this->gravity_ * this->segments_.densities_[seg].value() * this->perforation_segment_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const double pressure_cell = this->getPerfCellPressure(fs).value();
const int cell_idx = this->well_cells_[perf];
const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = int_quantities.fluidState();
// pressure difference between the segment and the perforation
const double perf_seg_press_diff = this->gravity_ * this->segments_.densities_[seg].value() * this->segments_.perforation_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const double pressure_cell = this->getPerfCellPressure(fs).value();
// calculating the b for the connection
std::vector<double> b_perf(this->num_components_);
for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
// calculating the b for the connection
std::vector<double> b_perf(this->num_components_);
for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
}
const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
b_perf[comp_idx] = fs.invB(phase).value();
}
const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
b_perf[comp_idx] = fs.invB(phase).value();
}
// the pressure difference between the connection and BHP
const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
const double pressure_diff = pressure_cell - h_perf;
// the pressure difference between the connection and BHP
const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
const double pressure_diff = pressure_cell - h_perf;
// do not take into consideration the crossflow here.
if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
deferred_logger.debug("CROSSFLOW_IPR",
"cross flow found when updateIPR for well " + this->name());
}
// do not take into consideration the crossflow here.
if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
deferred_logger.debug("CROSSFLOW_IPR",
"cross flow found when updateIPR for well " + this->name());
}
// the well index associated with the connection
const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
// the well index associated with the connection
const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> ipr_b_perf(this->ipr_b_.size());
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
ipr_b_perf[comp_idx] += tw_mob;
}
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> ipr_b_perf(this->ipr_b_.size());
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
ipr_b_perf[comp_idx] += tw_mob;
}
// we need to handle the rs and rv when both oil and gas are present
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const double rs = (fs.Rs()).value();
const double rv = (fs.Rv()).value();
// we need to handle the rs and rv when both oil and gas are present
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const double rs = (fs.Rs()).value();
const double rv = (fs.Rv()).value();
const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
ipr_a_perf[gas_comp_idx] += dis_gas_a;
ipr_a_perf[oil_comp_idx] += vap_oil_a;
ipr_a_perf[gas_comp_idx] += dis_gas_a;
ipr_a_perf[oil_comp_idx] += vap_oil_a;
const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
ipr_b_perf[gas_comp_idx] += dis_gas_b;
ipr_b_perf[oil_comp_idx] += vap_oil_b;
}
ipr_b_perf[gas_comp_idx] += dis_gas_b;
ipr_b_perf[oil_comp_idx] += vap_oil_b;
}
for (size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
}
for (size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
}
}
}
}
@ -1605,7 +1606,7 @@ namespace Opm
auto& perf_data = ws.perf_data;
auto& perf_rates = perf_data.phase_rates;
auto& perf_press_state = perf_data.pressure;
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> mob(this->num_components_, 0.0);
@ -1687,14 +1688,14 @@ namespace Opm
for (int seg = 0; seg < nseg; ++seg) {
const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = intQuants.fluidState();
// pressure difference between the segment and the perforation
const EvalWell perf_seg_press_diff = this->gravity_ * this->segments_.densities_[seg] * this->perforation_segment_depth_diffs_[perf];
const EvalWell perf_seg_press_diff = this->gravity_ * this->segments_.densities_[seg] * this->segments_.perforation_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
@ -1888,7 +1889,7 @@ namespace Opm
double max_pressure = 0.0;
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = int_quants.fluidState();
@ -1916,7 +1917,7 @@ namespace Opm
for (int seg = 0; seg < nseg; ++seg) {
// calculating the perforation rate for each perforation that belongs to this segment
const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg));
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations_[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<Scalar> mob(this->num_components_, 0.0);