major cleanup with Zach and Zizhou

polyfem · Oct 13, 2023 · cfc4092 · cfc4092
1 parent d77199d
commit cfc4092
Show file tree

Hide file tree

Showing 21 changed files with 83 additions and 192 deletions.
diff --git a/src/polysolve/nonlinear/BFGS.cpp b/src/polysolve/nonlinear/BFGS.cpp
@@ -7,9 +7,9 @@ namespace polysolve::nonlinear
 
     BFGS::BFGS(const json &solver_params,
                const json &linear_solver_params,
-               const double dt, const double characteristic_length,
+               const double characteristic_length,
                spdlog::logger &logger)
-        : Superclass(solver_params, dt, characteristic_length, logger)
+        : Superclass(solver_params, characteristic_length, logger)
     {
         linear_solver = polysolve::linear::Solver::create(
             linear_solver_params["solver"], linear_solver_params["precond"]);
@@ -34,6 +34,8 @@ namespace polysolve::nonlinear
         if (this->descent_strategy == 1)
             this->descent_strategy++;
 
+        assert(m_prev_x.size() > 0);
+        reset_history(m_prev_x.size());
         assert(this->descent_strategy <= 2);
     }
 
@@ -54,7 +56,7 @@ namespace polysolve::nonlinear
         this->descent_strategy = 2;
     }
 
-    bool BFGS::compute_update_direction(
+    void BFGS::compute_update_direction(
         Problem &objFunc,
         const TVector &x,
         const TVector &grad,
@@ -79,7 +81,8 @@ namespace polysolve::nonlinear
                 // warn if using gradient descent
                 m_logger.warn("Unable to factorize Hessian: \"{}\"; reverting to {}",
                               err.what(), this->descent_strategy_name());
-                return compute_update_direction(objFunc, x, grad, direction);
+                compute_update_direction(objFunc, x, grad, direction);
+                return;
             }
 
             TVector y = grad - m_prev_grad;
@@ -94,28 +97,5 @@ namespace polysolve::nonlinear
 
         m_prev_x = x;
         m_prev_grad = grad;
-
-        if (std::isnan(direction.squaredNorm()))
-        {
-            reset_history(x.size());
-            increase_descent_strategy();
-            m_logger.log(
-                this->descent_strategy == 2 ? spdlog::level::warn : spdlog::level::debug,
-                "nan in direction {} (||∇f||={}); reverting to {}",
-                direction.dot(grad), this->descent_strategy_name());
-            return compute_update_direction(objFunc, x, grad, direction);
-        }
-        else if (grad.squaredNorm() != 0 && direction.dot(grad) >= 0)
-        {
-            reset_history(x.size());
-            increase_descent_strategy();
-            m_logger.log(
-                this->descent_strategy == 2 ? spdlog::level::warn : spdlog::level::debug,
-                "BFGS direction is not a descent direction (Δx⋅g={}≥0); reverting to {}",
-                direction.dot(grad), this->descent_strategy_name());
-            return compute_update_direction(objFunc, x, grad, direction);
-        }
-
-        return true;
     }
 } // namespace polysolve::nonlinear
diff --git a/src/polysolve/nonlinear/BFGS.hpp b/src/polysolve/nonlinear/BFGS.hpp
@@ -20,7 +20,7 @@ namespace polysolve::nonlinear
 
         BFGS(const json &solver_params,
              const json &linear_solver_params,
-             const double dt, const double characteristic_length,
+             const double characteristic_length,
              spdlog::logger &logger);
 
         std::string name() const override { return "BFGS"; }
@@ -43,7 +43,7 @@ namespace polysolve::nonlinear
 
         void reset_history(const int ndof);
 
-        virtual bool compute_update_direction(
+        virtual void compute_update_direction(
             Problem &objFunc,
             const TVector &x,
             const TVector &grad,

diff --git a/src/polysolve/nonlinear/BoxConstraintSolver.cpp b/src/polysolve/nonlinear/BoxConstraintSolver.cpp
@@ -10,17 +10,16 @@ namespace polysolve::nonlinear
     std::unique_ptr<Solver> BoxConstraintSolver::create(const std::string &solver,
                                                         const json &solver_params,
                                                         const json &linear_solver_params,
-                                                        const double dt,
                                                         const double characteristic_length,
                                                         spdlog::logger &logger)
     {
         if (solver == "LBFGSB" || solver == "L-BFGS-B")
         {
-            return std::make_unique<LBFGSB>(solver_params, dt, characteristic_length, logger);
+            return std::make_unique<LBFGSB>(solver_params, characteristic_length, logger);
         }
         else if (solver == "MMA")
         {
-            return std::make_unique<MMA>(solver_params, dt, characteristic_length, logger);
+            return std::make_unique<MMA>(solver_params, characteristic_length, logger);
         }
         throw std::runtime_error("Unrecognized solver type: " + solver);
     }
@@ -32,10 +31,9 @@ namespace polysolve::nonlinear
     }
 
     BoxConstraintSolver::BoxConstraintSolver(const json &solver_params,
-                                             const double dt,
                                              const double characteristic_length,
                                              spdlog::logger &logger)
-        : Superclass(solver_params, dt, characteristic_length, logger)
+        : Superclass(solver_params, characteristic_length, logger)
     {
         if (solver_params["max_change"].is_number())
             max_change_val_ = solver_params["max_change"];

diff --git a/src/polysolve/nonlinear/BoxConstraintSolver.hpp b/src/polysolve/nonlinear/BoxConstraintSolver.hpp
@@ -16,7 +16,6 @@ namespace polysolve::nonlinear
         static std::unique_ptr<Solver> create(const std::string &solver,
                                               const json &solver_params,
                                               const json &linear_solver_params,
-                                              const double dt,
                                               const double characteristic_length,
                                               spdlog::logger &logger);
         // List available solvers
@@ -27,7 +26,6 @@ namespace polysolve::nonlinear
         using typename Superclass::TVector;
 
         BoxConstraintSolver(const json &solver_params,
-                            const double dt,
                             const double characteristic_length,
                             spdlog::logger &logger);
 

diff --git a/src/polysolve/nonlinear/DenseNewton.cpp b/src/polysolve/nonlinear/DenseNewton.cpp
@@ -6,11 +6,10 @@ namespace polysolve::nonlinear
     DenseNewton::DenseNewton(
         const json &solver_params,
         const json &linear_solver_params,
-        const double dt, const double characteristic_length,
+        const double characteristic_length,
         spdlog::logger &logger)
         : Superclass(solver_params,
                      linear_solver_params,
-                     dt,
                      characteristic_length,
                      logger)
     {

diff --git a/src/polysolve/nonlinear/DenseNewton.hpp b/src/polysolve/nonlinear/DenseNewton.hpp
@@ -16,7 +16,7 @@ namespace polysolve::nonlinear
 
         DenseNewton(const json &solver_params,
                     const json &linear_solver_params,
-                    const double dt, const double characteristic_length,
+                    const double characteristic_length,
                     spdlog::logger &logger);
 
         std::string name() const override { return "DenseNewton"; }

diff --git a/src/polysolve/nonlinear/GradientDescent.cpp b/src/polysolve/nonlinear/GradientDescent.cpp
@@ -4,10 +4,9 @@ namespace polysolve::nonlinear
 {
 
     GradientDescent::GradientDescent(const json &solver_params_,
-                                     const double dt,
                                      const double characteristic_length,
                                      spdlog::logger &logger)
-        : Superclass(solver_params_, dt, characteristic_length, logger)
+        : Superclass(solver_params_, characteristic_length, logger)
     {
     }
 
@@ -33,15 +32,13 @@ namespace polysolve::nonlinear
         this->descent_strategy = 1;
     }
 
-    bool GradientDescent::compute_update_direction(
+    void GradientDescent::compute_update_direction(
         Problem &objFunc,
         const TVector &x,
         const TVector &grad,
         TVector &direction)
     {
         direction = -grad;
-
-        return true;
     }
 
 } // namespace polysolve::nonlinear
diff --git a/src/polysolve/nonlinear/GradientDescent.hpp b/src/polysolve/nonlinear/GradientDescent.hpp
@@ -15,7 +15,6 @@ namespace polysolve::nonlinear
         using typename Superclass::TVector;
 
         GradientDescent(const json &solver_params_,
-                        const double dt,
                         const double characteristic_length,
                         spdlog::logger &logger);
 
@@ -31,7 +30,7 @@ namespace polysolve::nonlinear
     protected:
         void reset(const int ndof) override;
 
-        virtual bool compute_update_direction(
+        virtual void compute_update_direction(
             Problem &objFunc,
             const TVector &x,
             const TVector &grad,

diff --git a/src/polysolve/nonlinear/LBFGS.cpp b/src/polysolve/nonlinear/LBFGS.cpp
@@ -5,11 +5,9 @@
 namespace polysolve::nonlinear
 {
     LBFGS::LBFGS(const json &solver_params,
-                 const double dt,
                  const double characteristic_length,
                  spdlog::logger &logger)
         : Superclass(solver_params,
-                     dt,
                      characteristic_length,
                      logger)
     {
@@ -49,7 +47,7 @@ namespace polysolve::nonlinear
         this->descent_strategy = 2;
     }
 
-    bool LBFGS::compute_update_direction(
+    void LBFGS::compute_update_direction(
         Problem &objFunc,
         const TVector &x,
         const TVector &grad,
@@ -75,7 +73,5 @@ namespace polysolve::nonlinear
 
         m_prev_x = x;
         m_prev_grad = grad;
-
-        return true;
     }
 } // namespace polysolve::nonlinear
diff --git a/src/polysolve/nonlinear/LBFGS.hpp b/src/polysolve/nonlinear/LBFGS.hpp
@@ -17,7 +17,6 @@ namespace polysolve::nonlinear
         using typename Superclass::TVector;
 
         LBFGS(const json &solver_params,
-              const double dt,
               const double characteristic_length,
               spdlog::logger &logger);
 
@@ -33,7 +32,7 @@ namespace polysolve::nonlinear
 
         void reset(const int ndof) override;
 
-        bool compute_update_direction(
+        void compute_update_direction(
             Problem &objFunc,
             const TVector &x,
             const TVector &grad,

diff --git a/src/polysolve/nonlinear/LBFGSB.cpp b/src/polysolve/nonlinear/LBFGSB.cpp
@@ -6,10 +6,9 @@
 namespace polysolve::nonlinear
 {
     LBFGSB::LBFGSB(const json &solver_params,
-                   const double dt,
                    const double characteristic_length,
                    spdlog::logger &logger)
-        : Superclass(solver_params, dt, characteristic_length, logger)
+        : Superclass(solver_params, characteristic_length, logger)
     {
         m_history_size = solver_params.value("history_size", 6);
     }
@@ -54,7 +53,7 @@ namespace polysolve::nonlinear
         this->descent_strategy = 1;
     }
 
-    bool LBFGSB::compute_update_direction(
+    void LBFGSB::compute_update_direction(
         Problem &objFunc,
         const TVector &x,
         const TVector &grad,
@@ -105,6 +104,7 @@ namespace polysolve::nonlinear
         //     m_logger.debug("direc: {}", direction.transpose());
         // }
 
+        // maybe remove me?
         if (std::isnan(direction.squaredNorm()))
         {
             reset_history(x.size());
@@ -128,7 +128,5 @@ namespace polysolve::nonlinear
 
         m_prev_x = x;
         m_prev_grad = grad;
-
-        return true;
     }
 } // namespace polysolve::nonlinear
diff --git a/src/polysolve/nonlinear/LBFGSB.hpp b/src/polysolve/nonlinear/LBFGSB.hpp
@@ -17,7 +17,6 @@ namespace polysolve::nonlinear
         using typename Superclass::TVector;
 
         LBFGSB(const json &solver_params,
-               const double dt,
                const double characteristic_length,
                spdlog::logger &logger);
 
@@ -50,7 +49,7 @@ namespace polysolve::nonlinear
 
         void reset_history(const int ndof);
 
-        virtual bool compute_update_direction(
+        virtual void compute_update_direction(
             Problem &objFunc,
             const TVector &x,
             const TVector &grad,

diff --git a/src/polysolve/nonlinear/MMA.cpp b/src/polysolve/nonlinear/MMA.cpp
@@ -5,10 +5,9 @@
 namespace polysolve::nonlinear
 {
     MMA::MMA(const json &solver_params,
-             const double dt,
              const double characteristic_length,
              spdlog::logger &logger)
-        : Superclass(solver_params, dt, characteristic_length, logger)
+        : Superclass(solver_params, characteristic_length, logger)
     {
         this->m_line_search = nullptr;
     }
@@ -36,7 +35,7 @@ namespace polysolve::nonlinear
         this->descent_strategy++;
     }
 
-    bool MMA::compute_update_direction(
+    void MMA::compute_update_direction(
         Problem &objFunc,
         const TVector &x,
         const TVector &grad,
@@ -67,6 +66,7 @@ namespace polysolve::nonlinear
         mma->Update(y.data(), grad.data(), g.data(), dg.data(), lower_bound.data(), upper_bound.data());
         direction = y - x;
 
+        // maybe remove me
         if (std::isnan(direction.squaredNorm()))
         {
             log_and_throw_error(m_logger, "nan in direction.");
@@ -76,7 +76,5 @@ namespace polysolve::nonlinear
         //     polyfem::logger().error("Direction is not a descent direction, stop.");
         //     throw std::runtime_error("Direction is not a descent direction, stop.");
         // }
-
-        return true;
     }
 } // namespace polysolve::nonlinear
diff --git a/src/polysolve/nonlinear/MMA.hpp b/src/polysolve/nonlinear/MMA.hpp
@@ -13,7 +13,6 @@ namespace polysolve::nonlinear
         using typename Superclass::TVector;
 
         MMA(const json &solver_params,
-            const double dt,
             const double characteristic_length,
             spdlog::logger &logger);
 
@@ -37,7 +36,7 @@ namespace polysolve::nonlinear
 
         void reset(const int ndof) override;
 
-        virtual bool compute_update_direction(
+        virtual void compute_update_direction(
             Problem &objFunc,
             const TVector &x,
             const TVector &grad,