From ec3b50d092f822d10befb60efabdd93cd0ff74e0 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 26 Apr 2016 13:20:53 +0200 Subject: [PATCH] adding updateWellState to MultisegmentWells --- opm/autodiff/BlackoilMultiSegmentModel.hpp | 1 + .../BlackoilMultiSegmentModel_impl.hpp | 64 +----------- opm/autodiff/MultisegmentWells.cpp | 2 + opm/autodiff/MultisegmentWells.hpp | 14 +++ opm/autodiff/MultisegmentWells_impl.hpp | 97 +++++++++++++++++++ 5 files changed, 116 insertions(+), 62 deletions(-) create mode 100644 opm/autodiff/MultisegmentWells_impl.hpp diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index 45fa399d8..86c037f76 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -172,6 +172,7 @@ namespace Opm { void updateWellControls(WellState& xw) const; + // TODO: kept for now. to be removed soon. void updateWellState(const V& dwells, WellState& well_state); diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index 9a6977c9a..af11a8c32 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -1116,66 +1116,7 @@ namespace Opm { BlackoilMultiSegmentModel::updateWellState(const V& dwells, WellState& well_state) { - - if (!wellsMultiSegment().empty()) - { - const int np = numPhases(); - const int nw = wellsMultiSegment().size(); - const int nseg_total = well_state.numSegments(); - - // Extract parts of dwells corresponding to each part. - int varstart = 0; - const V dsegqs = subset(dwells, Span(np * nseg_total, 1, varstart)); - varstart += dsegqs.size(); - const V dsegp = subset(dwells, Span(nseg_total, 1, varstart)); - varstart += dsegp.size(); - assert(varstart == dwells.size()); - const double dpmaxrel = dpMaxRel(); - - - // segment phase rates update - // in dwells, the phase rates are ordered by phase. - // while in WellStateMultiSegment, the phase rates are ordered by segments - const DataBlock wsr = Eigen::Map(dsegqs.data(), np, nseg_total).transpose(); - const V dwsr = Eigen::Map(wsr.data(), nseg_total * np); - const V wsr_old = Eigen::Map(&well_state.segPhaseRates()[0], nseg_total * np); - const V sr = wsr_old - dwsr; - std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin()); - - - // segment pressure updates - const V segp_old = Eigen::Map(&well_state.segPress()[0], nseg_total, 1); - // TODO: applying the pressure change limiter to all the segments, not sure if it is the correct thing to do - const V dsegp_limited = sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel); - const V segp = segp_old - dsegp_limited; - std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin()); - - // update the well rates and bhps, which are not anymore primary vabriables. - // they are updated directly from the updated segment phase rates and segment pressures. - - // Bhp update. - V bhp = V::Zero(nw); - V wr = V::Zero(nw * np); - // it is better to use subset - - int start_segment = 0; - for (int w = 0; w < nw; ++w) { - bhp[w] = well_state.segPress()[start_segment]; - // insert can be faster - for (int p = 0; p < np; ++p) { - wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment]; - } - - const int nseg = wellsMultiSegment()[w]->numberOfSegments(); - start_segment += nseg; - } - - assert(start_segment == nseg_total); - std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); - std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); - - // TODO: handling the THP control related. - } + msWells().updateWellState(dwells, fluid_.numPhases(), dpMaxRel(), well_state); } @@ -1487,8 +1428,7 @@ namespace Opm { ADB::V total_residual_v = total_residual.value(); const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); assert(dx.size() == total_residual_v.size()); - // asImpl().updateWellState(dx.array(), well_state); - updateWellState(dx.array(), well_state); + asImpl().updateWellState(dx.array(), well_state); updateWellControls(well_state); } } while (it < 15); diff --git a/opm/autodiff/MultisegmentWells.cpp b/opm/autodiff/MultisegmentWells.cpp index e7e2c094d..f30cefbc0 100644 --- a/opm/autodiff/MultisegmentWells.cpp +++ b/opm/autodiff/MultisegmentWells.cpp @@ -18,6 +18,8 @@ along with OPM. If not, see . */ +#include + #include namespace Opm { diff --git a/opm/autodiff/MultisegmentWells.hpp b/opm/autodiff/MultisegmentWells.hpp index 7c90722d0..4382fb1f8 100644 --- a/opm/autodiff/MultisegmentWells.hpp +++ b/opm/autodiff/MultisegmentWells.hpp @@ -30,9 +30,12 @@ #include #include +#include + #include + namespace Opm { @@ -116,6 +119,15 @@ namespace Opm { Vector& segVDt() { return segvdt_; }; + + template + void + updateWellState(const Vector& dwells, + const int np, + const double dpmaxrel, + WellState& well_state) const; + + protected: // TODO: probably a wells_active_ will be required here. const std::vector wells_multisegment_; @@ -175,4 +187,6 @@ namespace Opm { } // namespace Opm +#include "MultisegmentWells_impl.hpp" + #endif // OPM_MULTISEGMENTWELLS_HEADER_INCLUDED diff --git a/opm/autodiff/MultisegmentWells_impl.hpp b/opm/autodiff/MultisegmentWells_impl.hpp new file mode 100644 index 000000000..a72571200 --- /dev/null +++ b/opm/autodiff/MultisegmentWells_impl.hpp @@ -0,0 +1,97 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED +#define OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED + + + +namespace Opm +{ + + template + void + MultisegmentWells:: + updateWellState(const Vector& dwells, + const int np, + const double dpmaxrel, + WellState& well_state) const + { + if (!wells().empty()) + { + const int nw = wells().size(); + const int nseg_total = nseg_total_; + + // Extract parts of dwells corresponding to each part. + int varstart = 0; + const Vector dsegqs = subset(dwells, Span(np * nseg_total, 1, varstart)); + varstart += dsegqs.size(); + const Vector dsegp = subset(dwells, Span(nseg_total, 1, varstart)); + varstart += dsegp.size(); + assert(varstart == dwells.size()); + + + // segment phase rates update + // in dwells, the phase rates are ordered by phase. + // while in WellStateMultiSegment, the phase rates are ordered by segments + const DataBlock wsr = Eigen::Map(dsegqs.data(), np, nseg_total).transpose(); + const Vector dwsr = Eigen::Map(wsr.data(), nseg_total * np); + const Vector wsr_old = Eigen::Map(&well_state.segPhaseRates()[0], nseg_total * np); + const Vector sr = wsr_old - dwsr; + std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin()); + + + // segment pressure updates + const Vector segp_old = Eigen::Map(&well_state.segPress()[0], nseg_total, 1); + // TODO: applying the pressure change limiter to all the segments, not sure if it is the correct thing to do + const Vector dsegp_limited = sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel); + const Vector segp = segp_old - dsegp_limited; + std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin()); + + // update the well rates and bhps, which are not anymore primary vabriables. + // they are updated directly from the updated segment phase rates and segment pressures. + + // Bhp update. + Vector bhp = Vector::Zero(nw); + Vector wr = Vector::Zero(nw * np); + // it is better to use subset + + int start_segment = 0; + for (int w = 0; w < nw; ++w) { + bhp[w] = well_state.segPress()[start_segment]; + // insert can be faster + for (int p = 0; p < np; ++p) { + wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment]; + } + + const int nseg = wells()[w]->numberOfSegments(); + start_segment += nseg; + } + + assert(start_segment == nseg_total); + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + + // TODO: handling the THP control related. + } + } + +} +#endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED