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 <logar...@gmail.com> 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); }