23 #ifndef OPM_WELLHELPERS_HEADER_INCLUDED
24 #define OPM_WELLHELPERS_HEADER_INCLUDED
26 #include <opm/common/OpmLog/OpmLog.hpp>
28 #include <opm/grid/cpgrid/GridHelpers.hpp>
30 #include <opm/simulators/wells/ParallelWellInfo.hpp>
32 #include <dune/istl/bcrsmatrix.hh>
33 #include <dune/common/dynmatrix.hh>
34 #include <dune/common/parallel/mpihelper.hh>
56 template<
typename Scalar>
60 using Block = Dune::DynamicMatrix<Scalar>;
61 using Matrix = Dune::BCRSMatrix<Block>;
64 : B_(&B), parallel_well_info_(¶llel_well_info)
68 template<
class X,
class Y>
69 void mv (
const X& x, Y& y)
const
71 #if !defined(NDEBUG) && HAVE_MPI
74 int cstring_size = parallel_well_info_->
name().size()+1;
75 std::vector<int> sizes(parallel_well_info_->communication().size());
76 parallel_well_info_->communication().allgather(&cstring_size, 1, sizes.data());
77 std::vector<int> offsets(sizes.size()+1, 0);
78 std::partial_sum(sizes.begin(), sizes.end(), offsets.begin() + 1);
79 std::vector<char> cstrings(offsets[sizes.size()]);
80 bool consistentWells =
true;
81 char* send =
const_cast<char*
>(parallel_well_info_->
name().c_str());
82 parallel_well_info_->communication().allgatherv(send, cstring_size,
83 cstrings.data(), sizes.data(),
85 for(std::size_t i = 0; i < sizes.size(); ++i)
87 std::string name(cstrings.data()+offsets[i]);
88 if (name != parallel_well_info_->
name())
90 if (parallel_well_info_->communication().rank() == 0)
93 std::string msg = std::string(
"Fatal Error: Not all ranks are computing for the same well")
94 +
" well should be " + parallel_well_info_->
name() +
" but is "
98 consistentWells =
false;
102 parallel_well_info_->communication().barrier();
104 if (!consistentWells)
106 MPI_Abort(MPI_COMM_WORLD, 1);
111 if (this->parallel_well_info_->communication().size() > 1)
118 using YField =
typename Y::block_type::value_type;
119 assert(y.size() == 1);
120 this->parallel_well_info_->communication().template allreduce<std::plus<YField>>(y[0].container().data(),
121 y[0].container().size());
126 template<
class X,
class Y>
127 void mmv (
const X& x, Y& y)
const
129 if (this->parallel_well_info_->communication().size() == 1)
149 double computeHydrostaticCorrection(
const double well_ref_depth,
const double vfp_ref_depth,
150 const double rho,
const double gravity) {
151 const double dh = vfp_ref_depth - well_ref_depth;
152 const double dp = rho * gravity * dh;
160 template<
typename Scalar,
typename Comm>
161 void sumDistributedWellEntries(Dune::DynamicMatrix<Scalar>& mat, Dune::DynamicVector<Scalar>& vec,
167 if (comm.size() == 1)
171 std::vector<Scalar> allEntries;
172 allEntries.reserve(mat.N()*mat.M()+vec.size());
173 for(
const auto& row: mat)
175 allEntries.insert(allEntries.end(), row.begin(), row.end());
177 allEntries.insert(allEntries.end(), vec.begin(), vec.end());
178 comm.sum(allEntries.data(), allEntries.size());
179 auto pos = allEntries.begin();
180 auto cols = mat.cols();
183 std::copy(pos, pos + cols, &(row[0]));
186 assert(std::size_t(allEntries.end() - pos) == vec.size());
187 std::copy(pos, allEntries.end(), &(vec[0]));
192 template <
int dim,
class C2F,
class FC>
193 std::array<double, dim>
194 getCubeDim(
const C2F& c2f,
195 FC begin_face_centroids,
198 std::array< std::vector<double>, dim > X;
200 const std::vector<double>::size_type
201 nf = std::distance(c2f[cell].begin(),
204 for (
int d = 0; d < dim; ++d) {
209 typedef typename C2F::row_type::const_iterator FI;
211 for (FI f = c2f[cell].begin(), e = c2f[cell].end(); f != e; ++f) {
212 using Opm::UgGridHelpers::increment;
213 using Opm::UgGridHelpers::getCoordinate;
215 const FC& fc = increment(begin_face_centroids, *f, dim);
217 for (
int d = 0; d < dim; ++d) {
218 X[d].push_back(getCoordinate(fc, d));
222 std::array<double, dim> cube;
223 for (
int d = 0; d < dim; ++d) {
224 typedef std::vector<double>::iterator VI;
225 typedef std::pair<VI,VI> PVI;
227 const PVI m = std::minmax_element(X[d].begin(), X[d].end());
229 cube[d] = *m.second - *m.first;
Class encapsulating some information about parallel wells.
Definition: ParallelWellInfo.hpp:252
const std::string & name() const
Name of the well.
Definition: ParallelWellInfo.hpp:347
A wrapper around the B matrix for distributed wells.
Definition: WellHelpers.hpp:58
void mv(const X &x, Y &y) const
y = A x
Definition: WellHelpers.hpp:69
void mmv(const X &x, Y &y) const
y = A x
Definition: WellHelpers.hpp:127
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27