NonlinearSolverEbos: put stabilizeNonlinearUpdate in compile unit

opm/simulators/flow/NonlinearSolverEbos.cpp

@@ -20,7 +20,13 @@
 #include <config.h>
 #include <opm/simulators/flow/NonlinearSolverEbos.hpp>
 
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+
+#include <opm/common/ErrorMacros.hpp>
+
 #include <cmath>
+#include <stdexcept>
 
 namespace Opm {
 namespace detail {
@@ -59,5 +65,56 @@ void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
     oscillate = (oscillatePhase > 1);
 }
 
+template <class BVector>
+void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
+                              const double omega,
+                              NonlinearRelaxType relaxType)
+{
+    // The dxOld is updated with dx.
+    // If omega is equal to 1., no relaxtion will be appiled.
+
+    BVector tempDxOld = dxOld;
+    dxOld = dx;
+
+    switch (relaxType) {
+    case NonlinearRelaxType::Dampen: {
+        if (omega == 1.) {
+            return;
+        }
+        for (auto& d : dx) {
+            d *= omega;
+        }
+        return;
+    }
+    case NonlinearRelaxType::SOR: {
+        if (omega == 1.) {
+            return;
+        }
+        auto i = dx.size();
+        for (i = 0; i < dx.size(); ++i) {
+            dx[i] *= omega;
+            tempDxOld[i] *= (1.-omega);
+            dx[i] += tempDxOld[i];
+        }
+        return;
+    }
+    default:
+        OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
+    }
+
+    return;
+}
+
+template<int Size> using BV = Dune::BlockVector<Dune::FieldVector<double,Size>>;
+#define INSTANCE(Size) \
+    template void stabilizeNonlinearUpdate<BV<Size>>(BV<Size>&, BV<Size>&, \
+                                                  const double, NonlinearRelaxType);
+INSTANCE(1)
+INSTANCE(2)
+INSTANCE(3)
+INSTANCE(4)
+INSTANCE(5)
+INSTANCE(6)
+
 } // namespace detail
 } // namespace Opm

opm/simulators/flow/NonlinearSolverEbos.hpp

@@ -80,15 +80,6 @@ struct NewtonRelaxationType<TypeTag, TTag::FlowNonLinearSolver> {
 
 namespace Opm {
 
-namespace detail {
-
-/// Detect oscillation or stagnation in a given residual history.
-void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
-                        const int it, const int numPhases, const double relaxRelTol,
-                        bool& oscillate, bool& stagnate);
-
-}
-
 class WellState;
 
 // Available relaxation scheme types.
@@ -97,6 +88,21 @@ enum class NonlinearRelaxType {
     SOR
 };
 
+namespace detail {
+
+/// Detect oscillation or stagnation in a given residual history.
+void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
+                        const int it, const int numPhases, const double relaxRelTol,
+                        bool& oscillate, bool& stagnate);
+
+/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
+/// Implemention for Dune block vectors.
+template <class BVector>
+void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
+                              const double omega, NonlinearRelaxType relaxType);
+
+}
+
     /// A nonlinear solver class suitable for general fully-implicit models,
     /// as well as pressure, transport and sequential models.
     template <class TypeTag, class PhysicalModel>
@@ -298,40 +304,7 @@ enum class NonlinearRelaxType {
         template <class BVector>
         void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
         {
-            // The dxOld is updated with dx.
+            detail::stabilizeNonlinearUpdate(dx, dxOld, omega, this->relaxType());
-            // If omega is equal to 1., no relaxtion will be appiled.
-
-            BVector tempDxOld = dxOld;
-            dxOld = dx;
-
-            switch (relaxType()) {
-            case NonlinearRelaxType::Dampen: {
-                if (omega == 1.) {
-                    return;
-                }
-                auto i = dx.size();
-                for (i = 0; i < dx.size(); ++i) {
-                    dx[i] *= omega;
-                }
-                return;
-            }
-            case NonlinearRelaxType::SOR: {
-                if (omega == 1.) {
-                    return;
-                }
-                auto i = dx.size();
-                for (i = 0; i < dx.size(); ++i) {
-                    dx[i] *= omega;
-                    tempDxOld[i] *= (1.-omega);
-                    dx[i] += tempDxOld[i];
-                }
-                return;
-            }
-            default:
-                OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
-            }
-
-            return;
         }
 
         /// The greatest relaxation factor (i.e. smallest factor) allowed.