PETSc will automatically avoid recomputing preconditioners if you call
solve multiple times with the same matrix and different right-hand-sides,
but we don't currently get this benefit because in PetscLinearSolver we
re-initialize the solver's matrix on each call to solve.
I tested out adding solve_with_same_matrix (see code snippet below) and
this works as expected. (This can verified by checking the output of
-log_summary from PETSc.) To use this you have to call
PetscLinearSolver::init(matrix) before you call solve_with_same_matrix.
I was wondering if there would be interest in adding something like this to
libMesh? If so, it would be nice if we made a unified LinearSolver
interface for this behavior.
David
-----------------------------------------------
template <typename T>
std::pair<unsigned int, Real>
PetscLinearSolver<T>::solve_with_same_matrix (NumericVector<T> &solution_in,
NumericVector<T> &rhs_in,
const double tol,
const unsigned int m_its)
{
START_LOG("solve_with_same_matrix()", "PetscLinearSolver");
libmesh_assert(this->initialized());
// Make sure the data passed in are really of Petsc types
PetscVector<T>* solution = cast_ptr<PetscVector<T>*>(&solution_in);
PetscVector<T>* rhs = cast_ptr<PetscVector<T>*>(&rhs_in);
PetscErrorCode ierr=0;
PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
PetscReal final_resid=0.;
// Close the vectors in case this wasn't already done.
solution->close ();
rhs->close ();
// Set the tolerances for the iterative solver. Use the user-supplied
// tolerance for the relative residual & leave the others at default
values.
ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
PETSC_DEFAULT, max_its);
LIBMESH_CHKERRABORT(ierr);
ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
LIBMESH_CHKERRABORT(ierr);
// Get the number of iterations required for convergence
ierr = KSPGetIterationNumber (_ksp, &its);
LIBMESH_CHKERRABORT(ierr);
// Get the norm of the final residual to return to the user.
ierr = KSPGetResidualNorm (_ksp, &final_resid);
LIBMESH_CHKERRABORT(ierr);
STOP_LOG("solve_with_same_matrix()", "PetscLinearSolver");
// return the # of its. and the final residual norm.
return std::make_pair(its, final_resid);
}
------------------------------------------------------------------------------
BPM Camp - Free Virtual Workshop May 6th at 10am PDT/1PM EDT
Develop your own process in accordance with the BPMN 2 standard
Learn Process modeling best practices with Bonita BPM through live exercises
http://www.bonitasoft.com/be-part-of-it/events/bpm-camp-virtual- event?utm_
source=Sourceforge_BPM_Camp_5_6_15&utm_medium=email&utm_campaign=VA_SF
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
