Skip to content

Commit

Permalink
Use unghosted vectors for some DiagonalMatrices.
Browse files Browse the repository at this point in the history
  • Loading branch information
drwells committed Aug 17, 2021
1 parent 8723405 commit 44909f6
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 21 deletions.
4 changes: 4 additions & 0 deletions doc/news/changes/incompatibilities/20210817DavidWells
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
Changed: The DiagonalMatrix objects set up by MatrixFreeOperators::Base and
it sinheriting classes no longer contain ghost entries.
<br>
(David Wells, 2021/08/17)
61 changes: 40 additions & 21 deletions include/deal.II/matrix_free/operators.h
Original file line number Diff line number Diff line change
Expand Up @@ -1864,8 +1864,8 @@ namespace MatrixFreeOperators
VectorType &inverse_diagonal_vector =
this->inverse_diagonal_entries->get_vector();
VectorType &diagonal_vector = this->diagonal_entries->get_vector();
this->initialize_dof_vector(inverse_diagonal_vector);
this->initialize_dof_vector(diagonal_vector);
this->initialize_dof_vector(inverse_diagonal_vector, false);
this->initialize_dof_vector(diagonal_vector, false);

// Set up the action of the mass matrix in a way that's compatible with
// MatrixFreeTools::compute_diagonal:
Expand All @@ -1886,18 +1886,24 @@ namespace MatrixFreeOperators
VectorizedArrayType> &)>
diagonal_evaluation_f(diagonal_evaluation);

// we need to have properly ghosted vectors for compute_diagonal()
VectorType temp;
this->initialize_dof_vector(temp);
Assert(this->selected_rows.size() > 0, ExcInternalError());
for (unsigned int block_n = 0; block_n < this->selected_rows.size();
++block_n)
MatrixFreeTools::compute_diagonal(*this->data,
BlockHelper::subblock(diagonal_vector,
block_n),
BlockHelper::subblock(temp, block_n),
diagonal_evaluation_f,
this->selected_rows[block_n]);

// Constrained entries will create zeros on the main diagonal, which we
// don't want
this->set_constrained_entries_to_one(diagonal_vector);
this->set_constrained_entries_to_one(temp);
diagonal_vector = temp;
// We want to have unghosted vectors in the diagonal matrices - verify that
// the MatrixFree object (or anything else) didn't change the partitioner
Assert(diagonal_vector.get_partitioner()->n_ghost_indices() == 0,
ExcInternalError());

inverse_diagonal_vector = diagonal_vector;

Expand All @@ -1910,8 +1916,8 @@ namespace MatrixFreeOperators
1. / inverse_diagonal_vector.local_element(i);
}

inverse_diagonal_vector.update_ghost_values();
diagonal_vector.update_ghost_values();
diagonal_vector.set_ghost_state(true);
inverse_diagonal_vector.set_ghost_state(true);
}


Expand Down Expand Up @@ -1944,13 +1950,22 @@ namespace MatrixFreeOperators
VectorType &inverse_lumped_diagonal_vector =
inverse_lumped_diagonal_entries->get_vector();
VectorType &lumped_diagonal_vector = lumped_diagonal_entries->get_vector();
this->initialize_dof_vector(inverse_lumped_diagonal_vector);
this->initialize_dof_vector(lumped_diagonal_vector);

// Re-use the inverse_lumped_diagonal_vector as the vector of 1s
inverse_lumped_diagonal_vector = Number(1.);
apply_add(lumped_diagonal_vector, inverse_lumped_diagonal_vector);
this->set_constrained_entries_to_one(lumped_diagonal_vector);
this->initialize_dof_vector(inverse_lumped_diagonal_vector, false);
this->initialize_dof_vector(lumped_diagonal_vector, false);

// we need properly ghosted vectors for apply_add()
VectorType ones;
this->initialize_dof_vector(ones);
ones = Number(1.);
VectorType temp;
this->initialize_dof_vector(temp);
apply_add(temp, ones);
this->set_constrained_entries_to_one(temp);
lumped_diagonal_vector = temp;
// We want to have unghosted vectors in the diagonal matrices - verify that
// the MatrixFree object (or anything else) didn't change the partitioner
Assert(lumped_diagonal_vector.get_partitioner()->n_ghost_indices() == 0,
ExcInternalError());

const size_type locally_owned_size =
inverse_lumped_diagonal_vector.locally_owned_size();
Expand All @@ -1967,8 +1982,8 @@ namespace MatrixFreeOperators
Number(1.) / lumped_diagonal_vector.local_element(i);
}

inverse_lumped_diagonal_vector.update_ghost_values();
lumped_diagonal_vector.update_ghost_values();
lumped_diagonal_vector.set_ghost_state(true);
inverse_lumped_diagonal_vector.set_ghost_state(true);
}


Expand Down Expand Up @@ -2187,14 +2202,18 @@ namespace MatrixFreeOperators
VectorType &inverse_diagonal_vector =
this->inverse_diagonal_entries->get_vector();
VectorType &diagonal_vector = this->diagonal_entries->get_vector();
this->initialize_dof_vector(inverse_diagonal_vector);
this->initialize_dof_vector(diagonal_vector);
this->initialize_dof_vector(inverse_diagonal_vector, false);
this->initialize_dof_vector(diagonal_vector, false);

this->data->cell_loop(&LaplaceOperator::local_diagonal_cell,
this,
diagonal_vector,
/*unused*/ diagonal_vector);
this->set_constrained_entries_to_one(diagonal_vector);
// We want to have unghosted vectors in the diagonal matrices - verify that
// the MatrixFree object (or anything else) didn't change the partitioner
Assert(diagonal_vector.get_partitioner()->n_ghost_indices() == 0,
ExcInternalError());

inverse_diagonal_vector = diagonal_vector;

Expand All @@ -2207,8 +2226,8 @@ namespace MatrixFreeOperators
else
inverse_diagonal_vector.local_element(i) = 1.;

inverse_diagonal_vector.update_ghost_values();
diagonal_vector.update_ghost_values();
diagonal_vector.set_ghost_state(true);
inverse_diagonal_vector.set_ghost_state(true);
}


Expand Down

0 comments on commit 44909f6

Please sign in to comment.