]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement ConstraintMatrix::distribute() in a different way. Adjust testcase.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 23 May 2013 19:19:51 +0000 (19:19 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 23 May 2013 19:19:51 +0000 (19:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@29558 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/lac/constraint_matrix.cc
tests/mpi/petsc_distribute_02.cc
tests/mpi/petsc_distribute_02/ncpu_1/cmp/generic
tests/mpi/petsc_distribute_02/ncpu_10/cmp/generic
tests/mpi/petsc_distribute_02/ncpu_2/cmp/generic
tests/mpi/petsc_distribute_02/ncpu_4/cmp/generic

index 306691df8bed740a5ff5bf3d4f8e9c6029dc80bb..0addcc6ad2b20096ca69aaba2a22a0b076e18aae 100644 (file)
@@ -1710,61 +1710,72 @@ ConstraintMatrix::distribute (PETScWrappers::MPI::Vector &vec) const
   Assert (vec.has_ghost_elements() == false,
           ExcMessage ("This operation can only be performed on vectors without ghost elements."));
 
-  typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
-  ConstraintLine index_comparison;
-  index_comparison.line = vec.local_range().first;
-  const constraint_iterator begin_my_constraints =
-    Utilities::lower_bound (lines.begin(),lines.end(),index_comparison);
-
-  index_comparison.line = vec.local_range().second;
-  const constraint_iterator end_my_constraints
-    = Utilities::lower_bound(lines.begin(),lines.end(),index_comparison);
-
-  // all indices we need to read from
-  IndexSet my_indices (vec.size());
-
-  const std::pair<unsigned int, unsigned int>
-  local_range = vec.local_range();
+  IndexSet vec_owned_elements (vec.size());
+  vec_owned_elements.add_range(vec.local_range().first,
+                               vec.local_range().second);
 
-  my_indices.add_range (local_range.first, local_range.second);
+  typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
 
-  std::set<unsigned int> individual_indices;
-  for (constraint_iterator it = begin_my_constraints;
-       it != end_my_constraints; ++it)
+#ifdef DEBUG
+  // the algorithm below pushes contributions c_ij x_j to a constrained
+  // DoF x_j *from the side of the owner of x_j*, as opposed to pulling
+  // it from the owner of the target side x_i. this relies on the
+  // assumption that both target and source know about the constraint
+  // on x_i.
+  //
+  // the disadvantage is that it is incredibly difficult to find out what
+  // is happening if the assumption is not satisfied. to help debug
+  // problems of this kind, we do the following in a first step if
+  // debugging is enabled:
+  // - every processor who owns an x_j where c_ij!=0 sends a one
+  //   to the owner of x_i to add up
+  // - the owner of x_i knows how many nonzero entries exist, so can
+  //   verify that the correct number of ones has been added
+  // - if the sum is not correct, then apparently one of the owners
+  //   of the x_j's did not know to send its one and, consequently,
+  //   will also not know to send the correct c_ij*x_j later on,
+  //   leading to a bug
+  set_zero (vec);
+  for (constraint_iterator it = lines.begin();
+       it != lines.end(); ++it)
     for (unsigned int i=0; i<it->entries.size(); ++i)
-      if ((it->entries[i].first < local_range.first)
-          ||
-          (it->entries[i].first >= local_range.second))
-        individual_indices.insert (it->entries[i].first);
-
-  my_indices.add_indices (individual_indices.begin(),
-                          individual_indices.end());
-
-  IndexSet local_range_is (vec.size());
-  local_range_is.add_range(local_range.first, local_range.second);
-
+      if (vec_owned_elements.is_element (it->entries[i].first))
+        vec(it->line) += 1;
+  vec.compress (VectorOperation::add);
+
+  DeclException2 (ExcIncorrectConstraint, int, int,
+                  << "While distributing the constraint for DoF "
+                  << arg1 << ", it turns out that one of the processors "
+                  << "who own the " << arg2
+                  << " degrees of freedom that x_" << arg1
+                  << " is constrained against does not know about "
+                  << "the constraint on x_" << arg1
+                  << ". Did you not initialize the ConstraintMatrix "
+                  << "with the appropriate locally_relevant set so "
+                  << "that every processor who owns a DoF that constrains "
+                  << "another DoF also knows about this constraint?");
+  for (constraint_iterator it = lines.begin();
+       it != lines.end(); ++it)
+    if (vec_owned_elements.is_element(it->line))
+      Assert (vec(it->line) == it->entries.size(),
+              ExcIncorrectConstraint(it->line, it->entries.size()));
+#endif
 
-  // create a vector and import those indices
-  PETScWrappers::MPI::Vector ghost_vec (vec.get_mpi_communicator(),
-                                        local_range_is,
-                                        my_indices);
-  ghost_vec = vec;
-  ghost_vec.update_ghost_values();
+  // set the values of all constrained DoFs to zero, then add the
+  // contributions for each constrained DoF of all the ones we
+  // own locally
+  set_zero (vec);
 
-  // finally do the distribution on own constraints
-  for (constraint_iterator it = begin_my_constraints;
-       it != end_my_constraints; ++it)
-    {
-      // fill entry in line next_constraint.line by adding the different
-      // contributions
-      PetscScalar new_value = it->inhomogeneity;
-      for (unsigned int i=0; i<it->entries.size(); ++i)
-        new_value += (PetscScalar(ghost_vec(it->entries[i].first)) *
-                      it->entries[i].second);
-      vec(it->line) = new_value;
-    }
+  for (constraint_iterator it = lines.begin();
+       it != lines.end(); ++it)
+    for (unsigned int i=0; i<it->entries.size(); ++i)
+      if (vec_owned_elements.is_element(it->entries[i].first))
+        vec(it->line) += it->entries[i].second *
+                         vec(it->entries[i].first);
 
-  vec.compress (VectorOperation::insert);
+  // now compress to communicate the entries that we added to
+  // and that weren't to local processors to the owner
+  vec.compress (VectorOperation::add);
 }
 
 
index 2ec7f95ff33cebce7017937e15206a1663726d1a..6656008f9c7f1c63ab1bf2e33764d32e885d2b35 100644 (file)
 //
 // compared to the _01 test, here the ConstraintMatrix object acts on an index
 // set that only includes the locally owned vector elements, without any
-// overlap. this verifies that we really only need to know about the *targets*
-// of constraints locally, not the *sources*. (at the time of writing this
-// test, ConstraintMatrix::distribute for PETSc distributed vectors actively
-// imported the source vector elements for those target elements we locally
-// own.)
+// overlap. this means that the owners of a source DoF in a constraint
+// do not know about the constraint if the constrained DoF is owned by
+// another processor. This is not allowed, and should trigger an
+// exception.
 
 
 #include "../tests.h"
@@ -100,52 +99,22 @@ void test()
     }
   cm.close ();
 
-  // now distribute these constraints
-  cm.distribute (vec);
-
-  // verify correctness
-  vec.compress ();
-
-  if (myid != 0)
-    Assert (vec(vec.local_range().first+10) == vec.local_range().first-25,
-           ExcInternalError());
-
-  if (myid != n_processes-1)
-    Assert (vec(vec.local_range().first+90) == vec.local_range().first+105,
-           ExcInternalError());
-
-  for (unsigned int i=vec.local_range().first; i<vec.local_range().second; ++i)
+  // now distribute these constraints. this should trigger an exception
+  bool exception_caught = false;
+  deal_II_exceptions::disable_abort_on_exception();
+  try
     {
-      if ((i != vec.local_range().first+10)
-         &&
-         (i != vec.local_range().first+90))
-       {
-         double val = vec(i);
-         Assert (std::fabs(val - i) <= 1e-6, ExcInternalError());
-       }
+      cm.distribute (vec);
     }
+  catch (...)
+    {
+      exception_caught = true;
+    }
+  if (n_processes > 1)
+    Assert (exception_caught == true, ExcInternalError());
 
-  {
-    double exact_l1 = 0;
-
-    // add up original values of vector entries
-    for (unsigned int i=0; i<vec.size(); ++i)
-      exact_l1 += i;
-
-    // but then correct for the constrained values
-    for (unsigned int p=0; p<n_processes; ++p)
-      {
-       if (p != 0)
-         exact_l1 = exact_l1 - (p*100+10) + (p*100-25);
-       if (p != n_processes-1)
-         exact_l1 = exact_l1 - (p*100+90) + (p*100+105);
-      }
-
-    Assert (vec.l1_norm() == exact_l1, ExcInternalError());
-
-    if (myid == 0)
-      deallog << "Norm = " << vec.l1_norm() << std::endl;
-  }
+  if (myid == 0)
+    deallog << "OK" << std::endl;
 }
 
 
index 2547450a428ff282205eceb66c1788461106f3eb..be8d055f86cdd3679faaa1b80a853068ea304944 100644 (file)
@@ -1,2 +1,2 @@
 
-DEAL:0::Norm = 4950.00
+DEAL:0::OK
index 464154b458dda40b804bbb12697277b01a3a320a..be8d055f86cdd3679faaa1b80a853068ea304944 100644 (file)
@@ -1,2 +1,2 @@
 
-DEAL:0::Norm = 499320.
+DEAL:0::OK
index a5749c49bf9db8553f0522aef73fc7a5f05a0c2c..be8d055f86cdd3679faaa1b80a853068ea304944 100644 (file)
@@ -1,2 +1,2 @@
 
-DEAL:0::Norm = 19880.0
+DEAL:0::OK
index 848c4e7ef53168b113bb1054bb1b5e4e73d97131..be8d055f86cdd3679faaa1b80a853068ea304944 100644 (file)
@@ -1,2 +1,2 @@
 
-DEAL:0::Norm = 79740.0
+DEAL:0::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.