Dear Roy,
On Thu, 22 Jan 2009, Roy Stogner wrote:
Currently, ex10 on 2 CPUs dies after the first refinement, when
operator() fails at TransientSystem::old_local_solution.operator()(i),
where i is outside of the local dofs on that processor, because at the
time the _global_to_local_map is empty.
I think I understand why it's failing now: because although with the
patch we initialize old_local_solution as a ghosted vector in
TransientSystem::init_data(), those vectors are getting reinitialized
as distributed vectors in System::project_vector() - the code tests
size() against local_size() to determine whether a vector is serial or
distributed, and when it sees a ghost vector it incorrectly concludes
the latter.
I invented two new methods in NumericVector that might help you: One
is NumericVector::is_ghosted(), that you can use to query whether a
given vector does contain ghost dofs, and the other is
NumericVector::init(const NumericVector&), which will initialize the
vector using the same storage scheme (i.e. parallel/serial/ghosted) as
its argument. See attached patch (which unfortunately again includes
your patch, but you get quite am impression about what I changed if
you diff the two patches against each other).
The offending code is lines 100-118 of system_projection.C, which I
think is some of the least comprehensible code I've ever added to
libMesh. And off the top of my head I'm not sure what the right
solution is.
So there's no reason for me to be worried if I don't understand that
code, is there? (-:
Do my new methods help?
Best Regards,
Tim
--
Dr. Tim Kroeger
tim.kroe...@mevis.fraunhofer.de Phone +49-421-218-7710
tim.kroe...@cevis.uni-bremen.de Fax +49-421-218-4236
Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany
Index: include/numerics/petsc_vector.h
===================================================================
--- include/numerics/petsc_vector.h (Revision 3231)
+++ include/numerics/petsc_vector.h (Arbeitskopie)
@@ -154,6 +154,20 @@
const std::vector<unsigned int>& /*ghost*/,
const bool /*fast*/ = false);
+ /**
+ * Creates a vector that has the same dimension and storage type as
+ * \p other, including ghost dofs.
+ */
+ virtual void init (const NumericVector<T>& other,
+ const bool fast = false);
+
+ /**
+ * Will return \p true if the vector contains ghost dofs. Default
+ * implementation is for solver packages for which ghost dofs are
+ * not yet implemented in libMesh.
+ */
+ virtual bool is_ghosted (void) const;
+
// /**
// * Change the dimension to that of the
// * vector \p V. The same applies as for
@@ -654,6 +668,13 @@
int petsc_n=static_cast<int>(n);
int petsc_n_local=static_cast<int>(n_local);
int petsc_n_ghost=static_cast<int>(ghost.size());
+
+ // If the mesh is disjoint, the following assertion will fail.
+ // If the mesh is not disjoint, every processor will either have
+ // all the dofs, none of the dofs, or some non-zero dofs at the
+ // boundary between processors.
+ libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
+
int* petsc_ghost = ghost.empty() ? PETSC_NULL :
const_cast<int*>(reinterpret_cast<const int*>(&ghost[0]));
@@ -686,6 +707,48 @@
template <typename T>
inline
+void PetscVector<T>::init (const NumericVector<T>& other,
+ const bool fast)
+{
+ const PetscVector<T>* other_petsc = libmesh_cast_ptr<const
PetscVector<T>*>(&other);
+
+ /* Check whether \p other is ghosted. */
+ if(!other_petsc->is_ghosted())
+ {
+ /* Not ghosted, do the default. */
+ NumericVector<T>::init(other,fast);
+ }
+ else
+ {
+ /* Ghosted. Hence, we need to find the ghost dofs of \p
+ other. */
+ std::vector<unsigned int>
ghost_vector(other_petsc->_global_to_local_map.size(),0);
+ GlobalToLocalMap::const_iterator it =
other_petsc->_global_to_local_map.begin();
+ GlobalToLocalMap::const_iterator it_end =
other_petsc->_global_to_local_map.end();
+ for(;it!=it_end; ++it)
+ {
+ libmesh_assert(it->second<ghost_vector.size());
+ ghost_vector[it->second] = it->first;
+ }
+ this->init(other.size(),other.local_size(),ghost_vector,fast);
+ }
+}
+
+
+
+
+template <typename T>
+inline
+bool PetscVector<T>::is_ghosted (void) const
+{
+ return !this->_global_to_local_map.empty();
+}
+
+
+
+
+template <typename T>
+inline
void PetscVector<T>::close ()
{
libmesh_assert (this->initialized());
Index: include/numerics/numeric_vector.h
===================================================================
--- include/numerics/numeric_vector.h (Revision 3231)
+++ include/numerics/numeric_vector.h (Arbeitskopie)
@@ -161,7 +161,21 @@
const unsigned int /*n_local*/,
const std::vector<unsigned int>& /*ghost*/,
const bool /*fast*/ = false) {}
+
+ /**
+ * Creates a vector that has the same dimension and storage type as
+ * \p other, including ghost dofs.
+ */
+ virtual void init (const NumericVector<T>& other,
+ const bool fast = false);
+ /**
+ * Will return \p true if the vector contains ghost dofs. Default
+ * implementation is for solver packages for which ghost dofs are
+ * not yet implemented in libMesh.
+ */
+ virtual bool is_ghosted (void) const { return false; }
+
// /**
// * Change the dimension to that of the
// * vector \p V. The same applies as for
Index: include/base/dof_map.h
===================================================================
--- include/base/dof_map.h (Revision 3231)
+++ include/base/dof_map.h (Arbeitskopie)
@@ -582,7 +582,7 @@
* when projecting the solution onto the refined mesh, hence the parent's
* DOF indices need to be included in the \p _send_list.
*/
- void augment_send_list_for_projection(const MeshBase &);
+ // void augment_send_list_for_projection(const MeshBase &);
/**
* Fills the vector di with the global degree of freedom indices
@@ -681,12 +681,9 @@
* variable in the system.
* Starts at index next_free_dof, and increments it to
* the post-final index.
- * If build_send_list is true, builds the send list. If
- * false, clears and reserves the send list
*/
void distribute_local_dofs_var_major (unsigned int& next_free_dof,
- MeshBase& mesh,
- const bool build_send_list);
+ MeshBase& mesh);
/**
* Distributes the global degrees of freedom, for dofs on
@@ -701,8 +698,7 @@
* false, clears and reserves the send list
*/
void distribute_local_dofs_node_major (unsigned int& next_free_dof,
- MeshBase& mesh,
- const bool build_send_list);
+ MeshBase& mesh);
/**
* Adds entries to the \p _send_list vector corresponding to DoFs
Index: src/numerics/numeric_vector.C
===================================================================
--- src/numerics/numeric_vector.C (Revision 3231)
+++ src/numerics/numeric_vector.C (Arbeitskopie)
@@ -189,7 +189,19 @@
}
+
+/* Default implementation for solver packages for which ghosted
+ vectors are not yet implemented. */
template <class T>
+void NumericVector<T>::init (const NumericVector<T>& other,
+ const bool fast)
+{
+ this->init(other.size(),other.local_size(),fast);
+}
+
+
+
+template <class T>
Real NumericVector<T>::subset_l1_norm (const std::set<unsigned int> & indices)
{
NumericVector<T> & v = *this;
Index: src/numerics/petsc_vector.C
===================================================================
--- src/numerics/petsc_vector.C (Revision 3231)
+++ src/numerics/petsc_vector.C (Arbeitskopie)
@@ -438,9 +438,11 @@
PetscVector<T>&
PetscVector<T>::operator = (const PetscVector<T>& v)
{
+ if (this->initialized())
+ this->clear();
+
if (v.initialized())
{
- this->init (v.size(), v.local_size());
this->_global_to_local_map = v._global_to_local_map;
this->_is_closed = v._is_closed;
this->_is_initialized = v._is_initialized;
@@ -449,8 +451,12 @@
{
int ierr = 0;
+ ierr = VecDuplicate (v._vec, &this->_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
ierr = VecCopy (v._vec, this->_vec);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
}
@@ -514,7 +520,7 @@
PetscVector<T>* v_local = libmesh_cast_ptr<PetscVector<T>*>(&v_local_in);
libmesh_assert (v_local != NULL);
- libmesh_assert (v_local->local_size() == this->size());
+ libmesh_assert (v_local->size() == this->size());
int ierr = 0;
const int n = this->size();
@@ -573,19 +579,21 @@
PetscVector<T>* v_local = libmesh_cast_ptr<PetscVector<T>*>(&v_local_in);
libmesh_assert (v_local != NULL);
- libmesh_assert (v_local->local_size() == this->size());
+ libmesh_assert (v_local->size() == this->size());
libmesh_assert (send_list.size() <= v_local->size());
int ierr=0;
- const int n_sl = send_list.size();
+ const unsigned int n_sl = send_list.size();
IS is;
VecScatter scatter;
- std::vector<int> idx(n_sl);
+ std::vector<int> idx(n_sl + this->local_size());
- for (int i=0; i<n_sl; i++)
+ for (unsigned int i=0; i<n_sl; i++)
idx[i] = static_cast<int>(send_list[i]);
+ for (unsigned int i = 0; i != this->local_size(); ++i)
+ idx[n_sl+i] = i + this->first_local_index();
// Create the index set & scatter object
if (idx.empty())
@@ -640,7 +648,7 @@
const std::vector<unsigned int>& send_list)
{
// Only good for serial vectors.
- libmesh_assert (this->size() == this->local_size());
+ // libmesh_assert (this->size() == this->local_size());
libmesh_assert (last_local_idx > first_local_idx);
libmesh_assert (send_list.size() <= this->size());
libmesh_assert (last_local_idx < this->size());
Index: src/base/dof_map.C
===================================================================
--- src/base/dof_map.C (Revision 3231)
+++ src/base/dof_map.C (Arbeitskopie)
@@ -666,9 +666,9 @@
// Set temporary DOF indices on this processor
if (node_major_dofs)
- this->distribute_local_dofs_node_major (next_free_dof, mesh, false);
+ this->distribute_local_dofs_node_major (next_free_dof, mesh);
else
- this->distribute_local_dofs_var_major (next_free_dof, mesh, false);
+ this->distribute_local_dofs_var_major (next_free_dof, mesh);
// Get DOF counts on all processors
std::vector<unsigned int> dofs_on_proc(n_proc, 0);
@@ -692,9 +692,9 @@
// Set permanent DOF indices on this processor
if (node_major_dofs)
- this->distribute_local_dofs_node_major (next_free_dof, mesh, true);
+ this->distribute_local_dofs_node_major (next_free_dof, mesh);
else
- this->distribute_local_dofs_var_major (next_free_dof, mesh, true);
+ this->distribute_local_dofs_var_major (next_free_dof, mesh);
libmesh_assert(next_free_dof == _end_df[proc_id]);
@@ -723,6 +723,8 @@
STOP_LOG("distribute_dofs()", "DofMap");
+ _send_list.clear();
+
// Note that in the add_neighbors_to_send_list nodes on processor
// boundaries that are shared by multiple elements are added for
// each element.
@@ -734,14 +736,13 @@
void DofMap::distribute_local_dofs_node_major(unsigned int &next_free_dof,
- MeshBase& mesh,
- const bool build_send_list)
+ MeshBase& mesh)
{
const unsigned int sys_num = this->sys_number();
const unsigned int n_vars = this->n_variables();
- // Number of elements to add to send_list
- unsigned int send_list_size = 0;
+ // We now only add remote dofs to the _send_list
+ // unsigned int send_list_size = 0;
// _var_first_local_df does not work with node_major dofs
_var_first_local_df.resize(n_vars+1);
@@ -781,17 +782,6 @@
0,
next_free_dof);
next_free_dof += node->n_comp(sys_num,var);
-
- // If we've been requested to, add local dofs to the
- // _send_list
- if (build_send_list)
- for (unsigned int index=0;
index<node->n_comp(sys_num,var);
- index++)
- _send_list.push_back(node->dof_number(sys_num,
- var,
- index));
- else
- send_list_size += node->n_comp(sys_num,var);
}
}
}
@@ -810,17 +800,6 @@
next_free_dof);
next_free_dof += elem->n_comp(sys_num,var);
-
- // If we've been requested to, add local dofs to the
- // _send_list
- if (build_send_list)
- for (unsigned int index=0; index<elem->n_comp(sys_num,var);
- index++)
- _send_list.push_back(elem->dof_number(sys_num,
- var,
- index));
- else
- send_list_size += elem->n_comp(sys_num,var);
}
} // done looping over elements
@@ -845,25 +824,18 @@
}
}
#endif // DEBUG
-
- if (!build_send_list)
- {
- _send_list.clear();
- _send_list.reserve(send_list_size);
- }
}
void DofMap::distribute_local_dofs_var_major(unsigned int &next_free_dof,
- MeshBase& mesh,
- const bool build_send_list)
+ MeshBase& mesh)
{
const unsigned int sys_num = this->sys_number();
const unsigned int n_vars = this->n_variables();
- // Number of elements to add to send_list
- unsigned int send_list_size = 0;
+ // We now only add remote dofs to the _send_list
+ // unsigned int send_list_size = 0;
// We will cache the first local index for each variable
_var_first_local_df.clear();
@@ -909,17 +881,6 @@
0,
next_free_dof);
next_free_dof += node->n_comp(sys_num,var);
-
- // If we've been requested to, add local dofs to the
- // _send_list
- if (build_send_list)
- for (unsigned int index=0; index<node->n_comp(sys_num,var);
- index++)
- _send_list.push_back(node->dof_number(sys_num,
- var,
- index));
- else
- send_list_size += node->n_comp(sys_num,var);
}
}
@@ -935,17 +896,6 @@
next_free_dof);
next_free_dof += elem->n_comp(sys_num,var);
-
- // If we've been requested to, add local dofs to the
- // _send_list
- if (build_send_list)
- for (unsigned int index=0; index<elem->n_comp(sys_num,var);
- index++)
- _send_list.push_back(elem->dof_number(sys_num,
- var,
- index));
- else
- send_list_size += elem->n_comp(sys_num,var);
}
} // end loop on elements
_var_first_local_df.push_back(next_free_dof);
@@ -972,12 +922,6 @@
}
}
#endif // DEBUG
-
- if (!build_send_list)
- {
- _send_list.clear();
- _send_list.reserve(send_list_size);
- }
}
@@ -987,7 +931,6 @@
START_LOG("add_neighbors_to_send_list()", "DofMap");
//-------------------------------------------------------------------------
- // The _send_list now only contains entries from the local processor.
// We need to add the DOFs from elements that live on neighboring processors
// that are neighbors of the elements on the local processor
//-------------------------------------------------------------------------
@@ -1032,10 +975,12 @@
// Get the DOF indices for this neighboring element
this->dof_indices (family[i], di);
- // Insert the DOF indices into the send list
- _send_list.insert (_send_list.end(),
- di.begin(), di.end());
- }
+ // Insert the remote DOF indices into the send list
+ for (unsigned int j=0; j != di.size(); ++j)
+ if (di[j] < this->first_dof() ||
+ di[j] >= this->end_dof())
+ _send_list.push_back(di[j]);
+ }
}
}
@@ -1070,9 +1015,11 @@
// Get the DOF indices for this neighboring element
this->dof_indices (elem, di);
- // Insert the DOF indices into the send list
- _send_list.insert (_send_list.end(),
- di.begin(), di.end());
+ // Insert the remote DOF indices into the send list
+ for (unsigned int j=0; j != di.size(); ++j)
+ if (di[j] < this->first_dof() &&
+ di[j] >= this->end_dof())
+ _send_list.push_back(di[j]);
}
}
@@ -1559,6 +1506,7 @@
+/*
void DofMap::augment_send_list_for_projection(const MeshBase& mesh)
{
// Loop over the active local elements in the mesh.
@@ -1610,6 +1558,7 @@
// The send-list might need to be sorted again.
if (needs_sorting) this->sort_send_list ();
}
+*/
#endif // LIBMESH_ENABLE_AMR
Index: src/solvers/transient_system.C
===================================================================
--- src/solvers/transient_system.C (Revision 3231)
+++ src/solvers/transient_system.C (Arbeitskopie)
@@ -98,8 +98,10 @@
Base::init_data();
// Initialize the old & older solutions
- old_local_solution->init (this->n_dofs());
- older_local_solution->init (this->n_dofs());
+ old_local_solution->init (this->n_dofs(), this->n_local_dofs(),
+ this->get_dof_map().get_send_list());
+ older_local_solution->init (this->n_dofs(), this->n_local_dofs(),
+ this->get_dof_map().get_send_list());
}
Index: src/solvers/unsteady_solver.C
===================================================================
--- src/solvers/unsteady_solver.C (Revision 3231)
+++ src/solvers/unsteady_solver.C (Arbeitskopie)
@@ -58,7 +58,8 @@
first_solve = false;
}
- old_local_nonlinear_solution->init (_system.n_dofs());
+ old_local_nonlinear_solution->init (_system.n_dofs(), _system.n_local_dofs(),
+ _system.get_dof_map().get_send_list());
_system.get_vector("_old_nonlinear_solution").localize
(*old_local_nonlinear_solution,
Index: src/solvers/system.C
===================================================================
--- src/solvers/system.C (Revision 3231)
+++ src/solvers/system.C (Arbeitskopie)
@@ -185,7 +185,8 @@
solution->init (this->n_dofs(), this->n_local_dofs());
// Resize the current_local_solution for the current mesh
- current_local_solution->init (this->n_dofs());
+ current_local_solution->init (this->n_dofs(), this->n_local_dofs(),
+ _dof_map->get_send_list());
// from now on, no chance to add additional vectors
_can_add_vectors = false;
@@ -211,19 +212,22 @@
v->init (this->n_dofs(), this->n_local_dofs());
}
+ const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();
+
// Restrict the solution on the coarsened cells
if (_solution_projection)
{
this->project_vector (*solution);
current_local_solution->clear();
- current_local_solution->init(this->n_dofs());
- const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();
+ current_local_solution->init(this->n_dofs(),
+ this->n_local_dofs(), send_list);
solution->localize (*current_local_solution, send_list);
}
else
{
current_local_solution->clear();
- current_local_solution->init(this->n_dofs());
+ current_local_solution->init(this->n_dofs(),
+ this->n_local_dofs(), send_list);
}
#endif
}
@@ -257,7 +261,8 @@
solution->init (this->n_dofs(), this->n_local_dofs());
libmesh_assert (solution->size() == current_local_solution->size());
- libmesh_assert (solution->size() == current_local_solution->local_size());
+ // Not true with ghosted vectors
+ // libmesh_assert (solution->size() == current_local_solution->local_size());
const unsigned int first_local_dof = solution->first_local_index();
const unsigned int local_size = solution->local_size();
@@ -275,7 +280,7 @@
const std::vector<unsigned int>& send_list = _dof_map->get_send_list ();
// Check sizes
- libmesh_assert (current_local_solution->local_size() == solution->size());
+ libmesh_assert (current_local_solution->size() == solution->size());
// More processors than elements => empty send_list
// libmesh_assert (!send_list.empty());
libmesh_assert (send_list.size() <= solution->size());
@@ -298,8 +303,9 @@
Utility::iota (send_list.begin(), send_list.end(), 0);
// Check sizes
- libmesh_assert (current_local_solution->local_size() == solution->size());
libmesh_assert (current_local_solution->size() == solution->size());
+ // Not true with ghosted vectors
+ // libmesh_assert (current_local_solution->local_size() == solution->size());
libmesh_assert (!send_list.empty());
libmesh_assert (send_list.size() <= solution->size());
------------------------------------------------------------------------------
This SF.net email is sponsored by:
SourcForge Community
SourceForge wants to tell your story.
http://p.sf.net/sfu/sf-spreadtheword
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
