First steps

opm/simulators/wells/MultisegmentWellEquations.cpp

@@ -55,9 +55,10 @@ MultisegmentWellEquations(const MultisegmentWellGeneric<Scalar>& well)
 
 template<class Scalar, int numWellEq, int numEq>
 void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
-init(const int num_cells,
- const int numPerfs,
- const std::vector<int>& cells,
+init(const int num_cells_this_process, //num_cells on this process
+ [[maybe_unused]] const int num_perfs_this_process, //num_perfs_this_process on this process
+ const int num_perfs_whole_mswell,
+ const std::vector<int>& cells, // cells for the perforations of all processes
  const std::vector<std::vector<int>>& segment_inlets,
  const std::vector<std::vector<int>>& perforations)
 {
@@ -78,8 +79,8 @@ init(const int num_cells,
  }
  duneD_.setSize(well_.numberOfSegments(), well_.numberOfSegments(), nnz_d);
  }
- duneB_.setSize(well_.numberOfSegments(), num_cells, numPerfs);
- duneC_.setSize(well_.numberOfSegments(), num_cells, numPerfs);
+ duneB_.setSize(well_.numberOfSegments(), num_cells_this_process, num_perfs_whole_mswell); //means: sgs x num_cells_this_process entries and num_perfs_whole_mswell nonzero entries
+ duneC_.setSize(well_.numberOfSegments(), num_cells_this_process, num_perfs_whole_mswell);
 
  // we need to add the off diagonal ones
  for (auto row = duneD_.createbegin(),

opm/simulators/wells/MultisegmentWellEquations.hpp

@@ -22,6 +22,8 @@
 #ifndef OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED
 #define OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED
 
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+#include <opm/simulators/wells/WellHelpers.hpp>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
@@ -70,13 +72,15 @@ class MultisegmentWellEquations
  MultisegmentWellEquations(const MultisegmentWellGeneric<Scalar>& well);
 
  //! \brief Setup sparsity pattern for the matrices.
- //! \param num_cells Total number of cells
- //! \param numPerfs Number of perforations
+ //! \param num_cells_this_process Total number of cells on this process
+ //! \param num_perfs_this_process Number of perforations on this process
+ //! \param num_perfs_whole_mswell Number of perforations of the whole well
  //! \param cells Cell indices for perforations
  //! \param segment_inlets Cell indices for segment inlets
  //! \param segment_perforations Cell indices for segment perforations
- void init(const int num_cells,
- const int numPerfs,
+ void init(const int num_cells_this_process,
+ const int num_perfs_this_process,
+ const int num_perfs_whole_mswell,
  const std::vector<int>& cells,
  const std::vector<std::vector<int>>& segment_inlets,
  const std::vector<std::vector<int>>& segment_perforations);
@@ -122,6 +126,9 @@ class MultisegmentWellEquations
  const int seg_pressure_var_ind,
  const WellState<Scalar>& well_state) const;
 
+ //! \brief Sum with off-process contribution.
+ void sumDistributed(Parallel::Communication comm);
+
  //! \brief Returns a const reference to the residual.
  const BVectorWell& residual() const
  {

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -53,8 +53,9 @@ namespace Opm
 
 template<typename FluidSystem, typename Indices>
 MultisegmentWellEval<FluidSystem,Indices>::
-MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices>& baseif)
+MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices>& baseif, const ParallelWellInfo<Scalar>& pw_info)
  : MultisegmentWellGeneric<Scalar>(baseif)
+ , pw_info_(pw_info)
  , baseif_(baseif)
  , linSys_(*this)
  , primary_variables_(baseif)
@@ -69,8 +70,17 @@ void
 MultisegmentWellEval<FluidSystem,Indices>::
 initMatrixAndVectors(const int num_cells)
 {
- linSys_.init(num_cells, baseif_.numPerfs(),
- baseif_.cells(), segments_.inlets(),
+ int num_perfs_this_process = baseif_.numPerfs();
+ int num_perfs_whole_mswell = this->pw_info_.communication().sum(num_perfs_this_process);
+
+ std::vector<int> cells_for_perfs_whole_well(num_perfs_whole_mswell, 0.0);
+ for (int perf = 0; perf < num_perfs_this_process; perf++) {
+ cells_for_perfs_whole_well[this->pw_info_.getGlobalPerfContainerFactory().localToGlobal(perf)] = baseif_.cells()[perf];
+ }
+ this->pw_info_.communication().sum(cells_for_perfs_whole_well.data(), num_perfs_whole_mswell);
+
+ linSys_.init(num_cells, num_perfs_this_process, num_perfs_whole_mswell,
+ cells_for_perfs_whole_well, segments_.inlets(),
  segments_.perforations());
  primary_variables_.resize(this->numberOfSegments());
 }
@@ -433,6 +443,7 @@ MultisegmentWellEval<FluidSystem,Indices>::
 getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
  DeferredLogger& deferred_logger) const
 {
+ std::cout << "B_avg.size() = " << B_avg.size() << " and baseif_.numComponents() = " << baseif_.numComponents() << std::endl;
  assert(int(B_avg.size() ) == baseif_.numComponents());
  std::vector<Scalar> residuals(numWellEq + 1, 0.0);
 

opm/simulators/wells/MultisegmentWellEval.hpp

@@ -67,9 +67,10 @@ class MultisegmentWellEval : public MultisegmentWellGeneric<typename FluidSystem
  //! \brief Returns a const reference to equation system.
  const Equations& linSys() const
  { return linSys_; }
+ const ParallelWellInfo<Scalar>& pw_info_;
 
 protected:
- MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices>& baseif);
+ MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices>& baseif, const ParallelWellInfo<Scalar>& pw_info);
 
  void initMatrixAndVectors(const int num_cells);
 

opm/simulators/wells/MultisegmentWellGeneric.cpp

@@ -81,7 +81,8 @@ scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inle
  }
 
  std::vector<Scalar> rates;
- WellState<Scalar>::calculateSegmentRates(segment_inlets,
+ WellState<Scalar>::calculateSegmentRates(ws,
+ segment_inlets,
  segment_perforations,
  perforation_rates,
  num_single_phase, 0, rates);

opm/simulators/wells/MultisegmentWellSegments.cpp

@@ -62,7 +62,7 @@ MultisegmentWellSegments<FluidSystem,Indices>::
 MultisegmentWellSegments(const int numSegments,
  WellInterfaceGeneric<Scalar>& well)
  : perforations_(numSegments)
- , perforation_depth_diffs_(well.numPerfs(), 0.0)
+ , perforation_depth_diffs_(well.wellEcl().getConnections().size(), 0.0)
  , inlets_(well.wellEcl().getSegments().size())
  , depth_diffs_(numSegments, 0.0)
  , densities_(numSegments, 0.0)
@@ -84,8 +84,8 @@ MultisegmentWellSegments(const int numSegments,
  // well_ecl_ and wells struct different
  // the current implementation is a temporary solution for now, it should be corrected from the parser
  // side
- int i_perf_wells = 0;
- well.perfDepth().resize(well_.numPerfs(), 0.);
+ // The perfDepth vector contains all perforations across all processes of this well!
+ well.perfDepth().resize(completion_set.size(), 0.);
  const auto& segment_set = well_.wellEcl().getSegments();
  for (std::size_t perf = 0; perf < completion_set.size(); ++perf) {
  const Connection& connection = completion_set.get(perf);
@@ -98,11 +98,10 @@ MultisegmentWellSegments(const int numSegments,
  connection.getI() + 1, connection.getJ() + 1,
  connection.getK() + 1));
  }
- perforations_[segment_index].push_back(i_perf_wells);
- well.perfDepth()[i_perf_wells] = connection.depth();
+ perforations_[segment_index].push_back(perf);
+ well.perfDepth()[perf] = connection.depth();
  const Scalar segment_depth = segment_set[segment_index].depth();
- perforation_depth_diffs_[i_perf_wells] = well_.perfDepth()[i_perf_wells] - segment_depth;
- i_perf_wells++;
+ perforation_depth_diffs_[perf] = well_.perfDepth()[perf] - segment_depth;
  }
  }
 
@@ -345,6 +344,7 @@ MultisegmentWellSegments<FluidSystem,Indices>::
 getPressureDiffSegPerf(const int seg,
  const int perf) const
 {
+ // Attention, perf here has a global numbering!
  return well_.gravity() * densities_[seg].value() * perforation_depth_diffs_[perf];
 }