[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Tue, 29 Jan 2019 23:59:32 -0500 (EST) |
branch: integration-point-variables
commit 89ee4f3edd5d6c8b089d3f74e8e9b119483f885f
Author: Konstantinos Poulios <address@hidden>
Date: Wed Jan 30 05:56:19 2019 +0100
White space changes
---
src/getfem/getfem_model_solvers.h | 238 +++++++++++++++++++-------------------
src/getfem_model_solvers.cc | 24 ++--
2 files changed, 133 insertions(+), 129 deletions(-)
diff --git a/src/getfem/getfem_model_solvers.h
b/src/getfem/getfem_model_solvers.h
index 278e7ba..6aa3b5e 100644
--- a/src/getfem/getfem_model_solvers.h
+++ b/src/getfem/getfem_model_solvers.h
@@ -233,10 +233,10 @@ namespace getfem {
// 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"); }
@@ -245,7 +245,7 @@ namespace getfem {
newton_search_with_step_control() {}
};
-
+
struct quadratic_newton_line_search : public abstract_newton_line_search {
double R0_, R1_;
@@ -408,9 +408,11 @@ namespace getfem {
/*********************************************************************/
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) {
+ 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;
@@ -436,103 +438,103 @@ namespace getfem {
// 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), coeff2 = coeff * R(1.5);
-
- 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) || res2 >= res0 * coeff2) {
- dec /= R(2); res = res2; coeff2 *= R(1.5);
- } else {
- gmm::add(gmm::scaled(dr, dec), pb.state_vector());
- break;
- }
- }
- dec *= R(2);
-
- nit++;
- coeff = std::max(R(1.05), coeff*R(0.93));
- bool near_end = (iter.get_iteration() > iter.get_maxiter()/2);
- bool cut = (alpha < R(1)) && near_end;
- if ((res > minres && nit > 4) || cut) {
- stagn++;
-
- if ((stagn > 10 && alpha > alpha0 + R(5E-2)) || cut) {
- alpha = (alpha + alpha0) / R(2);
- if (near_end) alpha = R(1);
- 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;
+ / 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), coeff2 = coeff * R(1.5);
+
+ 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) || res2 >= res0 * coeff2) {
+ dec /= R(2); res = res2; coeff2 *= R(1.5);
+ } else {
+ gmm::add(gmm::scaled(dr, dec), pb.state_vector());
+ break;
+ }
+ }
+ dec *= R(2);
+
+ nit++;
+ coeff = std::max(R(1.05), coeff*R(0.93));
+ bool near_end = (iter.get_iteration() > iter.get_maxiter()/2);
+ bool cut = (alpha < R(1)) && near_end;
+ if ((res > minres && nit > 4) || cut) {
+ stagn++;
+
+ if ((stagn > 10 && alpha > alpha0 + R(5E-2)) || cut) {
+ alpha = (alpha + alpha0) / R(2);
+ if (near_end) alpha = R(1);
+ 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 return;
+ 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 return;
}
}
}
@@ -544,9 +546,11 @@ namespace getfem {
/* ***************************************************************** */
template <typename PB>
- void classical_Newton(PB &pb, gmm::iteration &iter,
- const abstract_linear_solver<typename PB::MATRIX,
- typename PB::VECTOR> &linear_solver) {
+ void classical_Newton
+ (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;
@@ -750,10 +754,10 @@ namespace getfem {
#elif GETFEM_PARA_LEVEL > 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER
if (md.is_symmetric())
return std::make_shared
- <linear_solver_distributed_mumps_sym<MATRIX, VECTOR>>();
+ <linear_solver_distributed_mumps_sym<MATRIX, VECTOR>>();
else
- return std::make_shared
- <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
+ return std::make_shared
+ <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
#else
size_type ndof = md.nb_dof(), max3d = 15000, dim = md.leading_dimension();
# ifdef GMM_USES_MUMPS
@@ -762,24 +766,24 @@ namespace getfem {
if ((ndof<300000 && dim<=2) || (ndof<max3d && dim<=3) || (ndof<1000)) {
# ifdef GMM_USES_MUMPS
if (md.is_symmetric())
- return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>();
+ return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>();
else
- return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
+ return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
# else
return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>();
# endif
}
else {
if (md.is_coercive())
- return std::make_shared
- <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
+ return std::make_shared
+ <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
else {
if (dim <= 2)
- return std::make_shared
- <linear_solver_gmres_preconditioned_ilut<MATRIX,VECTOR>>();
- else
- return std::make_shared
- <linear_solver_gmres_preconditioned_ilu<MATRIX,VECTOR>>();
+ return std::make_shared
+ <linear_solver_gmres_preconditioned_ilut<MATRIX,VECTOR>>();
+ else
+ return std::make_shared
+ <linear_solver_gmres_preconditioned_ilu<MATRIX,VECTOR>>();
}
}
#endif
@@ -800,7 +804,7 @@ namespace getfem {
return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>();
# else
return std::make_shared
- <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
+ <linear_solver_distributed_mumps<MATRIX, VECTOR>>();
# endif
#else
GMM_ASSERT1(false, "Mumps is not interfaced");
@@ -808,16 +812,16 @@ namespace getfem {
}
else if (bgeot::casecmp(name, "cg/ildlt") == 0)
return std::make_shared
- <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
+ <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>();
else if (bgeot::casecmp(name, "gmres/ilu") == 0)
return std::make_shared
- <linear_solver_gmres_preconditioned_ilu<MATRIX, VECTOR>>();
+ <linear_solver_gmres_preconditioned_ilu<MATRIX, VECTOR>>();
else if (bgeot::casecmp(name, "gmres/ilut") == 0)
return std::make_shared
- <linear_solver_gmres_preconditioned_ilut<MATRIX, VECTOR>>();
+ <linear_solver_gmres_preconditioned_ilut<MATRIX, VECTOR>>();
else if (bgeot::casecmp(name, "gmres/ilutp") == 0)
return std::make_shared
- <linear_solver_gmres_preconditioned_ilutp<MATRIX, VECTOR>>();
+ <linear_solver_gmres_preconditioned_ilutp<MATRIX, VECTOR>>();
else if (bgeot::casecmp(name, "auto") == 0)
return default_linear_solver<MATRIX, VECTOR>(md);
else
diff --git a/src/getfem_model_solvers.cc b/src/getfem_model_solvers.cc
index 00da7d2..2a45398 100644
--- a/src/getfem_model_solvers.cc
+++ b/src/getfem_model_solvers.cc
@@ -29,7 +29,7 @@ namespace getfem {
static rmodel_plsolver_type rdefault_linear_solver(const model &md) {
return default_linear_solver<model_real_sparse_matrix,
model_real_plain_vector>(md);
- }
+ }
static cmodel_plsolver_type cdefault_linear_solver(const model &md) {
return default_linear_solver<model_complex_sparse_matrix,
@@ -49,7 +49,7 @@ namespace getfem {
max_ratio_reached = false;
}
- double default_newton_line_search::next_try(void) {
+ double default_newton_line_search::next_try() {
alpha_old = alpha; ++it;
// alpha *= 0.5;
if (alpha >= 0.4) alpha *= 0.5; else alpha *= alpha_mult;
@@ -60,7 +60,7 @@ namespace getfem {
// cout << "r = " << r << " alpha = " << alpha_old << " count_pat = " <<
count_pat << endl;
if (!max_ratio_reached && r < first_res * alpha_max_ratio) {
alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
- it_max_ratio_reached = it; max_ratio_reached = true;
+ it_max_ratio_reached = it; max_ratio_reached = true;
}
if (max_ratio_reached &&
r < r_max_ratio_reached * 0.5 &&
@@ -71,9 +71,9 @@ namespace getfem {
if (count == 0 || r < conv_r)
{ conv_r = r; conv_alpha = alpha_old; count = 1; }
if (conv_r < first_res) ++count;
-
+
if (r < first_res * alpha_min_ratio)
- { count_pat = 0; return true; }
+ { count_pat = 0; return true; }
if (count>=5 || (alpha < alpha_min && max_ratio_reached) || alpha<1e-15) {
if (conv_r < first_res * 0.99) count_pat = 0;
if (/*gmm::random() * 50. < -log(conv_alpha)-4.0 ||*/ count_pat >= 3)
@@ -91,10 +91,10 @@ namespace getfem {
/* ***************************************************************** */
template <typename MATRIX, typename VECTOR, typename PLSOLVER>
- void compute_init_values(model &md, gmm::iteration &iter,
- PLSOLVER lsolver,
- abstract_newton_line_search &ls, const MATRIX &K,
- const VECTOR &rhs) {
+ void compute_init_values(model &md, gmm::iteration &iter,
+ PLSOLVER lsolver,
+ abstract_newton_line_search &ls, const MATRIX &K,
+ const VECTOR &rhs) {
VECTOR state(md.nb_dof());
md.from_variables(state);
@@ -103,7 +103,7 @@ namespace getfem {
scalar_type dt = md.get_time_step();
scalar_type ddt = md.get_init_time_step();
scalar_type t = md.get_time();
-
+
// Solve for ddt
md.set_time_step(ddt);
gmm::iteration iter1 = iter;
@@ -148,9 +148,9 @@ namespace getfem {
else {
model_pb<MATRIX, VECTOR> mdpb(md, ls, state, rhs, K);
if (dynamic_cast<newton_search_with_step_control *>(&ls))
- Newton_with_step_control(mdpb, iter, *lsolver);
+ Newton_with_step_control(mdpb, iter, *lsolver);
else
- classical_Newton(mdpb, iter, *lsolver);
+ classical_Newton(mdpb, iter, *lsolver);
}
md.to_variables(state); // copy the state vector into the model variables
}