ThunderEgg  1.0.0
CG.h
Go to the documentation of this file.
1 /***************************************************************************
2  * ThunderEgg, a library for solvers on adaptively refined block-structured
3  * Cartesian grids.
4  *
5  * Copyright (c) 2021 Scott Aiton
6  *
7  * This program is free software: you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation, either version 3 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program. If not, see <https://www.gnu.org/licenses/>.
19  ***************************************************************************/
20 
21 #ifndef THUNDEREGG_ITERATIVE_CG_H
22 #define THUNDEREGG_ITERATIVE_CG_H
23 
31 #include <ThunderEgg/Operator.h>
32 #include <ThunderEgg/Timer.h>
33 
34 namespace ThunderEgg::Iterative {
40 template<int D>
41 class CG : public Solver<D>
42 {
43 private:
47  int max_iterations = 1000;
51  double tolerance = 1e-12;
55  std::shared_ptr<Timer> timer = nullptr;
56 
57  void applyWithPreconditioner(const Operator<D>* M_l,
58  const Operator<D>& A,
59  const Operator<D>* M_r,
60  const Vector<D>& x,
61  Vector<D>& b) const
62  {
63  if (M_l == nullptr && M_r == nullptr) {
64  A.apply(x, b);
65  } else if (M_l == nullptr && M_r != nullptr) {
66  Vector<D> tmp = b.getZeroClone();
67  M_r->apply(x, tmp);
68  A.apply(tmp, b);
69  }
70  }
71 
72 public:
78  CG<D>* clone() const override { return new CG<D>(*this); }
86  void setMaxIterations(int max_iterations_in) { max_iterations = max_iterations_in; };
94  int getMaxIterations() const { return max_iterations; }
102  void setTolerance(double tolerance_in) { tolerance = tolerance_in; };
110  double getTolerance() const { return tolerance; }
116  void setTimer(std::shared_ptr<Timer> timer_in) { timer = timer_in; }
117 
123  std::shared_ptr<Timer> getTimer() const { return timer; }
124 
125 public:
126  int solve(const Operator<D>& A,
127  Vector<D>& x,
128  const Vector<D>& b,
129  const Operator<D>* Mr = nullptr,
130  bool output = false,
131  std::ostream& os = std::cout) const override
132  {
133  Vector<D> resid = b.getZeroClone();
134 
135  A.apply(x, resid);
136  resid.scaleThenAdd(-1, b);
137 
138  Vector<D> initial_guess = x;
139  x.set(0);
140 
141  double r0_norm = b.twoNorm();
142  Vector<D> p = resid;
143  Vector<D> ap = b.getZeroClone();
144 
145  double rho = resid.dot(resid);
146 
147  int num_its = 0;
148  if (r0_norm == 0) {
149  return num_its;
150  }
151  double residual = resid.twoNorm() / r0_norm;
152  if (output) {
153  char buf[100];
154  sprintf(buf, "%5d %16.8e\n", num_its, residual);
155  os << std::string(buf);
156  }
157  while (residual > tolerance && num_its < max_iterations) {
158  if (timer) {
159  timer->start("Iteration");
160  }
161 
162  if (rho == 0) {
163  throw BreakdownError("CG broke down, rho was 0 on iteration " + std::to_string(num_its));
164  }
165 
166  applyWithPreconditioner(nullptr, A, Mr, p, ap);
167  double alpha = rho / p.dot(ap);
168  x.addScaled(alpha, p);
169  resid.addScaled(-alpha, ap);
170 
171  double rho_new = resid.dot(resid);
172  double beta = rho_new / rho;
173  p.scaleThenAdd(beta, resid);
174 
175  num_its++;
176  rho = rho_new;
177  residual = resid.twoNorm() / r0_norm;
178 
179  if (output) {
180  char buf[100];
181  sprintf(buf, "%5d %16.8e\n", num_its, residual);
182  os << std::string(buf);
183  }
184  if (timer) {
185  timer->stop("Iteration");
186  }
187  }
188  if (Mr != nullptr) {
189  Mr->apply(x, resid);
190  x.copy(resid);
191  }
192  x.add(initial_guess);
193  return num_its;
194  }
195 };
196 } // namespace ThunderEgg::Iterative
197 extern template class ThunderEgg::Iterative::CG<2>;
198 extern template class ThunderEgg::Iterative::CG<3>;
199 #endif
ThunderEgg::Vector::scaleThenAdd
void scaleThenAdd(double alpha, const Vector< D > &b)
this = alpha * this + b
Definition: Vector.h:509
ThunderEgg::Iterative::BreakdownError
Breakdown exception for iterative methods.
Definition: BreakdownError.h:35
Operator.h
Operator class.
ThunderEgg::Iterative::CG::getMaxIterations
int getMaxIterations() const
Get the maximum number of iterations.
Definition: CG.h:94
ThunderEgg::Vector::set
void set(double alpha)
set all value in the vector
Definition: Vector.h:391
Solver.h
Solver class.
ThunderEgg::Vector::addScaled
void addScaled(double alpha, const Vector< D > &b)
this = this + alpha * b
Definition: Vector.h:483
ThunderEgg::Iterative::CG::solve
int solve(const Operator< D > &A, Vector< D > &x, const Vector< D > &b, const Operator< D > *Mr=nullptr, bool output=false, std::ostream &os=std::cout) const override
Perform an iterative solve.
Definition: CG.h:126
ThunderEgg::Iterative
Iterative Solvers.
Definition: BiCGStab.h:34
ThunderEgg::Iterative::CG::setTolerance
void setTolerance(double tolerance_in)
Set the stopping tolerance.
Definition: CG.h:102
ThunderEgg::Iterative::CG::getTimer
std::shared_ptr< Timer > getTimer() const
Get the Timer object.
Definition: CG.h:123
Timer.h
Timer class.
ThunderEgg::Vector::copy
void copy(const Vector< D > &b)
copy the values of the other vector
Definition: Vector.h:443
ThunderEgg::Operator::apply
virtual void apply(const Vector< D > &x, Vector< D > &b) const =0
Virtual function that base classes have to implement.
ThunderEgg::Vector::twoNorm
double twoNorm() const
get the l2norm
Definition: Vector.h:553
ThunderEgg::Vector::add
void add(const Vector< D > &b)
add the other vector to this vector
Definition: Vector.h:471
ThunderEgg::Iterative::CG
CG iterative solver.
Definition: CG.h:41
ThunderEgg::Iterative::CG::setTimer
void setTimer(std::shared_ptr< Timer > timer_in)
Set the Timer object.
Definition: CG.h:116
ThunderEgg::Operator
Base class for operators.
Definition: Operator.h:37
ThunderEgg::Iterative::Solver
Abstract interface for Iterative solvers.
Definition: Solver.h:39
ThunderEgg::Vector
Vector class for use in thunderegg.
Definition: Vector.h:42
ThunderEgg::Vector::dot
double dot(const Vector< D > &b) const
get the dot product
Definition: Vector.h:583
ThunderEgg::Vector::getZeroClone
Vector< D > getZeroClone() const
Get a vector of the same length initialized to zero.
Definition: Vector.h:601
ThunderEgg::Iterative::CG::getTolerance
double getTolerance() const
Get the stopping tolerance.
Definition: CG.h:110
ThunderEgg::Iterative::CG::setMaxIterations
void setMaxIterations(int max_iterations_in)
Set the maximum number of iterations.
Definition: CG.h:86
BreakdownError.h
BreakdownError struct.
ThunderEgg::Iterative::CG::clone
CG< D > * clone() const override
Clone this solver.
Definition: CG.h:78