Flow123d  master-f44eb46
linsys_BDDC.cc
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 linsys_BDDC.cc
15  * @brief Solver based on Multilevel BDDC - using corresponding class of OpenFTL package
16  * @author Jakub Sistek
17  */
18 
19 #include <mpi.h>
20 #include "config.h"
21 
22 
23 // need BDDCML wrapper
24 #ifdef FLOW123D_HAVE_BDDCML
25  #include <map>
26  #include "la/bddcml_wrapper.hh"
27 #endif // FLOW123D_HAVE_BDDCML
28 
29 #include "la/linsys.hh"
30 #include "la/linsys_BDDC.hh"
31 #include "system/sys_profiler.hh"
32 #include "input/input_exception.hh"
33 
34 
35 
36 namespace it = Input::Type;
37 
38 
40  return it::Record("Bddc", "BDDCML (Balancing Domain Decomposition by Constraints - Multi-Level) solver settings.")
42  .declare_key("r_tol", it::Double(0.0, 1.0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
43  "If not, we use the value 1.0e-7."),
44  "Residual tolerance relative to the initial error.")
45  .declare_key("max_it", it::Integer(0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
46  "If not, we use the value 1000."),
47  "Maximum number of outer iterations of the linear solver.")
48 
49  .declare_key("max_nondecr_it", it::Integer(0), it::Default("30"),
50  "Maximum number of iterations of the linear solver with non-decreasing residual.")
51  .declare_key("number_of_levels", it::Integer(0), it::Default("2"),
52  "Number of levels in the multilevel method (=2 for the standard BDDC).")
53  .declare_key("use_adaptive_bddc", it::Bool(), it::Default("false"),
54  "Use adaptive selection of constraints in BDDCML.")
55  .declare_key("bddcml_verbosity_level", it::Integer(0,2), it::Default("0"),
56  "Level of verbosity of the BDDCML library:\n\n - 0 - no output,\n - 1 - mild output,\n - 2 - detailed output.")
57  .close();
58 }
59 
60 
62 
63 
64 
66  const bool swap_sign )
67  : LinSys( rows_ds ),
68  swap_sign_(swap_sign)
69 {
70 #ifdef FLOW123D_HAVE_BDDCML
71  // from the point of view of assembly, BDDC linsys is in the ADD state
73 #else
74  throw Input::ExcInputMessage() << EI_Message("Unsupported solver BDDC. Compiled without support for the BDDCML solver.");
75 #endif // FLOW123D_HAVE_BDDCML
76 }
77 
78 
79 void LinSys_BDDC::set_tolerances(double r_tol, FMT_UNUSED double a_tol, FMT_UNUSED double d_tol, unsigned int max_it)
80 {
81  if (! in_rec_.is_empty()) {
82  // input record is set
83  r_tol_ = in_rec_.val<double>("r_tol", r_tol);
84  // BDDC does not use a_tol_
85  a_tol_ = 0.01 * r_tol_;
86  max_it_ = in_rec_.val<unsigned int>("max_it", max_it);\
87  }
88 }
89 
90 
92  const int nDim, const int numNodes, FMT_UNUSED const int numDofs,
93  const std::vector<int> & inet,
94  const std::vector<int> & nnet,
95  const std::vector<int> & nndf,
96  const std::vector<int> & isegn,
97  const std::vector<int> & isngn,
98  const std::vector<int> & isvgvn,
99  const std::vector<double> & xyz,
100  const std::vector<double> & element_permeability,
101  const int meshDim )
102 {
103 #ifdef FLOW123D_HAVE_BDDCML
104 
105  uint num_of_local_subdomains = 1;
106  bddcml_ = new Bddcml_( size_,
107  nndf.size(),
108  matrix_type,
109  rows_ds_->get_comm(),
110  num_of_local_subdomains);
111 
112  // prepare space for local solution
113  locSolution_.resize( nndf.size() );
114  chkerr( VecCreateSeqWithArray( PETSC_COMM_SELF, 1, nndf.size(), &(locSolution_[0]), &locSolVec_ ));
115 
116  // simply pass the data to BDDCML solver
117  isngn_.resize(isngn.size());
118  std::copy( isngn.begin(), isngn.end(), isngn_.begin() );
119  ASSERT_PERMANENT_EQ( numDofs, static_cast<int>(size_) ).error( "Global problem size mismatch!" );
120 
121  bddcml_ -> loadRawMesh( nDim, numNodes, inet, nnet, nndf, isegn, isngn, isvgvn, xyz, element_permeability, meshDim );
122 
123  // create a map for BDDCML to PETSc vector
124  PetscErrorCode ierr;
125 
126  // local index set
127  PetscInt numDofsSubInt = static_cast<int>( isngn_.size( ) );
129 
130  ISLocalToGlobalMapping subdomainMapping;
131  ierr = ISLocalToGlobalMappingCreate( comm_, 1, numDofsSubInt, &(idx[0]), PETSC_COPY_VALUES, &subdomainMapping ); CHKERRV( ierr );
132 
133  IS subdomainIndexSet;
134  IS from;
135  ierr = ISCreateStride( PETSC_COMM_SELF, numDofsSubInt, 0, 1, &subdomainIndexSet );
136  ierr = ISLocalToGlobalMappingApplyIS( subdomainMapping, subdomainIndexSet, &from );
137 
138 
139  // create scatter from parallel PETSc vector to local indices of subdomain
140  ierr = VecScatterCreate( solution_, from, locSolVec_, subdomainIndexSet, &VSpetscToSubScatter_ ); CHKERRV( ierr );
141  chkerr(ISDestroy( &subdomainIndexSet ));
142  chkerr(ISDestroy( &from ));
143 
144  //double * locSolVecArray;
145  //ierr = VecGetArray( locSolVec_, &locSolVecArray ); CHKERRV( ierr );
146  //std::copy( locSolution_.begin(), locSolution_.end(), locSolVecArray );
147  //ierr = VecRestoreArray( locSolVec_, &locSolVecArray ); CHKERRV( ierr );
148 
149  // scatter local solutions back to global one
150  VecScatterBegin( VSpetscToSubScatter_, locSolVec_, solution_, INSERT_VALUES, SCATTER_REVERSE );
151  VecScatterEnd( VSpetscToSubScatter_, locSolVec_, solution_, INSERT_VALUES, SCATTER_REVERSE );
152 #endif // FLOW123D_HAVE_BDDCML
153 }
154 
155 void LinSys_BDDC::mat_set_values( int nrow, int *rows, int ncol, int *cols, double *vals )
156 {
157 #ifdef FLOW123D_HAVE_BDDCML
158  namespace ublas = boost::numeric::ublas;
159 
160  std::vector< unsigned > myRows( nrow );
161  std::vector< unsigned > myCols( ncol );
162  ublas::matrix< double > mat( nrow, ncol );
163 
164  std::copy( &(rows[0]), &(rows[nrow]), myRows.begin() );
165  std::copy( &(cols[0]), &(cols[ncol]), myCols.begin() );
166 
167  for ( int i = 0; i < nrow; i++ ) {
168  for ( int j = 0; j < ncol; j++ ) {
169  mat( i, j ) = vals[i*ncol + j];
170  }
171  }
172  if (swap_sign_) {
173  mat = -mat;
174  }
175 
176  bddcml_ -> insertToMatrix( mat, myRows, myCols );
177 #endif // FLOW123D_HAVE_BDDCML
178 }
179 
180 void LinSys_BDDC::rhs_set_values( int nrow, int *rows, double *vals)
181 {
182 #ifdef FLOW123D_HAVE_BDDCML
183  namespace ublas = boost::numeric::ublas;
184 
185  std::vector< unsigned > myRows( nrow );
186  ublas::vector< double > vec( nrow );
187 
188  std::copy( &(rows[0]), &(rows[nrow]), myRows.begin() );
189 
190  for ( int i = 0; i < nrow; i++ ) {
191  vec( i ) = vals[i];
192  }
193  if (swap_sign_) {
194  vec = -vec;
195  }
196 
197  bddcml_ -> insertToRhs( vec, myRows );
198 #endif // FLOW123D_HAVE_BDDCML
199 }
200 
201 void LinSys_BDDC::diagonal_weights_set_value( int global_index, double value )
202 {
203 #ifdef FLOW123D_HAVE_BDDCML
204  bddcml_ -> insertToDiagonalWeights( global_index, value );
205 #endif // FLOW123D_HAVE_BDDCML
206 }
207 
209 {
210 #ifdef FLOW123D_HAVE_BDDCML
211  bddcml_ -> clearMatrix( );
212 #endif // FLOW123D_HAVE_BDDCML
213  return 0;
214 }
215 
217 {
218 #ifdef FLOW123D_HAVE_BDDCML
219  bddcml_ -> clearRhs( );
220 #endif // FLOW123D_HAVE_BDDCML
221  return 0;
222 }
223 
225 {
226 #ifdef FLOW123D_HAVE_BDDCML
227  bddcml_ -> finishMatAssembly( );
228 #endif // FLOW123D_HAVE_BDDCML
229 }
230 
231 void LinSys_BDDC::apply_constrains( double scalar )
232 {
233 #ifdef FLOW123D_HAVE_BDDCML
234  bddcml_ -> applyConstraints( constraints_, 1., scalar );
235 #endif // FLOW123D_HAVE_BDDCML
236 }
237 
238 LinSys::SolveInfo LinSys_BDDC::solve() // ! params are not currently used
239 {
240 #ifdef FLOW123D_HAVE_BDDCML
241  std::vector<int> * numSubAtLevels = NULL; //!< number of subdomains at levels
242  START_TIMER("BDDC linear solver");
243 
244  {
245 
246  START_TIMER("BDDC linear iteration");
248 
249 
250  LogOut().fmt("BDDCML converged reason: {} ( 0 means OK ) \n", bddcml_ -> giveConvergedReason() );
251  LogOut().fmt("BDDCML converged in {} iterations. \n", bddcml_ -> giveNumIterations() );
252  LogOut().fmt("BDDCML estimated condition number is {} \n", bddcml_ -> giveCondNumber() );
253  ADD_CALLS(bddcml_ -> giveNumIterations());
254  }
255 
256  // download local solution
257  bddcml_ -> giveSolution( isngn_, locSolution_ );
258 
259 
260  double * locSolVecArray;
261  VecGetArray( locSolVec_, &locSolVecArray );
262  std::copy( locSolution_.begin(), locSolution_.end(), locSolVecArray );
263  VecRestoreArray( locSolVec_, &locSolVecArray );
264 
265  // scatter local solutions back to global one
266  VecScatterBegin( VSpetscToSubScatter_, locSolVec_, solution_, INSERT_VALUES, SCATTER_REVERSE );
267  VecScatterEnd( VSpetscToSubScatter_, locSolVec_, solution_, INSERT_VALUES, SCATTER_REVERSE );
268 
269  // upper bound on the residual error
271 
272  return LinSys::SolveInfo(bddcml_ -> giveConvergedReason(), bddcml_ -> giveNumIterations());
273 #else
274  return LinSys::SolveInfo(0,0);
275 #endif // FLOW123D_HAVE_BDDCML
276 
277 }
278 
280 {
281  LinSys::set_from_input(in_rec);
282 
283  // BDDCML specific parameters
284  max_nondecr_it_ = in_rec.val<int>("max_nondecr_it");
285  number_of_levels_ = in_rec.val<int>("number_of_levels");
286  use_adaptive_bddc_ = in_rec.val<bool>("use_adaptive_bddc");
287  bddcml_verbosity_level_ = in_rec.val<int>("bddcml_verbosity_level");
288 }
289 
291 {
292 #ifdef FLOW123D_HAVE_BDDCML
293  isngn_.clear();
294  locSolution_.clear();
295 
296  chkerr(VecDestroy( &locSolVec_ ));
297 
298  chkerr(VecScatterDestroy( &VSpetscToSubScatter_ ));
299 
300  delete bddcml_;
301 #endif // FLOW123D_HAVE_BDDCML
302 }
303 
304 // construct global solution
305 //void LinSys_BDDC::gatherSolution_( )
306 //{
307 //#ifdef FLOW123D_HAVE_BDDCML
308 // int ierr;
309 //
310 // // merge solution on root
311 // int rank;
312 // MPI_Comm_rank( comm_, &rank );
313 // int nProc;
314 // MPI_Comm_size( comm_, &nProc );
315 //
316 // globalSolution_.resize( size_ );
317 // std::vector<double> locSolutionNeib;
318 // if ( rank == 0 ) {
319 // // merge my own data
320 // for ( unsigned int i = 0; i < isngn_.size(); i++ ) {
321 // int ind = isngn_[i];
322 // globalSolution_[ind] = locSolution_[i];
323 // }
324 // for ( int iProc = 1; iProc < nProc; iProc++ ) {
325 // // receive length
326 // int length;
327 // MPI_Status status;
328 // ierr = MPI_Recv( &length, 1, MPI_INT, iProc, iProc, comm_, &status );
329 //
330 // // receive indices
331 // std::vector<int> inds( length );
332 // ierr = MPI_Recv( &(inds[0]), length, MPI_INT, iProc, iProc, comm_, &status );
333 //
334 // // receive local solution
335 // locSolutionNeib.resize( length );
336 // ierr = MPI_Recv( &(locSolutionNeib[0]), length, MPI_DOUBLE, iProc, iProc, comm_, &status );
337 //
338 // // merge data
339 // for ( int i = 0; i < length; i++ ) {
340 // int ind = inds[i];
341 // globalSolution_[ind] = locSolutionNeib[i];
342 // }
343 // }
344 // }
345 // else {
346 // // send my solution to root
347 // int length = isngn_.size();
348 // ierr = MPI_Send( &length, 1, MPI_INT, 0, rank, comm_ );
349 // ierr = MPI_Send( &(isngn_[0]), length, MPI_INT, 0, rank, comm_ );
350 // ierr = MPI_Send( &(locSolution_[0]), length, MPI_DOUBLE, 0, rank, comm_ );
351 // }
352 // // broadcast global solution from root
353 // ierr = MPI_Bcast( &(globalSolution_[0]), globalSolution_.size(), MPI_DOUBLE, 0, comm_ );
354 //#endif // FLOW123D_HAVE_BDDCML
355 //}
356 
358 {
359  double bnorm=0.0;
360  VecNorm(locSolVec_, NORM_2, &bnorm);
361 
362  return max(a_tol_, r_tol_*bnorm);
363 }
364 
365 
366 void LinSys_BDDC::print_matrix(std::ostream& out)
367 {
368  int rank;
370 
371  if(rank == 0){
372  out << "zzz = [\n";
373 
374  bddcml_->writeMatrix(out);
375  out << "];\n"
376  << "zzz(:,1:2) = zzz(:,1:2) + ones(size(zzz),2);\n" // fix matlab indices (+1)
377  << "matrix_bddc = spconvert(zzz);\n";
378  }
379 }
#define ASSERT_PERMANENT_EQ(a, b)
Definition of comparative assert macro (EQual)
Definition: asserts.hh:329
MPI_Comm get_comm() const
Returns communicator.
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
bool is_empty() const
Definition: accessors.hh:365
const Ret val(const string &key) const
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default read_time(const std::string &description)
The factory function to make an default value that will be specified at the time when a key will be r...
Definition: type_record.hh:97
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Class for declaration of the integral input data.
Definition: type_base.hh:483
Record type proxy class.
Definition: type_record.hh:182
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
int bddcml_verbosity_level_
Definition: linsys_BDDC.hh:129
void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals) override
Definition: linsys_BDDC.cc:155
Vec locSolVec_
local solution PETSc vector - sequential
Definition: linsys_BDDC.hh:138
std::vector< int > isngn_
indices of subdomain nodes in global numbering
Definition: linsys_BDDC.hh:136
static const int registrar
Registrar of class to factory.
Definition: linsys_BDDC.hh:123
void print_matrix(std::ostream &out)
Definition: linsys_BDDC.cc:366
PetscErrorCode rhs_zero_entries() override
Definition: linsys_BDDC.cc:216
void apply_constrains(double scalar=1.) override
Definition: linsys_BDDC.cc:231
la::BddcmlWrapper Bddcml_
Definition: linsys_BDDC.hh:141
double get_solution_precision() override
Definition: linsys_BDDC.cc:357
void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override
Sets tolerances. Note that BDDC does not use a_tol.
Definition: linsys_BDDC.cc:79
int number_of_levels_
number of levels in the multilevel method
Definition: linsys_BDDC.hh:127
void finish_assembly() override
Definition: linsys_BDDC.cc:224
static const Input::Type::Record & get_input_type()
Definition: linsys_BDDC.cc:39
const bool swap_sign_
swap sign of matrix and rhs entries, e.g. to make the matrix SPD
Definition: linsys_BDDC.hh:134
LinSys::SolveInfo solve() override
Definition: linsys_BDDC.cc:238
bool use_adaptive_bddc_
should adaptive BDDC be used?
Definition: linsys_BDDC.hh:128
la::BddcmlWrapper::MatrixType BDDCMatrixType
Definition: linsys_BDDC.hh:44
PetscErrorCode mat_zero_entries() override
Definition: linsys_BDDC.cc:208
int max_nondecr_it_
parameters expected from input file:
Definition: linsys_BDDC.hh:126
void load_mesh(BDDCMatrixType matrix_type, const int nDim, const int numNodes, const int numDofs, 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_permeability, const int meshDim)
Definition: linsys_BDDC.cc:91
VecScatter VSpetscToSubScatter_
scatter from solution_ to locSolVec_
Definition: linsys_BDDC.hh:139
void diagonal_weights_set_value(int global_index, double value)
Definition: linsys_BDDC.cc:201
std::vector< double > locSolution_
subdomain solution
Definition: linsys_BDDC.hh:137
void rhs_set_values(int nrow, int *rows, double *vals) override
Definition: linsys_BDDC.cc:180
LinSys_BDDC(const Distribution *rows_ds, const bool swap_sign=false)
Definition: linsys_BDDC.cc:65
Bddcml_ * bddcml_
BDDCML wrapper.
Definition: linsys_BDDC.hh:142
void set_from_input(const Input::Record in_rec) override
Definition: linsys_BDDC.cc:279
unsigned int max_it_
maximum number of iterations of linear solver
Definition: linsys.hh:668
MPI_Comm comm_
Definition: linsys.hh:670
virtual void set_from_input(const Input::Record in_rec)
Definition: linsys.hh:641
Input::Record in_rec_
Definition: linsys.hh:697
ConstraintVec_ constraints_
Definition: linsys.hh:693
unsigned size_
global number of matrix rows, i.e. problem size
Definition: linsys.hh:674
double a_tol_
absolute tolerance of linear solver
Definition: linsys.hh:666
@ ADD
Definition: linsys.hh:98
SetValuesMode status_
Set value status of the linear system.
Definition: linsys.hh:671
const Distribution * rows_ds_
final distribution of rows of MH matrix
Definition: linsys.hh:676
Vec solution_
PETSc vector constructed with vb array.
Definition: linsys.hh:686
double residual_norm_
local solution array pointing into Vec solution_
Definition: linsys.hh:691
double r_tol_
relative tolerance of linear solver
Definition: linsys.hh:665
static Input::Type::Abstract & get_input_type()
Definition: linsys.cc:29
double normRhs()
Get norm of right-hand side.
void writeMatrix(std::ostream &out)
Outputs.
static constexpr bool value
Definition: json.hpp:87
Wrappers for linear systems based on MPIAIJ and MATIS format.
Solver based on Multilevel BDDC - using corresponding class of OpenFTL package.
#define LogOut()
Macro defining 'log' record of log.
Definition: logger.hh:281
unsigned int uint
#define MPI_COMM_WORLD
Definition: mpi.h:123
#define MPI_Comm_rank
Definition: mpi.h:236
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
#define FMT_UNUSED
Definition: posix.h:75
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
#define START_TIMER(tag)
Starts a timer with specified tag.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142