Merge pull request #2809 from GitPaean/fixing_ct_aquifer

fixing the initialization of aquifer pressure for CT aquifer
This commit is contained in:
Atgeirr Flø Rasmussen 2020-09-28 12:28:02 +02:00 committed by GitHub
commit 64f9f1431b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 50 additions and 75 deletions

View File

@ -82,7 +82,7 @@ protected:
auto globalCellIdx = ugrid.globalCell();
// We hack the cell depth values for now. We can actually get it from elementcontext pos
this->cell_depth_.resize(this->size(), aquct_data_.d0);
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), 0.0);
@ -165,7 +165,7 @@ protected:
inline Scalar dpai(int idx)
{
Scalar dp = this->pa0_
+ this->rhow_.at(idx).value() * this->gravity_() * (this->cell_depth_.at(idx) - aquct_data_.d0)
+ this->rhow_.at(idx).value() * this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth())
- this->pressure_previous_.at(idx);
return dp;
}
@ -207,7 +207,7 @@ protected:
int pvttableIdx = aquct_data_.pvttableID - 1;
this->rhow_.resize(this->size(), 0.);
if (!aquct_data_.p0.first) {
this->pa0_ = calculateReservoirEquilibrium();
this->pa0_ = this->calculateReservoirEquilibrium();
} else {
this->pa0_ = aquct_data_.p0.second;
}
@ -232,42 +232,9 @@ protected:
mu_w_ = mu_w_aquifer.value();
}
// This function is for calculating the aquifer properties from equilibrium state with the reservoir
// TODO: this function can be moved to the Inteface class, since it is the same for both Aquifer models
inline Scalar calculateReservoirEquilibrium() override
virtual Scalar aquiferDepth() const override
{
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(this->ebos_simulator_);
const auto& gridView = this->ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
int idx = this->cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState();
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
this->rhow_[idx] = fs.density(waterPhaseIdx);
pw_aquifer.push_back(
(water_pressure_reservoir
- this->rhow_[idx].value() * this->gravity_() * (this->cell_depth_[idx] - aquct_data_.d0))
* this->alphai_[idx]);
}
// We take the average of the calculated equilibrium pressures.
Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / pw_aquifer.size();
return aquifer_pres_avg;
return aquct_data_.d0;
}
}; // class AquiferCarterTracy
} // namespace Opm

View File

@ -84,7 +84,7 @@ protected:
// We hack the cell depth values for now. We can actually get it from elementcontext pos
this->cell_depth_.resize(this->size(), aqufetp_data_.d0);
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), 0.0);
@ -171,7 +171,7 @@ protected:
inline Eval dpai(int idx)
{
const Eval dp = aquifer_pressure_ - this->pressure_current_.at(idx)
+ this->rhow_[idx] * this->gravity_() * (this->cell_depth_[idx] - aqufetp_data_.d0);
+ this->rhow_[idx] * this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
return dp;
}
@ -204,47 +204,16 @@ protected:
}
if (!aqufetp_data_.p0.first) {
this->pa0_ = calculateReservoirEquilibrium();
this->pa0_ = this->calculateReservoirEquilibrium();
} else {
this->pa0_ = aqufetp_data_.p0.second;
}
aquifer_pressure_ = this->pa0_;
}
inline Scalar calculateReservoirEquilibrium() override
virtual Scalar aquiferDepth() const override
{
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(this->ebos_simulator_);
const auto& gridView = this->ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
int idx = this->cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState();
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
this->rhow_[idx] = fs.density(waterPhaseIdx);
pw_aquifer.push_back(
(water_pressure_reservoir
- this->rhow_[idx].value() * this->gravity_() * (this->cell_depth_[idx] - aqufetp_data_.d0))
* this->alphai_[idx]);
}
// We take the average of the calculated equilibrium pressures.
const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.);
const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha;
return aquifer_pres_avg;
return aqufetp_data_.d0;
}
}; // Class AquiferFetkovich
} // namespace Opm

View File

@ -261,7 +261,46 @@ protected:
virtual void calculateAquiferConstants() = 0;
virtual Scalar calculateReservoirEquilibrium() = 0;
virtual Scalar aquiferDepth() const = 0;
// This function is for calculating the aquifer properties from equilibrium state with the reservoir
virtual Scalar calculateReservoirEquilibrium()
{
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
std::vector<Scalar> pw_aquifer;
Scalar water_pressure_reservoir;
ElementContext elemCtx(this->ebos_simulator_);
const auto& gridView = this->ebos_simulator_.gridView();
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
size_t cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
int idx = this->cellToConnectionIdx_[cellIdx];
if (idx < 0)
continue;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq0.fluidState();
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
this->rhow_[idx] = fs.density(waterPhaseIdx);
pw_aquifer.push_back(
(water_pressure_reservoir
- this->rhow_[idx].value() * this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth()))
* this->alphai_[idx]);
}
// We take the average of the calculated equilibrium pressures.
const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.);
const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha;
return aquifer_pres_avg;
}
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
};
} // namespace Opm