Lgr wells

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -82,6 +82,23 @@
 
 #include <fmt/format.h>
 
+// #include <opm/input/eclipse/Schedule/WellTraj/RigEclipseWellLogExtractor.hpp>
+// #include <external/resinsight/ReservoirDataModel/RigWellLogExtractionTools.h>
+// #include <external/resinsight/ReservoirDataModel/RigWellPath.h>
+// #include <external/resinsight/ReservoirDataModel/cvfGeometryTools.h>
+// #include <external/resinsight/ReservoirDataModel/RigWellLogExtractor.h>
+// #include <external/resinsight/ReservoirDataModel/RigCellGeometryTools.h>
+// #include <external/resinsight/CommonCode/cvfStructGrid.h>
+// #include <external/resinsight/LibGeometry/cvfBoundingBox.h>
+// #include <opm/grid/common/CartesianIndexMapper.hpp>
+// #include <dune/geometry/referenceelements.hh>
+// #include <dune/grid/common/mcmgmapper.hh>
+// #include <dune/grid/common/gridenums.hh>
+
+#include <opm/input/eclipse/Schedule/WellTraj/RigEclipseWellLogExtractorGrid.hpp>
+#include <opm/input/eclipse/Schedule/WellTraj/RigEclipseWellLogExtractorGrid_impl.hpp>//hack
+#include "BoundingBoxTree.hpp"
+
 namespace Opm {
 
 template<class Scalar>
@@ -334,20 +351,29 @@ initializeWellProdIndCalculators()
 
 template<class Scalar>
 void BlackoilWellModelGeneric<Scalar>::
-initializeWellPerfData()
+initializeWellPerfData(Dune::CpGrid* cpgrid)
 {
     well_perf_data_.resize(wells_ecl_.size());
 
     this->conn_idx_map_.clear();
     this->conn_idx_map_.reserve(wells_ecl_.size());
 
     int well_index = 0;
-    for (const auto& well : wells_ecl_) {
+    for (auto& well : wells_ecl_) {
         int connection_index = 0;
 
         // INVALID_ECL_INDEX marks no above perf available
         int connection_index_above = ParallelWellInfo<Scalar>::INVALID_ECL_INDEX;
-
+	bool recalculated = false;
+	auto& connections = well.getConnections();
+	if( ! (cpgrid == NULL) && connections.hasTraj()){
+	  // recalculate connections
+	  external::cvf::ref<external::cvf::BoundingBoxTree> cellSearchTree;
+	  external::buildBoundingBoxTree(cellSearchTree, *cpgrid);
+
+	  connections.recomputeConnections(*cpgrid,cellSearchTree);
+	  recalculated = true;
+	} 
         well_perf_data_[well_index].clear();
         well_perf_data_[well_index].reserve(well.getConnections().size());
 
@@ -365,9 +391,13 @@ initializeWellPerfData()
         parallelWellInfo.beginReset();
 
         for (const auto& connection : well.getConnections()) {
-            const auto active_index =
-                this->compressedIndexForInterior(connection.global_index());
-
+	  int active_index;
+	  if(!recalculated){
+            active_index =
+	      this->compressedIndexForInterior(connection.global_index()); // completly strang for LGR
+	  } else {
+	    active_index = connection.global_index();
+	  }
             const auto connIsOpen =
                 connection.state() == Connection::State::OPEN;
 

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -42,6 +42,8 @@
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/WGState.hpp>
 
+#include <opm/grid/CpGrid.hpp>
+
 #include <cstddef>
 #include <functional>
 #include <map>
@@ -358,7 +360,7 @@ class BlackoilWellModelGeneric
     createLocalParallelWellInfo(const std::vector<Well>& wells);
 
     void initializeWellProdIndCalculators();
-    void initializeWellPerfData();
+  void initializeWellPerfData(Dune::CpGrid * grid = NULL);
 
     bool wasDynamicallyShutThisTimeStep(const int well_index) const;
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -247,6 +247,7 @@ namespace Opm {
         DeferredLogger local_deferredLogger{};
 
         const auto& comm = this->simulator_.vanguard().grid().comm();
+        auto& cpgrid = this->simulator_.vanguard().grid();
 
         // Wells_ecl_ holds this rank's wells, both open and stopped/shut.
         this->wells_ecl_ = this->getLocalWells(reportStepIdx);
@@ -258,7 +259,7 @@ namespace Opm {
         // scope a bit.
         OPM_BEGIN_PARALLEL_TRY_CATCH()
         {
-            this->initializeWellPerfData();
+            this->initializeWellPerfData(&cpgrid);
             this->initializeWellState(reportStepIdx);
             this->wbp_.initializeWBPCalculationService();
 

opm/simulators/wells/BoundingBoxTree.hpp

@@ -0,0 +1,75 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 - 2017 Statoil ASA.
+  Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2016 - 2018 IRIS 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/>.
+*/
+
+#ifndef OPM_BOUNDINGBOXTREE_INCLUDED
+#define OPM_BOUNDINGBOXTREE_INCLUDED
+#include <opm/input/eclipse/Schedule/WellTraj/RigEclipseWellLogExtractor.hpp>
+#include <external/resinsight/ReservoirDataModel/RigWellLogExtractionTools.h>
+#include <external/resinsight/ReservoirDataModel/RigWellPath.h>
+#include <external/resinsight/ReservoirDataModel/cvfGeometryTools.h>
+#include <external/resinsight/ReservoirDataModel/RigWellLogExtractor.h>
+#include <external/resinsight/ReservoirDataModel/RigCellGeometryTools.h>
+#include <external/resinsight/CommonCode/cvfStructGrid.h>
+#include <external/resinsight/LibGeometry/cvfBoundingBox.h>
+#include <opm/grid/common/CartesianIndexMapper.hpp>
+#include <dune/geometry/referenceelements.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/common/gridenums.hh>
+
+#include <opm/input/eclipse/Schedule/WellTraj/RigEclipseWellLogExtractorGrid.hpp>
+//#include <opm/input/eclipse/Schedule/WellTraj/RigEclipseWellLogExtractorGrid_impl.hpp>
+namespace external
+{
+void
+buildBoundingBoxTree(cvf::ref<cvf::BoundingBoxTree>& m_cellSearchTree, const ::Dune::CpGrid& grid)
+{
+    using Grid = ::Dune::CpGrid;
+    using GridView = typename Grid::LeafGridView;
+    const auto& gv = grid.leafGridView();
+    size_t cellCount = gv.size(0);
+    std::vector<size_t> cellIndicesForBoundingBoxes;
+    std::vector<cvf::BoundingBox> cellBoundingBoxes;
+
+    std::array<double, 3> cornerPointArray;
+    cvf::Vec3d cornerPoint;
+    using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+    ElementMapper mapper(gv, Dune::mcmgElementLayout()); // used id sets interally
+    for (const auto& element : Dune::elements(gv)) {
+        int index = mapper.index(element);
+        auto geom = element.geometry();
+        cvf::BoundingBox cellBB;
+        cvf::Vec3d cornerPoint;
+        // NB order should not matter when adding to bounding box: dune ordring and resinsight ordering is different
+        //  dune 0 1 2 3 4 5 6 7 is resinsight 0 1 3 2 4 5 7 6 (i think)
+        for (std::size_t l = 0; l < geom.corners(); l++) {
+            auto cornerPointArray = geom.corner(l);
+            cornerPoint = cvf::Vec3d(cornerPointArray[0], cornerPointArray[1], cornerPointArray[2]);
+            cellBB.add(cornerPoint);
+        }
+        cellIndicesForBoundingBoxes.emplace_back(index);
+        cellBoundingBoxes.emplace_back(cellBB);
+    }
+    m_cellSearchTree = new cvf::BoundingBoxTree;
+    m_cellSearchTree->buildTreeFromBoundingBoxes(cellBoundingBoxes, &cellIndicesForBoundingBoxes);
+}
+}
+#endif // OPM_BLACKOILWELLMODEL_WBP_HEADER_INCLUDED

opm/simulators/wells/WellGroupHelpers.cpp

@@ -92,7 +92,7 @@ namespace Opm {
                       const int phasePos,
                       const bool injector)
     {
-
+        OPM_TIMEFUNCTION();
         Scalar rate = 0.0;
         for (const std::string& groupName : group.groups()) {
             const auto& groupTmp = schedule.getGroup(groupName, reportStepIdx);