Flow123d  release_2.2.0-34-g18a8075
linsys_PETSC.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_PETSC.cc
15  * @brief Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construction
16  * @author Jakub Sistek
17  */
18 
19 // derived from base linsys
20 #include "la/linsys_PETSC.hh"
21 #include "petscvec.h"
22 #include "petscksp.h"
23 #include "petscmat.h"
24 #include "system/sys_profiler.hh"
25 #include "system/system.hh"
26 
27 
28 //#include <boost/bind.hpp>
29 
30 namespace it = Input::Type;
31 
33  return it::Record("Petsc", "Interface to PETSc solvers. Convergence criteria is:\n"
34  "```\n"
35  "norm( res_n ) < max( norm( res_0 ) * r_tol, a_tol )\n"
36  "```\n"
37  "where res_i is the residuum vector after i-th iteration of the solver and res_0 is an estimate of the norm of initial residual. "
38  "If the initial guess of the solution is provided (usually only for transient equations) the residual of this estimate is used, "
39  "otherwise the norm of preconditioned RHS is used. "
40  "The default norm is L2 norm of preconditioned residual: (($ P^{-1}(Ax-b)$)), usage of other norm may be prescribed using the 'option' key. "
41  "See also PETSc documentation for KSPSetNormType.")
43  .declare_key("r_tol", it::Double(0.0, 1.0), it::Default::read_time("Defalut value set by nonlinear solver or equation. If not we use value 1.0e-7."),
44  "Relative residual tolerance, (to initial error).")
45  .declare_key("a_tol", it::Double(0.0), it::Default::read_time("Defalut value set by nonlinear solver or equation. If not we use value 1.0e-11."),
46  "Absolute residual tolerance.")
47  .declare_key("max_it", it::Integer(0), it::Default::read_time("Defalut value set by nonlinear solver or equation. If not we use value 1000."),
48  "Maximum number of outer iterations of the linear solver.")
49  .declare_key("options", it::String(), it::Default("\"\""), "Options passed to PETSC before creating KSP instead of default setting.")
50  .close();
51 }
52 
53 
55 
56 
57 LinSys_PETSC::LinSys_PETSC( const Distribution * rows_ds, const std::string &params)
58  : LinSys( rows_ds ),
59  params_(params),
60  init_guess_nonzero(false),
61  matrix_(0)
62 {
63  // create PETSC vectors:
64  PetscErrorCode ierr;
65  // rhs
66  v_rhs_= new double[ rows_ds_->lsize() + 1 ];
67  ierr = VecCreateMPIWithArray( comm_, 1, rows_ds_->lsize(), PETSC_DECIDE, v_rhs_, &rhs_ ); CHKERRV( ierr );
68  ierr = VecZeroEntries( rhs_ ); CHKERRV( ierr );
69  VecDuplicate(rhs_, &residual_);
70 
71  matrix_ = NULL;
72  solution_precision_ = std::numeric_limits<double>::infinity();
73  matrix_changed_ = true;
74  rhs_changed_ = true;
75 }
76 
78  : LinSys(other), params_(other.params_), v_rhs_(NULL), solution_precision_(other.solution_precision_)
79 {
80  MatCopy(other.matrix_, matrix_, DIFFERENT_NONZERO_PATTERN);
81  VecCopy(other.rhs_, rhs_);
82  VecCopy(other.on_vec_, on_vec_);
83  VecCopy(other.off_vec_, off_vec_);
84 }
85 
86 void LinSys_PETSC::set_tolerances(double r_tol, double a_tol, unsigned int max_it)
87 {
88  if (! in_rec_.is_empty()) {
89  // input record is set
90  r_tol_ = in_rec_.val<double>("r_tol", r_tol);
91  a_tol_ = in_rec_.val<double>("a_tol", a_tol);
92  max_it_ = in_rec_.val<unsigned int>("max_it", max_it);
93  } else {
94  r_tol_ = r_tol;
95  a_tol_ = a_tol;
96  max_it_ = max_it;
97 
98  }
99 }
100 
101 
103 {
104  PetscErrorCode ierr;
105 
106  ierr = VecCreateMPI( comm_, rows_ds_->lsize(), PETSC_DECIDE, &(on_vec_) ); CHKERRV( ierr );
107  ierr = VecDuplicate( on_vec_, &(off_vec_) ); CHKERRV( ierr );
108  status_ = ALLOCATE;
109 }
110 
112 {
113  switch ( status_ ) {
114  case ALLOCATE:
115  this->preallocate_matrix( );
116  break;
117  case INSERT:
118  this->finish_assembly( MAT_FLUSH_ASSEMBLY );
119  break;
120  case ADD:
121  case DONE:
122  break;
123  default:
124  OLD_ASSERT( false, "Can not set values. Matrix is not preallocated.\n");
125  }
126  status_ = ADD;
127 }
128 
130 {
131  switch ( status_ ) {
132  case ALLOCATE:
133  this->preallocate_matrix();
134  break;
135  case ADD:
136  this->finish_assembly( MAT_FLUSH_ASSEMBLY );
137  break;
138  case INSERT:
139  case DONE:
140  break;
141  default:
142  OLD_ASSERT( false, "Can not set values. Matrix is not preallocated.\n");
143  }
144  status_ = INSERT;
145 }
146 
147 void LinSys_PETSC::mat_set_values( int nrow, int *rows, int ncol, int *cols, double *vals )
148 {
149  PetscErrorCode ierr;
150 
151  // here vals would need to be converted from double to PetscScalar if it was ever something else than double :-)
152  switch (status_) {
153  case INSERT:
154  case ADD:
155  ierr = MatSetValues(matrix_,nrow,rows,ncol,cols,vals,(InsertMode)status_); CHKERRV( ierr );
156  break;
157  case ALLOCATE:
158  this->preallocate_values(nrow,rows,ncol,cols);
159  break;
160  default: DebugOut() << "LS SetValues with non allowed insert mode.\n";
161  }
162 
163  matrix_changed_ = true;
164 }
165 
166 void LinSys_PETSC::rhs_set_values( int nrow, int *rows, double *vals )
167 {
168  PetscErrorCode ierr;
169 
170  switch (status_) {
171  case INSERT:
172  case ADD:
173  ierr = VecSetValues(rhs_,nrow,rows,vals,(InsertMode)status_); CHKERRV( ierr );
174  break;
175  case ALLOCATE:
176  break;
177  default: OLD_ASSERT(false, "LinSys's status disallow set values.\n");
178  }
179 
180  rhs_changed_ = true;
181 }
182 
183 void LinSys_PETSC::preallocate_values(int nrow,int *rows,int ncol,int *cols)
184 {
185  int i,j;
186  int col;
187  PetscInt row;
188 
189  for (i=0; i<nrow; i++) {
190  row=rows[i];
191  for(j=0; j<ncol; j++) {
192  col = cols[j];
193  if (rows_ds_->get_proc(row) == rows_ds_->get_proc(col))
194  VecSetValue(on_vec_,row,1.0,ADD_VALUES);
195  else
196  VecSetValue(off_vec_,row,1.0,ADD_VALUES);
197  }
198  }
199 }
200 
202 {
203  OLD_ASSERT(status_ == ALLOCATE, "Linear system has to be in ALLOCATE status.");
204 
205  PetscErrorCode ierr;
206  PetscInt *on_nz, *off_nz;
207  PetscScalar *on_array, *off_array;
208 
209  // assembly and get values from counting vectors, destroy them
210  VecAssemblyBegin(on_vec_);
211  VecAssemblyBegin(off_vec_);
212 
213  on_nz = new PetscInt[ rows_ds_->lsize() ];
214  off_nz = new PetscInt[ rows_ds_->lsize() ];
215 
216  VecAssemblyEnd(on_vec_);
217  VecAssemblyEnd(off_vec_);
218 
219  VecGetArray( on_vec_, &on_array );
220  VecGetArray( off_vec_, &off_array );
221 
222  for ( unsigned int i=0; i<rows_ds_->lsize(); i++ ) {
223  on_nz[i] = static_cast<PetscInt>( on_array[i]+0.1 ); // small fraction to ensure correct rounding
224  off_nz[i] = static_cast<PetscInt>( off_array[i]+0.1 );
225  }
226 
227  VecRestoreArray(on_vec_,&on_array);
228  VecRestoreArray(off_vec_,&off_array);
229  VecDestroy(&on_vec_);
230  VecDestroy(&off_vec_);
231 
232  // create PETSC matrix with preallocation
233  if (matrix_ != NULL)
234  {
235  ierr = MatDestroy(&matrix_); CHKERRV( ierr );
236  }
237  ierr = MatCreateAIJ(PETSC_COMM_WORLD, rows_ds_->lsize(), rows_ds_->lsize(), PETSC_DETERMINE, PETSC_DETERMINE,
238  0, on_nz, 0, off_nz, &matrix_); CHKERRV( ierr );
239 
240  if (symmetric_) MatSetOption(matrix_, MAT_SYMMETRIC, PETSC_TRUE);
241  MatSetOption(matrix_, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
242 
243  delete[] on_nz;
244  delete[] off_nz;
245 }
246 
248 {
249  MatAssemblyType assemblyType = MAT_FINAL_ASSEMBLY;
250  this->finish_assembly( assemblyType );
251 }
252 
253 void LinSys_PETSC::finish_assembly( MatAssemblyType assembly_type )
254 {
255  PetscErrorCode ierr;
256 
257  if (status_ == ALLOCATE) {
258  WarningOut() << "Finalizing linear system without setting values.\n";
259  this->preallocate_matrix();
260  }
261  ierr = MatAssemblyBegin(matrix_, assembly_type); CHKERRV( ierr );
262  ierr = VecAssemblyBegin(rhs_); CHKERRV( ierr );
263  ierr = MatAssemblyEnd(matrix_, assembly_type); CHKERRV( ierr );
264  ierr = VecAssemblyEnd(rhs_); CHKERRV( ierr );
265 
266  if (assembly_type == MAT_FINAL_ASSEMBLY) status_ = DONE;
267 
268  matrix_changed_ = true;
269  rhs_changed_ = true;
270 }
271 
272 void LinSys_PETSC::apply_constrains( double scalar )
273 {
274  PetscErrorCode ierr;
275 
276  // check that system matrix is assembled
277  OLD_ASSERT ( status_ == DONE, "System matrix and right-hand side are not assembled when applying constraints." );
278 
279  // number of constraints
280  PetscInt numConstraints = static_cast<PetscInt>( constraints_.size() );
281 
282  // Additional multiplier for numerical reasons (criterion to be established)
283  const PetscScalar diagScalar = static_cast<PetscScalar>( scalar );
284 
285  std::vector<PetscInt> globalDofs;
287 
288  // Constraint iterators
289  ConstraintVec_::const_iterator cIter = constraints_.begin( );
290  ConstraintVec_::const_iterator cEnd = constraints_.end( );
291  // collect global dof indices and the correpsonding values
292  for ( ; cIter != cEnd; ++cIter ) {
293  globalDofs.push_back( static_cast<PetscInt>( cIter -> first ) );
294  values.push_back( static_cast<PetscScalar>( cIter -> second ) * diagScalar );
295  }
296 
297  // prepare pointers to be passed to PETSc
298  PetscInt * globalDofPtr = this->makePetscPointer_( globalDofs );
299  PetscScalar * valuePtr = this->makePetscPointer_( values );
300 
301  // set matrix rows to zero
302  ierr = MatZeroRows( matrix_, numConstraints, globalDofPtr, diagScalar, PETSC_NULL, PETSC_NULL ); CHKERRV( ierr );
303  matrix_changed_ = true;
304 
305  // set RHS entries to values (crashes if called with NULL pointers)
306  if ( numConstraints ) {
307  ierr = VecSetValues( rhs_, numConstraints, globalDofPtr, valuePtr, INSERT_VALUES ); CHKERRV( ierr );
308  rhs_changed_ = true;
309  }
310 
311  // perform communication in the rhs vector
312  ierr = VecAssemblyBegin( rhs_ ); CHKERRV( ierr );
313  ierr = VecAssemblyEnd( rhs_ ); CHKERRV( ierr );
314 }
315 
316 
318 {
319  init_guess_nonzero = set_nonzero;
320 }
321 
322 
324 {
325 
326  const char *petsc_dflt_opt;
327  int nits;
328 
329  // -mat_no_inode ... inodes are usefull only for
330  // vector problems e.g. MH without Schur complement reduction
331  if (rows_ds_->np() > 1) {
332  // parallel setting
333  if (this->is_positive_definite())
334  //petsc_dflt_opt="-ksp_type cg -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type icc -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
335  petsc_dflt_opt="-ksp_type bcgs -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type ilu -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
336  else
337  petsc_dflt_opt="-ksp_type bcgs -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type ilu -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
338 
339  }
340  else {
341  // serial setting
342  if (this->is_positive_definite())
343  //petsc_dflt_opt="-ksp_type cg -pc_type icc -pc_factor_levels 3 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
344  petsc_dflt_opt="-ksp_type bcgs -pc_type ilu -pc_factor_levels 5 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
345  else
346  petsc_dflt_opt="-ksp_type bcgs -pc_type ilu -pc_factor_levels 5 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
347  }
348 
349  if (params_ == "") params_ = petsc_dflt_opt;
350  LogOut().fmt("inserting petsc options: {}\n",params_.c_str());
351  PetscOptionsInsertString(params_.c_str()); // overwrites previous options values
352 
353  MatSetOption( matrix_, MAT_USE_INODES, PETSC_FALSE );
354 
355  chkerr(KSPCreate( comm_, &system ));
356  chkerr(KSPSetOperators(system, matrix_, matrix_));
357 
358 
359  // TODO take care of tolerances - shall we support both input file and command line petsc setting
360  chkerr(KSPSetTolerances(system, r_tol_, a_tol_, PETSC_DEFAULT,PETSC_DEFAULT));
361  chkerr(KSPSetTolerances(system, r_tol_, a_tol_, PETSC_DEFAULT, max_it_));
362  KSPSetFromOptions(system);
363  // We set the KSP flag set_initial_guess_nonzero
364  // unless KSP type is preonly.
365  // In such case PETSc fails (version 3.4.1)
366  if (init_guess_nonzero)
367  {
368  KSPType type;
369  KSPGetType(system, &type);
370  if (strcmp(type, KSPPREONLY) != 0)
371  KSPSetInitialGuessNonzero(system, PETSC_TRUE);
372  }
373 
374  {
375  START_TIMER("PETSC linear solver");
376  START_TIMER("PETSC linear iteration");
377  chkerr(KSPSolve(system, rhs_, solution_ ));
378  KSPGetConvergedReason(system,&reason);
379  KSPGetIterationNumber(system,&nits);
380  ADD_CALLS(nits);
381  }
382  // substitute by PETSc call for residual
383  VecNorm(rhs_, NORM_2, &residual_norm_);
384 
385  LogOut().fmt("convergence reason {}, number of iterations is {}\n", reason, nits);
386 
387  // get residual norm
388  KSPGetResidualNorm(system, &solution_precision_);
389 
390  // TODO: I do not understand this
391  //Profiler::instance()->set_timer_subframes("SOLVING MH SYSTEM", nits);
392 
393  KSPDestroy(&system);
394 
395  return static_cast<int>(reason);
396 
397 }
398 
400 {
401  std::string matFileName = "flow123d_matrix.m";
402  std::string rhsFileName = "flow123d_rhs.m";
403  std::string solFileName = "flow123d_sol.m";
404 
405  PetscViewer myViewer;
406 
407  PetscViewerASCIIOpen(comm_,matFileName.c_str(),&myViewer);
408  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
409  MatView( matrix_, myViewer );
410  PetscViewerDestroy(&myViewer);
411 
412  PetscViewerASCIIOpen(comm_,rhsFileName.c_str(),&myViewer);
413  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
414  VecView( rhs_, myViewer );
415  PetscViewerDestroy(&myViewer);
416 
417  if ( solution_ != NULL ) {
418  PetscViewerASCIIOpen(comm_,solFileName.c_str(),&myViewer);
419  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
420  VecView( solution_, myViewer );
421  PetscViewerDestroy(&myViewer);
422  }
423 }
424 
426 {
427  PetscErrorCode ierr;
428 
429  if (matrix_ != NULL) { ierr = MatDestroy(&matrix_); CHKERRV( ierr ); }
430  ierr = VecDestroy(&rhs_); CHKERRV( ierr );
431 
432  if (v_rhs_ != NULL) delete[] v_rhs_;
433 }
434 
435 
436 
438 {
439  LinSys::set_from_input( in_rec );
440 
441  // PETSC specific parameters
442  // If parameters are specified in input file, they are used,
443  // otherwise keep settings provided in constructor of LinSys_PETSC.
444  std::string user_params = in_rec.val<string>("options");
445  if (user_params != "") params_ = user_params;
446 }
447 
448 
450 {
451  return solution_precision_;
452 }
453 
454 
456 {
457  MatMult(matrix_, solution_, residual_);
458  VecAXPY(residual_,-1.0, rhs_);
459  double residual_norm;
460  VecNorm(residual_, NORM_2, &residual_norm);
461  return residual_norm;
462 }
unsigned int get_proc(unsigned int idx) const
get processor of the given index
SetValuesMode status_
Set value status of the linear system.
Definition: linsys.hh:583
void apply_constrains(double scalar=1.) override
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:598
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
virtual void set_from_input(const Input::Record in_rec)
Definition: linsys.hh:554
bool matrix_changed_
true if the matrix was changed since the last solve
Definition: linsys.hh:595
double * v_rhs_
local RHS array pointing to Vec rhs_
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
void set_from_input(const Input::Record in_rec) override
int solve() override
bool symmetric_
Definition: linsys.hh:590
ConstraintVec_ constraints_
Definition: linsys.hh:604
Vec rhs_
PETSc vector constructed with vx array.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:147
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
void start_allocation() override
Class for declaration of the integral input data.
Definition: type_base.hh:489
void preallocate_values(int nrow, int *rows, int ncol, int *cols)
Input::Record in_rec_
Definition: linsys.hh:608
void start_add_assembly() override
#define LogOut()
Macro defining &#39;log&#39; record of log.
Definition: logger.hh:239
bool rhs_changed_
true if the right hand side was changed since the last solve
Definition: linsys.hh:596
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:303
double get_solution_precision() override
void start_insert_assembly() override
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:195
#define OLD_ASSERT(...)
Definition: global_defs.h:131
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:540
MPI_Comm comm_
Definition: linsys.hh:582
unsigned int max_it_
maximum number of iterations of linear solver
Definition: linsys.hh:580
T * makePetscPointer_(std::vector< T > &array)
LinSys_PETSC(const Distribution *rows_ds, const std::string &params="")
Definition: linsys_PETSC.cc:57
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
const Ret val(const string &key) const
#define START_TIMER(tag)
Starts a timer with specified tag.
const Distribution * rows_ds_
final distribution of rows of MH matrix
Definition: linsys.hh:588
double a_tol_
absolute tolerance of linear solver
Definition: linsys.hh:579
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:490
bool is_empty() const
Definition: accessors.hh:366
unsigned int np() const
get num of processors
double residual_norm_
local solution array pointing into Vec solution_
Definition: linsys.hh:602
Vec on_vec_
Vectors for counting non-zero entries in diagonal block.
void preallocate_matrix()
double compute_residual() override
std::string params_
command-line-like options for the PETSc solver
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
KSPConvergedReason reason
Vec off_vec_
Vectors for counting non-zero entries in off-diagonal block.
double r_tol_
relative tolerance of linear solver
Definition: linsys.hh:578
void finish_assembly() override
void view() override
void rhs_set_values(int nrow, int *rows, double *vals) override
Abstract linear system class.
bool is_positive_definite()
Definition: linsys.hh:510
#define WarningOut()
Macro defining &#39;warning&#39; record of log.
Definition: logger.hh:236
void set_initial_guess_nonzero(bool set_nonzero=true)
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
void set_tolerances(double r_tol, double a_tol, unsigned int max_it) override
Definition: linsys_PETSC.cc:86
Mat matrix_
Petsc matrix of the problem.
static const int registrar
Registrar of class to factory.
Record type proxy class.
Definition: type_record.hh:182
void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals) override
static Input::Type::Abstract & get_input_type()
Definition: linsys.cc:29
bool init_guess_nonzero
flag for starting from nonzero guess
#define DebugOut()
Macro defining &#39;debug&#39; record of log.
Definition: logger.hh:242
Class for declaration of the input data that are in string format.
Definition: type_base.hh:588
Vec solution_
PETSc vector constructed with vb array.
Definition: linsys.hh:598
double solution_precision_
unsigned int lsize(int proc) const
get local size