From: bangerth Date: Thu, 23 May 2013 19:19:51 +0000 (+0000) Subject: Implement ConstraintMatrix::distribute() in a different way. Adjust testcase. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a8bda5adcadd0f23700b42ddc3f3201bdf88acf7;p=dealii-svn.git Implement ConstraintMatrix::distribute() in a different way. Adjust testcase. git-svn-id: https://svn.dealii.org/trunk@29558 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/lac/constraint_matrix.cc b/deal.II/source/lac/constraint_matrix.cc index 306691df8b..0addcc6ad2 100644 --- a/deal.II/source/lac/constraint_matrix.cc +++ b/deal.II/source/lac/constraint_matrix.cc @@ -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::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 - 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::const_iterator constraint_iterator; - std::set 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; ientries.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; ientries.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; ientries.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); } diff --git a/tests/mpi/petsc_distribute_02.cc b/tests/mpi/petsc_distribute_02.cc index 2ec7f95ff3..6656008f9c 100644 --- a/tests/mpi/petsc_distribute_02.cc +++ b/tests/mpi/petsc_distribute_02.cc @@ -26,11 +26,10 @@ // // 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 1) + Assert (exception_caught == true, ExcInternalError()); - { - double exact_l1 = 0; - - // add up original values of vector entries - for (unsigned int i=0; i