Flow123d  release_3.0.0-1212-g8801db3
bddcml_wrapper.hpp
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 bddcml_wrapper.hpp
15  * @author Jakub Sistek
16  * @brief
17  */
18 
19 #ifndef la_bddcml_wrapper_h
20 #define la_bddcml_wrapper_h
21 //------------------------------------------------------------------------------
22 #include <boost/numeric/ublas/matrix.hpp>
23 #include <boost/numeric/ublas/vector.hpp>
24 #include <ostream>
25 #include <vector>
26 #include <set>
27 #include <cmath>
28 #include <la/matrix_coo.hpp>
29 
30 #include <system/global_defs.h>
31 #include <system/xio.h>
32 #include "system/sys_profiler.hh"
33 
34 extern "C" {
35  #include <bddcml_interface_c.h>
36 }
37 
38 //------------------------------------------------------------------------------
39 namespace la{
40  class BddcmlWrapper;
41  namespace ublas = boost::numeric::ublas;
42 
43  namespace detail_{
44 
45  template< typename DOFHOLDER, typename MAPTYPE >
46  int getLocalNumber( DOFHOLDER * np, MAPTYPE & map ) {
47 
48  unsigned globalIndex = np -> giveId( );
49 
50  typename MAPTYPE::iterator pos = map.find( globalIndex );
51  OLD_ASSERT( pos != map.end(),
52  "Cannot remap node index %d to local indices. \n ", globalIndex );
53  int localIndexInt = pos -> second;
54 
55  return localIndexInt;
56  }
57  }
58 
59 }
60 
61 //------------------------------------------------------------------------------
62 /** \brief Multilevel BDDC based linear system solver
63  * \details This class implements solution of a system of linear equations
64  * by multilevel BDDC implementation in BDDCML solver.
65  *
66  * The class uses assembly of the matrix into
67  * the coordinate matrix format implemented in matrix_coo.hpp file.
68  *
69  * prepareMatAssembly - (optional) - reserves memory for sparse coordinate matrix
70  * insertToMatrix - called many times, inserts values into coordinate matrix
71  * insertToRhs - called many times, inserts values into RHS.
72  * If matrix is created and assembled, structure of the vector is deduced from it.
73  * Otherwise, it is created as empty.
74  * applyConstraints - currently only passes list of boundary conditions into the BDDCML solver,
75  * rather then modifying the matrix and rhs
76  * - as a result, no nonsymmetry is forced to SPD and symmetric indefinite problems by constraints,
77  * and these can be accelerated and be more memory efficient by specifying proper matrixType
78  * solveSystem - solves the system
79  */
81 {
82 public:
84  GENERAL, //!< general (default),
85  SPD, //!< symmetric positive definite,
86  SYMMETRICGENERAL, //!< general symmetric ( e.g. saddle point ),
87  SPD_VIA_SYMMETRICGENERAL //!< (advanced) interface problem is symm. pos. def.,
88  //!< but subdomain problems are general symmetric ( e.g. saddle point )
89  };
90  typedef enum matrixTypeEnum MatrixType;
91 
92 private:
93  typedef ublas::matrix< double > SubMat_;
94  typedef ublas::vector< double > SubVec_;
97 
98 public:
99  //! Constructor with
100  //! total number of dofs,
101  //! subdomain number of dofs,
102  //! (optional) type of matrix,
103  //! (optional) communicator ( MPI_COMM_WORLD is the default ),
104  //! (optional) number of subdomains ( if not given, use number of processes ),
105  BddcmlWrapper( const unsigned numDofs,
106  const unsigned numDofsSub,
107  const MatrixType matrixType = GENERAL,
108  const MPI_Comm comm = MPI_COMM_WORLD,
109  int numSubLoc = 1 );
110 
111  //! Destructor frees the BDDCML structures
112  ~BddcmlWrapper();
113 
114  //! Load raw data about mesh
115  void loadRawMesh( const int nDim, const int numNodes, const int numDofs,
116  const std::vector<int> & inet,
117  const std::vector<int> & nnet,
118  const std::vector<int> & nndf,
119  const std::vector<int> & isegn,
120  const std::vector<int> & isngn,
121  const std::vector<int> & isvgvn,
122  const std::vector<double> & xyz,
123  const std::vector<double> & element_data,
124  const int meshDim = 0 );
125 
126  //! Prepare assembly of matrix - reserve space in underlying coordinate matrix object
127  void prepareMatAssembly( unsigned numElements, unsigned elMatSize );
128 
129  //! insert submatrix to system matrix - into coordinate matrix object
130  void insertToMatrix( const SubMat_ & eStiff,
131  const VecUint_ & rowIndices,
132  const VecUint_ & colIndices );
133 
134  //! Finalize assembly of matrix
135  void finishMatAssembly( );
136 
137  //! insert subvector to system rhs vector
138  void insertToRhs( const SubVec_ & eForce, const VecUint_ & dofIndices );
139 
140  //! Outputs
141  void writeMatrix( std::ostream & out ) { coo_.write( out ); }
142 
143  //! Compatibility call which does absolutely nothing
144  void finishAssembly( ) { }
145 
146  //! Apply constraints
147  //! Currently, scalar is not used in the function and it is kept only for compatibility with constraint functors.
148  //! Proper scalar is enforced within the BDDC solver.
149  void applyConstraints( ConstraintVec_ & constraints, const double factor, const double scalar );
150 
151  //! Set any DOF to zero
152  //!
153  //! \param[in] index Global index of DOF to fix
154  //! \param[in] scalar Value to put on diagonal (default 1.)
155  void fixDOF( const unsigned index, const double scalar = 1. );
156 
157  //! Insert a value to diagonal for creating weights in BDDC
158  //!
159  //! \param[in] index Global index of DOF
160  //! \param[in] value Value to put on diagonal
161  void insertToDiagonalWeights( const int & index, const double & value );
162 
163  //! Finalize assembly of diagonal
164  void finishDiagAssembly( );
165 
166  //! Solve the system
167  void solveSystem( double tol = 1.e-7, //!< tolerance on relative residual ||res||/||rhs||
168  int numLevels = 2, //!< number of levels
169  std::vector<int> * numSubAtLevels = NULL, //!< number of subdomains at levels
170  int verboseLevel = 0, //!< level of verbosity of BDDCML library
171  //!< ( 0 - only fatal errors reported,
172  //!< 1 - mild output,
173  //!< 2 - detailed output )
174  int maxIt = 1000, //!< maximum number of iterations
175  int ndecrMax = 30, //!< maximum number of iterations with non-decreasing residual
176  //!< ( used to stop diverging process )
177  bool use_adaptive = false //!< if adaptive constraints should be used by BDDCML
178  );
179 
180  //! Get norm of right-hand side
181  double normRhs( ) { return normRhs_; }
182 
183  //! Get norm of solution
184  double normSol( ) { return normSol_; }
185 
186  //! Give number of iteration
187  int giveNumIterations() { return numIter_ ; }
188 
189  //! Give reason of convergence
190  int giveConvergedReason() { return convergedReason_ ; }
191 
192  //! Give estimate of condition number
193  double giveCondNumber() { return condNumber_ ; }
194 
195  //! Fill matrix with zeros
196  void clearMatrix( ) {
197  // remember length
198  unsigned length = coo_.nnz( );
199  coo_.clear( );
200  coo_.prepareAssembly( length );
201  }
202  //! Fill RHS with zeros
203  void clearRhs( ) {
204  std::fill( rhsVec_.begin( ), rhsVec_.end( ), 0. );
205  normRhs_ = 0.;
206  }
207 
208  //! Blank arrays for boundary conditions
209  void clearBC( ) {
210  std::fill( ifix_.begin(), ifix_.end(), 0 );
211  std::fill( fixv_.begin(), fixv_.end(), 0. );
212  }
213 
214  //! Fill solution with zeros
215  void clearSol( )
216  { std::fill( sol_.begin( ), sol_.end( ), 0. );
217  normSol_ = 0.;
218  }
219 
220  //! Destroy matrix
221  void destroyMatrix( ) {
222  coo_.clear( );
223  }
224  //! Destroy RHS
225  void destroyRhs( ) { rhsVec_.clear( ); }
226  //! Destroy boundary conditions
227  void destroyBC( ) {
228  ifix_.clear( );
229  fixv_.clear( );
230  }
231  //! Destroy solution
232  void destroySol( ) { sol_.clear( ); }
233 
234  //! Give access to the solution vector
235  template<typename VEC1,typename VEC2>
236  void giveSolution( const VEC1 & dofIndices, VEC2 & result ) const;
237 
238 private:
239  const int numDofs_; //!< number of total dofs
240  const int numDofsSub_; //!< number of subdomain dofs
241  const MatrixType matrixType_; //!< type of matrix
242  const MPI_Comm comm_; //!< communicator
243  int numSubLoc_; //!< number of subdomains on process
244 
245  int nProc_; //!< number of processors in communicator
246  int rank_; //!< ID of processors in communicator
247 
248  int nDim_; //!< space dimension
249  int meshDim_; //!< mesh dimension, may differ from space dimension - e.g. 2 for shells in 3D
250 
251  int numElemSub_; //!< number of elements in subdomain
252  int numElem_; //!< total number of elements
253  int numNodesSub_; //!< number of nodes in subdomain
254  int numNodes_; //!< total number of nodes
255 
256  bool meshLoaded_; //! safety state variable to warn if solver is called without loading subdomain mesh
257  std::vector<int> inet_; //!< indices of nodes on elements (linearized connectivity) - numbering w.r.t. local subdomain numbering
258  std::vector<int> nnet_; //!< numbers of nodes on elements
259  //!< (determines chunks of nodes for each element in inet),
260  //!< e.g. [ 8, 8, 8, ... ] for linear hexahedra
261  std::vector<int> nndf_; //!< number of nodal degrees of freedom ( entry for each node ),
262  //!< e.g. [ 3, 3, 3, ... ] for 3D solids
263  std::vector<double> xyz_; //!< linearized array of coordinates - all x values, all y values, all z values ( if present )
264  //!< i.e. [ x_0 x_1 x_2 ... x_(numNodes - 1) y_0 y_1 y_2 ... y_(numNodes - 1) z_0 z_1 z_2 ... z_(numNodes - 1) ]
265  std::vector<double> element_data_; //!< array with data on elements, here used for permeability for each element
266 
267  std::vector<int> isngn_; //!< Indices of Subdomain Nodes in Global Numbering
268  std::vector<int> isvgvn_; //!< Indices of Subdomain Variables in Global Variable Numbering ( i.e. dof mapping )
269  std::vector<int> isegn_; //!< Indices of Subdomain Elements in Global Numbering
270 
271  typedef std::map<unsigned,unsigned> Global2LocalMap_; //! type for storage of global to local map
272  Global2LocalMap_ global2LocalDofMap_; //!< map from global dof indices to subdomain local indices
273 
274 
275  std::vector<int> ifix_; //!< indices of fixed boundary conditions - ( 0 for free variable, 1 for constrained dof )
276  std::vector<double> fixv_; //!< values of fixed boundary conditions at places of ifix_ - values outside fixed dofs are ignored
277 
278  la::MatrixCoo<int,double> coo_; //!< matrix in coordinate format (COO)
279  bool isMatAssembled_; //!< true if matrix is assembled
280  la::MatrixCoo<int,double> diagWeightsCoo_; //!< diagonal of subdomain matrix in sparse format and global numbers
281  bool isDiagAssembled_; //!< true if diagonal is assembled
282 
283  std::vector<double> rhsVec_; //!< vector with RHS values restricted to subdomain, i.e. values are repeated at shared nodes
284  double normRhs_; //!< norm of the global right-hand side
285 
286  std::vector<double> sol_; //!< distributed solution vector
287  double normSol_; //!< norm of the global solution
288 
289  int numIter_; //!< required number of iterations
290  int convergedReason_; //!< reason of convergence:
291  //!< 0 - converged relative residual, desirable state
292  //!< -1 - reached maximal number of iterations
293  //!< -2 - reached maximal number of iterations with non-decreasing residual
294  double condNumber_; //!< estimated condition number of the preconditioned operator
295  //!< only provided for SPD problems solved by PCG ( -1.0 if not meaningful )
296 
297 };
298 
299 //-----------------------------------------------------------------------------
300 
301 #include "bddcml_wrapper.ipp"
302 
303 //-----------------------------------------------------------------------------
304 
305 #endif
double giveCondNumber()
Give estimate of condition number.
void destroyRhs()
Destroy RHS.
double normSol()
Get norm of solution.
double normRhs()
Get norm of right-hand side.
std::vector< double > fixv_
values of fixed boundary conditions at places of ifix_ - values outside fixed dofs are ignored ...
int MPI_Comm
Definition: mpi.h:141
int giveConvergedReason()
Give reason of convergence.
int numIter_
required number of iterations
int nProc_
number of processors in communicator
int rank_
ID of processors in communicator.
std::vector< int > inet_
safety state variable to warn if solver is called without loading subdomain mesh
std::vector< int > nndf_
std::vector< double > xyz_
int numNodes_
total number of nodes
void destroySol()
Destroy solution.
Global2LocalMap_ global2LocalDofMap_
type for storage of global to local map
enum matrixTypeEnum MatrixType
la::MatrixCoo< int, double > coo_
matrix in coordinate format (COO)
double normSol_
norm of the global solution
int nDim_
space dimension
I/O functions with filename storing, able to track current line in opened file. All standard stdio fu...
const MatrixType matrixType_
type of matrix
static constexpr bool value
Definition: json.hpp:87
int numElem_
total number of elements
int numSubLoc_
number of subdomains on process
#define OLD_ASSERT(...)
Definition: global_defs.h:131
std::vector< double > sol_
distributed solution vector
Global macros to enhance readability and debugging, general constants.
std::vector< int > isvgvn_
Indices of Subdomain Variables in Global Variable Numbering ( i.e. dof mapping )
std::vector< int > isegn_
Indices of Subdomain Elements in Global Numbering.
int numElemSub_
number of elements in subdomain
void clearSol()
Fill solution with zeros.
int getLocalNumber(DOFHOLDER *np, MAPTYPE &map)
std::vector< double > rhsVec_
vector with RHS values restricted to subdomain, i.e. values are repeated at shared nodes ...
bool isMatAssembled_
true if matrix is assembled
std::vector< int > ifix_
indices of fixed boundary conditions - ( 0 for free variable, 1 for constrained dof ) ...
symmetric positive definite,
general symmetric ( e.g. saddle point ),
void clearMatrix()
Fill matrix with zeros.
bool isDiagAssembled_
true if diagonal is assembled
void destroyBC()
Destroy boundary conditions.
std::vector< int > isngn_
Indices of Subdomain Nodes in Global Numbering.
std::vector< unsigned > VecUint_
int meshDim_
mesh dimension, may differ from space dimension - e.g. 2 for shells in 3D
int numNodesSub_
number of nodes in subdomain
std::map< unsigned, unsigned > Global2LocalMap_
std::vector< int > nnet_
#define MPI_COMM_WORLD
Definition: mpi.h:123
int giveNumIterations()
Give number of iteration.
const MPI_Comm comm_
communicator
ublas::vector< double > SubVec_
void clearRhs()
Fill RHS with zeros.
std::vector< double > element_data_
array with data on elements, here used for permeability for each element
void finishAssembly()
Compatibility call which does absolutely nothing.
la::MatrixCoo< int, double > diagWeightsCoo_
diagonal of subdomain matrix in sparse format and global numbers
void clearBC()
Blank arrays for boundary conditions.
void writeMatrix(std::ostream &out)
Outputs.
const int numDofsSub_
number of subdomain dofs
void destroyMatrix()
Destroy matrix.
Multilevel BDDC based linear system solver.
const int numDofs_
number of total dofs
ublas::matrix< double > SubMat_
double normRhs_
norm of the global right-hand side
std::vector< std::pair< unsigned, double > > ConstraintVec_