[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Thu, 12 Oct 2017 07:36:12 -0400 (EDT) |
branch: devel-yves-ctrl-newton
commit d8c94f49695de6572103b63b02b92d37931e5622
Author: Yves Renard <address@hidden>
Date: Thu Oct 12 13:35:50 2017 +0200
A attempt for a Newton algorithm with step control : work in progress
---
interface/src/gf_model_get.cc | 3 +-
.../tests/matlab/demo_Nitsche_contact_by_hand.m | 5 +-
interface/tests/matlab/demo_plasticity.m | 4 +-
src/getfem/getfem_model_solvers.h | 179 ++++++++++++++++++++-
src/getfem_model_solvers.cc | 11 +-
src/gmm/gmm_blas.h | 64 ++++++++
tests/nonlinear_elastostatic.param | 2 +-
tests/plasticity.cc | 5 +-
8 files changed, 256 insertions(+), 17 deletions(-)
diff --git a/interface/src/gf_model_get.cc b/interface/src/gf_model_get.cc
index ba0d365..a90cdb8 100644
--- a/interface/src/gf_model_get.cc
+++ b/interface/src/gf_model_get.cc
@@ -476,7 +476,8 @@ void gf_model_get(getfemint::mexargs_in& m_in,
if (alpha_mult < scalar_type(0))
alpha_mult = 3.0/5.0;
- getfem::default_newton_line_search default_ls;
+ // getfem::default_newton_line_search default_ls;
+ getfem::newton_search_with_step_control default_ls;
getfem::simplest_newton_line_search simplest_ls(size_type(-1),
alpha_max_ratio,
alpha_min, alpha_mult,
alpha_threshold_res);
getfem::systematic_newton_line_search systematic_ls(size_type(-1),
alpha_min, alpha_mult);
diff --git a/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
b/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
index 8075a97..0572754 100644
--- a/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
+++ b/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
@@ -24,7 +24,7 @@ if (asize(1)) is_automatic = true; else is_automatic = false;
end
gf_workspace('clear all');
clear all;
-approximation_type = 2 % 0 : Augmentend Lagrangian
+approximation_type = 0 % 0 : Augmentend Lagrangian
% 1 : Nitsche (biased)
% 2 : Unbiased Nitsche method
@@ -240,7 +240,8 @@ NX = Nxy(zz)
max_res = 1E-8;
max_iter = 100;
- gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', niter, 'noisy',
'lsearch', 'simplest', 'alpha min', 0.8, 'lsolver', 'mumps');
+ % gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', niter, 'noisy',
'lsearch', 'simplest', 'alpha min', 0.8, 'lsolver', 'mumps');
+ gf_model_get(md, 'solve', 'max_res', 1E-9, 'max_iter', niter, 'noisy',
'lsolver', 'mumps');
U1 = gf_model_get(md, 'variable', 'u1');
UU1 = gf_model_get(md, 'variable', 'u1');
diff --git a/interface/tests/matlab/demo_plasticity.m
b/interface/tests/matlab/demo_plasticity.m
index 78ff3d1..e23e9fd 100644
--- a/interface/tests/matlab/demo_plasticity.m
+++ b/interface/tests/matlab/demo_plasticity.m
@@ -187,8 +187,8 @@ for step=1:size(t,2),
end;
% Solve the system
- get(md, 'solve', 'noisy', 'lsearch', 'simplest', 'alpha min', 0.8,
'max_iter', 100, 'max_res', 1e-6);
- % get(md, 'solve', 'noisy', 'max_iter', 80);
+ % get(md, 'solve', 'noisy', 'lsearch', 'simplest', 'alpha min', 0.8,
'max_iter', 100, 'max_res', 1e-6);
+ get(md, 'solve', 'noisy', 'max_iter', 80);
% Retrieve the solution U
U = get(md, 'variable', 'u');
diff --git a/src/getfem/getfem_model_solvers.h
b/src/getfem/getfem_model_solvers.h
index 6198072..c39f6d7 100644
--- a/src/getfem/getfem_model_solvers.h
+++ b/src/getfem/getfem_model_solvers.h
@@ -231,6 +231,21 @@ namespace getfem {
virtual ~abstract_newton_line_search() { }
};
+ // Dummy linesearch for the newton with step control
+ struct newton_search_with_step_control : public abstract_newton_line_search {
+
+ virtual void init_search(double /*r*/, size_t /*git*/, double /*R0*/ = 0.0)
+ { GMM_ASSERT1(false, "Not to be used"); }
+
+ virtual double next_try(void)
+ { GMM_ASSERT1(false, "Not to be used"); }
+
+ virtual bool is_converged(double /*r*/, double /*R1*/ = 0.0)
+ { GMM_ASSERT1(false, "Not to be used"); }
+
+ newton_search_with_step_control() {}
+ };
+
struct quadratic_newton_line_search : public abstract_newton_line_search {
double R0_, R1_;
@@ -368,11 +383,161 @@ namespace getfem {
};
+ /* ***************************************************************** */
+ /* Newton(-Raphson) algorithm with step control. */
+ /* ***************************************************************** */
+ /* Still solves a problem F(X) = 0 sarting at X0 but setting */
+ /* B0 = F(X0) and solving */
+ /* F(X) = (1-alpha)B0 (1) */
+ /* with alpha growing from 0 to 1. */
+ /* A very basic line search is applied. */
+ /* */
+ /* Step 0 : set alpha0 = 0, alpha = 1, X0 given and B0 = F(X0). */
+ /* Step 1 : Set Ri = F(Xi) - (1-alpha)B0 */
+ /* If ||Ri|| < rho, Xi+1 = Xi and go to step 2 */
+ /* Perform Newton step on problem (1) */
+ /* If the Newton do not converge (stagnation) */
+ /* alpha <- max(alpha0+1E-4, (alpha+alpha0)/2) */
+ /* Loop on step 1 with Xi unchanged */
+ /* Step 2 : if alpha=1 stop */
+ /* else alpha0 <- alpha, */
+ /* alpha <- min(1,alpha0+2*(alpha-alpha0)), */
+ /* Go to step 1 with Xi+1 */
+ /* being the result of Newton iteraitons of step1. */
+ /* */
+ /*********************************************************************/
+
+ template <typename PB>
+ void Newton_with_step_control(PB &pb, gmm::iteration &iter,
+ const abstract_linear_solver<typename PB::MATRIX,
+ typename PB::VECTOR> &linear_solver) {
+ typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T;
+ typedef typename gmm::number_traits<T>::magnitude_type R;
+ gmm::iteration iter_linsolv0 = iter;
+ iter_linsolv0.reduce_noisy();
+ iter_linsolv0.set_resmax(iter.get_resmax()/20.0);
+ iter_linsolv0.set_maxiter(10000); // arbitrary
+
+ pb.compute_residual();
+ R approx_eln = pb.approx_external_load_norm();
+ if (approx_eln == R(0)) approx_eln = R(1);
+
+ typename PB::VECTOR b0(gmm::vect_size(pb.residual()));
+ gmm::copy(pb.residual(), b0);
+ typename PB::VECTOR Xi(gmm::vect_size(pb.residual()));
+ gmm::copy(pb.state_vector(), Xi);
+
+ typename PB::VECTOR dr(gmm::vect_size(pb.residual()));
+ typename PB::VECTOR b(gmm::vect_size(pb.residual()));
+
+ R alpha0(0), alpha(1), res0(gmm::vect_norm1(b0)), minres(res0);
+ // const newton_search_with_step_control *ls
+ // = dynamic_cast<const newton_search_with_step_control *>(&(pb.ls));
+ // GMM_ASSERT1(ls, "Internal error");
+ size_type nit = 0, stagn = 0;
+ R coeff = R(2);
+
+ scalar_type crit = pb.residual_norm() / approx_eln;
+ if (iter.finished(crit)) return;
+ for(;;) {
+
+ crit = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha))
+ / approx_eln;
+ if (!iter.converged(crit)) {
+ gmm::iteration iter_linsolv = iter_linsolv0;
+ if (iter.get_noisy() > 1)
+ cout << "starting tangent matrix computation" << endl;
+
+ int is_singular = 1;
+ while (is_singular) { // Linear system solve
+ pb.compute_tangent_matrix();
+ gmm::clear(dr);
+ gmm::copy(gmm::scaled(pb.residual(), pb.scale_residual()), b);
+ gmm::add(gmm::scaled(b0,alpha-R(1)), b);
+ if (iter.get_noisy() > 1) cout << "starting linear solver" << endl;
+ iter_linsolv.init();
+ linear_solver(pb.tangent_matrix(), dr, b, iter_linsolv);
+ if (!iter_linsolv.converged()) {
+ is_singular++;
+ if (is_singular <= 4) {
+ if (iter.get_noisy())
+ cout << "Singular tangent matrix:"
+ " perturbation of the state vector." << endl;
+ pb.perturbation();
+ pb.compute_residual();
+ } else {
+ if (iter.get_noisy())
+ cout << "Singular tangent matrix: perturbation failed, "
+ << "aborting." << endl;
+ return;
+ }
+ }
+ else is_singular = 0;
+ }
+ if (iter.get_noisy() > 1) cout << "linear solver done" << endl;
+
+
+ gmm::add(dr, pb.state_vector());
+ pb.compute_residual();
+ R res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
+ R dec = R(1)/R(2);
+
+ while (dec > R(1E-5) && res >= res0 * coeff) {
+ gmm::add(gmm::scaled(dr, -dec), pb.state_vector());
+ pb.compute_residual();
+ R res2 = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
+ if (res2 < res*R(0.95)) {
+ dec /= R(2); res = res2;
+ } else {
+ gmm::add(gmm::scaled(dr, dec), pb.state_vector());
+ break;
+ }
+ }
+ dec *= R(2);
+
+ nit++;
+ coeff = std::max(R(1), coeff*R(0.93));
+ if (res > minres && nit > 4) {
+ stagn++;
+
+ if (stagn > 10 && alpha > alpha0 + 5E-2) {
+ alpha = (alpha + alpha0) / R(2);
+ gmm::copy(Xi, pb.state_vector());
+ pb.compute_residual();
+ res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha));
+ nit = 0;
+ stagn = 0; coeff = R(2);
+ }
+ }
+ if (res < minres || (alpha == R(1) && nit < 10)) stagn = 0;
+ res0 = res;
+ if (nit < 5) minres = res0; else minres = std::min(minres, res0);
+
+ if (iter.get_noisy())
+ cout << "step control [" << std::setw(8) << alpha0 << ","
+ << std::setw(8) << alpha << "," << std::setw(10) << dec << "]";
+ ++iter;
+ // crit = std::min(res / approx_eln,
+ // gmm::vect_norm1(dr) / std::max(1E-25,pb.state_norm()));
+ crit = res / approx_eln;
+ }
+
+ if (iter.finished(crit)) {
+ if (iter.converged() && alpha < R(1)) {
+ R a = alpha;
+ alpha = std::min(R(1), alpha*R(3) - alpha0*R(2));
+ alpha0 = a;
+ gmm::copy(pb.state_vector(), Xi);
+ nit = 0; stagn = 0; coeff = R(2);
+ } else if (alpha == 1) return;
+ }
+ }
+ }
/* ***************************************************************** */
- /* Newton algorithm. */
+ /* Classicel Newton(-Raphson) algorithm. */
/* ***************************************************************** */
template <typename PB>
@@ -387,12 +552,13 @@ namespace getfem {
iter_linsolv0.set_maxiter(10000); // arbitrary
pb.compute_residual();
+ R approx_eln = pb.approx_external_load_norm();
+ if (approx_eln == R(0)) approx_eln = R(1);
typename PB::VECTOR dr(gmm::vect_size(pb.residual()));
typename PB::VECTOR b(gmm::vect_size(pb.residual()));
- scalar_type crit = pb.residual_norm()
- / std::max(1E-25, pb.approx_external_load_norm());
+ scalar_type crit = pb.residual_norm() / approx_eln;
while (!iter.finished(crit)) {
gmm::iteration iter_linsolv = iter_linsolv0;
if (iter.get_noisy() > 1)
@@ -429,8 +595,7 @@ namespace getfem {
//executes a pb.compute_residual();
if (iter.get_noisy()) cout << "alpha = " << std::setw(6) << alpha << " ";
++iter;
- crit = std::min(pb.residual_norm()
- / std::max(1E-25, pb.approx_external_load_norm()),
+ crit = std::min(pb.residual_norm() / approx_eln,
gmm::vect_norm1(dr) / std::max(1E-25, pb.state_norm()));
}
}
@@ -477,7 +642,9 @@ namespace getfem {
gmm::add(gmm::scaled(V, ampl), state);
}
- const VECTOR &residual(void) { return rhs; }
+ const VECTOR &residual(void) const { return rhs; }
+ const VECTOR &state_vector(void) const { return state; }
+ VECTOR &state_vector(void) { return state; }
R state_norm(void) const
{ return gmm::vect_norm1(state); }
diff --git a/src/getfem_model_solvers.cc b/src/getfem_model_solvers.cc
index dae6bf5..00da7d2 100644
--- a/src/getfem_model_solvers.cc
+++ b/src/getfem_model_solvers.cc
@@ -147,7 +147,10 @@ namespace getfem {
}
else {
model_pb<MATRIX, VECTOR> mdpb(md, ls, state, rhs, K);
- classical_Newton(mdpb, iter, *lsolver);
+ if (dynamic_cast<newton_search_with_step_control *>(&ls))
+ Newton_with_step_control(mdpb, iter, *lsolver);
+ else
+ classical_Newton(mdpb, iter, *lsolver);
}
md.to_variables(state); // copy the state vector into the model variables
}
@@ -175,12 +178,14 @@ namespace getfem {
void standard_solve(model &md, gmm::iteration &iter,
cmodel_plsolver_type lsolver) {
- default_newton_line_search ls;
+ newton_search_with_step_control ls;
+ // default_newton_line_search ls;
standard_solve(md, iter, lsolver, ls);
}
void standard_solve(model &md, gmm::iteration &iter) {
- getfem::default_newton_line_search ls;
+ newton_search_with_step_control ls;
+ // default_newton_line_search ls;
if (md.is_complex())
standard_solve(md, iter, cdefault_linear_solver(md), ls);
else
diff --git a/src/gmm/gmm_blas.h b/src/gmm/gmm_blas.h
index b237355..e227ac8 100644
--- a/src/gmm/gmm_blas.h
+++ b/src/gmm/gmm_blas.h
@@ -651,6 +651,38 @@ namespace gmm {
return res;
}
+ /** 1-distance between two vectors */
+ template <typename V1, typename V2> inline
+ typename number_traits<typename linalg_traits<V1>::value_type>
+ ::magnitude_type
+ vect_dist1(const V1 &v1, const V2 &v2) { // not fully optimized
+ typedef typename linalg_traits<V1>::value_type T;
+ typedef typename number_traits<T>::magnitude_type R;
+ auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
+ auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
+ size_type k1(0), k2(0);
+ R res(0);
+ while (it1 != ite1 && it2 != ite2) {
+ size_type i1 = index_of_it(it1, k1,
+ typename linalg_traits<V1>::storage_type());
+ size_type i2 = index_of_it(it2, k2,
+ typename linalg_traits<V2>::storage_type());
+
+ if (i1 == i2) {
+ res += gmm::abs(*it2 - *it1); ++it1; ++k1; ++it2; ++k2;
+ }
+ else if (i1 < i2) {
+ res += gmm::abs(*it1); ++it1; ++k1;
+ }
+ else {
+ res += gmm::abs(*it2); ++it2; ++k2;
+ }
+ }
+ while (it1 != ite1) { res += gmm::abs(*it1); ++it1; }
+ while (it2 != ite2) { res += gmm::abs(*it2); ++it2; }
+ return res;
+ }
+
/* ******************************************************************** */
/* vector Infinity norm */
/* ******************************************************************** */
@@ -666,6 +698,38 @@ namespace gmm {
return res;
}
+ /** Infinity distance between two vectors */
+ template <typename V1, typename V2> inline
+ typename number_traits<typename linalg_traits<V1>::value_type>
+ ::magnitude_type
+ vect_distinf(const V1 &v1, const V2 &v2) { // not fully optimized
+ typedef typename linalg_traits<V1>::value_type T;
+ typedef typename number_traits<T>::magnitude_type R;
+ auto it1 = vect_const_begin(v1), ite1 = vect_const_end(v1);
+ auto it2 = vect_const_begin(v2), ite2 = vect_const_end(v2);
+ size_type k1(0), k2(0);
+ R res(0);
+ while (it1 != ite1 && it2 != ite2) {
+ size_type i1 = index_of_it(it1, k1,
+ typename linalg_traits<V1>::storage_type());
+ size_type i2 = index_of_it(it2, k2,
+ typename linalg_traits<V2>::storage_type());
+
+ if (i1 == i2) {
+ res = std::max(res, gmm::abs(*it2 - *it1)); ++it1; ++k1; ++it2; ++k2;
+ }
+ else if (i1 < i2) {
+ res = std::max(res, gmm::abs(*it1)); ++it1; ++k1;
+ }
+ else {
+ res = std::max(res, gmm::abs(*it2)); ++it2; ++k2;
+ }
+ }
+ while (it1 != ite1) { res = std::max(res, gmm::abs(*it1)); ++it1; }
+ while (it2 != ite2) { res = std::max(res, gmm::abs(*it2)); ++it2; }
+ return res;
+ }
+
/* ******************************************************************** */
/* matrix norm1 */
/* ******************************************************************** */
diff --git a/tests/nonlinear_elastostatic.param
b/tests/nonlinear_elastostatic.param
index ad1fb31..9914d71 100644
--- a/tests/nonlinear_elastostatic.param
+++ b/tests/nonlinear_elastostatic.param
@@ -81,7 +81,7 @@ INTEGRATION = 'IM_TETRAHEDRON(6)'
%INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(7), 3)';
RESIDUAL = 1E-6; % residu for iterative solvers.
-MAXITER = 100;
+MAXITER = 100000;
DIRICHLET_VERSION=0;
NBSTEP = 40;
%%%%% saving parameters %%%%%
diff --git a/tests/plasticity.cc b/tests/plasticity.cc
index fe7d16f..99eb31b 100644
--- a/tests/plasticity.cc
+++ b/tests/plasticity.cc
@@ -327,8 +327,9 @@ bool elastoplasticity_problem::solve(plain_vector &U) {
// Generic solve.
cout << "Number of variables : "
<< model.nb_dof() << endl;
-
- getfem::simplest_newton_line_search ls;
+
+ getfem::newton_search_with_step_control ls;
+ // getfem::simplest_newton_line_search ls;
gmm::iteration iter(residual, 2, 40000);
getfem::standard_solve(model, iter,
getfem::rselect_linear_solver(model, "superlu"), ls);