Add minIter to the NewtonRaphson settings (#327)

Author: tarun-mitruka

ikarus/finiteelements/mechanics/materials/vanishingstress.hh

@@ -14,6 +14,7 @@
 #include <ikarus/finiteelements/mechanics/materials/interface.hh>
 #include <ikarus/solver/nonlinearsolver/newtonraphson.hh>
 #include <ikarus/solver/nonlinearsolver/nonlinearsolverfactory.hh>
+#include <ikarus/utils/concepts.hh>
 #include <ikarus/utils/nonlinearoperator.hh>
 
 namespace Ikarus {
@@ -37,6 +38,7 @@ struct VanishingStress : public Material<VanishingStress<stressIndexPair, MI>>
       countDiagonalIndices(fixedPairs);                                ///< Number of fixed diagonal indices.
   static constexpr auto freeStrains = freeVoigtIndices.size();         ///< Number of free strains.
   using ScalarType                  = typename Underlying::ScalarType; ///< Scalar type.
+  static constexpr bool isAutoDiff  = Concepts::AutodiffScalar<ScalarType>;
 
   static constexpr auto strainTag              = Underlying::strainTag;            ///< Strain tag.
   static constexpr auto stressTag              = Underlying::stressTag;            ///< Stress tag.
@@ -185,10 +187,13 @@ private:
         E(indexPair[0], indexPair[1]) += ecomps(ri++);
       }
     };
+
+    int minIter = isAutoDiff ? 1 : 0;
     // THE CTAD is broken for designated initializers in clang 16, when we drop support this can be simplified
     NewtonRaphsonConfig<decltype(linearSolver), decltype(updateFunction)> nrs{
-        .parameters = {.tol = tol_, .maxIter = 100},
-          .linearSolver = linearSolver, .updateFunction = updateFunction
+        .parameters     = {.tol = tol_, .maxIter = 100, .minIter = minIter},
+        .linearSolver   = linearSolver,
+        .updateFunction = updateFunction
     };
 
     auto nr = createNonlinearSolver(std::move(nrs), nonOp);

ikarus/solver/nonlinearsolver/newtonraphson.hh

@@ -25,6 +25,7 @@ struct NRSettings
 {
   double tol{1e-8};
   int maxIter{20};
+  int minIter{0};
 };
 
 /**
@@ -131,7 +132,12 @@ public:
    * \brief Set up the solver with the given settings.
    * \param settings Newton-Raphson settings.
    */
-  void setup(const Settings& settings) { settings_ = settings; }
+  void setup(const Settings& settings) {
+    if (settings.minIter > settings.maxIter)
+      DUNE_THROW(Dune::InvalidStateException,
+                 "Minimum number of iterations cannot be greater than maximum number of iterations");
+    settings_ = settings;
+  }
 
 #ifndef DOXYGEN
   struct NoPredictor
@@ -164,7 +170,7 @@ public:
     int iter{0};
     if constexpr (isLinearSolver)
       linearSolver_.analyzePattern(Ax);
-    while (rNorm > settings_.tol && iter < settings_.maxIter) {
+    while ((rNorm > settings_.tol && iter < settings_.maxIter) or iter < settings_.minIter) {
       this->notify(NonLinearSolverMessages::ITERATION_STARTED);
       if constexpr (isLinearSolver) {
         linearSolver_.factorize(Ax);

ikarus/utils/concepts.hh

@@ -19,6 +19,8 @@
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
+#include <autodiff/forward/dual/dual.hpp>
+
 #include "ikarus/assembler/dirichletbcenforcement.hh"
 #include "ikarus/finiteelements/mechanics/materials/tags.hh"
 #include <ikarus/utils/traits.hh>
@@ -584,5 +586,26 @@ namespace Concepts {
     { g.grid() };
   };
 
+  namespace Impl {
+    template <typename T>
+    struct is_dual : std::false_type
+    {
+    };
+
+    // Specialization for Dual<T, U>: this will be true for Dual types
+    template <typename T, typename U>
+    struct is_dual<autodiff::detail::Dual<T, U>> : std::true_type
+    {
+    };
+  } // namespace Impl
+
+  /**
+  \concept AutodiffScalar
+  \brief Concept to check if the underlying scalar type is a dual type.
+  \tparam T The scalar type to be checked if it is a dual type.
+  */
+  template <typename T>
+  concept AutodiffScalar = Impl::is_dual<T>::value;
+
 } // namespace Concepts
 } // namespace Ikarus

tests/src/checkfebyautodiff.hh

@@ -2,6 +2,9 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
 
 #pragma once
+
+#include "testhelpers.hh"
+
 #include <dune/common/float_cmp.hh>
 #include <dune/common/test/testsuite.hh>
 
@@ -11,19 +14,17 @@
 #include <ikarus/finiteelements/fefactory.hh>
 #include <ikarus/utils/basis.hh>
 
-template <typename GridView, typename PreBasis, typename Skills, typename AffordanceColl>
-auto checkFESByAutoDiff(const GridView& gridView, const PreBasis& pb, Skills&& skills, AffordanceColl affordance) {
-  auto basis = Ikarus::makeBasis(gridView, pb);
-  Eigen::VectorXd d;
-  d.setRandom(basis.flat().dimension());
+template <typename GridView, typename BasisHandler, typename Skills, typename AffordanceColl, typename VectorType>
+auto checkFESByAutoDiffImpl(const GridView& gridView, const BasisHandler& basis, Skills&& skills,
+                            AffordanceColl affordance, VectorType& d, const std::string& messageIfFailed = "") {
   double lambda = 7.3;
-
-  auto fe  = Ikarus::makeFE(basis, std::forward<Skills>(skills));
-  using FE = decltype(fe);
+  auto fe       = Ikarus::makeFE(basis, std::forward<Skills>(skills));
+  using FE      = decltype(fe);
 
   auto req = typename FE::Requirement();
   req.insertGlobalSolution(d).insertParameter(lambda);
-  Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation");
+  Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation" +
+                    messageIfFailed);
   for (auto element : elements(gridView)) {
     auto localView = basis.flat().localView();
     localView.bind(element);
@@ -51,7 +52,7 @@ auto checkFESByAutoDiff(const GridView& gridView, const PreBasis& pb, Skills&& s
     calculateVector(fe, req, affordance.vectorAffordance(), R);
     calculateVector(feAutoDiff, req, affordance.vectorAffordance(), RAutoDiff);
 
-    t.check(K.isApprox(KAutoDiff, tol),
+    t.check(isApproxSame(K, KAutoDiff, tol),
             "Mismatch between the stiffness matrices obtained from explicit implementation and the one based on "
             "automatic differentiation for " +
                 feClassName)
@@ -60,7 +61,7 @@ auto checkFESByAutoDiff(const GridView& gridView, const PreBasis& pb, Skills&& s
         << KAutoDiff << "\nThe difference is\n"
         << (K - KAutoDiff);
 
-    t.check(R.isApprox(RAutoDiff, tol),
+    t.check(isApproxSame(R, RAutoDiff, tol),
             "Mismatch between the residual vectors obtained from explicit implementation and the one based on "
             "automatic differentiation for " +
                 feClassName)
@@ -76,3 +77,16 @@ auto checkFESByAutoDiff(const GridView& gridView, const PreBasis& pb, Skills&& s
 
   return t;
 }
+
+template <typename GridView, typename PreBasis, typename Skills, typename AffordanceColl>
+auto checkFESByAutoDiff(const GridView& gridView, const PreBasis& pb, Skills&& skills, AffordanceColl affordance,
+                        const std::string& testName = "") {
+  Dune::TestSuite t("AutoDiff Test" + testName);
+  auto basis = Ikarus::makeBasis(gridView, pb);
+  Eigen::VectorXd d;
+  d.setZero(basis.flat().dimension());
+  t.subTest(checkFESByAutoDiffImpl(gridView, basis, skills, affordance, d, " Zero Displacements"));
+  d.setRandom(basis.flat().dimension());
+  t.subTest(checkFESByAutoDiffImpl(gridView, basis, skills, affordance, d, " Non-zero Displacements"));
+  return t;
+}

tests/src/testnonlinearelasticity.hh

@@ -3,6 +3,7 @@
 
 #include <config.h>
 
+#include "checkfebyautodiff.hh"
 #include "resultcollection.hh"
 #include "testcommon.hh"
 
@@ -287,3 +288,34 @@ auto SingleElementTest(const Material& mat) {
   }
   return t;
 }
+
+template <int gridDim, typename MAT, typename TestSuiteType>
+void autoDiffTest(TestSuiteType& t, const MAT& mat, const std::string& testName = "") {
+  using namespace Ikarus;
+  using namespace Dune::Functions::BasisFactory;
+  auto vL = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
+    Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fExt;
+    fExt.setZero();
+    fExt[1] = 2 * lamb;
+    return fExt;
+  };
+
+  auto nBL = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
+    Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fExt;
+    fExt.setZero();
+    fExt[0] = lamb / 40;
+    return fExt;
+  };
+
+  auto grid     = createUGGridFromCorners<gridDim>(CornerDistortionFlag::randomlyDistorted);
+  auto gridView = grid->leafGridView();
+
+  /// We artificially apply a Neumann load on the complete boundary
+  Dune::BitSetVector<1> neumannVertices(gridView.size(gridDim), true);
+  BoundaryPatch neumannBoundary(gridView, neumannVertices);
+
+  t.subTest(checkFESByAutoDiff(
+      gridView, power<gridDim>(lagrange<1>()),
+      skills(Ikarus::nonLinearElastic(mat), volumeLoad<gridDim>(vL), neumannBoundaryLoad(&neumannBoundary, nBL)),
+      Ikarus::AffordanceCollections::elastoStatics, testName));
+}

tests/src/testnonlinearelasticityneohooke.cpp

@@ -11,15 +11,26 @@
 using Dune::TestSuite;
 
 int main(int argc, char** argv) {
+  using namespace Ikarus;
   Ikarus::init(argc, argv);
-  TestSuite t;
+  TestSuite t("NonLinearElastic + NeoHooke Test");
 
-  auto matParameter = Ikarus::toLamesFirstParameterAndShearModulus({.emodul = 1000, .nu = 0.3});
+  auto matParameter1 = toLamesFirstParameterAndShearModulus({.emodul = 1000, .nu = 0.3});
+  auto matParameter2 = toLamesFirstParameterAndShearModulus({.emodul = 1000, .nu = 0.0});
 
-  Ikarus::NeoHooke matNH(matParameter);
+  NeoHooke matNH1(matParameter1);
+  NeoHooke matNH2(matParameter2);
 
-  t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::Alu>(matNH));
-  t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::Yasp>(matNH));
-  t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::IgaSurfaceIn2D>(matNH));
+  auto planeStressMat1 = planeStress(matNH1, 1e-8);
+  auto planeStressMat2 = planeStress(matNH2, 1e-8);
+
+  t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::Alu>(matNH1));
+  t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::Yasp>(matNH1));
+  t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::IgaSurfaceIn2D>(matNH1));
+
+  autoDiffTest<2>(t, planeStressMat1, " nu != 0");
+  autoDiffTest<2>(t, planeStressMat2, " nu = 0");
+  autoDiffTest<3>(t, matNH1, " nu != 0");
+  autoDiffTest<3>(t, matNH2, " nu = 0");
   return t.exit();
 }

tests/src/testnonlinearelasticitysvk.cpp

@@ -3,7 +3,6 @@
 
 #include <config.h>
 
-#include "checkfebyautodiff.hh"
 #include "testcommon.hh"
 #include "testnonlinearelasticity.hh"
 
@@ -13,74 +12,32 @@ using Dune::TestSuite;
 
 int main(int argc, char** argv) {
   using namespace Ikarus;
-  using namespace Dune::Functions::BasisFactory;
   Ikarus::init(argc, argv);
-  TestSuite t;
+  TestSuite t("NonLinearElastic + StVenantKirchhoff Test");
 
   auto matParameter1 = toLamesFirstParameterAndShearModulus({.emodul = 1000, .nu = 0.3});
   auto matParameter2 = toLamesFirstParameterAndShearModulus({.emodul = 1000, .nu = 0.0});
 
   StVenantKirchhoff matSVK1(matParameter1);
   StVenantKirchhoff matSVK2(matParameter2);
-  auto planeStressMat = planeStress(matSVK2, 1e-8);
-  auto planeStrainMat = planeStrain(matSVK1);
+  auto planeStressMat1 = planeStress(matSVK1, 1e-8);
+  auto planeStressMat2 = planeStress(matSVK2, 1e-8);
+  auto planeStrainMat  = planeStrain(matSVK1);
 
   t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::Alu>(matSVK1));
   t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::Yasp>(matSVK1));
   t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::IgaSurfaceIn2D>(matSVK1));
-  t.subTest(GreenLagrangeStrainTest<2>(planeStressMat));
+  t.subTest(GreenLagrangeStrainTest<2>(planeStressMat1));
+  t.subTest(GreenLagrangeStrainTest<2>(planeStressMat2));
   t.subTest(GreenLagrangeStrainTest<2>(planeStrainMat));
-  t.subTest(GreenLagrangeStrainTest<3>(matSVK2));
-  t.subTest(SingleElementTest(planeStressMat));
-
-  auto vL = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
-    Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fExt;
-    fExt.setZero();
-    fExt[1] = 2 * lamb;
-    return fExt;
-  };
-
-  auto nBL = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
-    Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fExt;
-    fExt.setZero();
-    fExt[0] = lamb / 40;
-    return fExt;
-  };
-
-  {
-    auto grid     = createUGGridFromCorners<2>(CornerDistortionFlag::randomlyDistorted);
-    auto gridView = grid->leafGridView();
-    /// We artificially apply a Neumann load on the complete boundary
-    Dune::BitSetVector<1> neumannVertices(gridView.size(2), true);
-    BoundaryPatch neumannBoundary(gridView, neumannVertices);
-    t.subTest(checkFESByAutoDiff(
-        gridView, power<2>(lagrange<1>()),
-        skills(Ikarus::nonLinearElastic(planeStressMat), volumeLoad<2>(vL), neumannBoundaryLoad(&neumannBoundary, nBL)),
-        Ikarus::AffordanceCollections::elastoStatics));
-  }
-  {
-    auto grid     = createUGGridFromCorners<2>(CornerDistortionFlag::randomlyDistorted);
-    auto gridView = grid->leafGridView();
-    /// We artificially apply a Neumann load on the complete boundary
-    Dune::BitSetVector<1> neumannVertices(gridView.size(2), true);
-    BoundaryPatch neumannBoundary(gridView, neumannVertices);
-    t.subTest(checkFESByAutoDiff(
-        gridView, power<2>(lagrange<1>()),
-        skills(Ikarus::nonLinearElastic(planeStrainMat), volumeLoad<2>(vL), neumannBoundaryLoad(&neumannBoundary, nBL)),
-        Ikarus::AffordanceCollections::elastoStatics));
-  }
-
-  {
-    auto grid     = createUGGridFromCorners<3>(CornerDistortionFlag::randomlyDistorted);
-    auto gridView = grid->leafGridView();
-    /// We artificially apply a Neumann load on the complete boundary
-    Dune::BitSetVector<1> neumannVertices(gridView.size(3), true);
-    BoundaryPatch neumannBoundary(gridView, neumannVertices);
-    t.subTest(checkFESByAutoDiff(
-        gridView, power<3>(lagrange<1>()),
-        skills(Ikarus::nonLinearElastic(matSVK1), volumeLoad<3>(vL), neumannBoundaryLoad(&neumannBoundary, nBL)),
-        Ikarus::AffordanceCollections::elastoStatics));
-  }
+  t.subTest(GreenLagrangeStrainTest<3>(matSVK1));
+  t.subTest(SingleElementTest(planeStressMat2));
+
+  autoDiffTest<2>(t, planeStressMat1, " nu != 0");
+  autoDiffTest<2>(t, planeStressMat2, " nu = 0");
+  autoDiffTest<2>(t, planeStrainMat, " nu != 0");
+  autoDiffTest<3>(t, matSVK1, " nu != 0");
+  autoDiffTest<3>(t, matSVK2, " nu = 0");
 
   return t.exit();
 }

tests/src/testnonlinearoperator.cpp

@@ -23,6 +23,8 @@ template <typename SolutionType, typename SolutionTypeExpected, typename NewtonR
 auto checkNewtonRaphson(NewtonRaphson& nr, SolutionType& x, double tolerance, int maxIter, int iterExpected,
                         const SolutionTypeExpected& xExpected, const auto& x_Predictor) {
   TestSuite t("checkNewtonRaphson");
+  t.checkThrow<Dune::InvalidStateException>([&]() { nr.setup({tolerance, maxIter, maxIter + 1}); },
+                                            "NewtonRaphson setup should fail if minIter > maxIter");
   nr.setup({tolerance, maxIter});
   const auto solverInfo = nr.solve(x_Predictor);