On Tue, 9 Sep 2008, John Peterson wrote:
Did you send a patch to the list? I'm going back through my email but
not seeing it...
No, just to Tim (and Ben, since he'd mentioned having some time to
benchmark a hopefully similar problem). When the patch didn't show
any clear improvement, I didn't bother with sending it to the rest of
the list... but I might as well attach it now, for anyone who's
curious. It was my attempt to cause the project_vector localize()
calls to only transfer necessary data, not fill a whole serial vector.
Depending on what the patch changes, the rest of this email might be
irrelevant ;-)
No, I think the set() calls are a separate issue.
---
Roy
Index: src/solvers/system_projection.C
===================================================================
--- src/solvers/system_projection.C (revision 3010)
+++ src/solvers/system_projection.C (working copy)
@@ -70,6 +70,8 @@
*/
new_v.clear();
+#ifdef ENABLE_AMR
+
// Resize the new vector and get a serial version.
NumericVector<Number> *new_vector_ptr;
AutoPtr<NumericVector<Number> > new_vector_built;
@@ -91,6 +93,24 @@
// we need to localize.
else
{
+ // Build a send list for efficient localization
+ std::vector<unsigned int> send_list;
+ Threads::parallel_for (
+ ConstElemRange (this->get_mesh().active_local_elements_begin(),
+ this->get_mesh().active_local_elements_end()),
+ BuildProjectionList(*this, send_list));
+ // Sort the send list. After this duplicated
+ // elements will be adjacent in the vector
+ std::sort(send_list.begin(), send_list.end());
+
+ // Now use std::unique to remove duplicate entries
+ std::vector<unsigned int>::iterator new_end =
+ std::unique (send_list.begin(), send_list.end());
+
+ // Remove the end of the send_list. Use the "swap trick"
+ // from Effective STL
+ std::vector<unsigned int> (send_list.begin(), new_end).swap (send_list);
+
new_v.init (this->n_dofs(),
this->n_local_dofs());
new_vector_built = NumericVector<Number>::build();
@@ -99,7 +119,7 @@
local_old_vector = local_old_vector_built.get();
new_vector_ptr->init(this->n_dofs(), this->n_dofs());
local_old_vector->init(old_v.size(), old_v.size());
- old_v.localize(*local_old_vector);
+ old_v.localize(*local_old_vector, send_list);
local_old_vector->close();
old_vector_ptr = local_old_vector;
}
@@ -108,8 +128,6 @@
NumericVector<Number> &new_vector = *new_vector_ptr;
const NumericVector<Number> &old_vector = *old_vector_ptr;
-#ifdef ENABLE_AMR
-
Threads::parallel_for (ConstElemRange
(this->get_mesh().active_local_elements_begin(),
this->get_mesh().active_local_elements_end(),
1000),
@@ -122,6 +140,9 @@
// If the old vector was serial, we probably need to send our values
// to other processors
+ //
+ // FIXME: I'm not sure how to make a NumericVector do that without
+ // creating a temporary parallel vector to use localize! - RHS
if (old_v.size() == old_v.local_size())
{
AutoPtr<NumericVector<Number> > dist_v = NumericVector<Number>::build();
@@ -134,7 +155,7 @@
dist_v->close();
- dist_v->localize (new_v);
+ dist_v->localize (new_v, this->get_dof_map().get_send_list());
new_v.close();
}
// If the old vector was parallel, we need to update it
@@ -155,7 +176,7 @@
#else
// AMR is disabled: simply copy the vector
- new_vector = old_vector;
+ new_v = old_v;
#endif // #ifdef ENABLE_AMR
@@ -604,8 +625,65 @@
}
+#ifndef ENABLE_AMR
+void System::BuildProjectionList::operator()(const ConstElemRange &) const
+{
+ libmesh_error();
+#else
+void System::BuildProjectionList::operator()(const ConstElemRange &range) const
+{
+ // The number of variables in this system
+ const unsigned int n_variables = system.n_vars();
+ // The DofMap for this system
+ const DofMap& dof_map = system.get_dof_map();
+ // Loop over all the variables in the system
+ for (unsigned int var=0; var<n_variables; var++)
+ {
+ // The old global DOF indices
+ std::vector<unsigned int> di;
+
+ // Iterate over the elements in the range
+ for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it !=
range.end(); ++elem_it)
+ {
+ const Elem* elem = *elem_it;
+ const Elem* parent = elem->parent();
+
+ if (elem->refinement_flag() == Elem::JUST_REFINED)
+ {
+ libmesh_assert (parent != NULL);
+ unsigned int old_parent_level = parent->p_level();
+
+ if (elem->p_refinement_flag() == Elem::JUST_REFINED)
+ {
+ // We may have done p refinement or coarsening as well;
+ // if so then we need to reset the parent's p level
+ // so we can get the right DoFs from it
+ libmesh_assert(elem->p_level() > 0);
+ (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() -
1);
+ }
+ else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
+ {
+ (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() +
1);
+ }
+
+ dof_map.old_dof_indices (parent, di, var);
+ send_list.insert(send_list.end(), di.begin(), di.end());
+
+ // Fix up the parent's p level in case we changed it
+ (const_cast<Elem *>(parent))->hack_p_level(old_parent_level);
+ }
+ else
+ dof_map.old_dof_indices (elem, di, var);
+
+ send_list.insert(send_list.end(), di.begin(), di.end());
+ } // end elem loop
+ } // end variables loop
+#endif // ENABLE_AMR
+}
+
+
void System::ProjectSolution::operator()(const ConstElemRange &range) const
{
// We need data to project
Index: include/solvers/system.h
===================================================================
--- include/solvers/system.h (revision 3010)
+++ include/solvers/system.h (working copy)
@@ -737,7 +737,7 @@
/**
* This class implements projecting a vector from
* an old mesh to the newly refined mesh. This
- * may be exectued in parallel on multiple threads.
+ * may be executed in parallel on multiple threads.
*/
class ProjectVector
{
@@ -760,6 +760,29 @@
/**
+ * This class builds the send_list of old dof indices
+ * whose coefficients are needed to perform a projection.
+ * This may be executed in parallel on multiple threads.
+ */
+ class BuildProjectionList
+ {
+ private:
+ const System &system;
+ std::vector<unsigned int> &send_list;
+
+ public:
+ BuildProjectionList (const System &system_in,
+ std::vector<unsigned int> &send_list_in) :
+ system(system_in),
+ send_list(send_list_in)
+ {}
+
+ void operator()(const ConstElemRange &range) const;
+ };
+
+
+
+ /**
* This class implements projecting a vector from
* an old mesh to the newly refined mesh. This
* may be exectued in parallel on multiple threads.
-------------------------------------------------------------------------
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
