[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch devel-logari81-internal-variabl
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] [getfem-commits] branch devel-logari81-internal-variables updated: Refactoring of model solvers |
|
Date: |
Wed, 29 Jan 2020 00:41:21 -0500 |
This is an automated email from the git hooks/post-receive script.
logari81 pushed a commit to branch devel-logari81-internal-variables
in repository getfem.
The following commit(s) were added to
refs/heads/devel-logari81-internal-variables by this push:
new f85112c Refactoring of model solvers
f85112c is described below
commit f85112cc3c7947eab82ed34cac44b751474852ed
Author: Konstantinos Poulios <address@hidden>
AuthorDate: Wed Jan 29 06:41:07 2020 +0100
Refactoring of model solvers
- make the model_pb object responsible for the solution of the linearized
system
---
src/getfem/getfem_model_solvers.h | 181 ++++++--------------------------
src/getfem/getfem_models.h | 18 ++++
src/getfem_model_solvers.cc | 213 +++++++++++++++++++++++++++++++++-----
3 files changed, 236 insertions(+), 176 deletions(-)
diff --git a/src/getfem/getfem_model_solvers.h
b/src/getfem/getfem_model_solvers.h
index 517b3d3..b286116 100644
--- a/src/getfem/getfem_model_solvers.h
+++ b/src/getfem/getfem_model_solvers.h
@@ -79,8 +79,10 @@ namespace getfem {
template <typename MAT, typename VECT>
struct abstract_linear_solver {
+ typedef MAT MATRIX;
+ typedef VECT VECTOR;
virtual void operator ()(const MAT &, VECT &, const VECT &,
- gmm::iteration &) const = 0;
+ gmm::iteration &) const = 0;
virtual ~abstract_linear_solver() {}
};
@@ -386,7 +388,7 @@ namespace getfem {
/* ***************************************************************** */
/* Newton(-Raphson) algorithm with step control. */
/* ***************************************************************** */
- /* Still solves a problem F(X) = 0 sarting at X0 but setting */
+ /* Still solves a problem F(X) = 0 starting at X0 but setting */
/* B0 = F(X0) and solving */
/* F(X) = (1-alpha)B0 (1) */
/* with alpha growing from 0 to 1. */
@@ -403,15 +405,12 @@ namespace getfem {
/* 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. */
+ /* being the result of Newton iterations 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)
+ void Newton_with_step_control(PB &pb, gmm::iteration &iter)
{
typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T;
typedef typename gmm::number_traits<T>::magnitude_type R;
@@ -430,7 +429,6 @@ namespace getfem {
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
@@ -447,18 +445,18 @@ namespace getfem {
/ 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(pb.residual(), b);
- gmm::add(gmm::scaled(b0,alpha-R(1)), b);
- if (iter.get_noisy() > 1) cout << "starting linear solver" << endl;
+ pb.add_to_residual(b0, alpha-R(1)); // canceled at next
compute_residual
iter_linsolv.init();
- linear_solver(pb.tangent_matrix(), dr, b, iter_linsolv);
+ if (iter.get_noisy() > 1)
+ cout << "starting tangent matrix computation" << endl;
+ pb.compute_tangent_matrix();
+ if (iter.get_noisy() > 1)
+ cout << "starting linear solver" << endl;
+ pb.linear_solve(dr, iter_linsolv);
if (!iter_linsolv.converged()) {
is_singular++;
if (is_singular <= 4) {
@@ -466,7 +464,7 @@ namespace getfem {
cout << "Singular tangent matrix:"
" perturbation of the state vector." << endl;
pb.perturbation();
- pb.compute_residual();
+ pb.compute_residual(); // cancels add_to_residual
} else {
if (iter.get_noisy())
cout << "Singular tangent matrix: perturbation failed, "
@@ -478,9 +476,8 @@ namespace getfem {
}
if (iter.get_noisy() > 1) cout << "linear solver done" << endl;
-
gmm::add(dr, pb.state_vector());
- pb.compute_residual();
+ pb.compute_residual(); // cancels add_to_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);
@@ -542,14 +539,11 @@ namespace getfem {
/* ***************************************************************** */
- /* Classicel Newton(-Raphson) algorithm. */
+ /* Classical Newton(-Raphson) algorithm. */
/* ***************************************************************** */
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)
{
typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T;
typedef typename gmm::number_traits<T>::magnitude_type R;
@@ -563,22 +557,21 @@ namespace getfem {
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() / approx_eln;
while (!iter.finished(crit)) {
gmm::iteration iter_linsolv = iter_linsolv0;
- if (iter.get_noisy() > 1)
- cout << "starting computing tangent matrix" << endl;
int is_singular = 1;
while (is_singular) {
- pb.compute_tangent_matrix();
gmm::clear(dr);
- gmm::copy(pb.residual(), 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.get_noisy() > 1)
+ cout << "starting computing tangent matrix" << endl;
+ pb.compute_tangent_matrix();
+ if (iter.get_noisy() > 1)
+ cout << "starting linear solver" << endl;
+ pb.linear_solve(dr, iter_linsolv);
if (!iter_linsolv.converged()) {
is_singular++;
if (is_singular <= 4) {
@@ -608,133 +601,17 @@ namespace getfem {
}
- /* ***************************************************************** */
- /* Intermediary structure for Newton algorithms with getfem::model. */
- /* ***************************************************************** */
-
- #define TRACE_SOL 0
-
- template <typename MAT, typename VEC>
- struct model_pb {
-
- typedef MAT MATRIX;
- typedef VEC VECTOR;
- typedef typename gmm::linalg_traits<VECTOR>::value_type T;
- typedef typename gmm::number_traits<T>::magnitude_type R;
-
- model &md;
- abstract_newton_line_search &ls;
- VECTOR stateinit, &state;
- const VECTOR &rhs;
- const MATRIX &K;
-
- void compute_tangent_matrix() {
- md.to_variables(state);
- md.assembly(model::BUILD_MATRIX);
- }
-
- const MATRIX &tangent_matrix() { return K; }
-
- void compute_residual() {
- md.to_variables(state);
- md.assembly(model::BUILD_RHS);
- }
-
- void perturbation() {
- R res = gmm::vect_norm2(state), ampl = std::max(res*R(1E-20), R(1E-50));
- std::vector<R> V(gmm::vect_size(state));
- gmm::fill_random(V);
- gmm::add(gmm::scaled(V, ampl), state);
- }
-
- const VECTOR &residual() const { return rhs; }
- const VECTOR &state_vector() const { return state; }
- VECTOR &state_vector() { return state; }
-
- R state_norm() const { return gmm::vect_norm1(state); }
-
- R approx_external_load_norm() { return md.approx_external_load(); }
-
- R residual_norm() {
- return gmm::vect_norm1(rhs); // A norm1 seems to be better than a norm2
- } // at least for contact problems.
-
- R compute_res(bool comp = true) {
- if (comp) compute_residual();
- return residual_norm();
- }
-
- R line_search(VECTOR &dr, const gmm::iteration &iter) {
- size_type nit = 0;
- gmm::resize(stateinit, md.nb_dof());
- gmm::copy(state, stateinit);
- R alpha(1), res, /* res_init, */ R0;
-
- /* res_init = */ res = compute_res(false);
- // cout << "first residual = " << residual() << endl << endl;
- R0 = gmm::real(gmm::vect_sp(dr, rhs));
-
-#if TRACE_SOL
- static int trace_number = 0;
- int trace_iter = 0;
- {
- std::stringstream trace_name;
- trace_name << "line_search_state" << std::setfill('0')
- << std::setw(3) << trace_number << "_000_init";
- gmm::vecsave(trace_name.str(),stateinit);
- }
- trace_number++;
-#endif
-
- ls.init_search(res, iter.get_iteration(), R0);
- do {
- alpha = ls.next_try();
- gmm::add(stateinit, gmm::scaled(dr, alpha), state);
-#if TRACE_SOL
- {
- trace_iter++;
- std::stringstream trace_name;
- trace_name << "line_search_state" << std::setfill('0')
- << std::setw(3) << trace_number << "_"
- << std::setfill('0') << std::setw(3) << trace_iter;
- gmm::vecsave(trace_name.str(), state);
- }
-#endif
- res = compute_res();
- // cout << "residual = " << residual() << endl << endl;
- R0 = gmm::real(gmm::vect_sp(dr, rhs));
-
- ++ nit;
- } while (!ls.is_converged(res, R0));
-
- if (alpha != ls.converged_value()) {
- alpha = ls.converged_value();
- gmm::add(stateinit, gmm::scaled(dr, alpha), state);
- res = ls.converged_residual();
- compute_residual();
- }
-
- return alpha;
- }
-
- model_pb(model &m, abstract_newton_line_search &ls_, VECTOR &st,
- const VECTOR &rhs_, const MATRIX &K_)
- : md(m), ls(ls_), state(st), rhs(rhs_), K(K_) {}
-
- };
//---------------------------------------------------------------------
// Default linear solver.
//---------------------------------------------------------------------
- typedef abstract_linear_solver<model_real_sparse_matrix,
- model_real_plain_vector> rmodel_linear_solver;
- typedef std::shared_ptr<rmodel_linear_solver> rmodel_plsolver_type;
- typedef abstract_linear_solver<model_complex_sparse_matrix,
- model_complex_plain_vector>
- cmodel_linear_solver;
- typedef std::shared_ptr<cmodel_linear_solver> cmodel_plsolver_type;
-
+ typedef std::shared_ptr<abstract_linear_solver<model_real_sparse_matrix,
+ model_real_plain_vector> >
+ rmodel_plsolver_type;
+ typedef std::shared_ptr<abstract_linear_solver<model_complex_sparse_matrix,
+ model_complex_plain_vector> >
+ cmodel_plsolver_type;
template<typename MATRIX, typename VECTOR>
std::shared_ptr<abstract_linear_solver<MATRIX, VECTOR>>
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 8b8c7f2..3e0657d 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -917,6 +917,15 @@ namespace getfem {
return rrhs;
}
+ /** Gives write access to the right hand side of the tangent linear system.
+ Some solvers need to manipulate the model rhs directly so that for
+ example internal condensed variables can be treated properly. */
+ model_real_plain_vector &set_real_rhs() const {
+ GMM_ASSERT1(!complex_version, "This model is a complex one");
+ context_check(); if (act_size_to_be_done) actualize_sizes();
+ return rrhs;
+ }
+
/** Gives access to the part of the right hand side of a term of
a particular nonlinear brick. Does not account of the eventual time
dispatcher. An assembly of the rhs has to be done first.
@@ -946,6 +955,15 @@ namespace getfem {
return crhs;
}
+ /** Gives write access to the right hand side of the tangent linear system.
+ Some solvers need to manipulate the model rhs directly so that for
+ example internal condensed variables can be treated properly. */
+ model_complex_plain_vector &set_complex_rhs() const {
+ GMM_ASSERT1(complex_version, "This model is a real one");
+ context_check(); if (act_size_to_be_done) actualize_sizes();
+ return crhs;
+ }
+
/** Gives access to the part of the right hand side of a term of a
particular nonlinear brick. Does not account of the eventual time
dispatcher. An assembly of the rhs has to be done first.
diff --git a/src/getfem_model_solvers.cc b/src/getfem_model_solvers.cc
index 68c2184..0a5c10d 100644
--- a/src/getfem_model_solvers.cc
+++ b/src/getfem_model_solvers.cc
@@ -25,7 +25,177 @@
namespace getfem {
+ #define TRACE_SOL 0
+ /* ***************************************************************** */
+ /* Intermediary structure for Newton algorithms with getfem::model. */
+ /* ***************************************************************** */
+
+ template <typename PLSOLVER>
+ class pb_base {
+ public:
+ typedef typename PLSOLVER::element_type::MATRIX MATRIX;
+ typedef typename PLSOLVER::element_type::VECTOR VECTOR;
+ typedef typename gmm::linalg_traits<VECTOR>::value_type T;
+ typedef typename gmm::number_traits<T>::magnitude_type R;
+
+ protected:
+ PLSOLVER linear_solver;
+ const MATRIX &K;
+ VECTOR &rhs;
+ VECTOR state;
+
+ public:
+ virtual const VECTOR &residual() const { return rhs; }
+ // A norm1 seems to be better than a norm2 (at least for contact problems).
+ virtual R residual_norm() { return gmm::vect_norm1(residual()); }
+ virtual const VECTOR &state_vector() const { return state; }
+ virtual VECTOR &state_vector() { return state; }
+ virtual R state_norm() const { return gmm::vect_norm1(state_vector()); }
+
+ virtual void perturbation() {
+ R res = gmm::vect_norm2(state), ampl = std::max(res*R(1E-20), R(1E-50));
+ std::vector<R> V(gmm::vect_size(state));
+ gmm::fill_random(V);
+ gmm::add(gmm::scaled(V, ampl), state);
+ }
+
+ virtual void add_to_residual(VECTOR &extra_rhs, R mult=1.) {
+ if (mult == R(1)) gmm::add(extra_rhs, rhs);
+ else if (mult != R(0)) gmm::add(gmm::scaled(extra_rhs, mult), rhs);
+ }
+
+ virtual void linear_solve(VECTOR &dr, gmm::iteration &iter) {
+ (*linear_solver)(K, dr, rhs, iter);
+ }
+
+ pb_base(PLSOLVER linsolv, const MATRIX &K_, VECTOR &rhs_)
+ : linear_solver(linsolv), K(K_), rhs(rhs_), state(gmm::vect_size(rhs_))
+ {}
+ };
+
+ /* ***************************************************************** */
+ /* Linear model problem. */
+ /* ***************************************************************** */
+ template <typename PLSOLVER>
+ class lin_model_pb : pb_base<PLSOLVER> {
+ model &md;
+ public:
+ using typename pb_base<PLSOLVER>::VECTOR;
+ using typename pb_base<PLSOLVER>::R;
+ using pb_base<PLSOLVER>::state_vector;
+ using pb_base<PLSOLVER>::linear_solve;
+ void compute_all() { md.assembly(model::BUILD_ALL); }
+ lin_model_pb(model &, PLSOLVER) = delete;
+ };
+
+ template <>
+ lin_model_pb<rmodel_plsolver_type>::lin_model_pb
+ (model &md_, rmodel_plsolver_type linsolv)
+ : pb_base<rmodel_plsolver_type>
+ (linsolv, md_.real_tangent_matrix(), md_.set_real_rhs()),
+ md(md_)
+ { md.from_variables(state_vector()); }
+ template <>
+ lin_model_pb<cmodel_plsolver_type>::lin_model_pb
+ (model &md_, cmodel_plsolver_type linsolv)
+ : pb_base<cmodel_plsolver_type>
+ (linsolv, md_.complex_tangent_matrix(), md_.set_complex_rhs()),
+ md(md_)
+ { md.from_variables(state_vector()); }
+
+
+ /* ***************************************************************** */
+ /* Non-linear model problem. */
+ /* ***************************************************************** */
+ template <typename PLSOLVER>
+ class nonlin_model_pb : protected pb_base<PLSOLVER> {
+ public:
+ using typename pb_base<PLSOLVER>::VECTOR;
+ using typename pb_base<PLSOLVER>::R;
+ protected:
+ model &md;
+ abstract_newton_line_search &ls;
+ private:
+ VECTOR stateinit; // just a temp used in line_search, could also be mutable
+ public:
+ using pb_base<PLSOLVER>::residual;
+ using pb_base<PLSOLVER>::residual_norm;
+ using pb_base<PLSOLVER>::state_vector;
+ using pb_base<PLSOLVER>::state_norm;
+ using pb_base<PLSOLVER>::add_to_residual;
+ using pb_base<PLSOLVER>::perturbation;
+ using pb_base<PLSOLVER>::linear_solve;
+
+ virtual R approx_external_load_norm() { return md.approx_external_load(); }
+
+ virtual void compute_tangent_matrix() {
+ md.to_variables(state_vector());
+ md.assembly(model::BUILD_MATRIX);
+ }
+
+ virtual void compute_residual() {
+ md.to_variables(state_vector());
+ md.assembly(model::BUILD_RHS);
+ }
+
+ virtual R line_search(VECTOR &dr, const gmm::iteration &iter) {
+ size_type nit = 0;
+ gmm::resize(stateinit, gmm::vect_size(state_vector()));
+ gmm::copy(state_vector(), stateinit);
+ R alpha(1), res, /* res_init, */ R0;
+
+ /* res_init = */ res = residual_norm();
+ // cout << "first residual = " << residual() << endl << endl;
+ R0 = gmm::real(gmm::vect_sp(dr, residual()));
+ ls.init_search(res, iter.get_iteration(), R0);
+ do {
+ alpha = ls.next_try();
+ gmm::add(stateinit, gmm::scaled(dr, alpha), state_vector());
+
+ compute_residual();
+ res = residual_norm();
+ // cout << "residual = " << residual() << endl << endl;
+ R0 = gmm::real(gmm::vect_sp(dr, residual()));
+ ++ nit;
+ } while (!ls.is_converged(res, R0));
+
+ if (alpha != ls.converged_value()) {
+ alpha = ls.converged_value();
+ gmm::add(stateinit, gmm::scaled(dr, alpha), state_vector());
+ res = ls.converged_residual();
+ compute_residual();
+ }
+
+ return alpha;
+ }
+
+ nonlin_model_pb(model &md_, abstract_newton_line_search &ls_,
+ PLSOLVER linear_solver_) = delete;
+ };
+
+ template <>
+ nonlin_model_pb<rmodel_plsolver_type>::nonlin_model_pb
+ (model &md_, abstract_newton_line_search &ls_, rmodel_plsolver_type
linsolv)
+ : pb_base<rmodel_plsolver_type>
+ (linsolv, md_.real_tangent_matrix(), md_.set_real_rhs()),
+ md(md_), ls(ls_)
+ { md.from_variables(state_vector()); }
+ template <>
+ nonlin_model_pb<cmodel_plsolver_type>::nonlin_model_pb
+ (model &md_, abstract_newton_line_search &ls_, cmodel_plsolver_type
linsolv)
+ : pb_base<cmodel_plsolver_type>
+ (linsolv, md_.complex_tangent_matrix(), md_.set_complex_rhs()),
+ md(md_), ls(ls_)
+ { md.from_variables(state_vector()); }
+
+
+
+
+
+ /* ***************************************************************** */
+ /* Linear solvers. */
+ /* ***************************************************************** */
static rmodel_plsolver_type rdefault_linear_solver(const model &md) {
return default_linear_solver<model_real_sparse_matrix,
model_real_plain_vector>(md);
@@ -90,13 +260,12 @@ namespace getfem {
/* time integration schemes. */
/* ***************************************************************** */
- template <typename MATRIX, typename VECTOR, typename PLSOLVER>
+ template <typename PLSOLVER>
void compute_init_values(model &md, gmm::iteration &iter,
PLSOLVER lsolver,
- abstract_newton_line_search &ls, const MATRIX &K,
- const VECTOR &rhs) {
+ abstract_newton_line_search &ls) {
- VECTOR state(md.nb_dof());
+ typename PLSOLVER::element_type::VECTOR state(md.nb_dof());
md.from_variables(state);
md.cancel_init_step();
md.set_time_integration(2);
@@ -107,7 +276,7 @@ namespace getfem {
// Solve for ddt
md.set_time_step(ddt);
gmm::iteration iter1 = iter;
- standard_solve(md, iter1, lsolver, ls, K, rhs);
+ standard_solve(md, iter1, lsolver, ls);
md.copy_init_time_derivative();
// Restore the model state
@@ -121,19 +290,15 @@ namespace getfem {
/* Standard solve. */
/* ***************************************************************** */
- template <typename MATRIX, typename VECTOR, typename PLSOLVER>
+ template <typename PLSOLVER>
void standard_solve(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); // copy the model variables in the state vector
+ abstract_newton_line_search &ls) {
int time_integration = md.is_time_integration();
if (time_integration) {
if (time_integration == 1 && md.is_init_step()) {
- compute_init_values(md, iter, lsolver, ls, K, rhs);
+ compute_init_values(md, iter, lsolver, ls);
return;
}
md.set_time(md.get_time() + md.get_time_step());
@@ -141,31 +306,31 @@ namespace getfem {
}
if (md.is_linear()) {
- md.assembly(model::BUILD_ALL);
- (*lsolver)(K, state, rhs, iter);
- }
- else {
- model_pb<MATRIX, VECTOR> mdpb(md, ls, state, rhs, K);
+ lin_model_pb<PLSOLVER> mdpb(md, lsolver);
+ mdpb.compute_all();
+ mdpb.linear_solve(mdpb.state_vector(), iter);
+ md.to_variables(mdpb.state_vector()); // copy the state vector into the
model variables
+ } else {
+ std::unique_ptr<nonlin_model_pb<PLSOLVER>> mdpb
+ = std::make_unique<nonlin_model_pb<PLSOLVER>>(md, ls, lsolver);
if (dynamic_cast<newton_search_with_step_control *>(&ls))
- Newton_with_step_control(mdpb, iter, *lsolver);
+ Newton_with_step_control(*mdpb, iter);
else
- classical_Newton(mdpb, iter, *lsolver);
+ classical_Newton(*mdpb, iter);
+ md.to_variables(mdpb->state_vector()); // copy the state vector into the
model variables
}
- md.to_variables(state); // copy the state vector into the model variables
}
void standard_solve(model &md, gmm::iteration &iter,
rmodel_plsolver_type lsolver,
abstract_newton_line_search &ls) {
- standard_solve(md, iter, lsolver, ls, md.real_tangent_matrix(),
- md.real_rhs());
+ standard_solve<rmodel_plsolver_type>(md, iter, lsolver, ls);
}
void standard_solve(model &md, gmm::iteration &iter,
cmodel_plsolver_type lsolver,
abstract_newton_line_search &ls) {
- standard_solve(md, iter, lsolver, ls, md.complex_tangent_matrix(),
- md.complex_rhs());
+ standard_solve<cmodel_plsolver_type>(md, iter, lsolver, ls);
}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch devel-logari81-internal-variables updated: Refactoring of model solvers,
Konstantinos Poulios <=