diff --git a/src/solver/solver_cg.cc b/src/solver/solver_cg.cc index c6885a1..7da51cb 100644 --- a/src/solver/solver_cg.cc +++ b/src/solver/solver_cg.cc @@ -1,128 +1,128 @@ /** * @file solver_cg.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief Implementation of cg solver * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "solver/solver_cg.hh" #include "solver/solver_error.hh" #include #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template SolverCG::SolverCG(Cell_t& cell, Real tol, Uint maxiter, - bool verbose) - :Parent(cell, tol, maxiter, verbose), + bool verbose, Communicator comm) + :Parent(cell, tol, maxiter, verbose), comm{comm}, r_k{make_field("residual r_k", this->collection)}, p_k{make_field("search direction r_k", this->collection)}, Ap_k{make_field("Effect of tangent A*p_k", this->collection)} {} /* ---------------------------------------------------------------------- */ template void SolverCG::solve(const Field_t & rhs, Field_t & x_f) { - x_f.eigenvec() = this-> solve(rhs.eigenvec(), x_f.eigenvec()); + x_f.eigenvec() = this->solve(rhs.eigenvec(), x_f.eigenvec()); }; //----------------------------------------------------------------------------// template typename SolverCG::SolvVectorOut SolverCG::solve(const SolvVectorInC rhs, SolvVectorIn x_0) { // Following implementation of algorithm 5.2 in Nocedal's Numerical Optimization (p. 112) auto r = this->r_k.eigen(); auto p = this->p_k.eigen(); auto Ap = this->Ap_k.eigen(); typename Field_t::EigenMap x(x_0.data(), r.rows(), r.cols()); // initialisation of algo r = this->cell.directional_stiffness_with_copy(x); r -= typename Field_t::ConstEigenMap(rhs.data(), r.rows(), r.cols()); p = -r; this->converged = false; - Real rdr = (r*r).sum(); - Real rhs_norm2 = rhs.squaredNorm(); + Real rdr = this->comm.sum((r*r).sum()); + Real rhs_norm2 = this->comm.sum(rhs.squaredNorm()); Real tol2 = ipow(this->tol,2)*rhs_norm2; size_t count_width{}; // for output formatting in verbose case if (this->verbose) { count_width = size_t(std::log10(this->maxiter))+1; } for (Uint i = 0; i < this->maxiter && (rdr > tol2 || i == 0); ++i, ++this->counter) { Ap = this->cell.directional_stiffness_with_copy(p); - Real alpha = rdr/(p*Ap).sum(); + Real alpha = rdr/this->comm.sum((p*Ap).sum()); x += alpha * p; r += alpha * Ap; - Real new_rdr = (r*r).sum(); + Real new_rdr = this->comm.sum((r*r).sum()); Real beta = new_rdr/rdr; rdr = new_rdr; if (this->verbose) { std::cout << " at CG step " << std::setw(count_width) << i << ": |r|/|b| = " << std::setw(15) << std::sqrt(rdr/rhs_norm2) << ", cg_tol = " << this->tol << std::endl; } p = -r+beta*p; } if (rdr < tol2) { this->converged=true; } else { std::stringstream err {}; err << " After " << this->counter << " steps, the solver " << " FAILED with |r|/|b| = " << std::setw(15) << std::sqrt(rdr/rhs_norm2) << ", cg_tol = " << this->tol << std::endl; throw ConvergenceError("Conjugate gradient has not converged." + err.str()); } return x_0; } /* ---------------------------------------------------------------------- */ template typename SolverCG::Tg_req_t SolverCG::get_tangent_req() const { return tangent_requirement; } template class SolverCG; //template class SolverCG; template class SolverCG; } // muSpectre diff --git a/src/solver/solver_cg.hh b/src/solver/solver_cg.hh index fffd4d0..6980881 100644 --- a/src/solver/solver_cg.hh +++ b/src/solver/solver_cg.hh @@ -1,113 +1,116 @@ /** * @file solver_cg.hh * * @author Till Junge * * @date 20 Dec 2017 * * @brief class for a simple implementation of a conjugate gradient * solver. This follows algorithm 5.2 in Nocedal's Numerical * Optimization (p 112) * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef SOLVER_CG_H #define SOLVER_CG_H #include "solver/solver_base.hh" +#include "common/communicator.hh" #include "common/field.hh" #include namespace muSpectre { /** * implements the `muSpectre::SolverBase` interface using a * conjugate gradient solver. This particular class is useful for * trouble shooting, as it can be made very verbose, but for * production runs, it is probably better to use * `muSpectre::SolverCGEigen`. */ template class SolverCG: public SolverBase { public: using Parent = SolverBase; //!< base class //! Input vector for solvers using SolvVectorIn = typename Parent::SolvVectorIn; //! Input vector for solvers using SolvVectorInC = typename Parent::SolvVectorInC; //! Output vector for solvers using SolvVectorOut = typename Parent::SolvVectorOut; using Cell_t = typename Parent::Cell_t; //!< cell type using Ccoord = typename Parent::Ccoord; //!< cell coordinates type //! kind of tangent that is required using Tg_req_t = typename Parent::TangentRequirement; //! cg only needs to handle fields that look like strain and stress using Field_t = TensorField< typename Parent::Collection_t, Real, secondOrder, DimM>; //! conjugate gradient needs directional stiffness constexpr static Tg_req_t tangent_requirement{Tg_req_t::NeedEffect}; //! Default constructor SolverCG() = delete; //! Constructor with domain resolutions, etc, - SolverCG(Cell_t& cell, Real tol, Uint maxiter=0, bool verbose =false); + SolverCG(Cell_t& cell, Real tol, Uint maxiter=0, bool verbose=false, + Communicator comm=Communicator()); //! Copy constructor SolverCG(const SolverCG &other) = delete; //! Move constructor SolverCG(SolverCG &&other) = default; //! Destructor virtual ~SolverCG() = default; //! Copy assignment operator SolverCG& operator=(const SolverCG &other) = delete; //! Move assignment operator SolverCG& operator=(SolverCG &&other) = default; bool has_converged() const override final {return this->converged;} //! actual solver void solve(const Field_t & rhs, Field_t & x); // this simplistic implementation has no initialisation phase so the default is ok SolvVectorOut solve(const SolvVectorInC rhs, SolvVectorIn x_0) override final; std::string name() const override final {return "CG";} protected: + Communicator comm; //!< communicator //! returns `muSpectre::Tg_req_t::NeedEffect` Tg_req_t get_tangent_req() const override final; Field_t & r_k; //!< residual Field_t & p_k; //!< search direction Field_t & Ap_k; //!< effect of tangent on search direction bool converged{false}; //!< whether the solver has converged private: }; } // muSpectre #endif /* SOLVER_CG_H */