diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index c0def89e056b20c4e1134c81d0b82f64042df6ab..261d2434d61efef7792b0ed674ac1783dfd49e52 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -28,12 +28,17 @@
#include <memory>
#include <iostream>
#include <type_traits>
+#include <algorithm>
+#include <numeric>
#include <dune/common/timer.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpicommunication.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/std/type_traits.hh>
+#include <dune/common/indices.hh>
+#include <dune/common/hybridutilities.hh>
+
#include <dune/istl/bvector.hh>
#include <dune/istl/multitypeblockvector.hh>
@@ -69,6 +74,81 @@ using DetectPVSwitch = typename Assembler::GridVariables::VolumeVariables::Prima
template<class Assembler>
using GetPVSwitch = Dune::Std::detected_or<int, DetectPVSwitch, Assembler>;
+// helpers to implement max relative shift
+template<class C> using dynamicIndexAccess = decltype(std::declval<C>()[0]);
+template<class C> using staticIndexAccess = decltype(std::declval<C>()[Dune::Indices::_0]);
+template<class C> static constexpr auto hasDynamicIndexAccess = Dune::Std::is_detected<dynamicIndexAccess, C>{};
+template<class C> static constexpr auto hasStaticIndexAccess = Dune::Std::is_detected<staticIndexAccess, C>{};
+
+template<class V, class Scalar, class Reduce, class Transform>
+auto hybridInnerProduct(const V& v1, const V& v2, Scalar init, Reduce&& r, Transform&& t)
+-> std::enable_if_t<hasDynamicIndexAccess<V>(), Scalar>
+{
+ return std::inner_product(v1.begin(), v1.end(), v2.begin(), init, std::forward<Reduce>(r), std::forward<Transform>(t));
+}
+
+template<class V, class Scalar, class Reduce, class Transform>
+auto hybridInnerProduct(const V& v1, const V& v2, Scalar init, Reduce&& r, Transform&& t)
+-> std::enable_if_t<hasStaticIndexAccess<V>() && !hasDynamicIndexAccess<V>(), Scalar>
+{
+ using namespace Dune::Hybrid;
+ forEach(std::make_index_sequence<V::N()>{}, [&](auto i){
+ init = r(init, hybridInnerProduct(v1[Dune::index_constant<i>{}], v2[Dune::index_constant<i>{}], init, std::forward<Reduce>(r), std::forward<Transform>(t)));
+ });
+ return init;
+}
+
+// Maximum relative shift at a degree of freedom.
+// For (primary variables) values below 1.0 we use
+// an absolute shift.
+template<class Scalar, class V>
+auto maxRelativeShift(const V& v1, const V& v2)
+-> std::enable_if_t<Dune::IsNumber<V>::value, Scalar>
+{
+ using std::abs; using std::max;
+ return abs(v1 - v2)/max<Scalar>(1.0, abs(v1 + v2)*0.5);
+}
+
+// Maximum relative shift for generic vector types.
+// Recursively calls maxRelativeShift until Dune::IsNumber is true.
+template<class Scalar, class V>
+auto maxRelativeShift(const V& v1, const V& v2)
+-> std::enable_if_t<!Dune::IsNumber<V>::value, Scalar>
+{
+ return hybridInnerProduct(v1, v2, Scalar(0.0),
+ [](const auto& a, const auto& b){ using std::max; return max(a, b); },
+ [](const auto& a, const auto& b){ return maxRelativeShift<Scalar>(a, b); }
+ );
+}
+
+template<class To, class From>
+void assign(To& to, const From& from)
+{
+ if constexpr (std::is_assignable<To&, From>::value)
+ to = from;
+
+ else if constexpr (hasStaticIndexAccess<To>() && hasStaticIndexAccess<To>() && !hasDynamicIndexAccess<From>() && !hasDynamicIndexAccess<From>())
+ {
+ using namespace Dune::Hybrid;
+ forEach(std::make_index_sequence<To::N()>{}, [&](auto i){
+ assign(to[Dune::index_constant<i>{}], from[Dune::index_constant<i>{}]);
+ });
+ }
+
+ else if constexpr (hasDynamicIndexAccess<From>() && hasDynamicIndexAccess<From>())
+ for (decltype(to.size()) i = 0; i < to.size(); ++i)
+ assign(to[i], from[i]);
+
+ else
+ DUNE_THROW(Dune::Exception, "Values are not assignable to each other!");
+}
+
+template<class T, std::enable_if_t<Dune::IsNumber<std::decay_t<T>>::value, int> = 0>
+constexpr std::size_t blockSize() { return 1; }
+
+template<class T, std::enable_if_t<!Dune::IsNumber<std::decay_t<T>>::value, int> = 0>
+constexpr std::size_t blockSize() { return std::decay_t<T>::size(); }
+
} // end namespace Detail
/*!
@@ -949,44 +1029,14 @@ private:
virtual void newtonUpdateShift_(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
- shift_ = 0;
- newtonUpdateShiftImpl_(uLastIter, deltaU);
+ auto uNew = uLastIter;
+ uNew -= deltaU;
+ shift_ = Detail::maxRelativeShift<Scalar>(uLastIter, uNew);
if (comm_.size() > 1)
shift_ = comm_.max(shift_);
}
- template<class SolVec>
- void newtonUpdateShiftImpl_(const SolVec &uLastIter,
- const SolVec &deltaU)
- {
- for (int i = 0; i < int(uLastIter.size()); ++i)
- {
- auto uNewI = uLastIter[i];
- uNewI -= deltaU[i];
-
- Scalar shiftAtDof = relativeShiftAtDof_(uLastIter[i], uNewI);
- using std::max;
- shift_ = max(shift_, shiftAtDof);
- }
- }
-
- template<class ...Args>
- void newtonUpdateShiftImpl_(const Dune::MultiTypeBlockVector<Args...> &uLastIter,
- const Dune::MultiTypeBlockVector<Args...> &deltaU)
- {
- // There seems to be a bug in g++5 which which prevents compilation when
- // passing the call to the implementation directly to Dune::Hybrid::forEach.
- // We therefore store this call in a lambda and pass it to the for loop afterwards.
- auto doUpdate = [&](const auto subVectorIdx)
- {
- this->newtonUpdateShiftImpl_(uLastIter[subVectorIdx], deltaU[subVectorIdx]);
- };
-
- using namespace Dune::Hybrid;
- forEach(integralRange(Dune::Hybrid::size(uLastIter)), doUpdate);
- }
-
virtual void lineSearchUpdate_(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
@@ -1050,19 +1100,14 @@ private:
//! to this field vector type in Dune ISTL
//! Could be avoided for vectors that already have the right type using SFINAE
//! but it shouldn't impact performance too much
- constexpr auto blockSize = std::decay_t<decltype(b[0])>::dimension;
+ constexpr auto blockSize = Detail::blockSize<decltype(b[0])>();
using BlockType = Dune::FieldVector<Scalar, blockSize>;
Dune::BlockVector<BlockType> xTmp; xTmp.resize(b.size());
Dune::BlockVector<BlockType> bTmp(xTmp);
- for (unsigned int i = 0; i < b.size(); ++i)
- for (unsigned int j = 0; j < blockSize; ++j)
- bTmp[i][j] = b[i][j];
+ Detail::assign(bTmp, b);
const int converged = ls.solve(A, xTmp, bTmp);
-
- for (unsigned int i = 0; i < x.size(); ++i)
- for (unsigned int j = 0; j < blockSize; ++j)
- x[i][j] = xTmp[i][j];
+ Detail::assign(x, xTmp);
return converged;
}
@@ -1188,7 +1233,7 @@ private:
nextPriVars -= uDelta[i];
// add the current relative shift for this degree of freedom
- auto shift = relativeShiftAtDof_(currentPriVars, nextPriVars);
+ auto shift = Detail::maxRelativeShift<Scalar>(currentPriVars, nextPriVars);
distanceFromLastLinearization_[i] += shift;
}
}
@@ -1213,30 +1258,6 @@ private:
DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector");
}
- /*!
- * \brief Returns the maximum relative shift between two vectors of
- * primary variables.
- *
- * \param priVars1 The first vector of primary variables
- * \param priVars2 The second vector of primary variables
- */
- template<class PrimaryVariables>
- Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1,
- const PrimaryVariables &priVars2) const
- {
- Scalar result = 0.0;
- using std::abs;
- using std::max;
- // iterate over all primary variables
- for (int j = 0; j < PrimaryVariables::dimension; ++j) {
- Scalar eqErr = abs(priVars1[j] - priVars2[j]);
- eqErr /= max<Scalar>(1.0,abs(priVars1[j] + priVars2[j])/2);
-
- result = max(result, eqErr);
- }
- return result;
- }
-
//! The communication object
Communication comm_;
diff --git a/test/nonlinear/newton/test_newton.cc b/test/nonlinear/newton/test_newton.cc
index 5f37e2c3656bbdaa8f0574888023e0823139cd6b..7394386149898ec3212761a79b835f00c57b6195 100644
--- a/test/nonlinear/newton/test_newton.cc
+++ b/test/nonlinear/newton/test_newton.cc
@@ -7,6 +7,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>
#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/bvector.hh>
#include <dumux/nonlinear/newtonsolver.hh>
/*
@@ -27,40 +28,10 @@
namespace Dumux {
-class MockScalarVariables
-{
-public:
- static constexpr int dimension = 1;
- MockScalarVariables() : var_(0.0) {}
- explicit MockScalarVariables(double v) : var_(v) {}
- MockScalarVariables& operator-=(const MockScalarVariables& other) { var_ -= other.var_; return *this; }
- double& operator[] (int i) { return var_; }
- const double& operator[] (int i) const { return var_; }
-private:
- double var_;
-};
-
-class MockScalarResidual
-{
-public:
- MockScalarResidual() : res_(0.0) {}
- explicit MockScalarResidual(double r) : res_(r) {}
- MockScalarResidual& operator=(double r) { res_[0] = r; return *this; }
- MockScalarResidual& operator*=(double a) { res_[0] *= a; return *this; }
- MockScalarResidual& operator+=(const MockScalarResidual& other) { res_[0] += other.res_[0]; return *this; }
- MockScalarResidual& operator-=(const MockScalarResidual& other) { res_[0] -= other.res_[0]; return *this; }
- constexpr std::size_t size() const { return 1; }
- MockScalarVariables& operator[] (int i) { return res_; }
- const MockScalarVariables& operator[] (int i) const { return res_; }
- double two_norm2() const { return res_[0]*res_[0]; }
-private:
- MockScalarVariables res_;
-};
-
class MockScalarAssembler
{
public:
- using ResidualType = MockScalarResidual;
+ using ResidualType = Dune::BlockVector<double>;
using Scalar = double;
using JacobianMatrix = double;
@@ -74,13 +45,14 @@ public:
void assembleResidual(const ResidualType& sol)
{
- res_[0][0] = sol[0][0]*sol[0][0] - 5.0;
+ res_.resize(1);
+ res_[0] = sol[0]*sol[0] - 5.0;
}
void assembleJacobianAndResidual (const ResidualType& sol)
{
assembleResidual(sol);
- jac_ = 2.0*sol[0][0];
+ jac_ = 2.0*sol[0];
}
JacobianMatrix& jacobian() { return jac_; }
@@ -90,7 +62,7 @@ public:
double residualNorm(const ResidualType& sol)
{
assembleResidual(sol);
- return res_[0][0];
+ return res_[0];
}
void updateGridVariables(const ResidualType& sol) {}
@@ -108,7 +80,7 @@ public:
template<class Vector>
bool solve (const double& A, Vector& x, const Vector& b) const
{
- x[0][0] = b[0][0]/A;
+ x[0] = b[0]/A;
return true;
}
};
@@ -137,16 +109,17 @@ int main(int argc, char* argv[]) try
auto solver = std::make_shared<Solver>(assembler, linearSolver);
double initialGuess = 0.1;
- MockScalarResidual x(initialGuess);
+ Dune::BlockVector<double> x(1);
+ x = initialGuess;
std::cout << "Solving: x^2 - 5 = 0" << std::endl;
solver->solve(x);
- std::cout << "Solution: " << std::setprecision(15) << x[0][0]
+ std::cout << "Solution: " << std::setprecision(15) << x[0]
<< ", exact: " << std::sqrt(5.0)
- << ", error: " << std::abs(x[0][0]-std::sqrt(5.0))/std::sqrt(5.0)*100 << "%" << std::endl;
+ << ", error: " << std::abs(x[0]-std::sqrt(5.0))/std::sqrt(5.0)*100 << "%" << std::endl;
- if (Dune::FloatCmp::ne(x[0][0], std::sqrt(5.0), 1e-13))
- DUNE_THROW(Dune::Exception, "Didn't find correct root: " << std::setprecision(15) << x[0][0] << ", exact: " << std::sqrt(5.0));
+ if (Dune::FloatCmp::ne(x[0], std::sqrt(5.0), 1e-13))
+ DUNE_THROW(Dune::Exception, "Didn't find correct root: " << std::setprecision(15) << x[0] << ", exact: " << std::sqrt(5.0));
return 0;