On Wed, 21 Jan 2009, Roy Stogner wrote:
I've found one of the bugs,
Great!
but don't know the best way to fix it.
With the attached patch, ex9 (at least) fails because the PetscVector
operator= is capable of taking a ghosted vector and producing a
non-ghosted vector copy with a non-empty global_to_local_map, which
then fails in operator() when VecGhostGetLocalForm is called on it.
I see. My idea is, since the vector is anyway cleared and
re-initialized in operator=, why not use PETSc's VecDuplicate()
function. This seems to fix ex9, see attached patch.
(Note on the patch: I did "svn diff > patch" as usual, which means
that my patch includes your changes as well. You might find it easier
to perform my changes manually, it's just replacing the contents of
operator=. If you know some easy method to create incremental
patches, please let me know.)
Best Regards,
Tim
--
Dr. Tim Kroeger
[email protected] Phone +49-421-218-7710
[email protected] 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 3225)
+++ include/numerics/petsc_vector.h (Arbeitskopie)
@@ -654,6 +654,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]));
Index: include/base/dof_map.h
===================================================================
--- include/base/dof_map.h (Revision 3225)
+++ 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/petsc_vector.C
===================================================================
--- src/numerics/petsc_vector.C (Revision 3225)
+++ src/numerics/petsc_vector.C (Arbeitskopie)
@@ -440,18 +440,18 @@
{
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;
-
+ if (this->initialized())
+ this->clear();
if (v.size() != 0)
{
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);
- }
+ this->_global_to_local_map = v._global_to_local_map;
+ this->_is_closed = v._is_closed;
+ this->_is_initialized = v._is_initialized;
}
return *this;
@@ -514,7 +514,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 +573,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())
Index: src/base/dof_map.C
===================================================================
--- src/base/dof_map.C (Revision 3225)
+++ 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 3225)
+++ 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/system.C
===================================================================
--- src/solvers/system.C (Revision 3225)
+++ 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
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
