Flow123d  master-f44eb46
schur.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file schur.hh
15  * @brief Assembly explicit Schur complement for the given linear system.
16  * Provides method for resolution of the full original vector of unknowns.
17  */
18 
19 #ifndef LA_SCHUR_HH_
20 #define LA_SCHUR_HH_
21 
22 #include <petscmat.h> // for Mat, _p_Mat
23 #include "la/linsys_PETSC.hh" // for LinSys_PETSC
24 #include "petscistypes.h" // for IS, _p_IS
25 #include "petscvec.h" // for Vec, _p_Vec
26 
27 class Distribution;
28 class LinSys;
29 
30 namespace Input { class Record; }
31 
32 
33 /**
34  * @brief Schur complement class for a PETSC based linear system
35  *
36  * Linear system consists of Matrix, inversion matrix, RHS, solution and pointer to Complement.
37  * It provides methods for:
38  * - set complement system
39  * - create distribution of complement system
40  * - create inversion matrix of system
41  * - form Schur complement and rhs
42  * - solve and resolve system
43  *
44  * Usage example:
45  * @CODE
46  * SchurComplement * schur = new SchurComplement(is, &distr);
47  *
48  * ... // allocation and assembly of matrix
49  *
50  * LinSys_PETSC * ls = new LinSys_PETSC( schur->make_complement_distribution() );
51  * schur->set_complement( ls );
52  * schur->solve();
53  * @ENDCODE
54  *
55  * Input record is passed to the complement system.
56  */
57 
58 typedef enum SchurState {
59  created, // need creation of all temporal matrixes ...
60  formed, // Schur complement ready to solve
61  solved // Resolved original Schur system
63 
64 typedef class SchurComplement : public LinSys_PETSC {
65 public:
66  /**
67  * Constructor
68  *
69  * Gets linear system with original matrix A and creates its inversion (IA matrix)
70  *
71  * ia - PETSC indexset for the eliminated block.
72  * ib - PETSc indexset for the Schur complement, complementary by default.
73  */
74  SchurComplement(Distribution *ds, IS ia, IS ib = nullptr);
75 
76  /**
77  * Copy constructor.
78  */
80 
81  void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override;
82 
83  /**
84  * Sets specific parameters defined by user in input file and used to calculate. Call set_from_input of complement
85  */
86  void set_from_input(const Input::Record in_rec) override;
87 
88  /**
89  * Returns pointer to LinSys object representing the schur complement.
90  */
91  LinSys *get_system() const {return (Compl);}
92 
93  /**
94  * Returns distribution of the original system (solved by class SchurComplement).
95  */
96  Distribution *get_distribution() const {return (ds_);}
97 
98  /**
99  * Destructor. In particular it also delete complement linear system if it was passed in
100  * through the @p set_complement() method.
101  */
103 
104  /** Compute only right hand side.
105  * This is useful when you change only rhs of the original system.
106  * TODO: We should ask original system if the matrix has changed (using LazyDependency) and
107  * possibly call only form_rhs, then this can be protected
108  */
109  void form_rhs();
110  /// Set complement LinSys object.
111  void set_complement(LinSys_PETSC *ls);
112  /// get distribution of complement object if complement is defined
114  /// get precision of solving
115  double get_solution_precision() override;
116 
117  /**
118  * Solve the system.
119  *
120  * - If matrix and/or rhs is changed the Schur complement is formed.
121  * - The inner linear solver is called for the Schur complement
122  * - Resolve is called to reconstruct eliminated part of the solution vector.
123  */
124  LinSys::SolveInfo solve() override;
125 
126  /**
127  * Only resolve the system with current solution vector. This is necessary for nonlinear solvers.
128  * - If matrix and/or rhs is changed the Schur complement is formed.
129  * - Resolve is called to reconstruct eliminated part of the solution vector.
130  */
131  void resolve();
132 
133  /**
134  * The solve or resolve must be called prior to computing the residual.
135  */
136  double compute_residual() override;
137  int loc_size_A, loc_size_B; // loc size of the A and B block
138  IS IsA, IsB; // parallel index sets of the A and B block
139 
140 protected:
141  /// create IA matrix
143 
144  void form_schur();
145 
146 
147 
148  Mat A; // Submatrix of matrix_ contains only data given by IsA parallel index set
149  Mat IA; // Inverse of block A
150 
151  Mat B, Bt; // B and B' block (could be different from real B transpose)
152  Mat C; // Sub matrix.
153  Mat xA; // Bt*IA*B
154  Mat IAB; // reconstruction matrix IA * B
155 
156  Vec RHS1, RHS2; // A and B - part of the RHS
157  Vec Sol1, Sol2; // A and B part of solution
158  VecScatter rhs1sc, rhs2sc; // scatter to parts of rhs
159  VecScatter sol1sc, sol2sc; // scatter to parts of solution
160 
161  SchurState state; // object internal state
162  int orig_lsize; ///< Size of local vector part of original system
163 
164  LinSys_PETSC *Compl; // Schur complement system: (C - B' IA B) * Sol2 = (B' * IA * RHS1 - RHS2)
165 
166  Distribution *ds_; // Distribution of B block
168 
169 #endif /* LA_SCHUR_HH_ */
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
friend class SchurComplement
Definition: linsys.hh:91
SchurState state
Definition: schur.hh:161
Distribution * ds_
Definition: schur.hh:166
VecScatter rhs2sc
Definition: schur.hh:158
void resolve()
Definition: schur.cc:418
Distribution * make_complement_distribution()
get distribution of complement object if complement is defined
Definition: schur.cc:294
void set_from_input(const Input::Record in_rec) override
Definition: schur.cc:138
void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override
Definition: schur.cc:281
VecScatter sol2sc
Definition: schur.hh:159
double compute_residual() override
Definition: schur.cc:434
void set_complement(LinSys_PETSC *ls)
Set complement LinSys object.
Definition: schur.cc:287
void form_schur()
Definition: schur.cc:166
void form_rhs()
Definition: schur.cc:258
LinSys::SolveInfo solve() override
Definition: schur.cc:387
int loc_size_A
Definition: schur.hh:137
~SchurComplement()
Definition: schur.cc:446
Distribution * get_distribution() const
Definition: schur.hh:96
int orig_lsize
Size of local vector part of original system.
Definition: schur.hh:162
double get_solution_precision() override
get precision of solving
Definition: schur.cc:378
VecScatter rhs1sc
Definition: schur.hh:158
VecScatter sol1sc
Definition: schur.hh:159
void create_inversion_matrix()
create IA matrix
Definition: schur.cc:301
LinSys * get_system() const
Definition: schur.hh:91
int loc_size_B
Definition: schur.hh:137
LinSys_PETSC * Compl
Definition: schur.hh:164
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
Abstract linear system class.
Definition: balance.hh:40
SchurState
Schur complement class for a PETSC based linear system.
Definition: schur.hh:58
@ formed
Definition: schur.hh:60
@ created
Definition: schur.hh:59
@ solved
Definition: schur.hh:61
SchurComplement SchurComplement
Definition: linsys.hh:91