LookUp Indices in eclgenericproblem

This commit is contained in:
Antonella Ritorto 2023-10-12 16:50:57 +02:00
parent b60eb25285
commit 17c7e2c91f
4 changed files with 38 additions and 18 deletions

View File

@ -29,6 +29,7 @@
#include <dune/grid/common/gridview.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/LookUpData.hh>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>

View File

@ -50,6 +50,7 @@ namespace Opm {
class Deck;
class EclipseState;
class Schedule;
template<typename Grid, typename GridView> class LookUpData;
int eclPositionalParameter(Dune::ParameterTree& tree,
std::set<std::string>& seenParams,
@ -267,6 +268,9 @@ public:
int numPressurePointsEquil() const
{ return numPressurePointsEquil_; }
auto getLookUpData() const
{ return lookUpData_; }
bool operator==(const EclGenericProblem& rhs) const;
@ -370,7 +374,12 @@ protected:
// equilibration parameters
int numPressurePointsEquil_;
// To lookup origin cell indices
using Grid = std::remove_cv_t< typename std::remove_reference<decltype(gridView_.grid())>::type>;
using LookUpData = Opm::LookUpData<Grid,GridView>;
const LookUpData lookUpData_;
private:
template<class T>
void updateNum(const std::string& name, std::vector<T>& numbers, std::size_t num_regions);

View File

@ -89,6 +89,7 @@ EclGenericProblem(const EclipseState& eclState,
, schedule_(schedule)
, gridView_(gridView)
, mixControls_(schedule)
, lookUpData_(gridView)
{
}
@ -168,14 +169,15 @@ readRockParameters_(const std::vector<Scalar>& cellCenterDepths,
rockTableIdx_.resize(numElem);
const auto& num = eclState_.fieldProps().get_int(rock_config.rocknum_property());
for (std::size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) {
rockTableIdx_[elemIdx] = num[elemIdx] - 1;
rockTableIdx_[elemIdx] = num[this->lookUpData_.getOriginIndex(elemIdx)] - 1;
auto fmtError =
[&num,elemIdx,&ijkIndex,&rock_config](const char* type, std::size_t size)
[&num,elemIdx,&ijkIndex,&rock_config, this](const char* type, std::size_t size)
{
return fmt::format("{} table index {} for elem {} read from {}"
" is is out of bounds for number of tables {}",
type, num[elemIdx], ijkIndex(elemIdx),
rock_config.rocknum_property(), size);
type, num[this->lookUpData_.getOriginIndex(elemIdx)],
ijkIndex(this-> lookUpData_.getOriginIndex(elemIdx)),
rock_config.rocknum_property(), size);
};
if (!rockCompPoroMult_.empty() &&
rockTableIdx_[elemIdx] >= rockCompPoroMult_.size()) {
@ -210,7 +212,8 @@ readRockParameters_(const std::vector<Scalar>& cellCenterDepths,
if (!rockTableIdx_.empty()) {
tableIdx = rockTableIdx_[elemIdx];
}
overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true);
overburdenPressure_[elemIdx] =
overburdenTables[tableIdx].eval(cellCenterDepths[this->lookUpData_.getOriginIndex(elemIdx)], /*extrapolation=*/true);
}
}
}
@ -317,7 +320,7 @@ rockCompressibility(unsigned globalSpaceIdx) const
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty()) {
tableIdx = this->rockTableIdx_[globalSpaceIdx];
tableIdx = this->rockTableIdx_[globalSpaceIdx];
}
return this->rockParams_[tableIdx].compressibility;
}
@ -331,7 +334,7 @@ rockReferencePressure(unsigned globalSpaceIdx) const
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty()) {
tableIdx = this->rockTableIdx_[globalSpaceIdx];
tableIdx = this->rockTableIdx_[globalSpaceIdx];
}
return this->rockParams_[tableIdx].referencePressure;
}
@ -353,7 +356,7 @@ rockFraction(unsigned elementIdx, unsigned timeIdx) const
// geometric volume of the element. Note that it can
// be larger than 1.0 if porevolume multipliers are used
// to for instance implement larger boundary cells
Scalar porosity = poro[elementIdx];
Scalar porosity = poro[this->lookUpData_.getOriginIndex(elementIdx)];
return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
}
@ -370,10 +373,10 @@ updateNum(const std::string& name, std::vector<T>& numbers, std::size_t num_regi
unsigned numElems = gridView_.size(/*codim=*/0);
numbers.resize(numElems);
for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
if (numData[elemIdx] > (int)num_regions) {
if (numData[this->lookUpData_.getOriginIndex(elemIdx)] > (int)num_regions) {
throw std::runtime_error("Values larger than maximum number of regions " + std::to_string(num_regions) + " provided in " + name);
} else if (numData[elemIdx] > 0) {
numbers[elemIdx] = static_cast<T>(numData[elemIdx]) - 1;
} else if (numData[this->lookUpData_.getOriginIndex(elemIdx)] > 0) {
numbers[elemIdx] = static_cast<T>(numData[this->lookUpData_.getOriginIndex(elemIdx)]) - 1;
} else {
throw std::runtime_error("zero or negative values provided for region array: " + name);
}

View File

@ -202,6 +202,7 @@ public:
using EclGenericProblem<GridView,FluidSystem,Scalar>::rockCompressibility;
using EclGenericProblem<GridView,FluidSystem,Scalar>::rockReferencePressure;
using EclGenericProblem<GridView,FluidSystem,Scalar>::porosity;
using EclGenericProblem<GridView,FluidSystem,Scalar>::getLookUpData;
/*!
* \copydoc FvBaseProblem::registerParameters
@ -2091,13 +2092,19 @@ protected:
this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
if constexpr (enablePolymer)
this->polymer_.concentration[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);
this->polymer_.concentration[elemIdx] =
eclWriter_->eclOutputModule().getPolymerConcentration(getLookUpData().getOriginIndex(elemIdx));
if constexpr (enableMICP){
this->micp_.microbialConcentration[elemIdx] = eclWriter_->eclOutputModule().getMicrobialConcentration(elemIdx);
this->micp_.oxygenConcentration[elemIdx] = eclWriter_->eclOutputModule().getOxygenConcentration(elemIdx);
this->micp_.ureaConcentration[elemIdx] = eclWriter_->eclOutputModule().getUreaConcentration(elemIdx);
this->micp_.biofilmConcentration[elemIdx] = eclWriter_->eclOutputModule().getBiofilmConcentration(elemIdx);
this->micp_.calciteConcentration[elemIdx] = eclWriter_->eclOutputModule().getCalciteConcentration(elemIdx);
this->micp_.microbialConcentration[elemIdx] =
eclWriter_->eclOutputModule().getMicrobialConcentration(getLookUpData().getOriginIndex(elemIdx));
this->micp_.oxygenConcentration[elemIdx] =
eclWriter_->eclOutputModule().getOxygenConcentration(getLookUpData().getOriginIndex(elemIdx));
this->micp_.ureaConcentration[elemIdx] =
eclWriter_->eclOutputModule().getUreaConcentration(getLookUpData().getOriginIndex(elemIdx));
this->micp_.biofilmConcentration[elemIdx] =
eclWriter_->eclOutputModule().getBiofilmConcentration(getLookUpData().getOriginIndex(elemIdx));
this->micp_.calciteConcentration[elemIdx]
= eclWriter_->eclOutputModule().getCalciteConcentration(getLookUpData().getOriginIndex(elemIdx));
}
// if we need to restart for polymer molecular weight simulation, we need to add related here
}