AICN cubic Newton algorithm (#165)

CHANGELOG.md

@@ -16,6 +16,7 @@
 6. `B31OS` and `EB31OS` associated transformations `B3DOSL`, `B3DOSC`; sections, `Fibre3DOS`, `Cell3DOS`; material wrappers `OS146`, `OS146S`
 7. elemental damping using Lee's model
 8. support Lode angle in CDPM2 [#163](https://github.com/TLCFEM/suanPan/pull/163)
+9. add `AICN` cubic Newton solver [#165](https://github.com/TLCFEM/suanPan/pull/165)
 
 ## version 3.1
 

Enhancement/suanPan.sublime-completions

@@ -5413,11 +5413,17 @@
       "trigger":"LAPACK"
     },
     {
-      "contents":"",
+      "contents":"solver LBFGS ${1:(1)} ${2:[2]}\n# (1) int, unique solver tag\n# [2] int, maximum number of history, default: 20",
       "details":"",
       "kind":"type",
       "trigger":"LBFGS"
     },
+    {
+      "contents":"solver AICN ${1:(1)} ${2:[2]}\n# (1) int, unique solver tag\n# [2] double, approx. length, default: 1.0",
+      "details":"AICN Solver",
+      "kind":"type",
+      "trigger":"AICN"
+    },
     {
       "contents":"",
       "details":"",

Enhancement/suanPan.sublime-syntax

@@ -32,7 +32,7 @@ contexts:
     - match: '\b(?i)(Hann|Hamming|Blackman|BlackmanNuttall|BlackmanHarris|FlatTop)\b'
       scope: string
       # other
-    - match: '\b(?i)((Abs|Rel)(Error|IncreEnergy|Residual)|ArcLength|Logic(AND|OR|XOR)|OALTS|Bathe(TwoStep|Explicit)|BFGS|Buckle|Combine|Constant|Converger|Decay|(Implicit|Explicit)Dynamic|FEAST|FixedNumber|FixedLength[2|3]D|Min(imum)?Gap[2|3]D|Max(imum)?Gap[2|3]D|Frequency|Generali(s|z)edAlpha(Explicit)|GSSSS(U|V)0|GSSSSOptimal|LeeNewmark(Full)?|Linear|(Max|Min)Displacement|MaxHistory|(Max|Min)Resistance|Modulated|MPC|MPDC|m?Newton|NZStrongMotion|Optimization|OutputType|ParticleCollision[2|3]D|LJPotential2D|Ramm|(Rayleigh|Nonviscous)?Newmark|Tchamwa|(Rel|Abs)(Incre)?(Disp|Acc)|(Finite)RigidWall(Penalty|Multiplier)|(Finite)RestitutionWall(Penalty)|Static|StrainEnergyEvolution|Tabular(Spline)|Sine|Cosine|WilsonPenzienNewmark|LeeElementalNewmark)\b'
+    - match: '\b(?i)((Abs|Rel)(Error|IncreEnergy|Residual)|ArcLength|Logic(AND|OR|XOR)|OALTS|AICN|Bathe(TwoStep|Explicit)|BFGS|Buckle|Combine|Constant|Converger|Decay|(Implicit|Explicit)Dynamic|FEAST|FixedNumber|FixedLength[2|3]D|Min(imum)?Gap[2|3]D|Max(imum)?Gap[2|3]D|Frequency|Generali(s|z)edAlpha(Explicit)|GSSSS(U|V)0|GSSSSOptimal|LeeNewmark(Full)?|Linear|(Max|Min)Displacement|MaxHistory|(Max|Min)Resistance|Modulated|MPC|MPDC|m?Newton|NZStrongMotion|Optimization|OutputType|ParticleCollision[2|3]D|LJPotential2D|Ramm|(Rayleigh|Nonviscous)?Newmark|Tchamwa|(Rel|Abs)(Incre)?(Disp|Acc)|(Finite)RigidWall(Penalty|Multiplier)|(Finite)RestitutionWall(Penalty)|Static|StrainEnergyEvolution|Tabular(Spline)|Sine|Cosine|WilsonPenzienNewmark|LeeElementalNewmark)\b'
       scope: string
       # output type
     - match: '\b(?i)(A[1-6]?|AR[1-3]?|AT|AXIAL|CSE|DAMAGE|DC|DF[1-6]?|DM[1-3]?|DT|E|E11|E12|E13|E22|E23|E33|ED|EE|EE11|EE12|EE13|EE22|EE23|EE33|EEEQ|EEP|EEP1|EEP2|EEP3|EEQ|EINT|EP[1-3]?|ES|HIST|HYDRO|IF[1-6]?|IM[1-3]{1}|K|KAPPAC|KAPPAP|KAPPAT|KE|LITR|M|MISES|MOMENT|MOMENTUM((R?)[XYZ])?|NL|NMISES|PE|PE11|PE12|PE13|PE22|PE23|PE33|PEEQ|PEP[1-3]?|PP|REBARE|REBARS|RESULTANT|RF[1-6]?|RM[1-3]?|RT|S|S11|S12|S13|S22|S23|S33|SD|SE|SHEAR|SINT|SINV|SP[1-3]?|SS|TORSION|TRESC|TSE|U[1-6]?|UR[1-3]?|UT|V[1-6]?|VD|VF|VR[1-3]?|VS|VT|GDF|YF|VE|BEAM[SE])\b'

Example/Solver/AICN.supan

@@ -0,0 +1,52 @@
+# A TEST MODEL FOR AICN SOLVER
+
+node 1 0 0
+node 2 1 0
+node 3 1 1
+node 4 0 1
+
+material BilinearJ2 2 5000 .2 4 .05
+material PlaneStress 1 2
+
+element CP4 1 1 2 3 4 1 1
+
+fix 1 1 1
+fix 2 2 1 2
+
+cload 1 0 10 1 3 4
+cload 2 0 20 2 3 4
+
+step static 1 1
+solver AICN 1 0.2
+set fixed_step_size 1
+set ini_step_size 1E-1
+set symm_mat 0
+
+converger RelIncreDisp 1 1E-11 20 1
+
+analyze
+
+peek element 1
+
+# Node 3:
+# Coordinate:
+#   1.0000e+00  1.0000e+00
+# Displacement:
+#   3.2825e-01 -1.1472e-02
+# Resistance:
+#   1.0000e+01  2.0000e+01
+# 
+# Node 4:
+# Coordinate:
+#   0.0000e+00  1.0000e+00
+# Displacement:
+#   4.2981e-01  2.7330e-01
+# Resistance:
+#   1.0000e+01  2.0000e+01
+peek node 3 4
+
+peek solver 1
+
+reset
+clear
+exit

Example/Solver/MPDC.supan

@@ -14,6 +14,7 @@ fix2 1 P 1
 displacement 1 0 1 2 3 ! wrong node
 
 step static 1 1
+solver Newton 1
 set fixed_step_size 1
 
 converger RelIncreDisp 1 1E-12 10 1

Solver/Newton.cpp

@@ -85,11 +85,13 @@ int Newton::analyze() {
         // call solver
         t_clock.tic();
 
-        auto flag = G->solve(samurai, G->get_force_residual());
+        const auto residual = G->get_force_residual();
+
+        auto flag = G->solve(samurai, residual);
 
         suanpan_assert([&] {
             if(!samurai.is_finite()) {
-                suanpan_fatal("Infinite number detected.\n");
+                suanpan_fatal("Non-finite number detected.\n");
                 flag = SUANPAN_FAIL;
             }
         });
@@ -100,6 +102,8 @@ int Newton::analyze() {
             return flag;
         }
 
+        if(const auto amp = amplification(samurai, residual); amp > 0.) samurai *= amp;
+
         // deal with mpc
         if(const auto n_size = W->get_size(); 0 != W->get_mpc()) {
             auto& border = W->get_auxiliary_stiffness();
@@ -152,3 +156,19 @@ int Newton::analyze() {
 void Newton::print() {
     suanpan_info("A solver based on Newton-Raphson method{}", initial_stiffness ? " using initial stiffness for each substep.\n" : ".\n");
 }
+
+double AICN::amplification(const vec& x, const vec& r) const {
+    const auto hessian_norm = dot(x, r);
+    if(hessian_norm <= 0.) return 0.;
+    const auto root_norm = l_est * std::sqrt(hessian_norm);
+    const auto amp = (std::sqrt(1. + 2. * root_norm) - 1.) / root_norm;
+    return std::isfinite(amp) ? amp : 0.;
+}
+
+AICN::AICN(const unsigned T, const double L)
+    : Newton(T)
+    , l_est(fabs(L)) {}
+
+void AICN::print() {
+    suanpan_info("A solver based on cubic Newton method. arXiv:2211.00140\n");
+}

Solver/Newton.h

@@ -30,9 +30,11 @@
 
 #include <Solver/Solver.h>
 
-class Newton final : public Solver {
+class Newton : public Solver {
     const bool initial_stiffness;
 
+    [[nodiscard]] virtual double amplification(const vec&, const vec&) const { return 0.; }
+
 public:
     explicit Newton(unsigned = 0, bool = false);
 
@@ -41,6 +43,17 @@ class Newton final : public Solver {
     void print() override;
 };
 
+class AICN : public Newton {
+    const double l_est;
+
+    [[nodiscard]] double amplification(const vec&, const vec&) const override;
+
+public:
+    explicit AICN(unsigned = 0, double = 1.);
+
+    void print() override;
+};
+
 #endif
 
 //! @}

Solver/SolverParser.cpp

@@ -347,6 +347,15 @@ int create_new_solver(const shared_ptr<DomainBase>& domain, istringstream& comma
     auto code = 0;
     if(is_equal(solver_type, "Newton")) { if(domain->insert(make_shared<Newton>(tag))) code = 1; }
     else if(is_equal(solver_type, "modifiedNewton") || is_equal(solver_type, "mNewton")) { if(domain->insert(make_shared<Newton>(tag, true))) code = 1; }
+    else if(is_equal(solver_type, "AICN")) {
+        auto length = 1.;
+        if(!command.eof() && !get_input(command, length)) {
+            suanpan_error("A valid length is required.\n");
+            return SUANPAN_SUCCESS;
+        }
+
+        if(domain->insert(make_shared<AICN>(tag, length))) code = 1;
+    }
     else if(is_equal(solver_type, "BFGS")) { if(domain->insert(make_shared<BFGS>(tag))) code = 1; }
     else if(is_equal(solver_type, "LBFGS")) {
         auto max_history = 20;

Step/ArcLength.cpp

@@ -41,7 +41,7 @@ int ArcLength::initialize() {
     modifier->set_domain(t_domain);
 
     // solver
-    if(nullptr != solver && !dynamic_cast<Ramm*>(solver.get())) solver = nullptr;
+    if(nullptr != solver && nullptr == std::dynamic_pointer_cast<Ramm>(solver)) solver = nullptr;
     if(nullptr == solver) solver = make_shared<Ramm>();
     solver->set_converger(tester);
     solver->set_integrator(modifier);

Step/Dynamic.cpp

@@ -63,7 +63,7 @@ int Dynamic::initialize() {
 
     // solver
     // avoid arc length solver
-    if(nullptr != solver) if(dynamic_cast<Ramm*>(solver.get())) solver = nullptr;
+    if(nullptr != solver) if(std::dynamic_pointer_cast<Ramm>(solver)) solver = nullptr;
     // automatically enable displacement controlled solver
     if(nullptr == solver) {
         auto flag = false;
@@ -75,7 +75,7 @@ int Dynamic::initialize() {
         flag ? solver = make_shared<MPDC>() : solver = make_shared<Newton>();
     }
 
-    if(dynamic_cast<BFGS*>(solver.get()) && dynamic_cast<LeeNewmarkBase*>(modifier.get())) {
+    if(std::dynamic_pointer_cast<BFGS>(solver) && std::dynamic_pointer_cast<LeeNewmarkBase>(modifier)) {
         suanpan_error("BFGS solver is not supported by Lee's damping model.\n");
         return SUANPAN_FAIL;
     }

Step/Static.cpp

@@ -46,14 +46,16 @@ int Static::initialize() {
 
     // solver
     // automatically enable displacement controlled solver
-    if(nullptr == solver) {
-        auto flag = false;
-        for(const auto& I : t_domain->get_load_pool())
-            if(I->if_displacement_control() && I->get_start_step() == get_tag()) {
-                flag = true;
-                break;
-            }
-        flag ? solver = make_shared<MPDC>() : solver = make_shared<Newton>();
+    auto flag = false;
+    for(const auto& I : t_domain->get_load_pool())
+        if(I->if_displacement_control() && I->get_start_step() == get_tag()) {
+            flag = true;
+            break;
+        }
+    if(nullptr == solver) flag ? solver = make_shared<MPDC>() : solver = make_shared<Newton>();
+    else if(flag && nullptr == std::dynamic_pointer_cast<MPDC>(solver)) {
+        suanpan_warning("Wrong solver assigned, using MPDC instead.\n");
+        solver = make_shared<MPDC>();
     }
 
     solver->set_integrator(modifier);