Expand Support for Geomechanical Workflows With Fracturing

opm/models/io/vtkmultiwriter.hh

@@ -156,6 +156,11 @@ public:
      * If the grid changes between two calls of beginWrite(), this
      * method _must_ be called before the second beginWrite()!
      */
+
+    void gridViewChanged(const GridView& gridView){
+        gridView_ = gridView;
+        gridChanged();
+    }
     void gridChanged()
     {
         elementMapper_.update(gridView_);
@@ -516,7 +521,8 @@ private:
         managedVectorBuffers_.clear();
     }
 
-    const GridView gridView_;
+    //const GridView gridView_;
+    GridView gridView_;
     ElementMapper elementMapper_;
     VertexMapper vertexMapper_;
 

opm/simulators/flow/MechContainer.cpp

@@ -167,7 +167,7 @@ outputRestart(data::Solution& sol)
         DataEntry{"MECHPOTF", UnitSystem::measure::pressure, &potentialForce_},
         DataEntry{"PRESPOTF", UnitSystem::measure::pressure, &potentialPressForce_},
         DataEntry{"STRAIN",   UnitSystem::measure::identity, &strain_},
-        DataEntry{"STRESS",   UnitSystem::measure::length,   &stress_},
+        DataEntry{"STRESS",   UnitSystem::measure::pressure,   &stress_},
         DataEntry{"TEMPPOTF", UnitSystem::measure::pressure, &potentialTempForce_},
     };
 

opm/simulators/linalg/ISTLSolver.hpp

@@ -337,9 +337,12 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         void initPrepare(const Matrix& M, Vector& b)
         {
             const bool firstcall = (matrix_ == nullptr);
-
+            const bool matrix_changed =  &M != matrix_;
+            if(matrix_changed){
+              force_recreate_ = true;
+            }
             // update matrix entries for solvers.
-            if (firstcall) {
+            if (firstcall || matrix_changed) {
                 // model will not change the matrix object. Hence simply store a pointer
                 // to the original one with a deleter that does nothing.
                 // Outch! We need to be able to scale the linear system! Hence const_cast
@@ -506,14 +509,23 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
                 flexibleSolver_[activeSolverNum_].pre_->update();
             }
         }
-
-
+    public:
+      void setForceReCreate(){
+        force_recreate_ = false;
+      }
+
+    protected:
         /// Return true if we should (re)create the whole solver,
         /// instead of just calling update() on the preconditioner.
-        bool shouldCreateSolver() const
+      bool shouldCreateSolver() //const
         {
             // Decide if we should recreate the solver or just do
             // a minimal preconditioner update.
+            if(force_recreate_){
+              force_recreate_ = false;
+              return true;
+            }
+
             if (flexibleSolver_.empty()) {
                 return true;
             }
@@ -667,6 +679,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
 
         std::shared_ptr< CommunicationType > comm_;
         std::unique_ptr<ElementChunksType> element_chunks_;
+        bool force_recreate_ = false;
     }; // end ISTLSolver
 
 } // namespace Opm

opm/simulators/utils/readDeck.cpp

@@ -182,7 +182,7 @@ namespace {
                                         const Opm::EclipseState&             eclipseState,
                                         const Opm::ParseContext&             parseContext,
                                         const bool                           lowActionParsingStrictness,
-                                        const bool                           keepKeywords,
+                                        const bool                           ,//keepKeywords,
                                         std::shared_ptr<Opm::Python>         python,
                                         std::shared_ptr<Opm::Schedule>&      schedule,
                                         std::unique_ptr<Opm::UDQState>&      udqState,
@@ -194,7 +194,7 @@ namespace {
         if (schedule == nullptr) {
             schedule = std::make_shared<Opm::Schedule>
                 (deck, eclipseState, parseContext,
-                 errorGuard, std::move(python), lowActionParsingStrictness, slaveMode, keepKeywords);
+                 errorGuard, std::move(python), lowActionParsingStrictness, slaveMode, true);
         }
         udqState = std::make_unique<Opm::UDQState>
             ((*schedule)[0].udq().params().undefinedValue());

opm/simulators/wells/ParallelWellInfo.cpp

@@ -806,7 +806,8 @@ checkAllConnectionsFound()
     }
     if (missingCells && pwinfo_.isOwner())
     {
-        OPM_THROW(std::runtime_error, msg);
+      std::cout << "Warning " << msg << std::endl;
+      //         OPM_THROW(std::runtime_error, msg);
     }
     return !missingCells;
 }

opm/simulators/wells/RuntimePerforation.hpp

@@ -35,6 +35,12 @@ struct RuntimePerforation
 
     /// Depth at which the new connection is created.
     double depth{};
+
+    /// Segment number (1-based), if applicable (MS well).
+    int segment{};
+
+    /// MD start/end along branch if applicable (MS well).
+    std::optional<std::pair<double, double>> perf_range{};
 };
 
 } // namespace Opm

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -700,21 +700,25 @@ void WellInterfaceGeneric<Scalar, IndexTraits>::addPerforations(const std::vecto
             this->well_index_[ind] += static_cast<Scalar>(perf.ctf);
         }
         else {
-            this->well_cells_.push_back(perf.cell);
-            this->well_index_.push_back(static_cast<Scalar>(perf.ctf));
-            this->perf_depth_.push_back(static_cast<Scalar>(perf.depth));
-
-            // Not strictly needed.
-            const double nan = std::nan("1");
-            this->perf_rep_radius_.push_back(nan);
-            this->perf_length_.push_back(nan);
-            this->bore_diameters_.push_back(nan);
-
-            // For now use the saturation table for the first cell.
-            this->saturation_table_number_
-                .push_back(this->saturation_table_number_.front());
-
-            ++this->number_of_local_perforations_;
+            std::cout << "Perforation to cell " << perf.cell
+                      << " does not exist in well " << this->name()
+                      << ", skip it. TimeStep cuts do to error in fracture solve?." << std::endl;
+            // assert(false);
+            // this->well_cells_.push_back(perf.cell);
+            // this->well_index_.push_back(static_cast<Scalar>(perf.ctf));
+            // this->perf_depth_.push_back(static_cast<Scalar>(perf.depth));
+
+            // // Not strictly needed.
+            // const double nan = std::nan("1");
+            // this->perf_rep_radius_.push_back(nan);
+            // this->perf_length_.push_back(nan);
+            // this->bore_diameters_.push_back(nan);
+
+            // // For now use the saturation table for the first cell.
+            // this->saturation_table_number_
+            //     .push_back(this->saturation_table_number_.front());
+
+            // ++this->number_of_local_perforations_;
         }
     }
 }