-
Notifications
You must be signed in to change notification settings - Fork 707
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Make get_cell_range_category() more consistent with get_cell_category() #14015
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,246 @@ | ||
// --------------------------------------------------------------------- | ||
// | ||
// Copyright (C) 2022 by the deal.II authors | ||
// | ||
// This file is part of the deal.II library. | ||
// | ||
// The deal.II library is free software; you can use it, redistribute | ||
// it, and/or modify it under the terms of the GNU Lesser General | ||
// Public License as published by the Free Software Foundation; either | ||
// version 2.1 of the License, or (at your option) any later version. | ||
// The full text of the license can be found in the file LICENSE.md at | ||
// the top level directory of deal.II. | ||
// | ||
// --------------------------------------------------------------------- | ||
|
||
|
||
|
||
// tests MatrixFree::get_cell_category() and | ||
// MatrixFree::get_cell_range_category() | ||
|
||
#include <deal.II/distributed/shared_tria.h> | ||
|
||
#include <deal.II/dofs/dof_handler.h> | ||
|
||
#include <deal.II/fe/fe_nothing.h> | ||
#include <deal.II/fe/fe_q.h> | ||
|
||
#include <deal.II/grid/filtered_iterator.h> | ||
#include <deal.II/grid/grid_generator.h> | ||
|
||
#include <deal.II/matrix_free/fe_evaluation.h> | ||
#include <deal.II/matrix_free/matrix_free.h> | ||
|
||
#include "../tests.h" | ||
|
||
template <int dim> | ||
void | ||
test() | ||
{ | ||
using Number = double; | ||
using VectorizedArrayType = VectorizedArray<Number>; | ||
|
||
Triangulation<dim> tria; | ||
|
||
GridGenerator::hyper_cube(tria); | ||
tria.refine_global(3); | ||
|
||
// caterorization - not strict | ||
{ | ||
DoFHandler<dim> dof_handler(tria); | ||
dof_handler.distribute_dofs(FE_Q<dim>(1)); | ||
|
||
typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data; | ||
data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none; | ||
data.mapping_update_flags = update_quadrature_points; | ||
|
||
std::vector<unsigned int> cell_vectorization_category( | ||
tria.n_active_cells()); | ||
for (unsigned int i = 0; i < cell_vectorization_category.size(); ++i) | ||
cell_vectorization_category[i] = i / 7; | ||
|
||
data.cell_vectorization_category = cell_vectorization_category; | ||
data.cell_vectorization_categories_strict = false; | ||
|
||
MatrixFree<dim, Number, VectorizedArrayType> matrix_free; | ||
matrix_free.reinit(MappingQ1<dim>(), | ||
dof_handler, | ||
AffineConstraints<Number>(), | ||
QGauss<dim>(2), | ||
data); | ||
|
||
double dummy = 0; | ||
|
||
matrix_free.template cell_loop<double, double>( | ||
[&](const auto &matrix_free, auto &, const auto &, const auto range) { | ||
unsigned int max_category_range = 0; | ||
|
||
for (unsigned int cell = range.first; cell < range.second; ++cell) | ||
{ | ||
unsigned int max_category = 0; | ||
for (unsigned int v = 0; | ||
v < matrix_free.n_active_entries_per_cell_batch(cell); | ||
++v) | ||
{ | ||
const auto category = cell_vectorization_category | ||
[matrix_free.get_cell_iterator(cell, v)->active_cell_index()]; | ||
|
||
// note: here and below we are using std::cout to output | ||
// information that is useful for debugging; however, since the | ||
// output depends on the hardware and might change if we | ||
// internally modify the way we loop over cells, we don't write | ||
// the data to the log but instead only check that the | ||
// properties that are described in the documentation are | ||
// fulfilled | ||
Comment on lines
+88
to
+94
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @kronbichler What do you think? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Excellent. |
||
std::cout << category << " "; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Is this output on purpose? It won't be compared, and one will only see it when manually running the test (which might be ok). There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I would keep it this way. I introduced these for visualization purposes. It looks like this:
In the test, I am only checking the properties, since outputting such an output would require that we create an output for each vectorization type. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I agree that we should not print this with vectorization. Maybe add a short comment in the code that explains why we have the |
||
max_category = std::max(max_category, category); | ||
} | ||
|
||
std::cout << std::endl; | ||
|
||
AssertDimension(max_category, | ||
matrix_free.get_cell_category(cell)); // test | ||
|
||
max_category_range = std::max(max_category_range, max_category); | ||
} | ||
|
||
AssertDimension(max_category_range, | ||
matrix_free.get_cell_range_category(range)); // test | ||
std::cout << std::endl; | ||
}, | ||
dummy, | ||
dummy); | ||
} | ||
deallog << "OK" << std::endl; | ||
|
||
// caterorization - strict | ||
{ | ||
DoFHandler<dim> dof_handler(tria); | ||
dof_handler.distribute_dofs(FE_Q<dim>(1)); | ||
|
||
typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data; | ||
data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none; | ||
data.mapping_update_flags = update_quadrature_points; | ||
data.mapping_update_flags_boundary_faces = update_quadrature_points; | ||
data.mapping_update_flags_inner_faces = update_quadrature_points; | ||
data.hold_all_faces_to_owned_cells = true; | ||
|
||
std::vector<unsigned int> cell_vectorization_category( | ||
tria.n_active_cells()); | ||
for (unsigned int i = 0; i < cell_vectorization_category.size(); ++i) | ||
cell_vectorization_category[i] = i / 7; | ||
|
||
data.cell_vectorization_category = cell_vectorization_category; | ||
data.cell_vectorization_categories_strict = true; | ||
|
||
MatrixFree<dim, Number, VectorizedArrayType> matrix_free; | ||
matrix_free.reinit(MappingQ1<dim>(), | ||
dof_handler, | ||
AffineConstraints<Number>(), | ||
QGauss<dim>(2), | ||
data); | ||
|
||
double dummy = 0; | ||
|
||
matrix_free.template cell_loop<double, double>( | ||
[&](const auto &matrix_free, auto &, const auto &, const auto range) { | ||
unsigned int max_category_range = 0; | ||
|
||
for (unsigned int cell = range.first; cell < range.second; ++cell) | ||
{ | ||
for (unsigned int v = 0; | ||
v < matrix_free.n_active_entries_per_cell_batch(cell); | ||
++v) | ||
{ | ||
const auto category = cell_vectorization_category | ||
[matrix_free.get_cell_iterator(cell, v)->active_cell_index()]; | ||
std::cout << category << " "; | ||
AssertDimension(category, | ||
matrix_free.get_cell_category(cell)); // test | ||
} | ||
|
||
std::cout << std::endl; | ||
|
||
max_category_range = | ||
std::max(max_category_range, matrix_free.get_cell_category(cell)); | ||
} | ||
|
||
AssertDimension(max_category_range, | ||
matrix_free.get_cell_range_category(range)); // test | ||
std::cout << std::endl; | ||
}, | ||
dummy, | ||
dummy); | ||
} | ||
deallog << "OK" << std::endl; | ||
|
||
|
||
// hp | ||
{ | ||
DoFHandler<dim> dof_handler(tria); | ||
|
||
unsigned int counter = 0; | ||
for (const auto &cell : dof_handler.active_cell_iterators()) | ||
if (cell->is_locally_owned()) | ||
cell->set_active_fe_index((counter++) / 7); | ||
|
||
hp::FECollection<dim> fe_collection; | ||
|
||
for (unsigned int i = 0; i < (tria.n_cells() + 6) / 7; ++i) | ||
fe_collection.push_back(FE_Q<dim>(1)); | ||
|
||
dof_handler.distribute_dofs(fe_collection); | ||
|
||
typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data; | ||
data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none; | ||
data.mapping_update_flags = update_quadrature_points; | ||
data.mapping_update_flags_boundary_faces = update_quadrature_points; | ||
data.mapping_update_flags_inner_faces = update_quadrature_points; | ||
data.hold_all_faces_to_owned_cells = true; | ||
|
||
MatrixFree<dim, Number, VectorizedArrayType> matrix_free; | ||
matrix_free.reinit(MappingQ1<dim>(), | ||
dof_handler, | ||
AffineConstraints<Number>(), | ||
QGauss<dim>(2), | ||
data); | ||
|
||
|
||
double dummy = 0; | ||
|
||
matrix_free.template cell_loop<double, double>( | ||
[&](const auto &matrix_free, auto &, const auto &, const auto range) { | ||
unsigned int max_category_range = 0; | ||
|
||
for (unsigned int cell = range.first; cell < range.second; ++cell) | ||
{ | ||
for (unsigned int v = 0; | ||
v < matrix_free.n_active_entries_per_cell_batch(cell); | ||
++v) | ||
{ | ||
const auto category = | ||
matrix_free.get_cell_iterator(cell, v)->active_fe_index(); | ||
std::cout << category << " "; | ||
AssertDimension(category, | ||
matrix_free.get_cell_category(cell)); // test | ||
AssertDimension( | ||
category, matrix_free.get_cell_range_category(range)); // test | ||
} | ||
|
||
std::cout << std::endl; | ||
} | ||
std::cout << std::endl; | ||
}, | ||
dummy, | ||
dummy); | ||
} | ||
deallog << "OK" << std::endl; | ||
} | ||
|
||
int | ||
main(int argc, char **argv) | ||
{ | ||
initlog(); | ||
|
||
test<2>(); | ||
} |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,4 @@ | ||
|
||
DEAL::OK | ||
DEAL::OK | ||
DEAL::OK |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@kronbichler There have been other places with references to
non-hp-DoFHandler
;)