From b540e1543d035ca4da784ad86e96ff64d24ce2f0 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 20 Oct 2015 17:01:03 -0400 Subject: [PATCH] update in/out GMG matrices - replace logic for in/out matrices for non-meshworker step-16/step-50 - avoid usage of get_refinement_edge_boundary_indices - test with/without meshworker in step-16-bdry1 --- examples/step-50/step-50.cc | 69 ++++++++++++---------------- tests/mpi/multigrid_adaptive.cc | 38 +++++++++++---- tests/multigrid/step-16-bdry1.cc | 61 +++++++++++++++++------- tests/multigrid/step-16-bdry1.output | 40 ++++++++++++++++ tests/multigrid/step-16.cc | 45 ++++++++++++------ 5 files changed, 175 insertions(+), 78 deletions(-) diff --git a/examples/step-50/step-50.cc b/examples/step-50/step-50.cc index 05964f5c92..46c1905c2d 100644 --- a/examples/step-50/step-50.cc +++ b/examples/step-50/step-50.cc @@ -103,8 +103,6 @@ #include #include #include -#include -#include // The last step is as in all // previous programs: @@ -113,11 +111,6 @@ namespace Step50 using namespace dealii; - - - - - // @sect3{The LaplaceProblem class template} // This main class is basically the same @@ -240,7 +233,7 @@ namespace Step50 const unsigned int) const { if (p.square() < 0.5*0.5) - return 20; + return 5; else return 1; } @@ -306,13 +299,6 @@ namespace Step50 deallog.depth_console(0); } - std::string id_to_string(const CellId &id) - { - std::ostringstream ss; - ss << id; - return ss.str(); - } - // @sect4{LaplaceProblem::setup_system} @@ -613,7 +599,7 @@ namespace Step50 // freedom at once can be done using // ConstraintMatrix::add_lines(): std::vector boundary_constraints (triangulation.n_global_levels()); - std::vector boundary_interface_constraints (triangulation.n_global_levels()); + ConstraintMatrix empty_constraints; for (unsigned int level=0; levellevel()); - + const unsigned int lvl = cell->level(); for (unsigned int i=0; ilevel(),local_dof_indices[i]) - || mg_constrained_dofs.at_refinement_edge(cell->level(),local_dof_indices[j])) */ - if ( !(interface_dofs_on_level.is_element(local_dof_indices[i])==true && - interface_dofs_on_level.is_element(local_dof_indices[j])==false)) - cell_matrix(i,j) = 0; - - boundary_interface_constraints[cell->level()] + if (interface_dofs_on_level.is_element(local_dof_indices[i]) // at_refinement_edge(i) + && + !interface_dofs_on_level.is_element(local_dof_indices[j]) // !at_refinement_edge(j) + && + ( + (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j]) + ) // ( !boundary(i) && !boundary(j) ) + || + ( + mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + local_dof_indices[i]==local_dof_indices[j] + ) // ( boundary(i) && boundary(j) && i==j ) + ) + ) + { + // do nothing, so add entries to interface matrix + } + else + { + cell_matrix(i,j) = 0; + } + + + empty_constraints .distribute_local_to_global (cell_matrix, local_dof_indices, mg_interface_matrices[cell->level()]); @@ -875,8 +872,6 @@ namespace Step50 SolverControl solver_control (500, 1e-8*system_rhs.l2_norm(), false); SolverGMRES cg (solver_control); - //solution = 0; - if (false) { // code to optionally compare to Trilinos ML @@ -927,7 +922,6 @@ namespace Step50 template void LaplaceProblem::refine_grid () { - Vector estimated_error_per_cell (triangulation.n_active_cells()); TrilinosWrappers::MPI::Vector temp_solution; @@ -945,8 +939,6 @@ namespace Step50 estimated_error_per_cell, 0.3, 0.0); - - triangulation.prepare_coarsening_and_refinement (); triangulation.execute_coarsening_and_refinement (); } @@ -1056,7 +1048,6 @@ namespace Step50 ? ")" : ", "); deallog << std::endl; - assemble_system (); assemble_multigrid (); diff --git a/tests/mpi/multigrid_adaptive.cc b/tests/mpi/multigrid_adaptive.cc index 0168acfc71..63dd185623 100644 --- a/tests/mpi/multigrid_adaptive.cc +++ b/tests/mpi/multigrid_adaptive.cc @@ -318,16 +318,12 @@ namespace Step50 = mg_constrained_dofs.get_refinement_edge_boundary_indices (); std::vector boundary_constraints (triangulation.n_levels()); - std::vector boundary_interface_constraints (triangulation.n_levels()); + ConstraintMatrix empty_constraints; for (unsigned int level=0; level::cell_iterator cell = mg_dof_handler.begin(), @@ -356,14 +352,36 @@ namespace Step50 .distribute_local_to_global (cell_matrix, local_dof_indices, mg_matrices[cell->level()]); + const unsigned int lvl = cell->level(); for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && - interface_dofs[cell->level()][local_dof_indices[j]]==false)) - cell_matrix(i,j) = 0; - - boundary_interface_constraints[cell->level()] + if (interface_dofs[lvl][local_dof_indices[i]] // at_refinement_edge(i) + && + !interface_dofs[lvl][local_dof_indices[j]] // !at_refinement_edge(j) + && + ( + (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j]) + ) // ( !boundary(i) && !boundary(j) ) + || + ( + mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + local_dof_indices[i]==local_dof_indices[j] + ) // ( boundary(i) && boundary(j) && i==j ) + ) + ) + { + // do nothing, so add entries to interface matrix + } + else + { + cell_matrix(i,j) = 0; + } + + empty_constraints .distribute_local_to_global (cell_matrix, local_dof_indices, mg_interface_matrices[cell->level()]); diff --git a/tests/multigrid/step-16-bdry1.cc b/tests/multigrid/step-16-bdry1.cc index 1306bd5756..26c3e20c2b 100644 --- a/tests/multigrid/step-16-bdry1.cc +++ b/tests/multigrid/step-16-bdry1.cc @@ -129,8 +129,8 @@ public: private: void setup_system (); void assemble_system (); - void assemble_multigrid (const bool &use_mw); - void solve (); + void assemble_multigrid (bool use_mw); + void solve (bool use_mw); void refine_grid (const std::string &reftype); void output_results (const unsigned int cycle) const; @@ -341,11 +341,15 @@ void LaplaceProblem::assemble_system () template -void LaplaceProblem::assemble_multigrid (const bool &use_mw) +void LaplaceProblem::assemble_multigrid (bool use_mw) { + deallog << "assemble_multigrid " << (use_mw ? "(mesh_worker)":"") << std::endl; + if (use_mw == true) { mg_matrices = 0.; + mg_interface_in = 0.; + mg_interface_out = 0.; MappingQ1 mapping; MeshWorker::IntegrationInfoBox info_box; @@ -392,20 +396,15 @@ void LaplaceProblem::assemble_multigrid (const bool &use_mw) std::vector > interface_dofs = mg_constrained_dofs.get_refinement_edge_indices (); - std::vector > boundary_interface_dofs - = mg_constrained_dofs.get_refinement_edge_boundary_indices (); std::vector boundary_constraints (triangulation.n_levels()); - std::vector boundary_interface_constraints (triangulation.n_levels()); + ConstraintMatrix empty_constraints; + for (unsigned int level=0; level::level_cell_iterator cell = mg_dof_handler.begin_mg(), @@ -434,13 +433,37 @@ void LaplaceProblem::assemble_multigrid (const bool &use_mw) local_dof_indices, mg_matrices[cell->level()]); + const unsigned int lvl = cell->level(); + for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && - interface_dofs[cell->level()][local_dof_indices[j]]==false)) - cell_matrix(i,j) = 0; + if (interface_dofs[lvl][local_dof_indices[i]] // at_refinement_edge(i) + && + !interface_dofs[lvl][local_dof_indices[j]] // !at_refinement_edge(j) + && + ( + (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j]) + ) // ( !boundary(i) && !boundary(j) ) + || + ( + mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + local_dof_indices[i]==local_dof_indices[j] + ) // ( boundary(i) && boundary(j) && i==j ) + ) + ) + { + // do nothing, so add entries to interface matrix + } + else + { + cell_matrix(i,j) = 0; + } + - boundary_interface_constraints[cell->level()] + empty_constraints .distribute_local_to_global (cell_matrix, local_dof_indices, mg_interface_in[cell->level()]); @@ -450,7 +473,7 @@ void LaplaceProblem::assemble_multigrid (const bool &use_mw) template -void LaplaceProblem::solve () +void LaplaceProblem::solve (bool use_mw) { MGTransferPrebuilt > mg_transfer(hanging_node_constraints, mg_constrained_dofs); mg_transfer.build_matrices(mg_dof_handler); @@ -470,6 +493,8 @@ void LaplaceProblem::solve () mg::Matrix<> mg_matrix(mg_matrices); mg::Matrix<> mg_interface_up(mg_interface_in); mg::Matrix<> mg_interface_down(mg_interface_in); + //if (use_mw) + //mg_interface_down.initialize(mg_interface_out); Multigrid > mg(mg_dof_handler, mg_matrix, @@ -588,9 +613,13 @@ void LaplaceProblem::run () deallog << std::endl; assemble_system (); + assemble_multigrid (false); + + solve (false); + assemble_multigrid (true); + solve (true); - solve (); // output_results (cycle); } } diff --git a/tests/multigrid/step-16-bdry1.output b/tests/multigrid/step-16-bdry1.output index 17d2f3bb22..9e63e59aab 100644 --- a/tests/multigrid/step-16-bdry1.output +++ b/tests/multigrid/step-16-bdry1.output @@ -3,6 +3,11 @@ DEAL::Cycle 0: DEAL:: Number of active cells: 4 DEAL::Number of degrees of freedom: 9 L0: 4 L1: 9 DEAL:: Number of degrees of freedom: 9 (by level: 4, 9) +DEAL::assemble_multigrid +DEAL:cg::Starting value 1.36931 +DEAL:cg::Convergence step 4 value 8.55703e-14 +DEAL:: 4 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 1.36931 DEAL:cg::Convergence step 4 value 8.55703e-14 DEAL:: 4 CG iterations needed to obtain convergence. @@ -10,6 +15,11 @@ DEAL::Cycle 1: DEAL:: Number of active cells: 16 DEAL::Number of degrees of freedom: 25 L0: 4 L1: 9 L2: 25 DEAL:: Number of degrees of freedom: 25 (by level: 4, 9, 25) +DEAL::assemble_multigrid +DEAL:cg::Starting value 0.843171 +DEAL:cg::Convergence step 6 value 9.78112e-15 +DEAL:: 6 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 0.843171 DEAL:cg::Convergence step 6 value 9.78112e-15 DEAL:: 6 CG iterations needed to obtain convergence. @@ -17,6 +27,11 @@ DEAL::Cycle 2: DEAL:: Number of active cells: 28 DEAL::Number of degrees of freedom: 41 L0: 4 L1: 9 L2: 25 L3: 25 DEAL:: Number of degrees of freedom: 41 (by level: 4, 9, 25, 25) +DEAL::assemble_multigrid +DEAL:cg::Starting value 0.737526 +DEAL:cg::Convergence step 8 value 8.81541e-15 +DEAL:: 8 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 0.737526 DEAL:cg::Convergence step 8 value 8.81541e-15 DEAL:: 8 CG iterations needed to obtain convergence. @@ -24,6 +39,11 @@ DEAL::Cycle 3: DEAL:: Number of active cells: 40 DEAL::Number of degrees of freedom: 57 L0: 4 L1: 9 L2: 25 L3: 25 L4: 25 DEAL:: Number of degrees of freedom: 57 (by level: 4, 9, 25, 25, 25) +DEAL::assemble_multigrid +DEAL:cg::Starting value 0.730417 +DEAL:cg::Convergence step 8 value 9.36194e-14 +DEAL:: 8 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 0.730417 DEAL:cg::Convergence step 8 value 9.36194e-14 DEAL:: 8 CG iterations needed to obtain convergence. @@ -31,6 +51,11 @@ DEAL::Cycle 4: DEAL:: Number of active cells: 52 DEAL::Number of degrees of freedom: 73 L0: 4 L1: 9 L2: 25 L3: 25 L4: 25 L5: 25 DEAL:: Number of degrees of freedom: 73 (by level: 4, 9, 25, 25, 25, 25) +DEAL::assemble_multigrid +DEAL:cg::Starting value 0.729970 +DEAL:cg::Convergence step 8 value 2.81558e-13 +DEAL:: 8 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 0.729970 DEAL:cg::Convergence step 8 value 2.81558e-13 DEAL:: 8 CG iterations needed to obtain convergence. @@ -38,6 +63,11 @@ DEAL::Cycle 5: DEAL:: Number of active cells: 160 DEAL::Number of degrees of freedom: 201 L0: 4 L1: 9 L2: 25 L3: 65 L4: 65 L5: 65 L6: 65 DEAL:: Number of degrees of freedom: 201 (by level: 4, 9, 25, 65, 65, 65, 65) +DEAL::assemble_multigrid +DEAL:cg::Starting value 0.492333 +DEAL:cg::Convergence step 9 value 5.07798e-13 +DEAL:: 9 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 0.492333 DEAL:cg::Convergence step 9 value 5.07798e-13 DEAL:: 9 CG iterations needed to obtain convergence. @@ -45,6 +75,11 @@ DEAL::Cycle 6: DEAL:: Number of active cells: 304 DEAL::Number of degrees of freedom: 357 L0: 4 L1: 9 L2: 25 L3: 81 L4: 81 L5: 81 L6: 81 L7: 153 DEAL:: Number of degrees of freedom: 357 (by level: 4, 9, 25, 81, 81, 81, 81, 153) +DEAL::assemble_multigrid +DEAL:cg::Starting value 0.404672 +DEAL:cg::Convergence step 8 value 5.14673e-13 +DEAL:: 8 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 0.404672 DEAL:cg::Convergence step 8 value 5.14673e-13 DEAL:: 8 CG iterations needed to obtain convergence. @@ -52,6 +87,11 @@ DEAL::Cycle 7: DEAL:: Number of active cells: 676 DEAL::Number of degrees of freedom: 761 L0: 4 L1: 9 L2: 25 L3: 81 L4: 81 L5: 81 L6: 121 L7: 209 L8: 465 DEAL:: Number of degrees of freedom: 761 (by level: 4, 9, 25, 81, 81, 81, 121, 209, 465) +DEAL::assemble_multigrid +DEAL:cg::Starting value 0.404553 +DEAL:cg::Convergence step 9 value 6.34913e-14 +DEAL:: 9 CG iterations needed to obtain convergence. +DEAL::assemble_multigrid (mesh_worker) DEAL:cg::Starting value 0.404553 DEAL:cg::Convergence step 9 value 6.34913e-14 DEAL:: 9 CG iterations needed to obtain convergence. diff --git a/tests/multigrid/step-16.cc b/tests/multigrid/step-16.cc index 2b0b81426e..ce75db645d 100644 --- a/tests/multigrid/step-16.cc +++ b/tests/multigrid/step-16.cc @@ -295,20 +295,14 @@ void LaplaceProblem::assemble_multigrid () std::vector > interface_dofs = mg_constrained_dofs.get_refinement_edge_indices (); - std::vector > boundary_interface_dofs - = mg_constrained_dofs.get_refinement_edge_boundary_indices (); std::vector boundary_constraints (triangulation.n_levels()); - std::vector boundary_interface_constraints (triangulation.n_levels()); + ConstraintMatrix empty_constraints; for (unsigned int level=0; level::cell_iterator cell = mg_dof_handler.begin(), @@ -382,13 +376,38 @@ void LaplaceProblem::assemble_multigrid () // (in the solve() // function) be able to just pass the // transpose matrix where necessary. + const unsigned int lvl = cell->level(); + for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && - interface_dofs[cell->level()][local_dof_indices[j]]==false)) - cell_matrix(i,j) = 0; - - boundary_interface_constraints[cell->level()] + if (interface_dofs[lvl][local_dof_indices[i]] // at_refinement_edge(i) + && + !interface_dofs[lvl][local_dof_indices[j]] // !at_refinement_edge(j) + && + ( + (!mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + !mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[j]) + ) // ( !boundary(i) && !boundary(j) ) + || + ( + mg_constrained_dofs.is_boundary_index(lvl, local_dof_indices[i]) + && + local_dof_indices[i]==local_dof_indices[j] + ) // ( boundary(i) && boundary(j) && i==j ) + ) + ) + { + // do nothing, so add entries to interface matrix + } + else + { + cell_matrix(i,j) = 0; + std::cout << i << " " << j << "\n"; + } + + + empty_constraints .distribute_local_to_global (cell_matrix, local_dof_indices, mg_interface_matrices[cell->level()]); @@ -584,7 +603,7 @@ void LaplaceProblem::run () assemble_multigrid (); solve (); -// output_results (cycle); + output_results (cycle); } } -- 2.39.5