On Wed, 15 Sep 2010, Jed Brown wrote:
On Wed, 15 Sep 2010 08:42:16 +0200 (CEST), Tim Kroeger
<[email protected]> wrote:
(One conjecture why SAME_NONZERO_PATTERN has been used here before is
that somebody thought it means that the system matrix and the
preconditioner matrix have the same nonzero pattern, but as far as I
understand now, this is just not true.)
It means that the nonzero pattern of the preconditioning matrix has not
changed since the last solve.
Yes, that's also what I understood now. Since at no point in the code
information about the nonzero pattern of the last solve is stored,
there is no reason for using SAME_NONZERO_PATTERN at all.
Another idea is to use VecCreateMPI(), giving only the global length
and using PETSC_DECIDE for the local length. This way, approximately
equal amounts of dofs are assigned to each processor, but will they
correspond to patches of the computational domain, or will they be
scattered all around? Well, since we are using a vector, this is of
course in the responsiblity of the user.
This chooses the sizes of each part, but you will still have to move
indices around.
This is done by VecScatter, isn't it? I think that the overhead of
doing this once before and after the solve will be small against the
gain in performance by having balanced vectors. Well, of course this
depends on a lot of things, but in the mean...
For a high-quality partition, you could extract the
submatrix with the naive (current) partition, then run a partitioner on
the submatrix, and finally extract indices using the resulting "good"
partition. Of course this may not be worthwhile if you will only solve
with it once.
I leave this for somebody else. (-:
Anyway, I have now attached the current state of my extension.
Questions are:
* (mainly to Jed:) Are the things in petsc_linear_solver.{h,C} now
correct, sensible, and reasonably efficient? If yes, I will do the
same things to the remaining overloaded solve() methods.
* (mainly ro Roy:) Do you like the API?
* (mainly to Roy:) Do you have some suggestion for a better name for
solve_only_on()? (Jed suggested set_active_domain().)
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/systems/system_subset_by_subdomain.h
===================================================================
--- include/systems/system_subset_by_subdomain.h (Revision 0)
+++ include/systems/system_subset_by_subdomain.h (Revision 0)
@@ -0,0 +1,173 @@
+// $Id: system.h 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+#ifndef __system_subset_by_subdomain_h__
+#define __system_subset_by_subdomain_h__
+
+// C++ includes
+#include <set>
+
+// Local Includes
+#include "system_subset.h"
+#include "id_types.h"
+
+namespace libMesh
+{
+
+// Forward Declarations
+
+/**
+ * This class represents a subset of the dofs of a \p System, selected
+ * by the \p subdomain_id and possible the variable numbers. The dofs
+ * in the subset will be sorted.
+ *
+ * @author Tim Kroeger, 2010.
+ */
+
+// ------------------------------------------------------------
+// SystemSubset class definition
+class SystemSubsetBySubdomain : public SystemSubset
+{
+public:
+
+ /**
+ * Subclass for user-specified selection of subdomain ids to be
+ * included in a \p SystemSubset.
+ */
+ class SubdomainSelection : public ReferenceCountedObject<SubdomainSelection>
+ {
+ public:
+
+ /**
+ * Constructor.
+ */
+ SubdomainSelection (void);
+
+ /**
+ * Destructor.
+ */
+ virtual ~SubdomainSelection (void);
+
+ /**
+ * Method that decides whether a given subdomain id is included in
+ * the subset or nor.
+ */
+ virtual bool operator()(const subdomain_id_type& subdomain_id)const=0;
+
+ private:
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure
+ * that the compiler doesn't implement a default.
+ */
+ SubdomainSelection(const SubdomainSelection&);
+
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure
+ * that the compiler doesn't implement a default.
+ */
+ SubdomainSelection& operator=(const SubdomainSelection&);
+
+ }; // subclass \p SubdomainSelection
+
+ /**
+ * Selection of subdomain ids by a list.
+ */
+ class SubdomainSelectionByList : public SubdomainSelection
+ {
+ public:
+ /**
+ * Constructor. Does not take a copy of the \p list, so make sure
+ * that the \p list does not go out of scope before \p *this does.
+ */
+ SubdomainSelectionByList (const std::set<subdomain_id_type>& list);
+
+ /**
+ * Method that decides whether a given subdomain id is included in
+ * the subset or nor.
+ */
+ virtual bool operator()(const subdomain_id_type& subdomain_id)const;
+
+ protected:
+ /**
+ * The actual list.
+ */
+ const std::set<subdomain_id_type>& _list;
+ };
+
+ /**
+ * Constructor. The subset will consist of those dofs which are
+ * associated to at least one mesh element that has a subdomain id
+ * contained in the \p subdomain_selection. If \p var_nums is not a
+ * \p NULL pointer, dofs that are associated to a variable number
+ * that is not contained in \p var_nums will not contain to the
+ * subset, no matter what elements they belong to.
+ */
+ SystemSubsetBySubdomain (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums = NULL);
+
+ /**
+ * Constructor. The subset will consist of those dofs which are
+ * associated to at least one mesh element that has a subdomain id
+ * contained in the set \p subdomain_ids. If \p var_nums is not a
+ * \p NULL pointer, dofs that are associated to a variable number
+ * that is not contained in \p var_nums will not contain to the
+ * subset, no matter what elements they belong to.
+ */
+ SystemSubsetBySubdomain (const System& system,
+ const std::set<subdomain_id_type>& subdomain_ids,
+ const std::set<unsigned int>* const var_nums = NULL);
+
+ /**
+ * Destructor.
+ */
+ virtual ~SystemSubsetBySubdomain ();
+
+ /**
+ * Method that returns the actual set of dofs that the subset
+ * consists of.
+ */
+ virtual const std::vector<unsigned int>& dof_ids(void)const;
+
+protected:
+
+ /**
+ * Initializes the class. Will be called by the constructors.
+ */
+ void init (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums);
+
+ /**
+ * The actual set of the dof ids.
+ */
+ std::vector<unsigned int> _dof_ids;
+
+}; // class SystemSubset
+
+} // namespace libMesh
+
+#endif
+
Index: include/systems/system.h
===================================================================
--- include/systems/system.h (Revision 3957)
+++ include/systems/system.h (Arbeitskopie)
@@ -56,6 +56,7 @@
template <typename T> class VectorValue;
typedef VectorValue<Number> NumberVectorValue;
typedef NumberVectorValue Gradient;
+class SystemSubset;
/**
* This is the base class for classes which contain
@@ -162,6 +163,13 @@
* This method is only implemented in some derived classes.
*/
virtual void assemble_residual_derivatives (const ParameterVector&
parameters);
+
+ /**
+ * After calling this method, any solve will be limited to the given
+ * subdomain. To disable this mode, call this method with \p subset
+ * being a \p NULL pointer.
+ */
+ virtual void solve_only_on (const SystemSubset* subset);
/**
* Solves the system. Must be overloaded in derived systems.
Index: include/systems/linear_implicit_system.h
===================================================================
--- include/systems/linear_implicit_system.h (Revision 3957)
+++ include/systems/linear_implicit_system.h (Arbeitskopie)
@@ -100,6 +100,13 @@
virtual void assemble () { ImplicitSystem::assemble(); }
/**
+ * After calling this method, any solve will be limited to the given
+ * subdomain. To disable this mode, call this method with \p subset
+ * being a \p NULL pointer.
+ */
+ virtual void solve_only_on (const SystemSubset* subset);
+
+ /**
* Assembles & solves the linear system A*x=b.
*/
virtual void solve ();
@@ -187,6 +194,11 @@
*/
ShellMatrix<Number>* _shell_matrix;
+ /**
+ * The current subset on which to solve (or \p NULL if none).
+ */
+ const SystemSubset* _subset;
+
};
Index: include/systems/system_subset.h
===================================================================
--- include/systems/system_subset.h (Revision 0)
+++ include/systems/system_subset.h (Revision 0)
@@ -0,0 +1,102 @@
+// $Id: system.h 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+#ifndef __system_subset_h__
+#define __system_subset_h__
+
+// C++ includes
+#include <vector>
+
+// Local Includes
+#include "reference_counted_object.h"
+
+namespace libMesh
+{
+
+// Forward Declarations
+class System;
+
+/**
+ * This is a base class for classes which represent subsets of the
+ * dofs of a \p System.
+ *
+ * @author Tim Kroeger, 2010.
+ */
+
+// ------------------------------------------------------------
+// SystemSubset class definition
+class SystemSubset : public ReferenceCountedObject<SystemSubset>
+{
+public:
+
+ /**
+ * Constructor.
+ */
+ SystemSubset (const System& system);
+
+public:
+
+ /**
+ * Destructor.
+ */
+ virtual ~SystemSubset (void);
+
+ /**
+ * Method that returns the actual set of dofs that the subset
+ * consists of.
+ */
+ virtual const std::vector<unsigned int>& dof_ids(void)const=0;
+
+ /**
+ * Returns the \p System to which we belong.
+ */
+ const System& get_system(void)const;
+
+protected:
+
+ /**
+ * A reference to the \p System we belong to.
+ */
+ const System& _system;
+
+private:
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure that
+ * the compiler doesn't implement a default.
+ */
+ SystemSubset(const SystemSubset&);
+
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure that
+ * the compiler doesn't implement a default.
+ */
+ SystemSubset& operator=(const SystemSubset&);
+
+}; // class SystemSubset
+
+} // namespace libMesh
+
+#endif
+
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::vector<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,35 @@
* Krylov subspace context
*/
KSP _ksp;
+
+ /**
+ * 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 the <em>local</em> dofs on which to
+ * solve. This index set is determined when a system is solved on a
+ * subset the first time. Before that point (or if no subset is
+ * selected), it will be a \p NULL pointer.
+ */
+ IS _solve_only_on_is_local;
+
+ /**
+ * Internal method that returns the size of \p _solve_only_on_is.
+ */
+ size_t _solve_only_on_is_size(void)const;
+
};
/*----------------------- functions ----------------------------------*/
template <typename T>
inline
-PetscLinearSolver<T>::PetscLinearSolver ()
+ PetscLinearSolver<T>::PetscLinearSolver ():
+ _solve_only_on_is(NULL),
+ _solve_only_on_is_local(NULL)
{
if (libMesh::n_processors() == 1)
this->_preconditioner_type = ILU_PRECOND;
@@ -277,6 +307,23 @@
}
+
+template <typename T>
+inline size_t
+PetscLinearSolver<T>::
+_solve_only_on_is_size(void)const
+{
+ libmesh_assert(_solve_only_on_is!=NULL);
+
+ PetscInt s;
+ int ierr = ISGetSize(_solve_only_on_is,&s);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ return static_cast<size_t>(s);
+}
+
+
+
} // namespace libMesh
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 <vector>
// Local includes
#include "libmesh_common.h"
@@ -117,6 +118,14 @@
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. The vector of
+ * dofs must not contain any duplicates.
+ */
+ virtual void solve_only_on (const std::vector<unsigned int>* const dofs);
+
+ /**
* This function calls the solver
* "_solver_type" preconditioned with the
* "_preconditioner_type" preconditioner. Note that this method
Index: src/systems/linear_implicit_system.C
===================================================================
--- src/systems/linear_implicit_system.C (Revision 3957)
+++ src/systems/linear_implicit_system.C (Arbeitskopie)
@@ -28,6 +28,7 @@
#include "numeric_vector.h" // for parameter sensitivity calcs
#include "parameter_vector.h"
#include "sparse_matrix.h" // for get_transpose
+#include "system_subset.h"
namespace libMesh
{
@@ -43,7 +44,8 @@
linear_solver (LinearSolver<Number>::build()),
_n_linear_iterations (0),
_final_linear_residual (1.e20),
- _shell_matrix(NULL)
+ _shell_matrix(NULL),
+ _subset(NULL)
{
}
@@ -61,6 +63,8 @@
{
// clear the linear solver
linear_solver->clear();
+
+ this->solve_only_on(NULL);
// clear the parent data
Parent::clear();
@@ -79,6 +83,17 @@
+void LinearImplicitSystem::solve_only_on (const SystemSubset* subset)
+{
+ _subset = subset;
+ if(subset!=NULL)
+ {
+ libmesh_assert(&subset->get_system()==this);
+ }
+}
+
+
+
void LinearImplicitSystem::solve ()
{
if (this->assemble_before_solve)
@@ -101,6 +116,11 @@
const unsigned int maxits =
es.parameters.get<unsigned int>("linear solver maximum iterations");
+ if(_subset!=NULL)
+ {
+ linear_solver->solve_only_on(&_subset->dof_ids());
+ }
+
// Solve the linear system. Several cases:
std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
if(_shell_matrix)
@@ -110,6 +130,11 @@
// 2.) No shell matrix, with or without user-supplied preconditioner
rval = linear_solver->solve (*matrix,
this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
+ if(_subset!=NULL)
+ {
+ linear_solver->solve_only_on(NULL);
+ }
+
// Store the number of linear iterations required to
// solve and the final residual.
_n_linear_iterations = rval.first;
Index: src/systems/system_subset.C
===================================================================
--- src/systems/system_subset.C (Revision 0)
+++ src/systems/system_subset.C (Revision 0)
@@ -0,0 +1,52 @@
+// $Id: system.C 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+// C++ includes
+
+
+// Local includes
+#include "system_subset.h"
+
+namespace libMesh
+{
+
+// ------------------------------------------------------------
+// SystemSubset implementation
+ SystemSubset::SystemSubset (const System& system):
+ _system(system)
+ {
+ }
+
+
+ SystemSubset::~SystemSubset (void)
+ {
+ }
+
+
+ const System&
+ SystemSubset::get_system(void)const
+ {
+ return _system;
+ }
+
+
+} // namespace libMesh
+
Index: src/systems/system_subset_by_subdomain.C
===================================================================
--- src/systems/system_subset_by_subdomain.C (Revision 0)
+++ src/systems/system_subset_by_subdomain.C (Revision 0)
@@ -0,0 +1,144 @@
+// $Id: system.C 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+// C++ includes
+
+
+// Local includes
+#include "system_subset_by_subdomain.h"
+#include "system.h"
+#include "dof_map.h"
+
+namespace libMesh
+{
+// ------------------------------------------------------------
+// SubdomainSelectio implementation
+ SystemSubsetBySubdomain::SubdomainSelection::
+ SubdomainSelection (void)
+ {
+ }
+
+ SystemSubsetBySubdomain::SubdomainSelection::
+ ~SubdomainSelection (void)
+ {
+ }
+
+ SystemSubsetBySubdomain::SubdomainSelectionByList::
+ SubdomainSelectionByList (const std::set<subdomain_id_type>& list):
+ _list(list)
+ {
+ }
+
+ bool
+ SystemSubsetBySubdomain::SubdomainSelectionByList::
+ operator()(const subdomain_id_type& subdomain_id)const
+ {
+ return _list.find(subdomain_id)!=_list.end();
+ }
+
+// ------------------------------------------------------------
+// SystemSubsetBySubdomain implementation
+
+ SystemSubsetBySubdomain::
+ SystemSubsetBySubdomain (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums):
+ SystemSubset(system),
+ _dof_ids()
+ {
+ this->init(system,subdomain_selection,var_nums);
+ }
+
+ SystemSubsetBySubdomain::
+ SystemSubsetBySubdomain (const System& system,
+ const std::set<subdomain_id_type>& subdomain_ids,
+ const std::set<unsigned int>* const var_nums):
+ SystemSubset(system),
+ _dof_ids()
+ {
+ SubdomainSelectionByList selection(subdomain_ids);
+ this->init(system,selection,var_nums);
+ }
+
+ SystemSubsetBySubdomain::
+ ~SystemSubsetBySubdomain (void)
+ {
+ }
+
+ const std::vector<unsigned int>&
+ SystemSubsetBySubdomain::
+ dof_ids(void)const
+ {
+ return _dof_ids;
+ }
+
+ void
+ SystemSubsetBySubdomain::
+ init (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums)
+ {
+ std::set<unsigned int> effective_var_nums;
+ if(var_nums!=NULL)
+ {
+ effective_var_nums = *var_nums;
+ }
+ else
+ {
+ for(unsigned int i=0; i<system.n_vars(); i++)
+ {
+ effective_var_nums.insert(i);
+ }
+ }
+
+ const DofMap & dof_map = system.get_dof_map();
+ std::vector<unsigned int> dof_indices;
+
+ const MeshBase& mesh = system.get_mesh();
+ MeshBase::const_element_iterator el =
mesh.active_elements_begin();
+ const MeshBase::const_element_iterator end_el =
mesh.active_elements_end();
+ for ( ; el != end_el; ++el)
+ {
+ const Elem* elem = *el;
+ if(subdomain_selection(elem->subdomain_id()))
+ {
+ std::set<unsigned int>::const_iterator it =
effective_var_nums.begin();
+ const std::set<unsigned int>::const_iterator itEnd =
effective_var_nums.end();
+ for (; it!=itEnd; ++it)
+ {
+ dof_map.dof_indices (elem, dof_indices, *it);
+ for(size_t i=0; i<dof_indices.size(); i++)
+ {
+ _dof_ids.push_back(dof_indices[i]);
+ }
+ }
+ }
+ }
+
+ /* Sort and unique the vector (using the same mechanism as in \p
+ DofMap::prepare_send_list()). */
+ std::sort(_dof_ids.begin(), _dof_ids.end());
+ std::vector<unsigned int>::iterator new_end = std::unique
(_dof_ids.begin(), _dof_ids.end());
+ std::vector<unsigned int> (_dof_ids.begin(), new_end).swap (_dof_ids);
+ }
+
+} // namespace libMesh
+
Index: src/systems/system.C
===================================================================
--- src/systems/system.C (Revision 3957)
+++ src/systems/system.C (Arbeitskopie)
@@ -367,6 +367,16 @@
+void System::solve_only_on (const SystemSubset* subset)
+{
+ if(subset!=NULL)
+ {
+ libmesh_not_implemented();
+ }
+}
+
+
+
void System::assemble ()
{
// Log how long the user's matrix assembly code takes
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;
@@ -294,7 +296,7 @@
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set operators. The input matrix works as the preconditioning matrix
- ierr = KSPSetOperators(_ksp, matrix->mat(),
matrix->mat(),SAME_NONZERO_PATTERN);
+ ierr = KSPSetOperators(_ksp, matrix->mat(),
matrix->mat(),DIFFERENT_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set user-specified solver and preconditioner types
@@ -360,9 +362,46 @@
+template <typename T>
+void
+PetscLinearSolver<T>::solve_only_on (const std::vector<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_local!=NULL)
+ {
+ ierr = ISDestroy(_solve_only_on_is_local);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ _solve_only_on_is_local = NULL;
+ }
+
+ if(dofs!=NULL)
+ {
+ PetscInt* petsc_dofs;
+ ierr = PetscMalloc(dofs->size()*sizeof(PetscInt),&petsc_dofs);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ for(size_t i=0; i<dofs->size(); i++)
+ {
+ petsc_dofs[i] = (*dofs)[i];
+ }
+
+ ierr =
ISCreateGeneralNC(libMesh::COMM_WORLD,dofs->size(),petsc_dofs,&_solve_only_on_is);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+}
+
+
+
template <typename T>
std::pair<unsigned int, Real>
PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in,
@@ -406,10 +445,15 @@
// 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(),
- SAME_NONZERO_PATTERN);
+ DIFFERENT_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the tolerances for the iterative solver. Use the user-supplied
@@ -431,6 +475,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 +524,77 @@
// 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;
+ Vec subrhs;
+ Vec subsolution;
+ VecScatter scatter;
+
+ // Set operators. Also restrict rhs and solution vector to
+ // subdomain if neccessary.
+ if(_solve_only_on_is!=NULL)
+ {
+ size_t is_size = this->_solve_only_on_is_size();
+
+ ierr = VecCreate(libMesh::COMM_WORLD,&subrhs);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetSizes(subrhs,PETSC_DECIDE,is_size);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetFromOptions(subrhs);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecCreate(libMesh::COMM_WORLD,&subsolution);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetSizes(subsolution,PETSC_DECIDE,is_size);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetFromOptions(subsolution);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecScatterCreate(rhs->vec(),_solve_only_on_is, subrhs,NULL,
&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);
+
+ /* Make the local index set if it doesn't exist yet. */
+ if(_solve_only_on_is_local==NULL)
+ {
+ PetscInt firstIndex = 0;
+ PetscInt lastIndex = 0;
+ ierr = VecGetOwnershipRange(subrhs,&firstIndex,&lastIndex);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr =
ISCreateStride(libMesh::COMM_WORLD,lastIndex-firstIndex,firstIndex,1,&_solve_only_on_is_local);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+
+ ierr = MatGetSubMatrix(matrix->mat(),
+ _solve_only_on_is_local,_solve_only_on_is,
+ PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = MatGetSubMatrix(precond->mat(),
+ _solve_only_on_is_local,_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 :
DIFFERENT_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 :
DIFFERENT_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 +603,16 @@
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)
+ {
+ ierr = KSPSolve (_ksp, subrhs, subsolution);
+ 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 +622,26 @@
ierr = KSPGetResidualNorm (_ksp, &final_resid);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ if(_solve_only_on_is!=NULL)
+ {
+ 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);
+ ierr = MatDestroy(submat);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = MatDestroy(subprecond);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+
#endif
STOP_LOG("solve()", "PetscLinearSolver");
@@ -574,7 +709,7 @@
// Set operators. The input matrix works as the preconditioning matrix
ierr = KSPSetOperators(_ksp, mat, mat,
- SAME_NONZERO_PATTERN);
+ DIFFERENT_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the tolerances for the iterative solver. Use the user-supplied
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::vector<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
