From 3b0bec52e0cccbc9ee2339af8a284d100ab1fa89 Mon Sep 17 00:00:00 2001 From: Theodore Chang Date: Mon, 23 Oct 2023 11:49:46 +0200 Subject: [PATCH 1/4] Add `AICN` solver --- Example/Solver/AICN.supan | 50 +++++++++++++++++++++++++++++++++++++++ Example/Solver/MPDC.supan | 1 + Solver/Newton.cpp | 24 +++++++++++++++++-- Solver/Newton.h | 15 +++++++++++- Solver/SolverParser.cpp | 9 +++++++ Step/Static.cpp | 18 +++++++------- 6 files changed, 106 insertions(+), 11 deletions(-) create mode 100644 Example/Solver/AICN.supan diff --git a/Example/Solver/AICN.supan b/Example/Solver/AICN.supan new file mode 100644 index 000000000..1c82294d2 --- /dev/null +++ b/Example/Solver/AICN.supan @@ -0,0 +1,50 @@ +# 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 + +reset +clear +exit \ No newline at end of file diff --git a/Example/Solver/MPDC.supan b/Example/Solver/MPDC.supan index cf2885be3..84eca181c 100644 --- a/Example/Solver/MPDC.supan +++ b/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 diff --git a/Solver/Newton.cpp b/Solver/Newton.cpp index 8fad73b8a..271055737 100644 --- a/Solver/Newton.cpp +++ b/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"); +} diff --git a/Solver/Newton.h b/Solver/Newton.h index 0088d70f1..c7500b685 100644 --- a/Solver/Newton.h +++ b/Solver/Newton.h @@ -30,9 +30,11 @@ #include -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 //! @} diff --git a/Solver/SolverParser.cpp b/Solver/SolverParser.cpp index 6f3ec0def..083c5a8a0 100644 --- a/Solver/SolverParser.cpp +++ b/Solver/SolverParser.cpp @@ -347,6 +347,15 @@ int create_new_solver(const shared_ptr& domain, istringstream& comma auto code = 0; if(is_equal(solver_type, "Newton")) { if(domain->insert(make_shared(tag))) code = 1; } else if(is_equal(solver_type, "modifiedNewton") || is_equal(solver_type, "mNewton")) { if(domain->insert(make_shared(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(tag, length))) code = 1; + } else if(is_equal(solver_type, "BFGS")) { if(domain->insert(make_shared(tag))) code = 1; } else if(is_equal(solver_type, "LBFGS")) { auto max_history = 20; diff --git a/Step/Static.cpp b/Step/Static.cpp index b37e74ce5..83a0e69af 100644 --- a/Step/Static.cpp +++ b/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() : solver = make_shared(); + 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() : solver = make_shared(); + else if(flag && nullptr == std::dynamic_pointer_cast(solver)) { + suanpan_warning("Wrong solver assigned, using MPDC instead.\n"); + solver = make_shared(); } solver->set_integrator(modifier); From ab9f20b49a334717e582f5bccfc00b5d0209af13 Mon Sep 17 00:00:00 2001 From: Theodore Chang Date: Mon, 23 Oct 2023 12:46:19 +0200 Subject: [PATCH 2/4] Add coverage --- Example/Solver/AICN.supan | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Example/Solver/AICN.supan b/Example/Solver/AICN.supan index 1c82294d2..f4d4861ef 100644 --- a/Example/Solver/AICN.supan +++ b/Example/Solver/AICN.supan @@ -45,6 +45,8 @@ peek element 1 # 1.0000e+01 2.0000e+01 peek node 3 4 +peek solver 1 + reset clear exit \ No newline at end of file From 9536fa1c43fe11305302222b50410927ade48c10 Mon Sep 17 00:00:00 2001 From: Theodore Chang Date: Mon, 23 Oct 2023 14:12:35 +0200 Subject: [PATCH 3/4] Use std::dynamic_pointer_cast --- Step/ArcLength.cpp | 2 +- Step/Dynamic.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Step/ArcLength.cpp b/Step/ArcLength.cpp index 3c24b47be..c6bc6a201 100644 --- a/Step/ArcLength.cpp +++ b/Step/ArcLength.cpp @@ -41,7 +41,7 @@ int ArcLength::initialize() { modifier->set_domain(t_domain); // solver - if(nullptr != solver && !dynamic_cast(solver.get())) solver = nullptr; + if(nullptr != solver && nullptr == std::dynamic_pointer_cast(solver)) solver = nullptr; if(nullptr == solver) solver = make_shared(); solver->set_converger(tester); solver->set_integrator(modifier); diff --git a/Step/Dynamic.cpp b/Step/Dynamic.cpp index a3381bd5c..e8c0a4eb9 100644 --- a/Step/Dynamic.cpp +++ b/Step/Dynamic.cpp @@ -63,7 +63,7 @@ int Dynamic::initialize() { // solver // avoid arc length solver - if(nullptr != solver) if(dynamic_cast(solver.get())) solver = nullptr; + if(nullptr != solver) if(std::dynamic_pointer_cast(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() : solver = make_shared(); } - if(dynamic_cast(solver.get()) && dynamic_cast(modifier.get())) { + if(std::dynamic_pointer_cast(solver) && std::dynamic_pointer_cast(modifier)) { suanpan_error("BFGS solver is not supported by Lee's damping model.\n"); return SUANPAN_FAIL; } From 4414eb2aaa6039c1e536ebb4d52c88c2173dc59b Mon Sep 17 00:00:00 2001 From: Theodore Chang Date: Mon, 23 Oct 2023 14:23:57 +0200 Subject: [PATCH 4/4] Update documentation, keyword, autocomplete --- CHANGELOG.md | 1 + Enhancement/suanPan.sublime-completions | 8 +++++++- Enhancement/suanPan.sublime-syntax | 2 +- 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 884662356..2426d20b2 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Enhancement/suanPan.sublime-completions b/Enhancement/suanPan.sublime-completions index 266a0c743..6bf5764fa 100644 --- a/Enhancement/suanPan.sublime-completions +++ b/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":"", diff --git a/Enhancement/suanPan.sublime-syntax b/Enhancement/suanPan.sublime-syntax index f287a8b37..6492e03ff 100644 --- a/Enhancement/suanPan.sublime-syntax +++ b/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'