Chombo + EB  3.0
LinearSolver.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 // BVS, June 18, 2003
12 
13 // We can assume that template class T has null construction.
14 
15 #ifndef _LINEARSOLVER_H_
16 #define _LINEARSOLVER_H_
17 
18 #include "REAL.H"
19 #include "Box.H"
20 #include <cmath>
21 #include "NamespaceHeader.H"
22 
23 ///
24 /**
25  Operator class for Linear solver for solving L(phi) = rhs.
26  */
27 template <class T>
28 class LinearOp
29 {
30 public:
31  ///
32  virtual ~LinearOp()
33  {
34  }
35 
36  ///
37  /**
38  Say you are solving L(phi) = rhs. Make a_lhs = L(a_phi) - a_rhs. If a_homogeneous is true,
39  evaluate the operator using homogeneous boundary conditions.
40  */
41  virtual void residual( T& a_lhs, const T& a_phi, const T& a_rhs, bool a_homogeneous = false) = 0;
42 
43  ///
44  /**
45  Given the current state of the residual the correction, apply your preconditioner to a_cor.
46  */
47  virtual void preCond( T& a_cor, const T& a_residual) = 0;
48 
49  ///
50  /**
51  In the context of solving L(phi) = rhs, set a_lhs = L(a_phi). If a_homogeneous is true,
52  evaluate the operator using homogeneous boundary conditions.
53  */
54  virtual void applyOp( T& a_lhs, const T& a_phi, bool a_homogeneous = false) = 0;
55 
56  ///
57  /**
58  Creat data holder a_lhs that mirrors a_rhs. You do not need to copy the data of a_rhs,
59  just make a holder the same size.
60  */
61  virtual void create( T& a_lhs, const T& a_rhs) = 0;
62 
63 
64  ///
65  /**
66  Opposite of create -- perform any operations required before lhs goes
67  out of scope. In general, the only time this needs to be defined in
68  a derived class is if the create() function called new. Otherwise, the
69  default no-op function is sufficient.
70  */
71  virtual void clear(T& a_lhs)
72  {
73  }
74 
75 
76  ///
77  /**
78  Set a_lhs equal to a_rhs.
79  */
80  virtual void assign( T& a_lhs, const T& a_rhs) = 0;
81 
82  virtual void assignLocal(T& a_lhs, const T& a_rhs)
83  {
84  this->assign(a_lhs, a_rhs);
85  }
86  ///
87  /**
88  Compute and return the dot product of a_1 and a_2. In most contexts, this
89  means return the sum over all data points of a_1*a_2.
90  */
91  virtual Real dotProduct(const T& a_1, const T& a_2) = 0;
92  /* multiple dot products (for GMRES) */
93  virtual void mDotProduct(const T& a_1, const int a_sz, const T a_2[], Real a_mdots[])
94  {
95  for (int j=0; j<a_sz; j++)
96  {
97  a_mdots[j] = dotProduct(a_1, a_2[j]);
98  }
99  }
100 
101  ///
102  /**
103  Increment by scaled amount (a_lhs += a_scale*a_x).
104  */
105  virtual void incr ( T& a_lhs, const T& a_x, Real a_scale) = 0;
106 
107  ///
108  /**
109  Set input to a scaled sum (a_lhs = a_a*a_x + a_b*a_y).
110  */
111  virtual void axby( T& a_lhs, const T& a_x, const T& a_y, Real a_a, Real a_b) = 0;
112 
113  ///
114  /**
115  Multiply the input by a given scale (a_lhs *= a_scale).
116  */
117  virtual void scale( T& a_lhs, const Real& a_scale) = 0;
118 
119  ///
120  /**
121  Return the norm of a_rhs.
122  a_ord == 0 max norm, a_ord == 1 sum(abs(a_rhs)), else, L(a_ord) norm.
123  */
124  virtual Real norm(const T& a_rhs, int a_ord) = 0;
125 
126  ///
127  /**
128  Return dx at this level of refinement
129  */
130  virtual Real dx() const
131  {
132  MayDay::Warning(" calling dx on base class\n");
133  return 0.;
134  }
135  ///
136  /**
137  Set a_lhs to zero.
138  */
139  virtual void setToZero(T& a_lhs) = 0;
140 
141  ///
142  /**
143  Debugging aid for solvers. Print out a "T" to a file named "filename"
144  default implementation is to print out a message saying "LinearOp::write not implemented"
145  */
146  virtual void write(const T* a, const char* filename)
147  {
148  MayDay::Warning("LinearOp::write not implemented");
149  }
150 
151 };
152 
153 ///
154 /**
155  Generic linear solver template. BiCGStab and others are built on top of this.
156 */
157 template <class T>
159 {
160 public:
161  virtual ~LinearSolver()
162  {
163  }
164 
165  ///
166  /**
167  reset whether the solver is homogeneous.
168  */
169  virtual void setHomogeneous(bool a_homogeneous) = 0;
170 
171  ///
172  /**
173  Define the operator and whether it is a homogeneous solver or not.
174  The LinearSolver does not take over ownership of this a_operator object.
175  It does not call delete on it when the LinearSolver is deleted. It is
176  meant to be like a late-binding reference. If you created a_operator
177  with new, you should call delete on it after LinearSolver is deleted
178  if you want to avoid memory leaks.
179 
180  */
181  virtual void define(LinearOp<T>* a_operator, bool a_homogeneous = false) = 0;
182 
183  ///
184  /**
185  Solve L(phi) = rhs (phi = L^-1 (rhs)).
186  */
187  virtual void solve(T& a_phi, const T& a_rhs) = 0;
188 
189 
190  /// Set a convergence metric, along with solver tolerance, if desired
191  /**
192  Default implementation does nothing, since there are probably
193  cases (liked direct solves), where this has no real meaning.
194  */
195  virtual void setConvergenceMetrics(Real a_metric, Real a_tolerance)
196  {
197  }
198 };
199 
200 #include "NamespaceFooter.H"
201 #endif /*_LINEARSOLVER_H_*/
Definition: LinearSolver.H:28
virtual void write(const T *a, const char *filename)
Definition: LinearSolver.H:146
virtual void scale(T &a_lhs, const Real &a_scale)=0
virtual void applyOp(T &a_lhs, const T &a_phi, bool a_homogeneous=false)=0
virtual void preCond(T &a_cor, const T &a_residual)=0
virtual void setConvergenceMetrics(Real a_metric, Real a_tolerance)
Set a convergence metric, along with solver tolerance, if desired.
Definition: LinearSolver.H:195
virtual Real norm(const T &a_rhs, int a_ord)=0
virtual Real dx() const
Definition: LinearSolver.H:130
virtual void setToZero(T &a_lhs)=0
virtual void solve(T &a_phi, const T &a_rhs)=0
virtual void axby(T &a_lhs, const T &a_x, const T &a_y, Real a_a, Real a_b)=0
double Real
Definition: REAL.H:33
virtual void incr(T &a_lhs, const T &a_x, Real a_scale)=0
virtual void define(LinearOp< T > *a_operator, bool a_homogeneous=false)=0
virtual void clear(T &a_lhs)
Definition: LinearSolver.H:71
virtual void setHomogeneous(bool a_homogeneous)=0
static void Warning(const char *const a_msg=m_nullString)
Print out message to cerr and continue.
virtual void mDotProduct(const T &a_1, const int a_sz, const T a_2[], Real a_mdots[])
Definition: LinearSolver.H:93
virtual void create(T &a_lhs, const T &a_rhs)=0
Definition: LinearSolver.H:158
virtual void assignLocal(T &a_lhs, const T &a_rhs)
Definition: LinearSolver.H:82
virtual ~LinearOp()
Definition: LinearSolver.H:32
virtual void assign(T &a_lhs, const T &a_rhs)=0
virtual ~LinearSolver()
Definition: LinearSolver.H:161
virtual Real dotProduct(const T &a_1, const T &a_2)=0
virtual void residual(T &a_lhs, const T &a_phi, const T &a_rhs, bool a_homogeneous=false)=0