OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5309 from akva2/aquifers_template_scalar

Browse Source 
Aquifer: use Scalar type
This commit is contained in:
Bård Skaflestad 2024-04-24 16:50:31 +02:00 committed by GitHub
commit 1dfdae3892
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
6 changed files with 59 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
opm/simulators/aquifers/AquiferAnalytical.hpp
View File

@ -54,6 +54,7 @@ class AquiferAnalytical : public AquiferInterface<TypeTag>
{
public:
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
@ -70,9 +71,8 @@ public:
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
static constexpr int numEq = BlackoilIndices::numEq;
using Scalar = double;
using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/numEq>;
using FluidState = BlackOilFluidState<Eval,
FluidSystem,
@ -99,12 +99,12 @@ public:
virtual ~AquiferAnalytical()
{}
void computeFaceAreaFraction(const std::vector<double>& total_face_area) override
void computeFaceAreaFraction(const std::vector<Scalar>& total_face_area) override
{
assert (total_face_area.size() >= static_cast<std::vector<double>::size_type>(this->aquiferID()));
assert (total_face_area.size() >= static_cast<typename std::vector<Scalar>::size_type>(this->aquiferID()));
const auto tfa = total_face_area[this->aquiferID() - 1];
const auto eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
const auto eps_sqrt = std::sqrt(std::numeric_limits<Scalar>::epsilon());
if (tfa < eps_sqrt) {
this->alphai_.assign(this->size(), Scalar{0});
@ -122,7 +122,7 @@ public:
this->area_fraction_ = this->totalFaceArea() / tfa;
}
double totalFaceArea() const override
Scalar totalFaceArea() const override
{
return this->total_face_area_;
}

4
opm/simulators/aquifers/AquiferCarterTracy.hpp
View File

@ -260,7 +260,7 @@ protected:
private:
Scalar timeConstantCO2Store() const
{
const auto press = this->aquct_data_.initial_pressure.value();
const Scalar press = this->aquct_data_.initial_pressure.value();
const auto temp = this->reservoirTemperatureCO2Store();
auto waterViscosity = Scalar { 0 };
@ -287,7 +287,7 @@ private:
Scalar waterDensityCO2Store() const
{
const auto press = this->aquct_data_.initial_pressure.value();
const Scalar press = this->aquct_data_.initial_pressure.value();
const auto temp = this->reservoirTemperatureCO2Store();
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {

25
opm/simulators/aquifers/AquiferConstantFlux.hpp
View File

@ -43,9 +43,10 @@ public:
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
static constexpr int numEq = BlackoilIndices::numEq;
using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/numEq>;
AquiferConstantFlux(const std::vector<Aquancon::AquancCell>& connections,
const Simulator& simulator,
@ -68,15 +69,15 @@ public:
virtual ~AquiferConstantFlux() = default;
void computeFaceAreaFraction(const std::vector<double>& total_face_area) override
void computeFaceAreaFraction(const std::vector<Scalar>& total_face_area) override
{
assert (total_face_area.size() >= static_cast<std::vector<double>::size_type>(this->aquiferID()));
assert (total_face_area.size() >= static_cast<typename std::vector<Scalar>::size_type>(this->aquiferID()));
this->area_fraction_ = this->totalFaceArea()
/ total_face_area[this->aquiferID() - 1];
}
double totalFaceArea() const override
Scalar totalFaceArea() const override
{
return this->total_face_area_;
}
@ -163,12 +164,12 @@ private:
SingleAquiferFlux aquifer_data_;
std::vector<Eval> connection_flux_{};
std::vector<int> cellToConnectionIdx_{};
double flux_rate_{};
double cumulative_flux_{};
double total_face_area_{0.0};
double area_fraction_{1.0};
Scalar flux_rate_{};
Scalar cumulative_flux_{};
Scalar total_face_area_{0.0};
Scalar area_fraction_{1.0};
double initializeConnections()
Scalar initializeConnections()
{
auto connected_face_area = 0.0;
@ -196,7 +197,7 @@ private:
return connected_face_area;
}
double computeFaceAreaFraction(const double connected_face_area) const
Scalar computeFaceAreaFraction(const Scalar connected_face_area) const
{
const auto tot_face_area = this->simulator_.vanguard()
.grid().comm().sum(connected_face_area);
@ -215,11 +216,11 @@ private:
return FluidSystem::waterCompIdx;
}
double totalFluxRate() const
Scalar totalFluxRate() const
{
return std::accumulate(this->connection_flux_.begin(),
this->connection_flux_.end(), 0.0,
[](const double rate, const auto& q)
[](const Scalar rate, const auto& q)
{
return rate + getValue(q);
});

5
opm/simulators/aquifers/AquiferInterface.hpp
View File

@ -37,6 +37,7 @@ public:
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
// Constructor
AquiferInterface(int aqID,
@ -58,8 +59,8 @@ public:
virtual data::AquiferData aquiferData() const = 0;
virtual void computeFaceAreaFraction(const std::vector<double>& total_face_area) = 0;
virtual double totalFaceArea() const = 0;
virtual void computeFaceAreaFraction(const std::vector<Scalar>& total_face_area) = 0;
virtual Scalar totalFaceArea() const = 0;
template <class Context>
void addToSource(RateVector& rates,

61
opm/simulators/aquifers/AquiferNumerical.hpp
View File

@ -50,12 +50,13 @@ public:
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
enum { dimWorld = GridView::dimensionworld };
enum { numPhases = FluidSystem::numPhases };
static constexpr int numEq = BlackoilIndices::numEq;
using Eval = DenseAd::Evaluation<double, numEq>;
using Eval = DenseAd::Evaluation<Scalar, numEq>;
using Toolbox = MathToolbox<Eval>;
using typename AquiferInterface<TypeTag>::RateVector;
@ -111,7 +112,10 @@ public:
if (const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
aqData != nullptr)
{
this->init_pressure_ = aqData->initPressure;
this->init_pressure_.resize(aqData->initPressure.size());
std::copy(aqData->initPressure.begin(),
aqData->initPressure.end(),
this->init_pressure_.begin());
}
this->solution_set_from_restart_ = true;
@ -136,7 +140,10 @@ public:
data.volume = this->cumulative_flux_;
auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
aquNum->initPressure = this->init_pressure_;
aquNum->initPressure.resize(this->init_pressure_.size());
std::copy(this->init_pressure_.begin(),
this->init_pressure_.end(),
aquNum->initPressure.begin());
return data;
}
@ -152,10 +159,10 @@ public:
this->cumulative_flux_ = 0.;
}
void computeFaceAreaFraction(const std::vector<double>& /*total_face_area*/) override
void computeFaceAreaFraction(const std::vector<Scalar>& /*total_face_area*/) override
{}
double totalFaceArea() const override
Scalar totalFaceArea() const override
{
return 1.0;
}
@ -177,7 +184,7 @@ public:
this->pressure_ == rhs.pressure_;
}
double cumulativeFlux() const
Scalar cumulativeFlux() const
{
return this->cumulative_flux_;
}
@ -205,16 +212,16 @@ private:
elemIt->partitionType() == Dune::InteriorEntity;
}
double calculateAquiferPressure() const
Scalar calculateAquiferPressure() const
{
auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
auto capture = std::vector<Scalar>(this->init_pressure_.size(), 0.0);
return this->calculateAquiferPressure(capture);
}
double calculateAquiferPressure(std::vector<double>& cell_pressure) const
Scalar calculateAquiferPressure(std::vector<Scalar>& cell_pressure) const
{
double sum_pressure_watervolume = 0.;
double sum_watervolume = 0.;
Scalar sum_pressure_watervolume = 0.;
Scalar sum_watervolume = 0.;
ElementContext elem_ctx(this->simulator_);
const auto& gridView = this->simulator_.gridView();
@ -236,12 +243,12 @@ private:
// TODO: the porosity of the cells are still wrong for numerical aquifer cells
// Because the dofVolume still based on the grid information.
// The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
const double water_saturation = fs.saturation(this->phaseIdx_()).value();
const double porosity = iq0.porosity().value();
const double volume = elem_ctx.dofTotalVolume(0, 0);
const Scalar water_saturation = fs.saturation(this->phaseIdx_()).value();
const Scalar porosity = iq0.porosity().value();
const Scalar volume = elem_ctx.dofTotalVolume(0, 0);
// TODO: not sure we should use water pressure here
const double water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
const double water_volume = volume * porosity * water_saturation;
const Scalar water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
const Scalar water_volume = volume * porosity * water_saturation;
sum_pressure_watervolume += water_volume * water_pressure_reservoir;
sum_watervolume += water_volume;
@ -260,16 +267,16 @@ private:
}
template <class ElemCtx>
double getWaterFlux(const ElemCtx& elem_ctx, unsigned face_idx) const
Scalar getWaterFlux(const ElemCtx& elem_ctx, unsigned face_idx) const
{
const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
const Scalar water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
return water_flux;
}
double calculateAquiferFluxRate() const
Scalar calculateAquiferFluxRate() const
{
double aquifer_flux = 0.0;
Scalar aquifer_flux = 0.0;
if (! this->connects_to_reservoir_) {
return aquifer_flux;
@ -312,11 +319,11 @@ private:
elem_ctx.updateAllIntensiveQuantities();
elem_ctx.updateAllExtensiveQuantities();
const double water_flux = getWaterFlux(elem_ctx,face_idx);
const Scalar water_flux = getWaterFlux(elem_ctx,face_idx);
const std::size_t up_id = water_flux >= 0.0 ? i : j;
const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
const double invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
const double face_area = face.area();
const Scalar invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
const Scalar face_area = face.area();
aquifer_flux += water_flux * invB * face_area;
}
@ -327,10 +334,10 @@ private:
return aquifer_flux;
}
double flux_rate_; // aquifer influx rate
double cumulative_flux_; // cumulative aquifer influx
std::vector<double> init_pressure_{};
double pressure_; // aquifer pressure
Scalar flux_rate_; // aquifer influx rate
Scalar cumulative_flux_; // cumulative aquifer influx
std::vector<Scalar> init_pressure_{};
Scalar pressure_; // aquifer pressure
bool solution_set_from_restart_ {false};
bool connects_to_reservoir_ {false};

2
opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
View File

@ -362,7 +362,7 @@ void BlackoilAquiferModel<TypeTag>::computeConnectionAreaFraction() const
maxAquID = this->simulator_.vanguard().grid().comm().max(maxAquID);
auto totalConnArea = std::vector<double>(maxAquID, 0.0);
auto totalConnArea = std::vector<Scalar>(maxAquID, 0.0);
for (const auto& aquifer : this->aquifers) {
totalConnArea[aquifer->aquiferID() - 1] += aquifer->totalFaceArea();
}