]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add tests
authorTimo Heister <timo.heister@gmail.com>
Tue, 21 Apr 2015 16:59:05 +0000 (12:59 -0400)
committerTimo Heister <timo.heister@gmail.com>
Tue, 21 Apr 2015 18:21:01 +0000 (14:21 -0400)
tests/mpi/interpolate_to_different_mesh_01.cc [new file with mode: 0644]
tests/mpi/interpolate_to_different_mesh_01.with_trilinos=true.mpirun=3.output [new file with mode: 0644]
tests/mpi/interpolate_to_different_mesh_02.cc [new file with mode: 0644]
tests/mpi/interpolate_to_different_mesh_02.with_trilinos=true.mpirun=3.output [new file with mode: 0644]

diff --git a/tests/mpi/interpolate_to_different_mesh_01.cc b/tests/mpi/interpolate_to_different_mesh_01.cc
new file mode 100644 (file)
index 0000000..8d82629
--- /dev/null
@@ -0,0 +1,197 @@
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <fstream>
+#include <deal.II/base/conditional_ostream.h>
+#include <cmath>
+#include <deal.II/lac/generic_linear_algebra.h>
+
+#include "../tests.h"
+
+// test VectorTools::interpolate_to_different_mesh in parallel
+
+using namespace dealii;
+
+namespace LA
+{
+  using namespace dealii::LinearAlgebraTrilinos;
+}
+
+template <int dim>
+class SomeFunction : public Function<dim>
+{
+public:
+  double value (const Point<dim> &p,
+                const unsigned int) const
+  {
+    return 1+p(0)*p(0);
+  }
+};
+
+template<int dim>
+void setup(DoFHandler<dim> & dh,
+          FE_Q<dim> & fe,
+          LA::MPI::Vector & vec,
+          LA::MPI::Vector & lr_vec)
+{
+  dh.distribute_dofs (fe);
+  vec.reinit(dh.locally_owned_dofs(), MPI_COMM_WORLD);
+  IndexSet locally_relevant;
+  DoFTools::extract_locally_relevant_dofs (dh, locally_relevant);
+  lr_vec.reinit(locally_relevant, MPI_COMM_WORLD);  
+}
+
+template<int dim>
+void output(DoFHandler<dim> & dh, LA::MPI::Vector & v, unsigned int loop, const std::string filename_)
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+  DataOut<dim> data_out;
+  data_out.add_data_vector (dh, v, "1");
+  data_out.build_patches (1);
+  std::ostringstream filename;
+  filename << filename_
+          << Utilities::int_to_string (loop, 2)
+          << "."
+          << Utilities::int_to_string (myid,2)
+          << ".vtu";
+       
+  std::ofstream output (filename.str().c_str());
+  data_out.write_vtu (output);
+  if (myid)
+    {
+      std::vector<std::string> filenames;
+      for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
+       filenames.push_back (filename_ +
+                            Utilities::int_to_string (loop, 2) +
+                            "." +
+                            Utilities::int_to_string(i, 2) +
+                            ".vtu");
+      const std::string
+       pvtu_master_filename = (filename_ +
+                               Utilities::int_to_string (loop, 2) +
+                               ".pvtu");
+      std::ofstream pvtu_master (pvtu_master_filename.c_str());
+      data_out.write_pvtu_record (pvtu_master, filenames);
+    }
+}
+
+template<int dim>
+void test()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD,
+                                              dealii::Triangulation<dim,dim>::none,
+                                              parallel::distributed::Triangulation<dim>::no_automatic_repartitioning);
+
+  GridGenerator::hyper_cube(tr);
+  tr.refine_global(2);
+  parallel::distributed::Triangulation<dim> tr2(MPI_COMM_WORLD,
+                                              dealii::Triangulation<dim,dim>::none,
+                                              parallel::distributed::Triangulation<dim>::no_automatic_repartitioning);
+
+  GridGenerator::hyper_cube(tr2);
+  tr2.refine_global(2);
+  
+  FE_Q<dim> fe(1);
+  DoFHandler<dim> dh(tr);
+  DoFHandler<dim> dh2(tr2);
+
+  SomeFunction<dim> func;
+  
+
+  for (unsigned int loop = 0; loop < 5; ++loop)
+    {
+      // randomly refine:
+      std::vector<bool> r_flags(tr.n_active_cells()*dim, false);
+      std::vector<bool> c_flags(tr.n_active_cells(), false);
+
+      for (unsigned int i=0;i<c_flags.size();++i)
+       {
+         int roll = Testing::rand()%4;
+         if (roll >= 2)
+           {
+             for (unsigned int j=0;j<dim;++j)
+               r_flags[i*dim+j] = true;
+           }
+         else if (roll == 1)
+           c_flags[i] = true;
+       }
+
+      tr.load_coarsen_flags(c_flags);
+      tr.load_refine_flags(r_flags);
+
+      tr.execute_coarsening_and_refinement ();
+      deallog << "locally owned cells: " << tr.n_locally_owned_active_cells()
+             << " / "
+             << tr.n_global_active_cells()
+             << std::endl;
+
+
+      LA::MPI::Vector vec1;
+      LA::MPI::Vector lr_vec1;
+      setup(dh, fe, vec1, lr_vec1);
+
+      LA::MPI::Vector vec2;
+      LA::MPI::Vector lr_vec2;
+      setup(dh2, fe, vec2, lr_vec2);
+
+      // interpolate func on old mesh:
+      VectorTools::interpolate(dh2, func, vec2);
+      lr_vec2 = vec2;
+
+      // interpolate from vec2 to vec1
+      VectorTools::interpolate_to_different_mesh(dh2, lr_vec2, dh, vec1);
+      lr_vec1 = vec1;
+
+      {
+       Vector<double> local_errors (tr.n_active_cells());
+       VectorTools::integrate_difference       (dh, lr_vec1, func, local_errors,
+                                                QGauss<dim>(3), VectorTools::L2_norm);
+       double total_local_error = local_errors.l2_norm();
+       const double total_global_error
+         = std::sqrt (Utilities::MPI::sum (total_local_error * total_local_error, MPI_COMM_WORLD));
+       if (myid == 0)
+         deallog << "err: " << total_global_error << std::endl;
+      }
+      
+      //output(dh, lr_vec1, loop, "solutionA-");
+      //output(dh2, lr_vec2, loop, "solutionB-");
+            
+      // also update tr2 to be the same as tr1
+      tr2.load_coarsen_flags(c_flags);
+      tr2.load_refine_flags(r_flags);
+      tr2.execute_coarsening_and_refinement ();
+
+      // repartition both in the same way
+      tr.repartition();
+      tr2.repartition();
+
+      // checks:
+      const unsigned int checksum = tr.get_checksum ();
+      const unsigned int checksum2 = tr2.get_checksum ();
+      if (myid == 0)
+       deallog << "Checksum: "
+               << checksum << " " << checksum2
+               << std::endl;
+    }
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll log;
+  test<2>();
+}
diff --git a/tests/mpi/interpolate_to_different_mesh_01.with_trilinos=true.mpirun=3.output b/tests/mpi/interpolate_to_different_mesh_01.with_trilinos=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..60a0396
--- /dev/null
@@ -0,0 +1,30 @@
+
+DEAL:0::locally owned cells: 0 / 49
+DEAL:0::err: 0.0114109
+DEAL:0::Checksum: 946342831 946342831
+DEAL:0::locally owned cells: 57 / 148
+DEAL:0::err: 0.00598894
+DEAL:0::Checksum: 585311429 585311429
+DEAL:0::locally owned cells: 148 / 415
+DEAL:0::err: 0.00203731
+DEAL:0::Checksum: 2102618372 2102618372
+DEAL:0::locally owned cells: 378 / 1162
+DEAL:0::err: 0.000842838
+DEAL:0::Checksum: 1959145913 1959145913
+DEAL:0::locally owned cells: 1146 / 3505
+DEAL:0::err: 0.000338503
+DEAL:0::Checksum: 1046055222 1046055222
+
+DEAL:1::locally owned cells: 0 / 49
+DEAL:1::locally owned cells: 51 / 148
+DEAL:1::locally owned cells: 140 / 415
+DEAL:1::locally owned cells: 395 / 1162
+DEAL:1::locally owned cells: 1172 / 3505
+
+
+DEAL:2::locally owned cells: 49 / 49
+DEAL:2::locally owned cells: 40 / 148
+DEAL:2::locally owned cells: 127 / 415
+DEAL:2::locally owned cells: 389 / 1162
+DEAL:2::locally owned cells: 1187 / 3505
+
diff --git a/tests/mpi/interpolate_to_different_mesh_02.cc b/tests/mpi/interpolate_to_different_mesh_02.cc
new file mode 100644 (file)
index 0000000..e91d5ef
--- /dev/null
@@ -0,0 +1,301 @@
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/distributed/solution_transfer.h>
+#include <fstream>
+#include <deal.II/base/conditional_ostream.h>
+#include <cmath>
+#include <deal.II/lac/generic_linear_algebra.h>
+#include <deal.II/lac/sparsity_tools.h>
+
+#include "../tests.h"
+
+// test VectorTools::interpolate_to_different_mesh in parallel
+// this is a slightly modified version from the example by Sam Cox from the mailing
+
+using namespace dealii;
+
+namespace LA
+{
+  using namespace dealii::LinearAlgebraTrilinos;
+}
+
+template <int dim>
+class SeventhProblem
+{
+public:
+    SeventhProblem (unsigned int prob_number);
+    ~SeventhProblem ();
+    void run (unsigned int cycle);
+private:
+    void setup_system ();
+    void setup_second_system ();
+    void assemble_system ();
+    void solve ();
+    MPI_Comm mpi_communicator;
+    parallel::distributed::Triangulation<2>::Settings settings;
+    parallel::distributed::Triangulation<dim> triangulation;
+    DoFHandler<dim> dof_handler;
+    FE_Q<dim> fe;
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+    ConstraintMatrix constraints;
+    TrilinosWrappers::SparseMatrix system_matrix;
+    TrilinosWrappers::MPI::Vector locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector interpolated_locally_relevant_solution;
+    TrilinosWrappers::MPI::Vector system_rhs;
+    parallel::distributed::Triangulation<dim> second_triangulation;
+    DoFHandler<dim> second_dof_handler;
+    FE_Q<dim> second_fe;
+    IndexSet second_locally_owned_dofs;
+    IndexSet second_locally_relevant_dofs;
+    TrilinosWrappers::MPI::Vector second_locally_relevant_solution;
+    ConditionalOStream pcout;
+    unsigned int prob_number;
+};
+template <int dim>
+SeventhProblem<dim>::SeventhProblem (unsigned int prob_number)
+:
+mpi_communicator (MPI_COMM_WORLD),
+settings (parallel::distributed::Triangulation<2,2>::no_automatic_repartitioning),
+triangulation (mpi_communicator,
+               typename Triangulation<dim>::MeshSmoothing
+               (Triangulation<dim>::smoothing_on_refinement |
+                Triangulation<dim>::smoothing_on_coarsening),
+               settings),
+dof_handler (triangulation),
+fe (2),
+second_triangulation (mpi_communicator,
+                      typename Triangulation<dim>::MeshSmoothing
+                      (Triangulation<dim>::smoothing_on_refinement |
+                       Triangulation<dim>::smoothing_on_coarsening),
+                      settings),
+second_dof_handler (second_triangulation),
+second_fe (2),
+pcout (deallog.get_file_stream(),
+       (Utilities::MPI::this_mpi_process(mpi_communicator)
+        == 0)),
+prob_number(prob_number)
+{}
+
+template <int dim>
+SeventhProblem<dim>::~SeventhProblem ()
+{
+    dof_handler.clear ();
+    second_dof_handler.clear ();
+}
+
+template <int dim>
+void SeventhProblem<dim>::setup_system ()
+{
+    dof_handler.distribute_dofs (fe);
+    locally_owned_dofs = dof_handler.locally_owned_dofs ();
+    DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                             locally_relevant_dofs);
+    locally_relevant_solution.reinit (locally_owned_dofs,
+                                      locally_relevant_dofs, mpi_communicator);
+    system_rhs.reinit (locally_owned_dofs, mpi_communicator);
+    system_rhs = 0;
+    constraints.clear ();
+    constraints.reinit (locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                              0,
+                                              ZeroFunction<dim>(),
+                                              constraints);
+    constraints.close ();
+    CompressedSimpleSparsityPattern csp (locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                     constraints, false);
+    SparsityTools::distribute_sparsity_pattern (csp,
+                                                dof_handler.n_locally_owned_dofs_per_processor(),
+                                                mpi_communicator,
+                                                locally_relevant_dofs);
+    system_matrix.reinit (locally_owned_dofs,
+                          locally_owned_dofs,
+                          csp,
+                          mpi_communicator);
+}
+
+template <int dim>
+void SeventhProblem<dim>::setup_second_system ()
+{
+    second_dof_handler.distribute_dofs (fe);
+    second_locally_owned_dofs = second_dof_handler.locally_owned_dofs ();
+    DoFTools::extract_locally_relevant_dofs (second_dof_handler,
+                                             second_locally_relevant_dofs);
+    second_locally_relevant_solution.reinit (second_locally_owned_dofs,
+                                             second_locally_relevant_dofs, mpi_communicator);
+    interpolated_locally_relevant_solution.reinit(second_locally_owned_dofs,
+                                                  second_locally_relevant_dofs, mpi_communicator);
+}
+
+template <int dim>
+void SeventhProblem<dim>::assemble_system ()
+{
+    const QGauss<dim> quadrature_formula(3);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values | update_gradients |
+                             update_quadrature_points |
+                             update_JxW_values);
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points = quadrature_formula.size();
+    FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double> cell_rhs (dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+        if (cell->is_locally_owned())
+        {
+            cell_matrix = 0;
+            cell_rhs = 0;
+            fe_values.reinit (cell);
+            for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+            {
+                const double
+                rhs_value
+                = (fe_values.quadrature_point(q_point)[1]
+                   >
+                   0.5+0.25*std::sin(4.0 * numbers::PI *
+                                     fe_values.quadrature_point(q_point)[0])
+                   ? 1 : -1);
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                {
+                    for (unsigned int j=0; j<dofs_per_cell; ++j)
+                        cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+                                             fe_values.shape_grad(j,q_point) *
+                                             fe_values.JxW(q_point));
+                    cell_rhs(i) += (rhs_value *
+                                    fe_values.shape_value(i,q_point) *
+                                    fe_values.JxW(q_point));
+                }
+            }
+            cell->get_dof_indices (local_dof_indices);
+            constraints.distribute_local_to_global (cell_matrix,
+                                                    cell_rhs,
+                                                    local_dof_indices,
+                                                    system_matrix,
+                                                    system_rhs);
+        }
+    system_matrix.compress (VectorOperation::add);
+    system_rhs.compress (VectorOperation::add);
+}
+template <int dim>
+void SeventhProblem<dim>::solve ()
+{
+    LA::MPI::Vector
+    completely_distributed_solution (locally_owned_dofs, mpi_communicator);
+    SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
+    TrilinosWrappers::SolverCG solver(solver_control, mpi_communicator);
+    TrilinosWrappers::PreconditionAMG preconditioner;
+    TrilinosWrappers::PreconditionAMG::AdditionalData data;
+    preconditioner.initialize(system_matrix, data);
+    solver.solve (system_matrix, completely_distributed_solution, system_rhs,
+                  preconditioner);
+    pcout << " Solved in " << solver_control.last_step()
+    << " iterations." << std::endl;
+    constraints.distribute (completely_distributed_solution);
+    locally_relevant_solution = completely_distributed_solution;
+}
+
+template <int dim>
+void SeventhProblem<dim>::run (unsigned int cycle)
+{
+    if (cycle == 0)
+    {
+        GridGenerator::hyper_cube (triangulation);
+        triangulation.refine_global (1);
+        GridGenerator::hyper_cube (second_triangulation);
+        second_triangulation.refine_global (1);
+        setup_system();
+    }
+    else
+    {
+
+        Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+        KellyErrorEstimator<dim>::estimate (dof_handler,
+                                            QGauss<dim-1>(3),
+                                            typename FunctionMap<dim>::type(),
+                                            locally_relevant_solution,
+                                            estimated_error_per_cell);
+        
+        parallel::distributed::GridRefinement::
+        refine_and_coarsen_fixed_number (triangulation,
+                                         estimated_error_per_cell,
+                                         0.5, 0.3);
+       std::vector<bool> r_flags;
+       std::vector<bool> c_flags;
+        triangulation.prepare_coarsening_and_refinement();
+        triangulation.save_refine_flags(r_flags);
+        triangulation.save_coarsen_flags(c_flags);        
+        
+        triangulation.execute_coarsening_and_refinement ();
+        
+        setup_system();
+        pcout << " Number of active cells: "
+        << triangulation.n_global_active_cells()
+        << std::endl
+        << " Number of degrees of freedom: "
+        << dof_handler.n_dofs()
+        << std::endl;
+        assemble_system ();
+        solve ();
+        
+        setup_second_system();
+        second_locally_relevant_solution = locally_relevant_solution;
+
+        VectorTools::interpolate_to_different_mesh(dof_handler, locally_relevant_solution,
+                                                   second_dof_handler, interpolated_locally_relevant_solution);
+        second_triangulation.load_coarsen_flags(c_flags);
+        second_triangulation.load_refine_flags(r_flags);
+        second_triangulation.execute_coarsening_and_refinement();
+    }
+}
+
+void seventh_grid()
+{
+    
+  ConditionalOStream           pcout(deallog.get_file_stream(),
+                                       (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)
+                                        == 0));
+    
+    pcout << "7th Starting" << std::endl;
+    SeventhProblem<2> lap(1);
+    const unsigned int n_cycles = 5;
+    for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
+    {
+        pcout << "Cycle " << cycle << ':' << std::endl;
+        lap.run(cycle);
+    }
+    pcout << "OK" << std::endl;
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll log;
+  deallog.depth_file(0);
+  
+  seventh_grid ();
+}
diff --git a/tests/mpi/interpolate_to_different_mesh_02.with_trilinos=true.mpirun=3.output b/tests/mpi/interpolate_to_different_mesh_02.with_trilinos=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..6246f25
--- /dev/null
@@ -0,0 +1,24 @@
+
+7th Starting
+Cycle 0:
+Cycle 1:
+ Number of active cells: 4
+ Number of degrees of freedom: 25
+ Solved in 1 iterations.
+Cycle 2:
+ Number of active cells: 16
+ Number of degrees of freedom: 81
+ Solved in 1 iterations.
+Cycle 3:
+ Number of active cells: 40
+ Number of degrees of freedom: 197
+ Solved in 1 iterations.
+Cycle 4:
+ Number of active cells: 124
+ Number of degrees of freedom: 581
+ Solved in 1 iterations.
+OK
+
+
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.