OEMSolver.h 7.86 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file OEMSolver.h */

/**
 * \defgroup Solver Solver module
 * @{ <img src="solver.png"> @}
 * 
 * \brief
 * Contains all classes needed for solving linear and non linear equation
 * systems.
 */

#ifndef AMDIS_OEM_SOLVER_H
#define AMDIS_OEM_SOLVER_H

#include "Global.h"
35
36
37
38
39
40
41
#include "AMDiS_fwd.h"
#include "Parameters.h"
#include "SolverMatrix.h"
#include "DOFVector.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "ITL_Preconditioner.h"
42
43
44
45
46
47
48
49
50
51
52
53

namespace AMDiS {

  /**
   * \ingroup Solver
   * 
   *\brief
   * Solver for linear equation systems.
   */
  class OEMSolver
  {
  public:
54
55
56
57

    typedef DOFMatrix::base_matrix_type   matrix_type;
    typedef DOFMatrix::value_type         value_type;

58
59
60
    /** \brief
     * The constructor reads needed parameters and sets solvers \ref name.
     */
61
62
63
64
65
66
67
68
69
70
71
    OEMSolver(std::string str) 
      : name(str),
	tolerance(0),
	relative(0),
	max_iter(1000),
	info(0),
	residual(0),
	print_cycle(100),
	leftPrecon(NULL),
	rightPrecon(NULL)
    {}
72
73

    /** \brief
74
     *      
75
     */
76
    virtual ~OEMSolver() 
77
    {
78
79
80
      if (leftPrecon) delete leftPrecon;
      if (rightPrecon) delete rightPrecon;
    }
81

82
83
84
    void initParameters()
    {
      FUNCNAME("OEMSolver::initParameters()");
85
      
86
87
88
89
90
91
      GET_PARAMETER(0, name + "->tolerance", "%f", &tolerance);
      GET_PARAMETER(0, name + "->relative tolerance", "%f", &relative);
      GET_PARAMETER(0, name + "->max iteration", "%d", &max_iter);
      GET_PARAMETER(0, name + "->print cycle", "%d", &print_cycle);
      GET_PARAMETER(0, name + "->info", "%d", &info);
    }
92

93
94
95
96
97
98
99
100
101
102
    /** Set left Preconditioner
     *
     *  It is assumed that the solver owns the preconditioner from this call.
     *  That means that the preconditioner is deleted by the solver when not needed.
     */
    void setLeftPrecon(ITL_BasePreconditioner* p) 
    {
	if (leftPrecon) delete leftPrecon;
	leftPrecon= p;
    }
103

104
105
106
107
108
109
110
111
112
113
    /** Set right Preconditioner
     *
     *  It is assumed that the solver owns the preconditioner from this call.
     *  That means that the preconditioner is deleted by the solver when not needed.
     */
    void setRightPrecon(ITL_BasePreconditioner* p) 
    {
	if (rightPrecon) delete rightPrecon;
	rightPrecon= p;
    }
114

115
116
117
118
119
    /// Linear System to be solved in the derived class
    virtual int solveSystem(const DOFMatrix::base_matrix_type&      A,
			    mtl::dense_vector<value_type>&          x,
			    const mtl::dense_vector<value_type>&    b) = 0;
    
120

121
122
123
124
125
126
127
128
129
    /// Solve a linear system for a scalar problem.
    int solveSystem(const SolverMatrix<DOFMatrix>& A,
		    DOFVector<double>&             x, 
		    DOFVector<double>&             b) 
    {
      mtl::dense_vector<value_type>   xx(x.getUsedSize(), &x[0]),
	                              bb(b.getUsedSize(), &b[0]);
      return solveSystem(A.getMatrix(), xx, bb);
    }
130

131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
    /// Solve a linear system for a vectorial problem.
    int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
		    SystemVector&                            x, 
		    SystemVector&                            b) 
    {
      int ns=  x.getSize(),  // Number of systems.
          size= x.getUsedSize(); // Size of all DOFVectors

      // Copy vectors
      mtl::dense_vector<value_type>   xx(size), bb(size);
      for (int i = 0, counter = 0; i < ns; i++) {
	  DOFVector<double>::Iterator it(b.getDOFVector(i), USED_DOFS);	  
	  for (it.reset(); !it.end(); ++it)	
	      bb[counter++] = *it;	
      }

      // Solver on DOFVector for single system
      int r= solveSystem(A.getMatrix(), xx, bb);

      for (int i = 0, counter = 0; i < ns; i++) {
	  DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);	  
	  for (it.reset(); !it.end(); ++it)	
	      *it = xx[counter++];	
      }
      return r;
156
    }
157
158
159
160
161
162
163
164
165
166

    // ===== getting-methods ======================================================

    /** \name getting methods
     * \{ 
     */

    /** \brief
     * Returns solvers \ref name.
     */
167
    inline const std::string& getName() { 
168
      return name; 
169
    }
170
171
172
173
174
175

    /** \brief
     * Returns \ref tolerance
     */
    inline double getTolerance() {
      return tolerance;
176
    }  
177
178
179
180
181
182

    /** \brief
     * Returns \ref max_iter
     */
    inline int getMaxIterations() {
      return max_iter;
183
    }
184
185
186
187
188
189

    /** \brief
     * Returns \ref residual
     */
    inline double getResidual() {
      return residual;
190
    }
191
192
193
194
195
196
197
198
199
200
201
202

    /** \} */ 

    // ===== setting-methods ======================================================

    /** \name setting methods
     * \{ 
     */

    /** \brief
     * Sets \ref tolerance
     */
203
    inline void setTolerance(double tol) {
204
      tolerance = tol;
205
    }
206
207
208
209

    /** \brief
     * Sets \ref relative
     */
210
    inline void setRelative(bool rel) {
211
      relative = rel;
212
    }
213
214
215
216

    /** \brief
     * Sets \ref max_iter
     */
217
    inline void setMaxIterations(int i) {
218
      max_iter = i;
219
    }
220
221
222
223

    /** \brief
     * Sets \ref info
     */
224
    inline void setInfo(int i) {
225
      info = i;
226
    }
227
228
229
230
231
232
233

    /** \} */

  protected:
    /** \brief
     * solvers name.
     */
234
    std::string name;
235
236

    /** \brief
237
     * Solver tolerance |r|. Set in OEMSolver's constructor. 
238
239
240
241
     */
    double tolerance;

    /** \brief
242
     * Relative solver tolerance |r|/|r0|. Set in OEMSolver's constructor. 
243
     */
244
    double relative;
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260

    /** \brief
     * maximal number of iterations. Set in OEMSolver's constructor.
     */
    int max_iter;

    /** \brief
     * info level during solving the system. Set in OEMSolver's constructor.
     */
    int info;

    /** \brief
     * current residual.
     */
    double residual;

261
262
    /// Print cycle, after how many iterations the residuum norm is logged.
    int print_cycle;
263

264
    ITL_BasePreconditioner *leftPrecon;
265

266
    ITL_BasePreconditioner *rightPrecon;
267
268
269
270
271
272
273
274
275
276
277
278
  };

  // ============================================================================
  // ===== class OEMSolverCreator ===============================================
  // ============================================================================

  /** 
   * \ingroup Solver
   *
   * \brief
   * Interface for creators of concrete OEMSolvers. 
   */
279
  class OEMSolverCreator : public CreatorInterface<OEMSolver>
280
281
  { 
  public:
282
    virtual ~OEMSolverCreator() {}
283
284
285
286

    /** \brief
     * Sets \ref problem
     */
287
288
    void setName(std::string str) { 
      name = str; 
289
    }
290
291
292
293
294
295

  protected:
    /** \brief
     * Pointer to the problem the solver will belong to. Needed as parameter
     * when constructing an OEMSolver
     */
296
    std::string name;
297
298
299
  };
}

300
#include "ITL_Solver.h"
301
302

#endif // AMDIS_OEM_SOLVER_H