Consistent auto-stepping

Enhancement/suanPan.sublime-completions

@@ -6427,7 +6427,7 @@
       "trigger":"RambergOsgood"
     },
     {
-      "contents":"solver Ramm (1) [2] [3] [4] [5]\n# (1) int, unique solver tag\n# [2] double, arc length, default: 0.1\n# [3] bool string, fixed arc length switch, default: false\n# [4] double, minimum arc length, default: 0.\n# [5] double, maximum arc length, default: 0.",
+      "contents":"solver Ramm (1)\n# (1) int, unique solver tag",
       "details":"Ramm's arc-length algorithm",
       "kind":"type",
       "trigger":"Ramm"

Example/Solver/Newton.supan

@@ -18,7 +18,7 @@ set ini_step_size .1
 
 cload 1 0 100 2 2
 
-converger RelIncreDisp 1 1E-10 10 1
+converger RelIncreDisp 1 1E-10 4 1
 
 analyze
 

Example/Solver/Ramm.supan

@@ -43,21 +43,24 @@ fix 2 2 1 17
 recorder 1 hdf5 Node U 11
 
 step arclength 1 11 2 -20
-solver Ramm 1 .1 false .05 .2
+solver Ramm 1
+set ini_step_size .1
+set min_step_size .05
+set max_step_size .2
 
 criterion MinDisplacement 1 11 2 -5.5
 
-converger RelIncreDisp 1 1E-10 4 1
+converger RelIncreDisp 1 1E-10 6 1
 
 analyze
 
 # Node 11:
 # Coordinate:
 #   2.0000e+00  8.0000e+00
 # Displacement:
-#   5.9432e+00 -5.5015e+00  6.7582e-01
+#   5.9392e+00 -5.5060e+00  6.7793e-01
 # Resistance:
-#  -2.7777e-08 -7.7092e+00  2.2737e-12
+#  -2.6251e-08 -9.1406e+00  3.5243e-12
 peek node 11
 
 # save recorder 1

Example/Solver/SHALLOW.ARC.supan

@@ -56,16 +56,17 @@ set symm_mat false
 
 criterion MinDisplacement 1 7 2 -1.8
 
-converger RelIncreDisp 1 1E-4 4 1
+converger RelIncreDisp 1 1E-10 6 1
 
 analyze
 
 # Node 7:
-# 4.0223 0.8396
+# Coordinate:
+#   4.0223e+00  8.3955e-01
 # Displacement:
-# 0.0168 -1.8050 0.1679
+#   2.8077e-02 -1.8912e+00  1.5941e-01
 # Resistance:
-# -1.4197e-10 -1.8196e+02 -2.2510e-11
+#  -1.2903e-10 -4.1933e+02 -1.3188e-11
 peek node 7
 
 # save recorder 1

Solver/BFGS.cpp

@@ -36,7 +36,7 @@ int BFGS::analyze() {
     const auto max_iteration = C->get_max_iteration();
 
     // iteration counter
-    unsigned counter = 0;
+    auto counter = 0u;
 
     // lambda alias
     auto& aux_lambda = W->modify_auxiliary_lambda();
@@ -59,6 +59,8 @@ int BFGS::analyze() {
     };
 
     while(true) {
+        set_step_amplifier(sqrt(max_iteration / (counter + 1.)));
+
         // update for nodes and elements
         if(SUANPAN_SUCCESS != G->update_trial_status()) return SUANPAN_FAIL;
         // process modifiers

Solver/MPDC.cpp

@@ -35,9 +35,11 @@ int MPDC::analyze() {
     wall_clock t_clock;
 
     // iteration counter
-    unsigned counter = 0;
+    auto counter = 0u;
 
     while(true) {
+        set_step_amplifier(sqrt(max_iteration / (counter + 1.)));
+
         // update for nodes and elements
         t_clock.tic();
         if(SUANPAN_SUCCESS != G->update_trial_status()) return SUANPAN_FAIL;

Solver/Newton.cpp

@@ -36,7 +36,7 @@ int Newton::analyze() {
     const auto max_iteration = C->get_max_iteration();
 
     // iteration counter
-    unsigned counter = 0;
+    auto counter = 0u;
 
     vec samurai, pre_samurai;
 
@@ -45,6 +45,8 @@ int Newton::analyze() {
     wall_clock t_clock;
 
     while(true) {
+        set_step_amplifier(sqrt(max_iteration / (counter + 1.)));
+
         // update for nodes and elements
         t_clock.tic();
         if(SUANPAN_SUCCESS != G->update_trial_status()) return SUANPAN_FAIL;

Solver/Ramm.cpp

@@ -21,18 +21,6 @@
 #include <Domain/Factory.hpp>
 #include <Solver/Integrator/Integrator.h>
 
-void Ramm::bound_arc_length() {
-    if(min_arc_length > 0. && arc_length < min_arc_length) arc_length = min_arc_length;
-    if(max_arc_length > 0. && arc_length > max_arc_length) arc_length = max_arc_length;
-}
-
-Ramm::Ramm(const unsigned T, const double L, const bool F, const double MINL, const double MAXL)
-    : Solver(T)
-    , arc_length(L)
-    , fixed_arc_length(F)
-    , min_arc_length(std::max(0., MINL))
-    , max_arc_length(std::max(0., MAXL)) {}
-
 int Ramm::analyze() {
     auto& C = get_converger();
     auto& G = get_integrator();
@@ -47,9 +35,11 @@ int Ramm::analyze() {
     vec samurai, disp_a, disp_ref;
 
     // iteration counter
-    unsigned counter = 0;
+    auto counter = 0u;
 
     while(true) {
+        set_step_amplifier(sqrt(max_iteration / (counter + 1.)));
+
         // update for nodes and elements
         if(SUANPAN_SUCCESS != G->update_trial_status()) return SUANPAN_FAIL;
         // process modifiers
@@ -95,21 +85,9 @@ int Ramm::analyze() {
         G->erase_machine_error(samurai);
 
         // exit if converged
-        if(C->is_converged(counter)) {
-            if(!fixed_arc_length) {
-                arc_length *= sqrt(max_iteration / static_cast<double>(counter));
-                bound_arc_length();
-            }
-            return SUANPAN_SUCCESS;
-        }
+        if(C->is_converged(counter)) return SUANPAN_SUCCESS;
         // exit if maximum iteration is hit
-        if(++counter > max_iteration) {
-            if(!fixed_arc_length) {
-                arc_length *= .5;
-                bound_arc_length();
-            }
-            return SUANPAN_FAIL;
-        }
+        if(++counter > max_iteration) return SUANPAN_FAIL;
 
         // update trial displacement
         G->update_from_ninja();
@@ -124,6 +102,8 @@ int Ramm::analyze() {
     }
 }
 
+void Ramm::set_step_size(const double S) { arc_length = S; }
+
 void Ramm::print() {
     suanpan_info("A solver using Ramm's arc length method.\n");
 }

Solver/Ramm.h

@@ -31,22 +31,15 @@
 #include <Solver/Solver.h>
 
 class Ramm final : public Solver {
-    double arc_length;
-    const bool fixed_arc_length;
-    const double min_arc_length, max_arc_length;
-
-    void bound_arc_length();
+    double arc_length = 1.;
 
 public:
-    explicit Ramm(unsigned = 0, // tag
-                  double = .1,  // arc length
-                  bool = false, // if fix arc length
-                  double = 0.,  // minimum arc length
-                  double = 0.   // maximum arc length
-    );
+    using Solver::Solver;
 
     int analyze() override;
 
+    void set_step_size(double) override;
+
     void print() override;
 };
 

Solver/Solver.cpp

@@ -38,6 +38,10 @@ int Solver::initialize() {
     return SUANPAN_SUCCESS;
 }
 
+void Solver::set_step_amplifier(const double S) { step_amplifier = std::max(1., S); }
+
+double Solver::get_step_amplifier() const { return step_amplifier; }
+
 void Solver::set_converger(const shared_ptr<Converger>& C) { converger = C; }
 
 const shared_ptr<Converger>& Solver::get_converger() const { return converger; }

Solver/Solver.h

@@ -39,6 +39,8 @@ class Solver : public Tag {
     shared_ptr<Converger> converger = nullptr;
     shared_ptr<Integrator> modifier = nullptr;
 
+    double step_amplifier = 1.0;
+
 public:
     explicit Solver(unsigned = 0);
     Solver(const Solver&) = default;
@@ -51,6 +53,11 @@ class Solver : public Tag {
 
     virtual int analyze() = 0;
 
+    virtual void set_step_size(double) {}
+
+    void set_step_amplifier(double);
+    [[nodiscard]] double get_step_amplifier() const;
+
     void set_converger(const shared_ptr<Converger>&);
     [[nodiscard]] const shared_ptr<Converger>& get_converger() const;
 

Solver/SolverParser.cpp

@@ -357,29 +357,6 @@ int create_new_solver(const shared_ptr<DomainBase>& domain, istringstream& comma
 
         if(domain->insert(make_shared<BFGS>(tag, max_history))) code = 1;
     }
-    else if(is_equal(solver_type, "Ramm")) {
-        auto arc_length = .1, min_arc_length = 0., max_arc_length = 0.;
-        string fixed_arc_length = "False";
-
-        if(!command.eof() && !get_input(command, arc_length)) {
-            suanpan_error("A valid arc length is required.\n");
-            return SUANPAN_SUCCESS;
-        }
-        if(!command.eof() && !get_input(command, fixed_arc_length)) {
-            suanpan_error("A valid fixed arc length switch is required.\n");
-            return SUANPAN_SUCCESS;
-        }
-        if(!command.eof() && !get_input(command, min_arc_length)) {
-            suanpan_error("A valid arc length limit is required.\n");
-            return SUANPAN_SUCCESS;
-        }
-        if(!command.eof() && !get_input(command, max_arc_length)) {
-            suanpan_error("A valid arc length limit is required.\n");
-            return SUANPAN_SUCCESS;
-        }
-
-        if(domain->insert(make_shared<Ramm>(tag, arc_length, is_true(fixed_arc_length), min_arc_length, max_arc_length))) code = 1;
-    }
     else if(is_equal(solver_type, "FEAST") || is_equal(solver_type, "QuadraticFEAST")) {
         unsigned eigen_number;
         if(!get_input(command, eigen_number)) {
@@ -402,6 +379,7 @@ int create_new_solver(const shared_ptr<DomainBase>& domain, istringstream& comma
         if(domain->insert(make_shared<FEAST>(tag, eigen_number, centre, radius, is_equal(solver_type, "QuadraticFEAST")))) code = 1;
     }
     else if(is_equal(solver_type, "DisplacementControl") || is_equal(solver_type, "MPDC")) { if(domain->insert(make_shared<MPDC>(tag))) code = 1; }
+    else if(is_equal(solver_type, "Ramm")) { if(domain->insert(make_shared<Ramm>(tag))) code = 1; }
     else
         suanpan_error("Cannot identify the solver type.\n");
 

Step/ArcLength.cpp

@@ -70,21 +70,38 @@ int ArcLength::analyze() {
     auto& S = get_solver();
     auto& G = get_integrator();
 
-    unsigned num_iteration = 0;
+    auto num_iteration = 0u;
+
+    auto arc_length = get_ini_step_size();
+    const auto min_arc_length = get_min_step_size();
+    const auto max_arc_length = get_max_step_size();
 
     while(true) {
         if(num_iteration++ > get_max_substep()) {
             suanpan_warning("The maximum sub-step number {} reached.\n", get_max_substep());
             return SUANPAN_FAIL;
         }
+        S->set_step_size(arc_length);
         if(auto code = S->analyze(); code == SUANPAN_SUCCESS) {
             G->stage_and_commit_status();
             G->record();
+            // adjust arc length, always increase
+            if(!is_fixed_step_size()) {
+                arc_length *= S->get_step_amplifier();
+                if(max_arc_length > 0. && arc_length > max_arc_length) arc_length = max_arc_length;
+            }
             // if exit is returned, the analysis shall be terminated
             code = G->process_criterion();
             if(SUANPAN_SUCCESS != code) return code;
         }
-        else if(code == SUANPAN_FAIL) G->reset_status();
+        else if(code == SUANPAN_FAIL) {
+            G->reset_status();
+            // no way to converge
+            if(is_fixed_step_size() || suanpan::approx_equal(arc_length, min_arc_length)) return SUANPAN_FAIL;
+            // decrease arc length
+            arc_length *= .5;
+            if(min_arc_length > 0. && arc_length < min_arc_length) arc_length = min_arc_length;
+        }
         else return SUANPAN_FAIL;
     }
 }

Step/Dynamic.cpp

@@ -102,7 +102,7 @@ int Dynamic::analyze() {
 
     while(true) {
         // check if the target time point is hit
-        if(remain_time <= 1E-7) return SUANPAN_SUCCESS;
+        if(remain_time <= 1E-10) return SUANPAN_SUCCESS;
         // check if the maximum substep number is hit
         if(++num_increment > get_max_substep()) {
             suanpan_warning("The maximum sub-step number {} reached.\n", get_max_substep());
@@ -120,7 +120,7 @@ int Dynamic::analyze() {
             G->record();
             if(G->allow_to_change_time_step()) {
                 if(!is_fixed_step_size() && ++num_converged_step > 5) {
-                    step_time = std::min(get_max_step_size(), step_time * time_step_amplification);
+                    step_time = std::min(get_max_step_size(), step_time * S->get_step_amplifier());
                     num_converged_step = 0;
                 }
                 // check if time overflows

Step/Dynamic.h

@@ -33,8 +33,6 @@
 enum class IntegratorType;
 
 class Dynamic final : public Step {
-    const double time_step_amplification = 1.2;
-
     const IntegratorType analysis_type;
 
 public:

Step/Static.cpp

@@ -77,7 +77,7 @@ int Static::analyze() {
 
     while(true) {
         // check if the target time point is hit
-        if(remain_time <= 1E-7) return SUANPAN_SUCCESS;
+        if(remain_time <= 1E-10) return SUANPAN_SUCCESS;
         // check if the maximum substep number is hit
         if(++num_increment > get_max_substep()) {
             suanpan_warning("The maximum sub-step number {} reached.\n", get_max_substep());
@@ -95,7 +95,7 @@ int Static::analyze() {
             // update time left which will be used in for example criterion
             set_time_left(remain_time -= step_time);
             if(!is_fixed_step_size() && ++num_converged_step > 5) {
-                step_time = std::min(get_max_step_size(), step_time * time_step_amplification);
+                step_time = std::min(get_max_step_size(), step_time * S->get_step_amplifier());
                 num_converged_step = 0;
             }
             // check if time overflows

Step/Static.h

@@ -34,8 +34,6 @@
 #include <Step/Step.h>
 
 class Static : public Step {
-    const double time_step_amplification = 1.2;
-
 public:
     explicit Static(unsigned = 0, // tag
                     double = 1.   // step time period