comments inline
On Mon, 13 Sep 2010 15:18:30 +0200 (CEST), Tim Kroeger
<[email protected]> wrote:
> Index: include/solvers/petsc_linear_solver.h
> ===================================================================
> --- include/solvers/petsc_linear_solver.h (Revision 3957)
> +++ include/solvers/petsc_linear_solver.h (Arbeitskopie)
> @@ -122,7 +122,15 @@
> * Initialize data structures if not done so already plus much more
> */
> void init (PetscMatrix<T>* matrix);
> +
> /**
> + * After calling this method, all successive solves will be limited
> + * to the given set of dofs. This mode can be disabled by calling
> + * this method with \p dofs being a \p NULL pointer.
> + */
> + virtual void solve_only_on (const std::set<unsigned int>* const dofs);
> +
> + /**
> * Call the Petsc solver. It calls the method below, using the
> * same matrix for the system and preconditioner matrices.
> */
> @@ -253,13 +261,36 @@
> * Krylov subspace context
> */
> KSP _ksp;
> +
> + /**
> + * Set of dofs to which any solve is limited (\p NULL means solve on
> + * all dofs).
> + */
> + const std::set<unsigned int>* _solve_only_on_dofs;
I don't think a set is the correct data structure because the order
actually matters (it affects memory performance and the strength of
preconditioners). I also don't see why it needs to be kept separate
from the IS, you're duplicating information, it uses more memory and
will have to be kept consistent.
> +
> + /**
> + * PETSc index set containing the dofs on which to solve (\p NULL
> + * means solve on all dofs).
> + */
> + IS _solve_only_on_is;
> +
> + /**
> + * PETSc index set containing dofs from 0 to n-1 where n is the
> + * number of dofs on which to solve (\p NULL means solve on all
> + * dofs).
> + */
> + IS _solve_only_on_is_seq;
You don't need this IS, see below.
> +
> };
>
>
> /*----------------------- functions ----------------------------------*/
> template <typename T>
> inline
> -PetscLinearSolver<T>::PetscLinearSolver ()
> + PetscLinearSolver<T>::PetscLinearSolver ():
> + _solve_only_on_dofs(NULL),
> + _solve_only_on_is(NULL),
> + _solve_only_on_is_seq(NULL)
> {
> if (libMesh::n_processors() == 1)
> this->_preconditioner_type = ILU_PRECOND;
> Index: include/solvers/linear_solver.h
> ===================================================================
> --- include/solvers/linear_solver.h (Revision 3957)
> +++ include/solvers/linear_solver.h (Arbeitskopie)
> @@ -24,6 +24,7 @@
>
>
> // C++ includes
> +#include <set>
>
> // Local includes
> #include "libmesh_common.h"
> @@ -117,6 +118,13 @@
> void attach_preconditioner(Preconditioner<T> * preconditioner);
>
> /**
> + * After calling this method, all successive solves will be limited
> + * to the given set of dofs. This mode can be disabled by calling
> + * this method with \p dofs being a \p NULL pointer.
> + */
> + virtual void solve_only_on (const std::set<unsigned int>* const dofs);
> +
> + /**
> * This function calls the solver
> * "_solver_type" preconditioned with the
> * "_preconditioner_type" preconditioner. Note that this method
> Index: src/solvers/petsc_linear_solver.C
> ===================================================================
> --- src/solvers/petsc_linear_solver.C (Revision 3957)
> +++ src/solvers/petsc_linear_solver.C (Arbeitskopie)
> @@ -98,6 +98,8 @@
> {
> if (this->initialized())
> {
> + this->solve_only_on(NULL);
> +
> this->_is_initialized = false;
>
> int ierr=0;
> @@ -360,9 +362,59 @@
>
>
>
> +template <typename T>
> +void
> +PetscLinearSolver<T>::solve_only_on (const std::set<unsigned int>* const
> dofs)
> +{
> + libmesh_assert(this->initialized());
>
> + int ierr=0;
>
> + if(_solve_only_on_is!=NULL)
> + {
> + ierr = ISDestroy(_solve_only_on_is);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + _solve_only_on_is = NULL;
> + }
>
> + if(_solve_only_on_is_seq!=NULL)
> + {
> + ierr = ISDestroy(_solve_only_on_is_seq);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + _solve_only_on_is_seq = NULL;
> + }
> +
> + _solve_only_on_dofs = dofs;
So the LinearSolver doesn't take ownership of this set? It seems
horrible to make the user responsible for that.
> +
> + if(dofs!=NULL)
> + {
> + PetscInt* petsc_dofs;
> + PetscInt* petsc_dofs_seq;
> + ierr = PetscMalloc(dofs->size()*sizeof(PetscInt),&petsc_dofs);
> + ierr = PetscMalloc(dofs->size()*sizeof(PetscInt),&petsc_dofs_seq);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + std::set<unsigned int>::const_iterator it = dofs->begin();
> + std::set<unsigned int>::const_iterator itEnd = dofs->end();
> + size_t count = 0;
> + while(it!=itEnd)
> + {
> + petsc_dofs[count] = *it;
> + petsc_dofs_seq[count] = count;
> + count++;
> + it++;
> + }
> + libmesh_assert(count=dofs->size());
> +
> + ierr =
> ISCreateGeneralNC(libMesh::COMM_WORLD,count,petsc_dofs,&_solve_only_on_is);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr =
> ISCreateGeneralNC(libMesh::COMM_WORLD,count,petsc_dofs_seq,&_solve_only_on_is_seq);
Could build this IS with
PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt
first,PetscInt step,IS *is);
but you don't need the IS at all, see below.
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> +}
> +
> +
> +
> template <typename T>
> std::pair<unsigned int, Real>
> PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in,
> @@ -406,6 +458,11 @@
>
> // 2.1.x & earlier style
> #if PETSC_VERSION_LESS_THAN(2,2,0)
> +
> + if(_solve_only_on_is!=NULL)
> + {
> + libmesh_not_implemented();
> + }
>
> // Set operators. The input matrix works as the preconditioning matrix
> ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
> @@ -431,6 +488,11 @@
> // 2.2.0
> #elif PETSC_VERSION_LESS_THAN(2,2,1)
>
> + if(_solve_only_on_is!=NULL)
> + {
> + libmesh_not_implemented();
> + }
> +
> // Set operators. The input matrix works as the preconditioning matrix
> //ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> // SAME_NONZERO_PATTERN);
> @@ -475,19 +537,34 @@
> // 2.2.1 & newer style
> #else
>
> - // Set operators. The input matrix works as the preconditioning matrix
> - if(!this->same_preconditioner)
> - {
> - ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> - SAME_NONZERO_PATTERN);
> - CHKERRABORT(libMesh::COMM_WORLD,ierr);
> - }
> + Mat submat = NULL;
> + Mat subprecond = NULL;
> +
> + // Set operators.
> + if(_solve_only_on_is!=NULL)
> + {
> + libmesh_assert(_solve_only_on_is_seq!=NULL);
> +
> + ierr = MatGetSubMatrix(matrix->mat(),
> + _solve_only_on_is,_solve_only_on_is,
> + PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = MatGetSubMatrix(precond->mat(),
> + _solve_only_on_is,_solve_only_on_is,
> + PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = KSPSetOperators(_ksp, submat, subprecond,
> + this->same_preconditioner ? SAME_PRECONDITIONER :
> SAME_NONZERO_PATTERN);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> else
> - {
> - ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> - SAME_PRECONDITIONER);
> - CHKERRABORT(libMesh::COMM_WORLD,ierr);
> - }
> + {
> + ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
> + this->same_preconditioner ? SAME_PRECONDITIONER :
> SAME_NONZERO_PATTERN);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
Use of this->same_preconditioner is an unrelated change, correct? Note
that there are more than two values in that enum. Also, you really have
a DIFFERENT_NONZERO_PATTERN if you just switched matrices to solve on a
different index set.
>
> // Set the tolerances for the iterative solver. Use the user-supplied
> // tolerance for the relative residual & leave the others at default
> values.
> @@ -496,8 +573,56 @@
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> // Solve the linear system
> - ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
> - CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + if(_solve_only_on_is!=NULL)
> + {
> + libmesh_assert(_solve_only_on_is_seq!=NULL);
> + libmesh_assert(_solve_only_on_dofs!=NULL);
> + Vec subrhs;
> + Vec subsolution;
> +
> + ierr =
> VecCreateSeq(libMesh::COMM_WORLD,_solve_only_on_dofs->size(),&subrhs);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr =
> VecCreateSeq(libMesh::COMM_WORLD,_solve_only_on_dofs->size(),&subsolution);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
Why are you specializing this to only work in serial?
> +
> + VecScatter scatter;
> + ierr = VecScatterCreate(rhs->vec(),_solve_only_on_is,
> + subrhs,_solve_only_on_is_seq,
> + &scatter);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
Can use NULL for the second IS, it means "scatter to the whole vector".
> +
> + ierr =
> VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr =
> VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr =
> VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr =
> VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = KSPSolve (_ksp, subrhs, subsolution);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr =
> VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr =
> VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = VecScatterDestroy(scatter);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> +
> + ierr = VecDestroy(subsolution);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = VecDestroy(subrhs);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> + else
> + {
> + ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
>
> // Get the number of iterations required for convergence
> ierr = KSPGetIterationNumber (_ksp, &its);
> @@ -507,6 +632,14 @@
> ierr = KSPGetResidualNorm (_ksp, &final_resid);
> CHKERRABORT(libMesh::COMM_WORLD,ierr);
>
> + if(_solve_only_on_is!=NULL)
> + {
> + ierr = MatDestroy(submat);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + ierr = MatDestroy(subprecond);
> + CHKERRABORT(libMesh::COMM_WORLD,ierr);
> + }
> +
> #endif
>
> STOP_LOG("solve()", "PetscLinearSolver");
> Index: src/solvers/linear_solver.C
> ===================================================================
> --- src/solvers/linear_solver.C (Revision 3957)
> +++ src/solvers/linear_solver.C (Arbeitskopie)
> @@ -115,7 +115,18 @@
> _preconditioner = preconditioner;
> }
>
> +template <typename T>
> +void
> +LinearSolver<T>::solve_only_on(const std::set<unsigned int>* const dofs)
> +{
> + if(dofs!=NULL)
> + {
> + libmesh_not_implemented();
> + }
> +}
If dofs == NULL, this function is a no-op?
>
> +
> +
> //------------------------------------------------------------------
> // Explicit instantiations
> template class LinearSolver<Number>;
------------------------------------------------------------------------------
Start uncovering the many advantages of virtual appliances
and start using them to simplify application deployment and
accelerate your shift to cloud computing
http://p.sf.net/sfu/novell-sfdev2dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
