FreeFOAM The Cross-Platform CFD Toolkit
lduMatrix.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Class
25  Foam::lduMatrix
26 
27 Description
28  lduMatrix is a general matrix class in which the coefficients are
29  stored as three arrays, one for the upper triangle, one for the
30  lower triangle and a third for the diagonal.
31 
32  Addressing arrays must be supplied for the upper and lower triangles.
33 
34  It might be better if this class were organised as a hierachy starting
35  from an empty matrix, then deriving diagonal, symmetric and asymmetric
36  matrices.
37 
38 SourceFiles
39  lduMatrixATmul.C
40  lduMatrix.C
41  lduMatrixTemplates.C
42  lduMatrixOperations.C
43  lduMatrixSolver.C
44  lduMatrixPreconditioner.C
45  lduMatrixTests.C
46  lduMatrixUpdateMatrixInterfaces.C
47 
48 \*---------------------------------------------------------------------------*/
49 
50 #ifndef lduMatrix_H
51 #define lduMatrix_H
52 
53 #include <OpenFOAM/lduMesh.H>
55 #include <OpenFOAM/FieldField.H>
57 #include <OpenFOAM/typeInfo.H>
58 #include <OpenFOAM/autoPtr.H>
60 
61 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
62 
63 namespace Foam
64 {
65 
66 // Forward declaration of friend functions and operators
67 
68 class lduMatrix;
69 Ostream& operator<<(Ostream&, const lduMatrix&);
70 
71 
72 /*---------------------------------------------------------------------------*\
73  Class lduMatrix Declaration
74 \*---------------------------------------------------------------------------*/
75 
76 class lduMatrix
77 {
78  // private data
79 
80  //- LDU mesh reference
81  const lduMesh& lduMesh_;
82 
83  //- Coefficients (not including interfaces)
84  scalarField *lowerPtr_, *diagPtr_, *upperPtr_;
85 
86 
87 public:
88 
89  //- Class returned by the solver, containing performance statistics
91  {
92  word solverName_;
93  word fieldName_;
94  scalar initialResidual_;
95  scalar finalResidual_;
96  label noIterations_;
97  bool converged_;
98  bool singular_;
99 
100 
101  public:
102 
103  // Constructors
104 
106  :
107  initialResidual_(0),
108  finalResidual_(0),
109  noIterations_(0),
110  converged_(false),
111  singular_(false)
112  {}
113 
114 
116  (
117  const word& solverName,
118  const word& fieldName,
119  const scalar iRes = 0,
120  const scalar fRes = 0,
121  const label nIter = 0,
122  const bool converged = false,
123  const bool singular = false
124  )
125  :
126  solverName_(solverName),
127  fieldName_(fieldName),
128  initialResidual_(iRes),
129  finalResidual_(fRes),
130  noIterations_(nIter),
131  converged_(converged),
132  singular_(singular)
133  {}
134 
135 
136  // Member functions
137 
138  //- Return solver name
139  const word& solverName() const
140  {
141  return solverName_;
142  }
143 
144  //- Return initial residual
145  scalar initialResidual() const
146  {
147  return initialResidual_;
148  }
149 
150  //- Return initial residual
151  scalar& initialResidual()
152  {
153  return initialResidual_;
154  }
155 
156 
157  //- Return final residual
158  scalar finalResidual() const
159  {
160  return finalResidual_;
161  }
162 
163  //- Return final residual
164  scalar& finalResidual()
165  {
166  return finalResidual_;
167  }
168 
169 
170  //- Return number of iterations
171  label nIterations() const
172  {
173  return noIterations_;
174  }
175 
176  //- Return number of iterations
177  label& nIterations()
178  {
179  return noIterations_;
180  }
181 
182 
183  //- Has the solver converged?
184  bool converged() const
185  {
186  return converged_;
187  }
188 
189  //- Is the matrix singular?
190  bool singular() const
191  {
192  return singular_;
193  }
194 
195  //- Convergence test
196  bool checkConvergence
197  (
198  const scalar tolerance,
199  const scalar relTolerance
200  );
201 
202  //- Singularity test
203  bool checkSingularity(const scalar residual);
204 
205  //- Print summary of solver performance
206  void print() const;
207  };
208 
209 
210  //- Abstract base-class for lduMatrix solvers
211  class solver
212  {
213  protected:
214 
215  // Protected data
216 
222 
223  //- dictionary of controls
225 
226  //- Maximum number of iterations in the solver
227  label maxIter_;
228 
229  //- Final convergence tolerance
230  scalar tolerance_;
231 
232  //- Convergence tolerance relative to the initial
233  scalar relTol_;
234 
235 
236  // Protected Member Functions
237 
238  //- Read the control parameters from the controlDict_
239  virtual void readControls();
240 
241 
242  public:
243 
244  //- Runtime type information
245  virtual const word& type() const = 0;
246 
247 
248  // Declare run-time constructor selection tables
249 
251  (
252  autoPtr,
253  solver,
254  symMatrix,
255  (
256  const word& fieldName,
257  const lduMatrix& matrix,
261  const dictionary& solverControls
262  ),
263  (
264  fieldName,
265  matrix,
266  interfaceBouCoeffs,
267  interfaceIntCoeffs,
268  interfaces,
269  solverControls
270  )
271  );
272 
274  (
275  autoPtr,
276  solver,
277  asymMatrix,
278  (
279  const word& fieldName,
280  const lduMatrix& matrix,
281  const FieldField<Field, scalar>& interfaceBouCoeffs,
282  const FieldField<Field, scalar>& interfaceIntCoeffs,
283  const lduInterfaceFieldPtrsList& interfaces,
284  const dictionary& solverControls
285  ),
286  (
287  fieldName,
288  matrix,
289  interfaceBouCoeffs,
290  interfaceIntCoeffs,
291  interfaces,
292  solverControls
293  )
294  );
295 
296 
297  // Constructors
298 
299  solver
300  (
301  const word& fieldName,
302  const lduMatrix& matrix,
303  const FieldField<Field, scalar>& interfaceBouCoeffs,
304  const FieldField<Field, scalar>& interfaceIntCoeffs,
305  const lduInterfaceFieldPtrsList& interfaces,
306  const dictionary& solverControls
307  );
308 
309  // Selectors
310 
311  //- Return a new solver
312  static autoPtr<solver> New
313  (
314  const word& fieldName,
315  const lduMatrix& matrix,
316  const FieldField<Field, scalar>& interfaceBouCoeffs,
317  const FieldField<Field, scalar>& interfaceIntCoeffs,
318  const lduInterfaceFieldPtrsList& interfaces,
319  const dictionary& solverControls
320  );
321 
322 
323 
324  // Destructor
325 
326  virtual ~solver()
327  {}
328 
329 
330  // Member functions
331 
332  // Access
333 
334  const word& fieldName() const
335  {
336  return fieldName_;
337  }
338 
339  const lduMatrix& matrix() const
340  {
341  return matrix_;
342  }
343 
345  {
346  return interfaceBouCoeffs_;
347  }
348 
350  {
351  return interfaceIntCoeffs_;
352  }
353 
355  {
356  return interfaces_;
357  }
358 
359 
360  //- Read and reset the solver parameters from the given stream
361  virtual void read(const dictionary&);
362 
363  virtual solverPerformance solve
364  (
365  scalarField& psi,
366  const scalarField& source,
367  const direction cmpt=0
368  ) const = 0;
369 
370  //- Return the matrix norm used to normalise the residual for the
371  // stopping criterion
372  scalar normFactor
373  (
374  const scalarField& psi,
375  const scalarField& source,
376  const scalarField& Apsi,
377  scalarField& tmpField
378  ) const;
379  };
380 
381 
382  //- Abstract base-class for lduMatrix smoothers
383  class smoother
384  {
385  protected:
386 
387  // Protected data
388 
394 
395 
396  public:
397 
398  //- Find the smoother name (directly or from a sub-dictionary)
399  static word getName(const dictionary&);
400 
401  //- Runtime type information
402  virtual const word& type() const = 0;
403 
404 
405  // Declare run-time constructor selection tables
406 
408  (
409  autoPtr,
410  smoother,
411  symMatrix,
412  (
413  const word& fieldName,
414  const lduMatrix& matrix,
415  const FieldField<Field, scalar>& interfaceBouCoeffs,
416  const FieldField<Field, scalar>& interfaceIntCoeffs,
417  const lduInterfaceFieldPtrsList& interfaces
418  ),
419  (
420  fieldName,
421  matrix,
422  interfaceBouCoeffs,
423  interfaceIntCoeffs,
424  interfaces
425  )
426  );
427 
429  (
430  autoPtr,
431  smoother,
432  asymMatrix,
433  (
434  const word& fieldName,
435  const lduMatrix& matrix,
436  const FieldField<Field, scalar>& interfaceBouCoeffs,
437  const FieldField<Field, scalar>& interfaceIntCoeffs,
438  const lduInterfaceFieldPtrsList& interfaces
439  ),
440  (
441  fieldName,
442  matrix,
443  interfaceBouCoeffs,
444  interfaceIntCoeffs,
445  interfaces
446  )
447  );
448 
449 
450  // Constructors
451 
452  smoother
453  (
454  const word& fieldName,
455  const lduMatrix& matrix,
456  const FieldField<Field, scalar>& interfaceBouCoeffs,
457  const FieldField<Field, scalar>& interfaceIntCoeffs,
458  const lduInterfaceFieldPtrsList& interfaces
459  );
460 
461 
462  // Selectors
463 
464  //- Return a new smoother
465  static autoPtr<smoother> New
466  (
467  const word& fieldName,
468  const lduMatrix& matrix,
469  const FieldField<Field, scalar>& interfaceBouCoeffs,
470  const FieldField<Field, scalar>& interfaceIntCoeffs,
471  const lduInterfaceFieldPtrsList& interfaces,
472  const dictionary& solverControls
473  );
474 
475 
476  // Destructor
477 
478  virtual ~smoother()
479  {}
480 
481 
482  // Member functions
483 
484  // Access
485 
486  const word& fieldName() const
487  {
488  return fieldName_;
489  }
490 
491  const lduMatrix& matrix() const
492  {
493  return matrix_;
494  }
495 
497  {
498  return interfaceBouCoeffs_;
499  }
500 
502  {
503  return interfaceIntCoeffs_;
504  }
505 
507  {
508  return interfaces_;
509  }
510 
511 
512  //- Smooth the solution for a given number of sweeps
513  virtual void smooth
514  (
515  scalarField& psi,
516  const scalarField& source,
517  const direction cmpt,
518  const label nSweeps
519  ) const = 0;
520  };
521 
522 
523  //- Abstract base-class for lduMatrix preconditioners
525  {
526  protected:
527 
528  // Protected data
529 
530  //- Reference to the base-solver this preconditioner is used with
531  const solver& solver_;
532 
533 
534  public:
535 
536  //- Find the preconditioner name (directly or from a sub-dictionary)
537  static word getName(const dictionary&);
538 
539  //- Runtime type information
540  virtual const word& type() const = 0;
541 
542 
543  // Declare run-time constructor selection tables
544 
546  (
547  autoPtr,
549  symMatrix,
550  (
551  const solver& sol,
552  const dictionary& solverControls
553  ),
554  (sol, solverControls)
555  );
556 
558  (
559  autoPtr,
561  asymMatrix,
562  (
563  const solver& sol,
564  const dictionary& solverControls
565  ),
566  (sol, solverControls)
567  );
568 
569 
570  // Constructors
571 
573  (
574  const solver& sol
575  )
576  :
577  solver_(sol)
578  {}
579 
580 
581  // Selectors
582 
583  //- Return a new preconditioner
585  (
586  const solver& sol,
587  const dictionary& solverControls
588  );
589 
590 
591  // Destructor
592 
593  virtual ~preconditioner()
594  {}
595 
596 
597  // Member functions
598 
599  //- Read and reset the preconditioner parameters
600  // from the given stream
601  virtual void read(const dictionary&)
602  {}
603 
604  //- Return wA the preconditioned form of residual rA
605  virtual void precondition
606  (
607  scalarField& wA,
608  const scalarField& rA,
609  const direction cmpt=0
610  ) const = 0;
611 
612  //- Return wT the transpose-matrix preconditioned form of
613  // residual rT.
614  // This is only required for preconditioning asymmetric matrices.
615  virtual void preconditionT
616  (
617  scalarField& wT,
618  const scalarField& rT,
619  const direction cmpt=0
620  ) const
621  {
623  (
624  type() +"::preconditionT"
625  "(scalarField& wT, const scalarField& rT, "
626  "const direction cmpt)"
627  );
628  }
629  };
630 
631 
632  // Static data
633 
634  // Declare name of the class and its debug switch
635  ClassName("lduMatrix");
636 
637  //- Large scalar for the use in solvers
638  static const scalar great_;
639 
640  //- Small scalar for the use in solvers
641  static const scalar small_;
642 
643 
644  // Constructors
645 
646  //- Construct given an LDU addressed mesh.
647  // The coefficients are initially empty for subsequent setting.
648  lduMatrix(const lduMesh&);
649 
650  //- Construct as copy
651  lduMatrix(const lduMatrix&);
652 
653  //- Construct as copy or re-use as specified.
654  lduMatrix(lduMatrix&, bool reUse);
655 
656  //- Construct given an LDU addressed mesh and an Istream
657  // from which the coefficients are read
658  lduMatrix(const lduMesh&, Istream&);
659 
660 
661  // Destructor
662 
663  ~lduMatrix();
664 
665 
666  // Member functions
667 
668  // Access to addressing
669 
670  //- Return the LDU mesh from which the addressing is obtained
671  const lduMesh& mesh() const
672  {
673  return lduMesh_;
674  }
675 
676  //- Return the LDU addressing
677  const lduAddressing& lduAddr() const
678  {
679  return lduMesh_.lduAddr();
680  }
681 
682  //- Return the patch evaluation schedule
683  const lduSchedule& patchSchedule() const
684  {
685  return lduAddr().patchSchedule();
686  }
687 
688 
689  // Access to coefficients
690 
691  scalarField& lower();
692  scalarField& diag();
693  scalarField& upper();
694 
695  const scalarField& lower() const;
696  const scalarField& diag() const;
697  const scalarField& upper() const;
698 
699  bool hasDiag() const
700  {
701  return (diagPtr_);
702  }
703 
704  bool hasUpper() const
705  {
706  return (upperPtr_);
707  }
708 
709  bool hasLower() const
710  {
711  return (lowerPtr_);
712  }
713 
714  bool diagonal() const
715  {
716  return (diagPtr_ && !lowerPtr_ && !upperPtr_);
717  }
718 
719  bool symmetric() const
720  {
721  return (diagPtr_ && (!lowerPtr_ && upperPtr_));
722  }
723 
724  bool asymmetric() const
725  {
726  return (diagPtr_ && lowerPtr_ && upperPtr_);
727  }
728 
729 
730  // operations
731 
732  void sumDiag();
733  void negSumDiag();
734 
735  void sumMagOffDiag(scalarField& sumOff) const;
736 
737  //- Matrix multiplication with updated interfaces.
738  void Amul
739  (
740  scalarField&,
741  const tmp<scalarField>&,
744  const direction cmpt
745  ) const;
746 
747  //- Matrix transpose multiplication with updated interfaces.
748  void Tmul
749  (
750  scalarField&,
751  const tmp<scalarField>&,
754  const direction cmpt
755  ) const;
756 
757 
758  //- Sum the coefficients on each row of the matrix
759  void sumA
760  (
761  scalarField&,
764  ) const;
765 
766 
767  void residual
768  (
769  scalarField& rA,
770  const scalarField& psi,
771  const scalarField& source,
772  const FieldField<Field, scalar>& interfaceBouCoeffs,
773  const lduInterfaceFieldPtrsList& interfaces,
774  const direction cmpt
775  ) const;
776 
778  (
779  const scalarField& psi,
780  const scalarField& source,
781  const FieldField<Field, scalar>& interfaceBouCoeffs,
782  const lduInterfaceFieldPtrsList& interfaces,
783  const direction cmpt
784  ) const;
785 
786 
787  //- Initialise the update of interfaced interfaces
788  // for matrix operations
790  (
791  const FieldField<Field, scalar>& interfaceCoeffs,
792  const lduInterfaceFieldPtrsList& interfaces,
793  const scalarField& psiif,
794  scalarField& result,
795  const direction cmpt
796  ) const;
797 
798  //- Update interfaced interfaces for matrix operations
800  (
801  const FieldField<Field, scalar>& interfaceCoeffs,
802  const lduInterfaceFieldPtrsList& interfaces,
803  const scalarField& psiif,
804  scalarField& result,
805  const direction cmpt
806  ) const;
807 
808 
809  template<class Type>
810  tmp<Field<Type> > H(const Field<Type>&) const;
811 
812  template<class Type>
813  tmp<Field<Type> > H(const tmp<Field<Type> >&) const;
814 
815  tmp<scalarField> H1() const;
816 
817  template<class Type>
818  tmp<Field<Type> > faceH(const Field<Type>&) const;
819 
820  template<class Type>
821  tmp<Field<Type> > faceH(const tmp<Field<Type> >&) const;
822 
823 
824  // Member operators
825 
826  void operator=(const lduMatrix&);
827 
828  void negate();
829 
830  void operator+=(const lduMatrix&);
831  void operator-=(const lduMatrix&);
832 
833  void operator*=(const scalarField&);
834  void operator*=(scalar);
835 
836 
837  // Ostream operator
838 
839  friend Ostream& operator<<(Ostream&, const lduMatrix&);
840 };
841 
842 
843 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
844 
845 } // End namespace Foam
846 
847 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
848 
849 #ifdef NoRepository
851 #endif
852 
853 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
854 
855 #endif
856 
857 // ************************ vim: set sw=4 sts=4 et: ************************ //