mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
wellhelpers: make templates private
to increase header encapsulation
This commit is contained in:
Arne Morten Kvarving
parent
9611edfee4
commit
69db0dd175
@ -96,6 +96,7 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
opm/simulators/wells/VFPProdProperties.cpp
|
||||
opm/simulators/wells/VFPInjProperties.cpp
|
||||
opm/simulators/wells/WellGroupHelpers.cpp
|
||||
opm/simulators/wells/WellHelpers.cpp
|
||||
opm/simulators/wells/WellInterfaceEval.cpp
|
||||
opm/simulators/wells/WellInterfaceFluidSystem.cpp
|
||||
opm/simulators/wells/WellInterfaceGeneric.cpp
|
||||
|
||||
240
opm/simulators/wells/WellHelpers.cpp
Normal file
240
opm/simulators/wells/WellHelpers.cpp
Normal file
@ -0,0 +1,240 @@
|
||||
/*
|
||||
Copyright 2016 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2016 Statoil ASA.
|
||||
Copyright 2020 OPM-OP AS.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
|
||||
#include <opm/simulators/wells/WellHelpers.hpp>
|
||||
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
|
||||
#include <opm/grid/cpgrid/GridHelpers.hpp>
|
||||
|
||||
#include <opm/simulators/wells/ParallelWellInfo.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
namespace wellhelpers {
|
||||
|
||||
template<typename Scalar>
|
||||
ParallelStandardWellB<Scalar>::
|
||||
ParallelStandardWellB(const Matrix& B, const ParallelWellInfo& parallel_well_info)
|
||||
: B_(&B), parallel_well_info_(¶llel_well_info)
|
||||
{}
|
||||
|
||||
template<typename Scalar>
|
||||
template<class X, class Y>
|
||||
void ParallelStandardWellB<Scalar>::
|
||||
mv (const X& x, Y& y) const
|
||||
{
|
||||
#if !defined(NDEBUG) && HAVE_MPI
|
||||
// We need to make sure that all ranks are actually computing
|
||||
// for the same well. Doing this by checking the name of the well.
|
||||
int cstring_size = parallel_well_info_->name().size()+1;
|
||||
std::vector<int> sizes(parallel_well_info_->communication().size());
|
||||
parallel_well_info_->communication().allgather(&cstring_size, 1, sizes.data());
|
||||
std::vector<int> offsets(sizes.size()+1, 0); //last entry will be accumulated size
|
||||
std::partial_sum(sizes.begin(), sizes.end(), offsets.begin() + 1);
|
||||
std::vector<char> cstrings(offsets[sizes.size()]);
|
||||
bool consistentWells = true;
|
||||
char* send = const_cast<char*>(parallel_well_info_->name().c_str());
|
||||
parallel_well_info_->communication().allgatherv(send, cstring_size,
|
||||
cstrings.data(), sizes.data(),
|
||||
offsets.data());
|
||||
for(std::size_t i = 0; i < sizes.size(); ++i)
|
||||
{
|
||||
std::string name(cstrings.data()+offsets[i]);
|
||||
if (name != parallel_well_info_->name())
|
||||
{
|
||||
if (parallel_well_info_->communication().rank() == 0)
|
||||
{
|
||||
//only one process per well logs, might not be 0 of MPI_COMM_WORLD, though
|
||||
std::string msg = std::string("Fatal Error: Not all ranks are computing for the same well")
|
||||
+ " well should be " + parallel_well_info_->name() + " but is "
|
||||
+ name;
|
||||
OpmLog::debug(msg);
|
||||
}
|
||||
consistentWells = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
parallel_well_info_->communication().barrier();
|
||||
// As not all processes are involved here we need to use MPI_Abort and hope MPI kills them all
|
||||
if (!consistentWells)
|
||||
{
|
||||
MPI_Abort(MPI_COMM_WORLD, 1);
|
||||
}
|
||||
#endif
|
||||
B_->mv(x, y);
|
||||
|
||||
if (this->parallel_well_info_->communication().size() > 1)
|
||||
{
|
||||
// Only do communication if we must.
|
||||
// The B matrix is basically a component-wise multiplication
|
||||
// with a vector followed by a parallel reduction. We do that
|
||||
// reduction to all ranks computing for the well to save the
|
||||
// broadcast when applying C^T.
|
||||
using YField = typename Y::block_type::value_type;
|
||||
assert(y.size() == 1);
|
||||
this->parallel_well_info_->communication().template allreduce<std::plus<YField>>(y[0].container().data(),
|
||||
y[0].container().size());
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Scalar>
|
||||
template<class X, class Y>
|
||||
void ParallelStandardWellB<Scalar>::
|
||||
mmv (const X& x, Y& y) const
|
||||
{
|
||||
if (this->parallel_well_info_->communication().size() == 1)
|
||||
{
|
||||
// Do the same thing as before. The else branch
|
||||
// produces different rounding errors and results
|
||||
// slightly different iteration counts / well curves
|
||||
B_->mmv(x, y);
|
||||
}
|
||||
else
|
||||
{
|
||||
Y temp(y);
|
||||
mv(x, temp); // includes parallel reduction
|
||||
y -= temp;
|
||||
}
|
||||
}
|
||||
|
||||
double computeHydrostaticCorrection(const double well_ref_depth, const double vfp_ref_depth,
|
||||
const double rho, const double gravity)
|
||||
{
|
||||
const double dh = vfp_ref_depth - well_ref_depth;
|
||||
const double dp = rho * gravity * dh;
|
||||
|
||||
return dp;
|
||||
}
|
||||
|
||||
template<typename Scalar, typename Comm>
|
||||
void sumDistributedWellEntries(Dune::DynamicMatrix<Scalar>& mat,
|
||||
Dune::DynamicVector<Scalar>& vec,
|
||||
const Comm& comm)
|
||||
{
|
||||
// DynamicMatrix does not use one contiguous array for storing the data
|
||||
// but a DynamicVector of DynamicVectors. Hence we need to copy the data
|
||||
// to contiguous memory for MPI.
|
||||
if (comm.size() == 1)
|
||||
{
|
||||
return;
|
||||
}
|
||||
std::vector<Scalar> allEntries;
|
||||
allEntries.reserve(mat.N()*mat.M()+vec.size());
|
||||
for(const auto& row: mat)
|
||||
{
|
||||
allEntries.insert(allEntries.end(), row.begin(), row.end());
|
||||
}
|
||||
allEntries.insert(allEntries.end(), vec.begin(), vec.end());
|
||||
comm.sum(allEntries.data(), allEntries.size());
|
||||
auto pos = allEntries.begin();
|
||||
auto cols = mat.cols();
|
||||
for(auto&& row: mat)
|
||||
{
|
||||
std::copy(pos, pos + cols, &(row[0]));
|
||||
pos += cols;
|
||||
}
|
||||
assert(std::size_t(allEntries.end() - pos) == vec.size());
|
||||
std::copy(pos, allEntries.end(), &(vec[0]));
|
||||
}
|
||||
|
||||
template <int dim, class C2F, class FC>
|
||||
std::array<double, dim>
|
||||
getCubeDim(const C2F& c2f,
|
||||
FC begin_face_centroids,
|
||||
int cell)
|
||||
{
|
||||
std::array<std::vector<double>, dim> X;
|
||||
{
|
||||
const std::vector<double>::size_type
|
||||
nf = std::distance(c2f[cell].begin(),
|
||||
c2f[cell].end ());
|
||||
|
||||
for (int d = 0; d < dim; ++d) {
|
||||
X[d].reserve(nf);
|
||||
}
|
||||
}
|
||||
|
||||
for (const auto& f : c2f[cell]) {
|
||||
using Opm::UgGridHelpers::increment;
|
||||
using Opm::UgGridHelpers::getCoordinate;
|
||||
|
||||
const FC& fc = increment(begin_face_centroids, f, dim);
|
||||
|
||||
for (int d = 0; d < dim; ++d) {
|
||||
X[d].push_back(getCoordinate(fc, d));
|
||||
}
|
||||
}
|
||||
|
||||
std::array<double, dim> cube;
|
||||
for (int d = 0; d < dim; ++d) {
|
||||
using VI = std::vector<double>::iterator;
|
||||
using PVI = std::pair<VI,VI>;
|
||||
|
||||
const PVI m = std::minmax_element(X[d].begin(), X[d].end());
|
||||
|
||||
cube[d] = *m.second - *m.first;
|
||||
}
|
||||
|
||||
return cube;
|
||||
}
|
||||
|
||||
template <class DenseMatrix>
|
||||
DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M)
|
||||
{
|
||||
DenseMatrix tmp{M.cols(), M.rows()};
|
||||
for (size_t i = 0; i < M.rows(); ++i) {
|
||||
for (size_t j = 0; j < M.cols(); ++j) {
|
||||
tmp[j][i] = M[i][j];
|
||||
}
|
||||
}
|
||||
return tmp;
|
||||
}
|
||||
|
||||
template class ParallelStandardWellB<double>;
|
||||
|
||||
template<int Dim> using Vec = Dune::BlockVector<Dune::FieldVector<double,Dim>>;
|
||||
using DynVec = Dune::BlockVector<Dune::DynamicVector<double>>;
|
||||
|
||||
#define INSTANCE(Dim) \
|
||||
template void ParallelStandardWellB<double>::mv<Vec<Dim>,DynVec>(const Vec<Dim>&,DynVec&) const; \
|
||||
template void ParallelStandardWellB<double>::mmv<Vec<Dim>,DynVec>(const Vec<Dim>&,DynVec&) const;
|
||||
|
||||
INSTANCE(1)
|
||||
INSTANCE(2)
|
||||
INSTANCE(3)
|
||||
INSTANCE(4)
|
||||
INSTANCE(5)
|
||||
INSTANCE(6)
|
||||
|
||||
using Comm = Parallel::Communication;
|
||||
template void sumDistributedWellEntries<double,Comm>(Dune::DynamicMatrix<double>& mat,
|
||||
Dune::DynamicVector<double>& vec,
|
||||
const Comm& comm);
|
||||
|
||||
using DMatrix = Dune::DynamicMatrix<double>;
|
||||
template DMatrix transposeDenseDynMatrix<DMatrix>(const DMatrix&);
|
||||
|
||||
} // namespace wellhelpers
|
||||
} // namespace Opm
|
||||
@ -23,22 +23,14 @@
|
||||
#ifndef OPM_WELLHELPERS_HEADER_INCLUDED
|
||||
#define OPM_WELLHELPERS_HEADER_INCLUDED
|
||||
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
|
||||
#include <opm/grid/cpgrid/GridHelpers.hpp>
|
||||
|
||||
#include <opm/simulators/wells/ParallelWellInfo.hpp>
|
||||
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/common/dynmatrix.hh>
|
||||
#include <dune/common/parallel/mpihelper.hh>
|
||||
|
||||
#include <vector>
|
||||
#include <array>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
|
||||
|
||||
class ParallelWellInfo;
|
||||
|
||||
namespace wellhelpers
|
||||
{
|
||||
@ -60,132 +52,30 @@ namespace Opm {
|
||||
using Block = Dune::DynamicMatrix<Scalar>;
|
||||
using Matrix = Dune::BCRSMatrix<Block>;
|
||||
|
||||
ParallelStandardWellB(const Matrix& B, const ParallelWellInfo& parallel_well_info)
|
||||
: B_(&B), parallel_well_info_(¶llel_well_info)
|
||||
{}
|
||||
ParallelStandardWellB(const Matrix& B, const ParallelWellInfo& parallel_well_info);
|
||||
|
||||
//! y = A x
|
||||
template<class X, class Y>
|
||||
void mv (const X& x, Y& y) const
|
||||
{
|
||||
#if !defined(NDEBUG) && HAVE_MPI
|
||||
// We need to make sure that all ranks are actually computing
|
||||
// for the same well. Doing this by checking the name of the well.
|
||||
int cstring_size = parallel_well_info_->name().size()+1;
|
||||
std::vector<int> sizes(parallel_well_info_->communication().size());
|
||||
parallel_well_info_->communication().allgather(&cstring_size, 1, sizes.data());
|
||||
std::vector<int> offsets(sizes.size()+1, 0); //last entry will be accumulated size
|
||||
std::partial_sum(sizes.begin(), sizes.end(), offsets.begin() + 1);
|
||||
std::vector<char> cstrings(offsets[sizes.size()]);
|
||||
bool consistentWells = true;
|
||||
char* send = const_cast<char*>(parallel_well_info_->name().c_str());
|
||||
parallel_well_info_->communication().allgatherv(send, cstring_size,
|
||||
cstrings.data(), sizes.data(),
|
||||
offsets.data());
|
||||
for(std::size_t i = 0; i < sizes.size(); ++i)
|
||||
{
|
||||
std::string name(cstrings.data()+offsets[i]);
|
||||
if (name != parallel_well_info_->name())
|
||||
{
|
||||
if (parallel_well_info_->communication().rank() == 0)
|
||||
{
|
||||
//only one process per well logs, might not be 0 of MPI_COMM_WORLD, though
|
||||
std::string msg = std::string("Fatal Error: Not all ranks are computing for the same well")
|
||||
+ " well should be " + parallel_well_info_->name() + " but is "
|
||||
+ name;
|
||||
OpmLog::debug(msg);
|
||||
}
|
||||
consistentWells = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
parallel_well_info_->communication().barrier();
|
||||
// As not all processes are involved here we need to use MPI_Abort and hope MPI kills them all
|
||||
if (!consistentWells)
|
||||
{
|
||||
MPI_Abort(MPI_COMM_WORLD, 1);
|
||||
}
|
||||
#endif
|
||||
B_->mv(x, y);
|
||||
|
||||
if (this->parallel_well_info_->communication().size() > 1)
|
||||
{
|
||||
// Only do communication if we must.
|
||||
// The B matrix is basically a component-wise multiplication
|
||||
// with a vector followed by a parallel reduction. We do that
|
||||
// reduction to all ranks computing for the well to save the
|
||||
// broadcast when applying C^T.
|
||||
using YField = typename Y::block_type::value_type;
|
||||
assert(y.size() == 1);
|
||||
this->parallel_well_info_->communication().template allreduce<std::plus<YField>>(y[0].container().data(),
|
||||
y[0].container().size());
|
||||
}
|
||||
}
|
||||
void mv (const X& x, Y& y) const;
|
||||
|
||||
//! y = A x
|
||||
template<class X, class Y>
|
||||
void mmv (const X& x, Y& y) const
|
||||
{
|
||||
if (this->parallel_well_info_->communication().size() == 1)
|
||||
{
|
||||
// Do the same thing as before. The else branch
|
||||
// produces different rounding errors and results
|
||||
// slightly different iteration counts / well curves
|
||||
B_->mmv(x, y);
|
||||
}
|
||||
else
|
||||
{
|
||||
Y temp(y);
|
||||
mv(x, temp); // includes parallel reduction
|
||||
y -= temp;
|
||||
}
|
||||
}
|
||||
void mmv (const X& x, Y& y) const;
|
||||
|
||||
private:
|
||||
const Matrix* B_;
|
||||
const ParallelWellInfo* parallel_well_info_;
|
||||
};
|
||||
|
||||
inline
|
||||
double computeHydrostaticCorrection(const double well_ref_depth, const double vfp_ref_depth,
|
||||
const double rho, const double gravity) {
|
||||
const double dh = vfp_ref_depth - well_ref_depth;
|
||||
const double dp = rho * gravity * dh;
|
||||
|
||||
return dp;
|
||||
}
|
||||
|
||||
double computeHydrostaticCorrection(const double well_ref_depth,
|
||||
const double vfp_ref_depth,
|
||||
const double rho, const double gravity);
|
||||
|
||||
|
||||
/// \brief Sums entries of the diagonal Matrix for distributed wells
|
||||
template<typename Scalar, typename Comm>
|
||||
void sumDistributedWellEntries(Dune::DynamicMatrix<Scalar>& mat, Dune::DynamicVector<Scalar>& vec,
|
||||
const Comm& comm)
|
||||
{
|
||||
// DynamicMatrix does not use one contiguous array for storing the data
|
||||
// but a DynamicVector of DynamicVectors. Hence we need to copy the data
|
||||
// to contiguous memory for MPI.
|
||||
if (comm.size() == 1)
|
||||
{
|
||||
return;
|
||||
}
|
||||
std::vector<Scalar> allEntries;
|
||||
allEntries.reserve(mat.N()*mat.M()+vec.size());
|
||||
for(const auto& row: mat)
|
||||
{
|
||||
allEntries.insert(allEntries.end(), row.begin(), row.end());
|
||||
}
|
||||
allEntries.insert(allEntries.end(), vec.begin(), vec.end());
|
||||
comm.sum(allEntries.data(), allEntries.size());
|
||||
auto pos = allEntries.begin();
|
||||
auto cols = mat.cols();
|
||||
for(auto&& row: mat)
|
||||
{
|
||||
std::copy(pos, pos + cols, &(row[0]));
|
||||
pos += cols;
|
||||
}
|
||||
assert(std::size_t(allEntries.end() - pos) == vec.size());
|
||||
std::copy(pos, allEntries.end(), &(vec[0]));
|
||||
}
|
||||
const Comm& comm);
|
||||
|
||||
|
||||
|
||||
@ -193,58 +83,12 @@ namespace Opm {
|
||||
std::array<double, dim>
|
||||
getCubeDim(const C2F& c2f,
|
||||
FC begin_face_centroids,
|
||||
int cell)
|
||||
{
|
||||
std::array< std::vector<double>, dim > X;
|
||||
{
|
||||
const std::vector<double>::size_type
|
||||
nf = std::distance(c2f[cell].begin(),
|
||||
c2f[cell].end ());
|
||||
int cell);
|
||||
|
||||
for (int d = 0; d < dim; ++d) {
|
||||
X[d].reserve(nf);
|
||||
}
|
||||
}
|
||||
|
||||
typedef typename C2F::row_type::const_iterator FI;
|
||||
|
||||
for (FI f = c2f[cell].begin(), e = c2f[cell].end(); f != e; ++f) {
|
||||
using Opm::UgGridHelpers::increment;
|
||||
using Opm::UgGridHelpers::getCoordinate;
|
||||
|
||||
const FC& fc = increment(begin_face_centroids, *f, dim);
|
||||
|
||||
for (int d = 0; d < dim; ++d) {
|
||||
X[d].push_back(getCoordinate(fc, d));
|
||||
}
|
||||
}
|
||||
|
||||
std::array<double, dim> cube;
|
||||
for (int d = 0; d < dim; ++d) {
|
||||
typedef std::vector<double>::iterator VI;
|
||||
typedef std::pair<VI,VI> PVI;
|
||||
|
||||
const PVI m = std::minmax_element(X[d].begin(), X[d].end());
|
||||
|
||||
cube[d] = *m.second - *m.first;
|
||||
}
|
||||
|
||||
return cube;
|
||||
}
|
||||
|
||||
// explicite transpose of dense matrix due to compilation problems
|
||||
// used for caclulating quasiimpes well weights
|
||||
// explicit transpose of a dense matrix due to compilation problems
|
||||
// used for calculating quasiimpes well weights
|
||||
template <class DenseMatrix>
|
||||
DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M)
|
||||
{
|
||||
DenseMatrix tmp{M.cols(), M.rows()};
|
||||
for (size_t i = 0; i < M.rows(); ++i) {
|
||||
for (size_t j = 0; j < M.cols(); ++j) {
|
||||
tmp[j][i] = M[i][j];
|
||||
}
|
||||
}
|
||||
return tmp;
|
||||
}
|
||||
DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M);
|
||||
|
||||
} // namespace wellhelpers
|
||||
}
|
||||
|
||||
@ -30,6 +30,7 @@
|
||||
|
||||
#include <opm/grid/UnstructuredGrid.h>
|
||||
#include <opm/grid/GridManager.hpp>
|
||||
#include <opm/grid/cpgrid/GridHelpers.hpp>
|
||||
|
||||
#include <opm/input/eclipse/Units/Units.hpp>
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user