From a97cf027a5ea681c2e7b480476f0a94e1696e674 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 20 Oct 2023 18:59:58 +0200 Subject: [PATCH] MGTransferGlobalCoarsening: test operator with ones on operator --- .../multigrid_a_01.cc | 30 ++++--- ...id_a_01.mpirun=1.with_trilinos=true.output | 18 +++++ .../multigrid_p_01.cc | 42 +++++++--- ...id_p_01.mpirun=1.with_trilinos=true.output | 78 +++++++++++++++++++ .../multigrid_util.h | 15 +++- 5 files changed, 161 insertions(+), 22 deletions(-) diff --git a/tests/multigrid-global-coarsening/multigrid_a_01.cc b/tests/multigrid-global-coarsening/multigrid_a_01.cc index f1df401b39..d7bf93e9c4 100644 --- a/tests/multigrid-global-coarsening/multigrid_a_01.cc +++ b/tests/multigrid-global-coarsening/multigrid_a_01.cc @@ -25,7 +25,8 @@ template void test(const unsigned int n_refinements, const unsigned int fe_degree_fine, - const bool do_simplex_mesh) + const bool do_simplex_mesh, + const bool ones_on_diagonal) { using VectorType = LinearAlgebra::distributed::Vector; @@ -89,7 +90,12 @@ test(const unsigned int n_refinements, constraint.close(); // set up operator - op.reinit(*mapping, dof_handler, *quad, constraint); + op.reinit(*mapping, + dof_handler, + *quad, + constraint, + numbers::invalid_unsigned_int, + ones_on_diagonal); } // set up transfer operator @@ -161,11 +167,17 @@ main(int argc, char **argv) deallog.precision(8); - for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) - for (unsigned int degree = 1; degree <= 4; ++degree) - test<2>(n_refinements, degree, false /*quadrilateral*/); - - for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) - for (unsigned int degree = 1; degree <= 2; ++degree) - test<2>(n_refinements, degree, true /*triangle*/); + for (const auto ones_on_diagonal : {false, true}) + { + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 1; degree <= 4; ++degree) + test<2>(n_refinements, + degree, + false /*quadrilateral*/, + ones_on_diagonal); + + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 1; degree <= 2; ++degree) + test<2>(n_refinements, degree, true /*triangle*/, ones_on_diagonal); + } } diff --git a/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output index 07875b2dbe..5a1eb699ca 100644 --- a/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output +++ b/tests/multigrid-global-coarsening/multigrid_a_01.mpirun=1.with_trilinos=true.output @@ -17,3 +17,21 @@ DEAL:0::2 1 3 tri 3 DEAL:0::2 2 3 tri 3 DEAL:0::2 1 4 tri 3 DEAL:0::2 2 4 tri 3 +DEAL:0::2 1 2 quad 2 +DEAL:0::2 2 2 quad 3 +DEAL:0::2 3 2 quad 3 +DEAL:0::2 4 2 quad 3 +DEAL:0::2 1 3 quad 3 +DEAL:0::2 2 3 quad 3 +DEAL:0::2 3 3 quad 3 +DEAL:0::2 4 3 quad 3 +DEAL:0::2 1 4 quad 3 +DEAL:0::2 2 4 quad 3 +DEAL:0::2 3 4 quad 3 +DEAL:0::2 4 4 quad 3 +DEAL:0::2 1 2 tri 3 +DEAL:0::2 2 2 tri 3 +DEAL:0::2 1 3 tri 3 +DEAL:0::2 2 3 tri 3 +DEAL:0::2 1 4 tri 3 +DEAL:0::2 2 4 tri 3 diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.cc b/tests/multigrid-global-coarsening/multigrid_p_01.cc index bf55c9ba8e..5ca1e98eeb 100644 --- a/tests/multigrid-global-coarsening/multigrid_p_01.cc +++ b/tests/multigrid-global-coarsening/multigrid_p_01.cc @@ -26,7 +26,8 @@ void test(const unsigned int n_refinements, const unsigned int fe_degree_fine, const bool do_simplex_mesh, - const unsigned int mesh_type) + const unsigned int mesh_type, + const bool ones_on_diagonal) { using VectorType = LinearAlgebra::distributed::Vector; @@ -117,7 +118,12 @@ test(const unsigned int n_refinements, constraint.close(); // set up operator - op.reinit(*mapping, dof_handler, *quad, constraint); + op.reinit(*mapping, + dof_handler, + *quad, + constraint, + numbers::invalid_unsigned_int, + ones_on_diagonal); } // set up transfer operator @@ -193,14 +199,26 @@ main(int argc, char **argv) deallog.precision(8); - for (unsigned int mesh_type = 0; mesh_type < 2; ++mesh_type) - { - for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) - for (unsigned int degree = 2; degree <= 4; ++degree) - test<2>(n_refinements, degree, false /*quadrilateral*/, mesh_type); - - for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) - for (unsigned int degree = 2; degree <= 2; ++degree) - test<2>(n_refinements, degree, true /*triangle*/, mesh_type); - } + for (const auto ones_on_diagonal : {false, true}) + for (unsigned int mesh_type = 0; mesh_type < 2; ++mesh_type) + { + for (unsigned int n_refinements = 2; n_refinements <= 4; + ++n_refinements) + for (unsigned int degree = 2; degree <= 4; ++degree) + test<2>(n_refinements, + degree, + false /*quadrilateral*/, + mesh_type, + ones_on_diagonal); + + + for (unsigned int n_refinements = 2; n_refinements <= 4; + ++n_refinements) + for (unsigned int degree = 2; degree <= 2; ++degree) + test<2>(n_refinements, + degree, + true /*triangle*/, + mesh_type, + ones_on_diagonal); + } } diff --git a/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output index f2980eab6f..4637dd6777 100644 --- a/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output +++ b/tests/multigrid-global-coarsening/multigrid_p_01.mpirun=1.with_trilinos=true.output @@ -77,3 +77,81 @@ DEAL:0::Norm interpolated solution on level 1: 0.41924192 DEAL:0::2 2 4 tri 4 DEAL:0::Norm interpolated solution on level 0: 0.40718520 DEAL:0::Norm interpolated solution on level 1: 0.80533381 +DEAL:0::2 2 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.33877875 +DEAL:0::Norm interpolated solution on level 1: 0.66014529 +DEAL:0::2 3 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.33872299 +DEAL:0::Norm interpolated solution on level 1: 0.99016130 +DEAL:0::2 4 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.33877890 +DEAL:0::Norm interpolated solution on level 1: 0.66010393 +DEAL:0::Norm interpolated solution on level 2: 1.3201957 +DEAL:0::2 2 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.66459543 +DEAL:0::Norm interpolated solution on level 1: 1.3203627 +DEAL:0::2 3 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.66458798 +DEAL:0::Norm interpolated solution on level 1: 1.9805368 +DEAL:0::2 4 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.66459725 +DEAL:0::Norm interpolated solution on level 1: 1.3203574 +DEAL:0::Norm interpolated solution on level 2: 2.6407133 +DEAL:0::2 2 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225826 +DEAL:0::Norm interpolated solution on level 1: 2.6407347 +DEAL:0::2 3 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225816 +DEAL:0::Norm interpolated solution on level 1: 3.9611011 +DEAL:0::2 4 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225830 +DEAL:0::Norm interpolated solution on level 1: 2.6407340 +DEAL:0::Norm interpolated solution on level 2: 5.2814679 +DEAL:0::2 2 2 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.33806168 +DEAL:0::Norm interpolated solution on level 1: 0.66011744 +DEAL:0::2 2 3 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.66445481 +DEAL:0::Norm interpolated solution on level 1: 1.3203595 +DEAL:0::2 2 4 tri 3 +DEAL:0::Norm interpolated solution on level 0: 1.3225585 +DEAL:0::Norm interpolated solution on level 1: 2.6407344 +DEAL:0::2 2 2 quad 2 +DEAL:0::Norm interpolated solution on level 0: 0.13157424 +DEAL:0::Norm interpolated solution on level 1: 0.24449632 +DEAL:0::2 3 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.13144458 +DEAL:0::Norm interpolated solution on level 1: 0.35587970 +DEAL:0::2 4 2 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.13139421 +DEAL:0::Norm interpolated solution on level 1: 0.24404614 +DEAL:0::Norm interpolated solution on level 2: 0.46558561 +DEAL:0::2 2 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.23622379 +DEAL:0::Norm interpolated solution on level 1: 0.45479680 +DEAL:0::2 3 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.23657515 +DEAL:0::Norm interpolated solution on level 1: 0.66773399 +DEAL:0::2 4 3 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.23652211 +DEAL:0::Norm interpolated solution on level 1: 0.45456995 +DEAL:0::Norm interpolated solution on level 2: 0.87733284 +DEAL:0::2 2 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.43119420 +DEAL:0::Norm interpolated solution on level 1: 0.85199621 +DEAL:0::2 3 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.43119824 +DEAL:0::Norm interpolated solution on level 1: 1.2615628 +DEAL:0::2 4 4 quad 3 +DEAL:0::Norm interpolated solution on level 0: 0.43118187 +DEAL:0::Norm interpolated solution on level 1: 0.85199072 +DEAL:0::Norm interpolated solution on level 2: 1.6669458 +DEAL:0::2 2 2 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.13575160 +DEAL:0::Norm interpolated solution on level 1: 0.24280231 +DEAL:0::2 2 3 tri 3 +DEAL:0::Norm interpolated solution on level 0: 0.21915908 +DEAL:0::Norm interpolated solution on level 1: 0.41924192 +DEAL:0::2 2 4 tri 4 +DEAL:0::Norm interpolated solution on level 0: 0.40718520 +DEAL:0::Norm interpolated solution on level 1: 0.80533381 diff --git a/tests/multigrid-global-coarsening/multigrid_util.h b/tests/multigrid-global-coarsening/multigrid_util.h index af2b57b6d8..a0a65e1db2 100644 --- a/tests/multigrid-global-coarsening/multigrid_util.h +++ b/tests/multigrid-global-coarsening/multigrid_util.h @@ -75,7 +75,8 @@ public: const DoFHandler &dof_handler, const Quadrature &quad, const AffineConstraints &constraints, - const unsigned int mg_level = numbers::invalid_unsigned_int) + const unsigned int mg_level = numbers::invalid_unsigned_int, + const bool ones_on_diagonal = false) { // Clear internal data structures (if operator is reused). this->system_matrix.clear(); @@ -92,6 +93,12 @@ public: data.mg_level = mg_level; matrix_free.reinit(mapping, dof_handler, constraints, quad, data); + + constrained_indices.clear(); + + if (ones_on_diagonal) + for (auto i : this->matrix_free.get_constrained_dofs()) + constrained_indices.push_back(i); } virtual types::global_dof_index @@ -122,6 +129,10 @@ public: { this->matrix_free.cell_loop( &Operator::do_cell_integral_range, this, dst, src, true); + + for (unsigned int i = 0; i < constrained_indices.size(); ++i) + dst.local_element(constrained_indices[i]) = + src.local_element(constrained_indices[i]); } void @@ -269,6 +280,8 @@ private: AffineConstraints constraints; mutable TrilinosWrappers::SparseMatrix system_matrix; + + std::vector constrained_indices; }; struct GMGParameters -- 2.39.5