6 #ifdef FLOW123D_HAVE_BDDCML
12 const unsigned numDofsSub,
13 const MatrixType matrixType,
16 : numDofs_( static_cast<int>( numDofs ) ),
17 numDofsSub_( static_cast<int>( numDofsSub ) ),
18 matrixType_( matrixType ),
20 numSubLoc_ ( numSubLoc )
28 rhsVec_.resize( numDofsSub_, 0. );
31 ifix_.resize( numDofsSub_, 0 );
32 fixv_.resize( numDofsSub_, 0. );
35 sol_.resize( numDofsSub_, 0. );
43 isMatAssembled_ =
false;
101 meshDim_ =
static_cast<int>( meshDim );
105 numElemSub_ = isegn.size( );
106 numNodesSub_ = isngn.size( );
109 numNodes_ = numNodes;
112 int numElemSubInt =
static_cast<int>( numElemSub_ );
114 &(procElemStarts[0]), 1,
MPI_INT, comm_ );
117 for (
unsigned i = 1; i < static_cast<unsigned>( nProc_ ); i++ )
118 procElemStarts[i] = procElemStarts[i-1] + procElemStarts[i];
120 for (
unsigned i = nProc_; i > 0; i-- )
121 procElemStarts[i] = procElemStarts[i-1];
123 procElemStarts[0] = 0;
126 numElem_ = procElemStarts[ nProc_ ];
130 .error(
"array inet size mismatch \n");
132 .error(
"array nndf content mismatch \n");
134 .error(
"array isvgvn size mismatch \n");
136 .error(
"array nnet size mismatch \n");
138 .error(
"array nndf size mismatch \n");
140 .error(
"array xyz size mismatch \n");
143 inet_.resize( inet.size() );
144 std::copy( inet.begin(), inet.end(), inet_.begin() );
145 nnet_.resize( nnet.size() );
146 std::copy( nnet.begin(), nnet.end(), nnet_.begin() );
147 nndf_.resize( nndf.size() );
148 std::copy( nndf.begin(), nndf.end(), nndf_.begin() );
149 isegn_.resize( isegn.size() );
150 std::copy( isegn.begin(), isegn.end(), isegn_.begin() );
151 isngn_.resize( isngn.size() );
152 std::copy( isngn.begin(), isngn.end(), isngn_.begin() );
153 isvgvn_.resize( isvgvn.size() );
154 std::copy( isvgvn.begin(), isvgvn.end(), isvgvn_.begin() );
155 xyz_.resize( xyz.size() );
156 std::copy( xyz.begin(), xyz.end(), xyz_.begin() );
157 element_data_.resize( element_data.size() );
158 std::copy( element_data.begin(), element_data.end(), element_data_.begin() );
161 for (
unsigned ind = 0; ind < isvgvn_.size(); ++ind ) {
162 global2LocalDofMap_.insert( std::make_pair(
static_cast<unsigned>( isvgvn_[ind] ), ind ) );
176 unsigned length = numElements * elMatSize;
177 coo_.prepareAssembly( length );
185 const VecUint_ & rowIndices,
186 const VecUint_ & colIndices )
189 for (
unsigned i = 0; i < rowIndices.size(); i++) {
190 for (
unsigned j = 0; j < colIndices.size(); j++) {
193 if (subMat(i,j) != 0.0)
194 coo_.insert( rowIndices[i], colIndices[j], subMat(i,j) );
199 isMatAssembled_ =
false;
210 if ( isMatAssembled_ )
return;
213 coo_.finishAssembly( );
216 isMatAssembled_ =
true;
224 const VecUint_ & dofIndices )
226 for (
unsigned i = 0; i < dofIndices.size( ); i ++ ) {
229 Global2LocalMap_::iterator pos = global2LocalDofMap_.find( dofIndices[i] );
234 if (pos == global2LocalDofMap_.end()) {
235 DebugOut().every_proc() <<
"RHS, Missing map to local for idx: " << dofIndices[i];
237 const unsigned indLoc = pos -> second;
239 rhsVec_[ indLoc ] += subVec[i];
250 const double factor,
FMT_UNUSED const double scalar )
253 ConstraintVec_::const_iterator cIter = constraints.begin( );
254 ConstraintVec_::const_iterator cEnd = constraints.end( );
256 for ( ; cIter != cEnd; ++cIter ) {
257 unsigned globalDof = cIter -> first;
258 double value = (cIter -> second );
261 Global2LocalMap_::iterator pos = global2LocalDofMap_.find( globalDof );
262 if ( pos != global2LocalDofMap_.end() ) {
264 unsigned indLoc = pos -> second;
267 fixv_[ indLoc ] = factor *
value;
279 thisConstraint.push_back( std::make_pair( index, 0. ) );
280 this -> applyConstraints( thisConstraint, 1., scalar );
281 thisConstraint.clear();
289 const double &
value )
292 diagWeightsCoo_.insert( index, index,
value );
295 isDiagAssembled_ =
false;
306 if ( isDiagAssembled_ )
return;
309 diagWeightsCoo_.finishAssembly( );
312 isDiagAssembled_ =
true;
321 int verboseLevel,
int maxIt,
int ndecrMax,
bool use_adaptive )
325 ASSERT_PERMANENT(meshLoaded_).error(
"Subdomain mesh was not loaded, perhaps missing call to loadMesh() member function. \n ");
334 if ( numSubAtLevels != NULL ) {
338 .error(
"Inconsistent size of numbers of subdomains at given levels\n");
340 .error(
"Inconsistent number of subdomains at first level\n");
342 std::copy( (*numSubAtLevels).begin(), (*numSubAtLevels).end(),
348 if ( numLevels == 2 ) {
349 numSubLev[0] = numSub;
352 else if ( numLevels > 2 ) {
353 double coarsening = pow( numSub, 1. /
double(numLevels - 1) );
355 numSubLev[0] = numSub;
356 for (
int i = 1; i <= numLevels-2; i++) {
357 int ir = numLevels - i - 1;
358 numSubLev[i] = (int) pow( coarsening, ir );
359 if ( numSubLev[i]%2 != 0 ) {
360 numSubLev[i] = numSubLev[i] + 1;
363 numSubLev[numLevels-1] = 1;
366 ASSERT_PERMANENT(
false ).error(
"Missing numbers of subdomains at levels. \n" );
373 int lnumSubLev = numSubLev.size();
376 int justDirectSolve = 0;
378 bddcml_init_c( &numLevels, &(numSubLev[0]), &lnumSubLev, &numSubLoc_, &commInt, &verboseLevel, &numBase, &justDirectSolve );
385 int linet = inet_.size();
386 int lnnet = nnet_.size();
387 int lnndf = nndf_.size();
388 int lisngn = isngn_.size();
389 int lisvgvn = isvgvn_.size();
390 int lisegn = isegn_.size();
392 int lxyz1 = numNodesSub_;
395 int lifix = ifix_.size();
396 int lfixv = fixv_.size();
398 int lrhsVec = rhsVec_.size();
399 int isRhsCompleteInt = 0;
407 bool onlyUpperTriangle;
408 if ( matrixType_ == GENERAL ) {
409 onlyUpperTriangle =
false;
412 onlyUpperTriangle =
true;
416 this -> finishMatAssembly();
417 coo_.extractArrays( i_sparse, j_sparse, a_sparse, onlyUpperTriangle );
425 for (
unsigned inz = 0; inz < a_sparse.size(); inz++ ) {
426 if ( i_sparse[inz] != indRow ) {
427 indRow = i_sparse[inz];
428 Global2LocalMap_::iterator pos = global2LocalDofMap_.find(
static_cast<unsigned> ( indRow ) );
430 .error(
"Cannot remap 'indRow' to local indices. \n");
431 indRowLoc =
static_cast<int> ( pos -> second );
433 if ( j_sparse[inz] != indCol ) {
434 indCol = j_sparse[inz];
435 Global2LocalMap_::iterator pos = global2LocalDofMap_.find(
static_cast<unsigned> ( indCol ) );
437 .error(
"Cannot remap 'indCol' to local indices. \n");
438 indColLoc =
static_cast<int> ( pos -> second );
441 i_sparse[ inz ] = indRowLoc;
442 j_sparse[ inz ] = indColLoc;
447 int isMatAssembledInt = 0;
449 int la = a_sparse.size();
453 .error(
"It appears that diagonal weights for BDDC are not loaded. This is currently mandatory. \n");
460 this -> finishDiagAssembly();
461 diagWeightsCoo_.extractArrays( i_diag_sparse, j_diag_sparse, diag_sparse );
462 diagWeightsCoo_.clear();
467 for (
unsigned inz = 0; inz < diag_sparse.size(); inz++ ) {
468 indRow = i_diag_sparse[inz];
469 Global2LocalMap_::iterator pos = global2LocalDofMap_.find(
static_cast<unsigned> ( indRow ) );
471 .error(
"Cannot remap 'indRow' to local indices. \n");
472 indRowLoc =
static_cast<int> ( pos -> second );
474 i_diag_sparse[ inz ] = indRowLoc;
477 int lsub_diagonal = numDofsSub_;
479 .error(
"Array length mismatch.\n");
483 for (
int i = 0; i < lsub_diagonal; i++ ){
484 indRowLoc = i_diag_sparse[i];
485 double value = diag_sparse[i];
486 sub_diagonal[indRowLoc] =
value;
488 ASSERT_PERMANENT( std::find(sub_diagonal.begin(), sub_diagonal.end(), -1. ) == sub_diagonal.end() )
489 .error(
"There are missing entries in the diagonal of weights. \n");
492 int numDofsInt =
static_cast<int> ( numDofs_ );
493 int numDofsSubInt =
static_cast<int> ( numDofsSub_ );
496 if ( matrixType_ == GENERAL ) matrixTypeInt = 0 ;
497 else if ( matrixType_ == SPD ) matrixTypeInt = 1 ;
498 else if ( matrixType_ == SYMMETRICGENERAL ) matrixTypeInt = 2 ;
499 else if ( matrixType_ == SPD_VIA_SYMMETRICGENERAL ) matrixTypeInt = 2 ;
503 int lsol = sol_.size();
507 int lUserConstraints1 = 0;
508 int lUserConstraints2 = 0;
510 int lelement_data1 = 1;
511 int lelement_data2 = element_data_.size();
516 int find_components_int = 1;
517 int use_dual_mesh_graph_int = 0;
518 int neighbouring = 0;
521 bddcml_upload_subdomain_data_c( &numElem_, &numNodes_, &numDofsInt, &nDim_, &meshDim_,
522 &iSub, &numElemSub_, &numNodesSub_, &numDofsSubInt,
523 &(inet_[0]),&linet, &(nnet_[0]),&lnnet, &(nndf_[0]),&lnndf,
524 &(isngn_[0]),&lisngn, &(isvgvn_[0]),&lisvgvn, &(isegn_[0]),&lisegn,
525 &(xyz_[0]),&lxyz1,&lxyz2,
526 &(ifix_[0]),&lifix, &(fixv_[0]),&lfixv,
527 &(rhsVec_[0]),&lrhsVec, &isRhsCompleteInt,
529 &matrixTypeInt, &(i_sparse[0]), &(j_sparse[0]), &(a_sparse[0]), &
la, &isMatAssembledInt,
530 &(userConstraints[0]),&lUserConstraints1,&lUserConstraints2,
531 &(element_data_[0]),&lelement_data1,&lelement_data2,
532 &(sub_diagonal[0]), &lsub_diagonal,
533 &find_components_int, &use_dual_mesh_graph_int, &neighbouring );
541 int use_defaults_int = 0;
542 int parallel_division_int = 0;
543 int use_arithmetic_int = 1;
544 int use_adaptive_int;
545 if ( use_adaptive ) {
546 use_adaptive_int = 1;
549 use_adaptive_int = 0;
551 int use_user_constraints_int = 0;
552 int weights_type = 3;
555 int use_corner_constraints_int = 1;
559 bddcml_setup_preconditioner_c( & matrixTypeInt, & use_defaults_int,
560 & parallel_division_int, & use_arithmetic_int,
561 & use_corner_constraints_int,
562 & use_adaptive_int, & use_user_constraints_int,
567 double normSquaredLoc = 0.;
569 bddcml_dotprod_subdomain_c( &iSub, &(rhsVec_[0]), &lrhsVec, &(rhsVec_[0]), &lrhsVec, &normSquaredLoc );
571 for (
double elem: rhsVec_ ) normSquaredLoc += elem*elem;
575 double normSquared = 0.;
577 normRhs_ = std::sqrt( normSquared );
581 if ( matrixType_ == SPD )
583 else if ( matrixType_ == SPD_VIA_SYMMETRICGENERAL)
588 int recycling_int = 0;
589 int max_number_of_stored_vectors = 100;
593 bddcml_solve_c( & commInt, & method, & tol, & maxIt, & ndecrMax,
594 & recycling_int, & max_number_of_stored_vectors,
595 & numIter_, & convergedReason_, & condNumber_);
604 bddcml_download_local_solution_c( &iSub, &(sol_[0]), &lsol );
609 bddcml_dotprod_subdomain_c( &iSub, &(sol_[0]), &lsol, &(sol_[0]), &lsol, &normSquaredLoc );
611 for (
double elem: sol_ ) normSquaredLoc += elem*elem;
617 normSol_ = std::sqrt( normSquared );
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
#define ASSERT_PERMANENT_GT(a, b)
Definition of comparative assert macro (Greater Than)
#define ASSERT_PERMANENT_EQ(a, b)
Definition of comparative assert macro (EQual)
std::vector< int > isvgvn_
Indices of Subdomain Variables in Global Variable Numbering ( i.e. dof mapping )
std::vector< double > sol_
distributed solution vector
la::MatrixCoo< int, double > coo_
matrix in coordinate format (COO)
void insertToMatrix(const SubMat_ &eStiff, const VecUint_ &rowIndices, const VecUint_ &colIndices)
insert submatrix to system matrix - into coordinate matrix object
std::vector< double > xyz_
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
std::vector< int > isngn_
Indices of Subdomain Nodes in Global Numbering.
void applyConstraints(ConstraintVec_ &constraints, const double factor, const double scalar)
void insertToDiagonalWeights(const int &index, const double &value)
std::vector< int > inet_
safety state variable to warn if solver is called without loading subdomain mesh
std::vector< double > rhsVec_
vector with RHS values restricted to subdomain, i.e. values are repeated at shared nodes
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.
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.
void fixDOF(const unsigned index, const double scalar=1.)
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.
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 finishDiagAssembly()
Finalize assembly of diagonal.
void clear()
destroy all contents
static constexpr bool value
#define DebugOut()
Macro defining 'debug' record of log.
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
#define MPI_Allgather(sendbuf, sendcount, sendtype, recvbuf, recvcount, recvtype, comm)
#define MPI_Comm_c2f(comm)
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.