Chombo + EB + MF  3.2
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  Opposite of create -- perform any operations required before lhs goes
66  out of scope. In general, the only time this needs to be defined in
67  a derived class is if the create() function called new. Otherwise, the
68  default no-op function is sufficient.
69  */
70  virtual void clear(T& a_lhs)
71  {
72  }
73 
74  ///
75  /**
76  Set a_lhs equal to a_rhs.
77  */
78  virtual void assign( T& a_lhs, const T& a_rhs) = 0;
79 
80  virtual void assignLocal(T& a_lhs, const T& a_rhs)
81  {
82  this->assign(a_lhs, a_rhs);
83  }
84  ///
85  /**
86  Compute and return the dot product of a_1 and a_2. In most contexts, this
87  means return the sum over all data points of a_1*a_2.
88  */
89  virtual Real dotProduct(const T& a_1, const T& a_2) = 0;
90  /* multiple dot products (for GMRES) */
91  virtual void mDotProduct(const T& a_1, const int a_sz, const T a_2[], Real a_mdots[])
92  {
93  for (int j=0; j<a_sz; j++)
94  {
95  a_mdots[j] = dotProduct(a_1, a_2[j]);
96  }
97  }
98 
99  ///
100  /**
101  Increment by scaled amount (a_lhs += a_scale*a_x).
102  */
103  virtual void incr ( T& a_lhs, const T& a_x, Real a_scale) = 0;
104 
105  ///
106  /**
107  Set input to a scaled sum (a_lhs = a_a*a_x + a_b*a_y).
108  */
109  virtual void axby( T& a_lhs, const T& a_x, const T& a_y, Real a_a, Real a_b) = 0;
110 
111  ///
112  /**
113  Multiply the input by a given scale (a_lhs *= a_scale).
114  */
115  virtual void scale( T& a_lhs, const Real& a_scale) = 0;
116 
117  ///
118  /**
119  Return the norm of a_rhs.
120  a_ord == 0 max norm, a_ord == 1 sum(abs(a_rhs)), else, L(a_ord) norm.
121  */
122  virtual Real norm(const T& a_rhs, int a_ord) = 0;
123 
124  ///
125  /**
126  Return dx at this level of refinement
127  */
128  virtual Real dx() const
129  {
130  MayDay::Warning(" calling dx on base class\n");
131  return 0.;
132  }
133  ///
134  /**
135  Set a_lhs to zero.
136  */
137  virtual void setToZero(T& a_lhs) = 0;
138 
139  ///
140  /**
141  Debugging aid for solvers. Print out a "T" to a file named "filename"
142  default implementation is to print out a message saying "LinearOp::write not implemented"
143  */
144  virtual void write(const T* a, const char* filename)
145  {
146  MayDay::Warning("LinearOp::write not implemented");
147  }
148 
149 };
150 
151 ///
152 /**
153  Generic linear solver template. BiCGStab and others are built on top of this.
154 */
155 template <class T>
157 {
158 public:
159  virtual ~LinearSolver()
160  {
161  }
162 
163  ///
164  /**
165  reset whether the solver is homogeneous.
166  */
167  virtual void setHomogeneous(bool a_homogeneous) = 0;
168 
169  ///
170  /**
171  Define the operator and whether it is a homogeneous solver or not.
172  The LinearSolver does not take over ownership of this a_operator object.
173  It does not call delete on it when the LinearSolver is deleted. It is
174  meant to be like a late-binding reference. If you created a_operator
175  with new, you should call delete on it after LinearSolver is deleted
176  if you want to avoid memory leaks.
177 
178  */
179  virtual void define(LinearOp<T>* a_operator, bool a_homogeneous = false) = 0;
180 
181  ///
182  /**
183  Solve L(phi) = rhs (phi = L^-1 (rhs)).
184  */
185  virtual void solve(T& a_phi, const T& a_rhs) = 0;
186 
187  /// Set a convergence metric, along with solver tolerance, if desired
188  /**
189  Default implementation does nothing, since there are probably
190  cases (liked direct solves), where this has no real meaning.
191  */
192  virtual void setConvergenceMetrics(Real a_metric, Real a_tolerance)
193  {
194  }
195 };
196 
197 #include "NamespaceFooter.H"
198 
199 #endif /*_LINEARSOLVER_H_*/
Definition: LinearSolver.H:28
virtual void write(const T *a, const char *filename)
Definition: LinearSolver.H:144
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:192
virtual Real norm(const T &a_rhs, int a_ord)=0
virtual void setToZero(T &a_lhs)=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 clear(T &a_lhs)
Definition: LinearSolver.H:70
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:91
virtual void create(T &a_lhs, const T &a_rhs)=0
Definition: LinearSolver.H:156
virtual void assignLocal(T &a_lhs, const T &a_rhs)
Definition: LinearSolver.H:80
virtual ~LinearOp()
Definition: LinearSolver.H:32
virtual void assign(T &a_lhs, const T &a_rhs)=0
virtual ~LinearSolver()
Definition: LinearSolver.H:159
virtual Real dotProduct(const T &a_1, const T &a_2)=0
virtual Real dx() const
Definition: LinearSolver.H:128
virtual void residual(T &a_lhs, const T &a_phi, const T &a_rhs, bool a_homogeneous=false)=0