33 return it::Record(
"Petsc",
"PETSc solver settings.\n It provides interface to various PETSc solvers. The convergence criteria is:\n"
35 "norm( res_i ) < max( norm( res_0 ) * r_tol, a_tol )\n"
37 "where ```res_i``` is the residuum vector after i-th iteration of the solver and ```res_0``` is the estimate of the norm of the 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 (($L_2$)) 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.")
44 "If not, we use the value 1.0e-7."),
45 "Residual tolerance relative to the initial error.")
47 "If not, we use the value 1.0e-11."),
48 "Absolute residual tolerance.")
50 "If not, we use the value 10000."),
51 "Tolerance for divergence.")
53 "If not, we use the value 1000."),
54 "Maximum number of outer iterations of the linear solver.")
56 "If the string is left empty (by default), the internal default options is used.")
67 init_guess_nonzero(false),
75 ierr = VecZeroEntries(
rhs_ ); CHKERRV( ierr );
85 :
LinSys(other), params_(other.params_), v_rhs_(NULL), solution_precision_(other.solution_precision_)
133 ASSERT_PERMANENT(
false).error(
"Can not set values. Matrix is not preallocated.\n");
151 ASSERT_PERMANENT(
false).error(
"Can not set values. Matrix is not preallocated.\n");
167 default:
DebugOut() <<
"LS SetValues with non allowed insert mode.\n";
180 ierr = VecSetValues(
rhs_,nrow,
rows,vals,(InsertMode)
status_); CHKERRV( ierr );
184 default:
ASSERT_PERMANENT(
false).error(
"LinSys's status disallow set values.\n");
196 for (i=0; i<nrow; i++) {
198 for(j=0; j<ncol; j++) {
201 VecSetValue(
on_vec_,row,1.0,ADD_VALUES);
203 VecSetValue(
off_vec_,row,1.0,ADD_VALUES);
213 PetscInt *on_nz, *off_nz;
214 PetscScalar *on_array, *off_array;
226 VecGetArray(
on_vec_, &on_array );
227 VecGetArray(
off_vec_, &off_array );
234 VecRestoreArray(
on_vec_,&on_array);
235 VecRestoreArray(
off_vec_,&off_array);
245 0, on_nz, 0, off_nz, &
matrix_); CHKERRV( ierr );
248 MatSetOption(
matrix_, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
253 MatSetOption(
matrix_, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
263 MatAssemblyType assemblyType = MAT_FINAL_ASSEMBLY;
272 WarningOut() <<
"Finalizing linear system without setting values.\n";
275 ierr = MatAssemblyBegin(
matrix_, assembly_type); CHKERRV( ierr );
276 ierr = VecAssemblyBegin(
rhs_); CHKERRV( ierr );
277 ierr = MatAssemblyEnd(
matrix_, assembly_type); CHKERRV( ierr );
278 ierr = VecAssemblyEnd(
rhs_); CHKERRV( ierr );
280 if (assembly_type == MAT_FINAL_ASSEMBLY)
status_ =
DONE;
296 ASSERT_EQ(
status_,
DONE).error(
"System matrix and right-hand side are not assembled when applying constraints." );
299 PetscInt numConstraints =
static_cast<PetscInt
>(
constraints_.size() );
302 const PetscScalar diagScalar =
static_cast<PetscScalar
>( scalar );
308 ConstraintVec_::const_iterator cIter =
constraints_.begin( );
309 ConstraintVec_::const_iterator cEnd =
constraints_.end( );
311 for ( ; cIter != cEnd; ++cIter ) {
312 globalDofs.push_back(
static_cast<PetscInt
>( cIter -> first ) );
313 values.push_back(
static_cast<PetscScalar
>( cIter -> second ) * diagScalar );
321 ierr = MatZeroRows(
matrix_, numConstraints, globalDofPtr, diagScalar, PETSC_NULL, PETSC_NULL ); CHKERRV( ierr );
325 if ( numConstraints ) {
326 ierr = VecSetValues(
rhs_, numConstraints, globalDofPtr, valuePtr, INSERT_VALUES ); CHKERRV( ierr );
331 ierr = VecAssemblyBegin(
rhs_ ); CHKERRV( ierr );
332 ierr = VecAssemblyEnd(
rhs_ ); CHKERRV( ierr );
345 const char *petsc_dflt_opt;
360 petsc_dflt_opt=
"-ksp_type cg -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_type asm -pc_asm_type basic -pc_asm_overlap 4 -sub_pc_type icc -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
363 petsc_dflt_opt=
"-ksp_type bcgs -ksp_diagonal_scale -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";
369 petsc_dflt_opt=
"-ksp_type cg -pc_type icc -pc_factor_levels 3 -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
372 petsc_dflt_opt=
"-ksp_type bcgs -pc_type ilu -pc_factor_levels 5 -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
376 LogOut().fmt(
"inserting petsc options: {}\n",
params_.c_str());
380 PetscOptionsInsertString(NULL,
params_.c_str());
382 MatSetOption(
matrix_, MAT_USE_INODES, PETSC_FALSE );
391 KSPSetFromOptions(
system);
398 KSPGetType(
system, &type);
399 if (strcmp(type, KSPPREONLY) != 0)
400 KSPSetInitialGuessNonzero(
system, PETSC_TRUE);
408 KSPGetIterationNumber(
system,&nits);
414 LogOut().fmt(
"convergence reason {}, number of iterations is {}\n",
reason, nits);
430 FilePath matFileName(text +
"_flow123d_matrix.m",FilePath::FileType::output_file);
431 FilePath rhsFileName(text +
"_flow123d_rhs.m",FilePath::FileType::output_file);
432 FilePath solFileName(text +
"_flow123d_sol.m",FilePath::FileType::output_file);
434 PetscViewer myViewer;
437 PetscViewerASCIIOpen(
comm_,((
string)matFileName).c_str(),&myViewer);
438 PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
440 PetscViewerDestroy(&myViewer);
443 WarningOut() <<
"PetscViewer: the matrix of LinSys is not set.\n";
445 if (
rhs_ != NULL ) {
446 PetscViewerASCIIOpen(
comm_,((
string)rhsFileName).c_str(),&myViewer);
447 PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
448 VecView(
rhs_, myViewer );
449 PetscViewerDestroy(&myViewer);
452 WarningOut() <<
"PetscViewer: the rhs of LinSys is not set.\n";
455 PetscViewerASCIIOpen(
comm_,((
string)solFileName).c_str(),&myViewer);
456 PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
458 PetscViewerDestroy(&myViewer);
461 WarningOut() <<
"PetscViewer: the solution of LinSys is not set.\n";
482 std::string user_params = in_rec.
val<
string>(
"options");
483 if (user_params !=
"")
params_ = user_params;
497 double residual_norm;
498 VecNorm(
residual_, NORM_2, &residual_norm);
499 return residual_norm;
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
unsigned int lsize(int proc) const
get local size
unsigned int get_proc(unsigned int idx) const
get processor of the given index
unsigned int size() const
get global size
unsigned int np() const
get num of processors
Dedicated class for storing path to input and output files.
void finish_assembly() override
void start_add_assembly() override
void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals) override
void apply_constrains(double scalar=1.) override
double compute_residual() override
Mat matrix_
Petsc matrix of the problem.
void preallocate_matrix()
LinSys::SolveInfo solve() override
double solution_precision_
void preallocate_values(int nrow, int *rows, int ncol, int *cols)
bool init_guess_nonzero
flag for starting from nonzero guess
void start_insert_assembly() override
void set_from_input(const Input::Record in_rec) override
static const Input::Type::Record & get_input_type()
void set_initial_guess_nonzero(bool set_nonzero=true)
std::string params_
command-line-like options for the PETSc solver
LinSys_PETSC(const Distribution *rows_ds, const std::string ¶ms="")
KSPConvergedReason reason
static const int registrar
Registrar of class to factory.
Vec rhs_
PETSc vector constructed with vx array.
Vec off_vec_
Vectors for counting non-zero entries in off-diagonal block.
Vec on_vec_
Vectors for counting non-zero entries in diagonal block.
T * makePetscPointer_(std::vector< T > &array)
void view(string text="") override
void rhs_set_values(int nrow, int *rows, double *vals) override
double * v_rhs_
local RHS array pointing to Vec rhs_
void start_allocation() override
double get_solution_precision() override
void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override
unsigned int max_it_
maximum number of iterations of linear solver
double d_tol_
tolerance for divergence of linear solver
virtual void set_from_input(const Input::Record in_rec)
bool is_positive_definite()
ConstraintVec_ constraints_
double a_tol_
absolute tolerance of linear solver
SetValuesMode status_
Set value status of the linear system.
const Distribution * rows_ds_
final distribution of rows of MH matrix
Vec solution_
PETSc vector constructed with vb array.
double residual_norm_
local solution array pointing into Vec solution_
double r_tol_
relative tolerance of linear solver
bool matrix_changed_
true if the matrix was changed since the last solve
static Input::Type::Abstract & get_input_type()
bool rhs_changed_
true if the right hand side was changed since the last solve
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
#define WarningOut()
Macro defining 'warning' record of log.
#define LogOut()
Macro defining 'log' record of log.
#define DebugOut()
Macro defining 'debug' record of log.
#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.