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