Hi,
I have had to make a few changes to the dof_map and petsc_matrix
files (both .h and .C) to accommodate some needs of my work. I am
solving a heat transfer problem, which involves conduction (solved
using finite elements) and internal cavity radiation due to heat
emission/absorption/reflection. Due to this, hundreds of nodes get
coupled with each other and the sparsity pattern has a large number of
non-zeros per row.
I am using petsc as the solver and for my application, and I have
found that providing the detailed sparsity pattern to the matrix
creation routine becomes important. Also, due to the way Petsc adds
new non-zeros in each row, I found it was important to specify all the
non-zero columns per row during the initialization (it improved the
performance of my code by an order of magnitude). Since the library
did not have an update-sparsity pattern for the petsc-matrix, I
implemented it and have attached it the following diffs.
The change in dof_map is not really necessary, but since I am
initializing a few matrices much deeper into the code (rather than
only the beginning), I made this change to store the sparsity_pattern
in the dof_map object.
Please use / ignore them as you see fit.
Regards,
Manav
Index: /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/
include/base/dof_map.h
===================================================================
--- /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/include/
base/dof_map.h (revision 3098)
+++ /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/include/
base/dof_map.h (working copy)
@@ -274,6 +274,16 @@
*/
void compute_sparsity (const MeshBase&);
+ /**
+ * @returns \p SparsityPattern::Graph object of this \p DofMap
+ */
+ const SparsityPattern::Graph& get_sparsity_pattern() const
+ {
+ libmesh_assert(this->_sp.get() != NULL);
+ return this->_sp->sparsity_pattern;
+ }
+
+
/**
* Returns a constant reference to the \p _send_list for this
processor. The
* \p _send_list contains the global indices of all the variables
in the
@@ -783,6 +793,11 @@
PeriodicBoundaries _periodic_boundaries;
#endif
+ /**
+ * sparsity pattern graph object for this dof map
+ */
+ AutoPtr<SparsityPattern::Build> _sp;
+
friend class SparsityPattern::Build;
};
Index: /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/src/
base/dof_map.C
===================================================================
--- /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/src/base/
dof_map.C (revision 3069)
+++ /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/src/base/
dof_map.C (working copy)
@@ -1116,26 +1116,37 @@
// Even better, if the full sparsity pattern is not needed then
// the number of nonzeros per row can be estimated from the
// sparsity patterns created on each thread.
- SparsityPattern::Build sp (mesh,
- *this,
- _dof_coupling,
- implicit_neighbor_dofs,
- need_full_sparsity_pattern);
+ // *** this has been changed to a member of the object since it
will be stored
+ // and used for matrices beyond this function ***
+ this->_sp.reset(new SparsityPattern::Build(mesh,
+ *this,
+ _dof_coupling,
+ implicit_neighbor_dofs,
+
need_full_sparsity_pattern));
Threads::parallel_reduce (ConstElemRange
(mesh.active_elements_begin(),
- mesh.active_elements_end()), sp);
+
mesh.active_elements_end()), *(this->_sp.get()));
#ifndef NDEBUG
// Avoid declaring these variables unless asserts are enabled.
const unsigned int proc_id = mesh.processor_id();
const unsigned int n_dofs_on_proc = this-
>n_dofs_on_processor(proc_id);
#endif
- libmesh_assert (sp.sparsity_pattern.size() == n_dofs_on_proc);
+ libmesh_assert (this->_sp->sparsity_pattern.size() ==
n_dofs_on_proc);
// steal the n_nz and n_oz arrays from sp -- it won't need them
any more,
// and this is more efficient than copying them.
- _n_nz.swap(sp.n_nz);
- _n_oz.swap(sp.n_oz);
+ //_n_nz.swap(sp.n_nz);
+ //_n_oz.swap(sp.n_oz);
+ // *** this has been changed to copy since the sparsity pattern is
being stored
+ // in this object. This is because the matrices will not
necessarily be initialized
+ // only once through this method, but at arbitrary times through
the sparsity_pattern***
+ _n_nz.resize(this->_sp->n_nz.size());
+ _n_oz.resize(this->_sp->n_oz.size());
+ for (unsigned int i=0; i< _n_nz.size(); i++)
+ _n_nz[i] = this->_sp->n_nz[i];
+ for (unsigned int i=0; i< _n_oz.size(); i++)
+ _n_oz[i] = this->_sp->n_oz[i];
STOP_LOG("compute_sparsity()", "DofMap");
@@ -1147,7 +1158,7 @@
end = _matrices.end();
for (; pos != end; ++pos)
- (*pos)->update_sparsity_pattern (sp.sparsity_pattern);
+ (*pos)->update_sparsity_pattern (this->_sp->sparsity_pattern);
}
Index: /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/
include/numerics/petsc_matrix.h
===================================================================
--- /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/include/
numerics/petsc_matrix.h (revision 3062)
+++ /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/include/
numerics/petsc_matrix.h (working copy)
@@ -32,9 +32,11 @@
+
+ /**
+ * The \p PetscMatrix needs the full sparsity pattern.
+ */
+ bool need_full_sparsity_pattern () const
+ { return true; }
+
+ /**
+ * Updates the matrix sparsity pattern.
+ */
+ virtual void update_sparsity_pattern (const SparsityPattern::Graph
&);
+
Index: /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/src/
numerics/petsc_matrix.C
===================================================================
--- /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/src/
numerics/petsc_matrix.C (revision 3062)
+++ /Users/manav/Documents/codes/fem/FESystem/contrib/libmesh/src/
numerics/petsc_matrix.C (working copy)
@@ -164,12 +164,106 @@
CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
+ // now update the sparsity pattern
+ const SparsityPattern::Graph& sparsity_pattern =
+ this->_dof_map->get_sparsity_pattern();
+ this->update_sparsity_pattern(sparsity_pattern);
+
this->zero();
}
+
template <typename T>
+void PetscMatrix<T>::update_sparsity_pattern
+(const SparsityPattern::Graph &sparsity_pattern)
+{
+ // make sure that the matrix has been initialized
+ assert (this->initialized());
+ unsigned int first_local_dof = this->_dof_map->first_dof();
+ unsigned int last_local_dof = this->_dof_map->last_dof();
+
+#ifdef DEBUG
+ assert ((last_local_dof - first_local_dof + 1) ==
sparsity_pattern.size());
+
+ // it is assumed that the values are sorted, but it can be checked
for debugging
+ for (unsigned int i=0; i < sparsity_pattern.size(); i++)
+ {
+ for (unsigned int j=1; j < sparsity_pattern[i].size(); j++)
+ assert (sparsity_pattern[i][j-1] < sparsity_pattern[i][j]);
+ }
+#endif
+
+
+ PetscErrorCode ierr=0;
+ std::vector<double> vals;
+ unsigned int row_id = 0;
+
+ ierr = MatSetOption (_mat, MAT_YES_NEW_NONZERO_LOCATIONS);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // iterate over all rows, and add the entries to the
+ for (unsigned int i=0; i < sparsity_pattern.size(); i++)
+ {
+ if (vals.size() != sparsity_pattern[i].size())
+ {
+ vals.resize(sparsity_pattern[i].size());
+ std::fill(vals.begin(), vals.end(), 1.0);
+ }
+
+ row_id = first_local_dof + i;
+
+ // now tell the matrix to set the values for this row
+ ierr = MatSetValues (_mat,
+ 1, (PetscInt*) &row_id,
+ sparsity_pattern[i].size(), (PetscInt*)
&sparsity_pattern[i][0],
+ (PetscScalar*) &vals[0], INSERT_VALUES);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+
+ ierr = MatSetOption (_mat, MAT_NEW_NONZERO_LOCATION_ERR);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = MatSetOption (_mat, MAT_KEEP_ZEROED_ROWS);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // finally, call the assemble functions
+ ierr = MatAssemblyBegin (_mat, MAT_FLUSH_ASSEMBLY);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = MatAssemblyEnd (_mat, MAT_FLUSH_ASSEMBLY);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ this->zero();
+}
+
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
