diff --git a/opm/simulators/aquifers/AquiferAnalytical.hpp b/opm/simulators/aquifers/AquiferAnalytical.hpp index 5b785ef41..bc73a15ca 100644 --- a/opm/simulators/aquifers/AquiferAnalytical.hpp +++ b/opm/simulators/aquifers/AquiferAnalytical.hpp @@ -54,6 +54,7 @@ class AquiferAnalytical : public AquiferInterface { public: using Simulator = GetPropType; + using Scalar = GetPropType; using ElementContext = GetPropType; using FluidSystem = GetPropType; using BlackoilIndices = GetPropType; @@ -70,9 +71,8 @@ public: enum { enableSaltPrecipitation = getPropValue() }; static constexpr int numEq = BlackoilIndices::numEq; - using Scalar = double; - using Eval = DenseAd::Evaluation; + using Eval = DenseAd::Evaluation; using FluidState = BlackOilFluidState& total_face_area) override + void computeFaceAreaFraction(const std::vector& total_face_area) override { - assert (total_face_area.size() >= static_cast::size_type>(this->aquiferID())); + assert (total_face_area.size() >= static_cast::size_type>(this->aquiferID())); const auto tfa = total_face_area[this->aquiferID() - 1]; - const auto eps_sqrt = std::sqrt(std::numeric_limits::epsilon()); + const auto eps_sqrt = std::sqrt(std::numeric_limits::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_; } diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index df3338a57..2e27732e0 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -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)) { diff --git a/opm/simulators/aquifers/AquiferConstantFlux.hpp b/opm/simulators/aquifers/AquiferConstantFlux.hpp index 4d0662e35..4f57ba90e 100644 --- a/opm/simulators/aquifers/AquiferConstantFlux.hpp +++ b/opm/simulators/aquifers/AquiferConstantFlux.hpp @@ -43,9 +43,10 @@ public: using ElementMapper = GetPropType; using FluidSystem = GetPropType; using BlackoilIndices = GetPropType; + using Scalar = GetPropType; static constexpr int numEq = BlackoilIndices::numEq; - using Eval = DenseAd::Evaluation; + using Eval = DenseAd::Evaluation; AquiferConstantFlux(const std::vector& connections, const Simulator& simulator, @@ -68,15 +69,15 @@ public: virtual ~AquiferConstantFlux() = default; - void computeFaceAreaFraction(const std::vector& total_face_area) override + void computeFaceAreaFraction(const std::vector& total_face_area) override { - assert (total_face_area.size() >= static_cast::size_type>(this->aquiferID())); + assert (total_face_area.size() >= static_cast::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 connection_flux_{}; std::vector 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); }); diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 3681cfe09..92f18abd9 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -37,6 +37,7 @@ public: using FluidSystem = GetPropType; using RateVector = GetPropType; using Simulator = GetPropType; + using Scalar = GetPropType; // Constructor AquiferInterface(int aqID, @@ -58,8 +59,8 @@ public: virtual data::AquiferData aquiferData() const = 0; - virtual void computeFaceAreaFraction(const std::vector& total_face_area) = 0; - virtual double totalFaceArea() const = 0; + virtual void computeFaceAreaFraction(const std::vector& total_face_area) = 0; + virtual Scalar totalFaceArea() const = 0; template void addToSource(RateVector& rates, diff --git a/opm/simulators/aquifers/AquiferNumerical.hpp b/opm/simulators/aquifers/AquiferNumerical.hpp index 54341e0ac..f186c6f8c 100644 --- a/opm/simulators/aquifers/AquiferNumerical.hpp +++ b/opm/simulators/aquifers/AquiferNumerical.hpp @@ -50,12 +50,13 @@ public: using IntensiveQuantities = GetPropType; using MaterialLaw = GetPropType; using Simulator = GetPropType; + using Scalar = GetPropType; enum { dimWorld = GridView::dimensionworld }; enum { numPhases = FluidSystem::numPhases }; static constexpr int numEq = BlackoilIndices::numEq; - using Eval = DenseAd::Evaluation; + using Eval = DenseAd::Evaluation; using Toolbox = MathToolbox; using typename AquiferInterface::RateVector; @@ -111,7 +112,10 @@ public: if (const auto* aqData = xaqPos->second.typeData.template get(); 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(); - 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& /*total_face_area*/) override + void computeFaceAreaFraction(const std::vector& /*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(this->init_pressure_.size(), 0.0); + auto capture = std::vector(this->init_pressure_.size(), 0.0); return this->calculateAquiferPressure(capture); } - double calculateAquiferPressure(std::vector& cell_pressure) const + Scalar calculateAquiferPressure(std::vector& 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 - 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 init_pressure_{}; - double pressure_; // aquifer pressure + Scalar flux_rate_; // aquifer influx rate + Scalar cumulative_flux_; // cumulative aquifer influx + std::vector init_pressure_{}; + Scalar pressure_; // aquifer pressure bool solution_set_from_restart_ {false}; bool connects_to_reservoir_ {false}; diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index 1f36476c9..3d944ea32 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -362,7 +362,7 @@ void BlackoilAquiferModel::computeConnectionAreaFraction() const maxAquID = this->simulator_.vanguard().grid().comm().max(maxAquID); - auto totalConnArea = std::vector(maxAquID, 0.0); + auto totalConnArea = std::vector(maxAquID, 0.0); for (const auto& aquifer : this->aquifers) { totalConnArea[aquifer->aquiferID() - 1] += aquifer->totalFaceArea(); }