On Mon, 13 Sep 2010, Jed Brown wrote:
On Mon, 13 Sep 2010 11:35:17 +0200 (CEST), Tim Kroeger
<[email protected]> wrote:
On Thu, 9 Sep 2010, Jed Brown wrote:
Do you want to assemble X, or are you really only working with X2? If
the former, you can MatGetSubMatrix (you just need an index set defining
the subdomain) the part you want and solve with that.
What will I have to do with the right hand side and the solution
vectors then? I assumed there would be something like
VecGetSubVector(), but that doesn't seem to be true. Or is there no
need to do anything at all?
PCFieldSplit does this for you, if you are extracting submatrices
yourself, then you would create a VecScatter to extract the part of the
vector you want (with parallel redistribution if desired).
I have now implemented it for one of the solve() methods of the
PetscLinearSolver class. Since my experience with PETSc is still low,
I'm somehow unsure whether the things I did are the correct ones.
Jed, can you have a look at it and let me know whether this looks
right?
If it does, I will continue doing the same thing for the remaining
solve() methods and then work on the System class following Roy's idea
with the subset class. I will not check in anything before I have
extended my application to successfully use this, but this may take
some time.
Best Regards,
Tim
--
Dr. Tim Kroeger
CeVis -- Center of Complex Systems and Visualization
University of Bremen [email protected]
Universitaetsallee 29 [email protected]
D-28359 Bremen Phone +49-421-218-7710
Germany Fax +49-421-218-4236
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;
+
+ /**
+ * 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;
+
};
/*----------------------- 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;
+
+ 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);
+ 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);
+ }
// 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);
+
+ VecScatter scatter;
+ ierr = VecScatterCreate(rhs->vec(),_solve_only_on_is,
+ subrhs,_solve_only_on_is_seq,
+ &scatter);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ 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();
+ }
+}
+
+
//------------------------------------------------------------------
// 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
