-// this tests the correctness of matrix free matrix-vector products by
-// comparing the result with a Trilinos sparse matrix assembled in the usual
-// way. The mesh is distributed among processors (hypercube) and has both
-// hanging nodes (by randomly refining some cells, so the mesh is going to be
-// different when run with different numbers of processors) and Dirichlet
-// boundary conditions
+// this tests the correctness of MPI-parallel matrix free matrix-vector
+// products by comparing the result with a Trilinos sparse matrix assembled in
+// the usual way. The mesh is distributed among processors (hypercube) and has
+// both hanging nodes (by randomly refining some cells, so the mesh is going
+// to be different when run with different numbers of processors) and
+// Dirichlet boundary conditions
//
// this test does not use multithreading within the MPI nodes.
#include "../tests.h"
-#include <deal.II/matrix_free/matrix_free.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
+#include "matrix_vector_mf.h"
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/function.h>
-#include <deal.II/lac/parallel_vector.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <iostream>
-template <int dim, int fe_degree, typename Number>
-void
-helmholtz_operator (const MatrixFree<dim,Number> &data,
- parallel::distributed::Vector<Number> &dst,
- const parallel::distributed::Vector<Number> &src,
- const std::pair<unsigned int,unsigned int> &cell_range)
-{
- FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
- const unsigned int n_q_points = fe_eval.n_q_points;
-
- for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
- {
- fe_eval.reinit (cell);
-
- fe_eval.read_dof_values (src);
- fe_eval.evaluate (true, true, false);
- for (unsigned int q=0; q<n_q_points; ++q)
- {
- fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
- fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
- }
- fe_eval.integrate (true,true);
- fe_eval.distribute_local_to_global (dst);
- }
-}
-
-
-
-template <int dim, int fe_degree, typename Number>
-class MatrixFreeTest
-{
-public:
- typedef VectorizedArray<Number> vector_t;
- static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
-
- MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
- data (data_in)
- {};
-
- void vmult (parallel::distributed::Vector<Number> &dst,
- const parallel::distributed::Vector<Number> &src) const
- {
- dst = 0;
- const std_cxx1x::function<void(const MatrixFree<dim,Number> &,
- parallel::distributed::Vector<Number> &,
- const parallel::distributed::Vector<Number> &,
- const std::pair<unsigned int,unsigned int> &)>
- wrap = helmholtz_operator<dim,fe_degree,Number>;
- data.cell_loop (wrap, dst, src);
- };
-
-private:
- const MatrixFree<dim,Number> &data;
-};
-
-
-
template <int dim, int fe_degree>
mf_data.reinit (dof, constraints, quad, data);
}
- MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+ MatrixFreeTest<dim,fe_degree,number,parallel::distributed::Vector<number> > mf (mf_data);
parallel::distributed::Vector<number> in, out, ref;
mf_data.initialize_dof_vector (in);
out.reinit (in);
// this tests the correctness of matrix free matrix-vector products when
// threads are additionally used to MPI by comparing the result with a
-// non-parallel version. Otherwise same problem as matrix_free_01.cc
+// non-parallel version. Otherwise same problem as matrix_vector_10.cc
#include "../tests.h"
-#include <deal.II/matrix_free/matrix_free.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
+#include "matrix_vector_mf.h"
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <iostream>
-template <int dim, int fe_degree, typename Number>
-void
-helmholtz_operator (const MatrixFree<dim,Number> &data,
- parallel::distributed::Vector<Number> &dst,
- const parallel::distributed::Vector<Number> &src,
- const std::pair<unsigned int,unsigned int> &cell_range)
-{
- FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
- const unsigned int n_q_points = fe_eval.n_q_points;
-
- for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
- {
- fe_eval.reinit (cell);
-
- fe_eval.read_dof_values (src);
- fe_eval.evaluate (true, true, false);
- for (unsigned int q=0; q<n_q_points; ++q)
- {
- fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
- fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
- }
- fe_eval.integrate (true,true);
- fe_eval.distribute_local_to_global (dst);
- }
-}
-
-
-
-template <int dim, int fe_degree, typename Number>
-class MatrixFreeTest
-{
-public:
- typedef VectorizedArray<Number> vector_t;
- static const std::size_t n_vectors = VectorizedArray<Number>::n_array_elements;
-
- MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
- data (data_in)
- {};
-
- void vmult (parallel::distributed::Vector<Number> &dst,
- const parallel::distributed::Vector<Number> &src) const
- {
- dst = 0;
- const std_cxx1x::function<void(const MatrixFree<dim,Number> &,
- parallel::distributed::Vector<Number> &,
- const parallel::distributed::Vector<Number> &,
- const std::pair<unsigned int,unsigned int> &)>
- wrap = helmholtz_operator<dim,fe_degree,Number>;
- data.cell_loop (wrap, dst, src);
- };
-
-private:
- const MatrixFree<dim,Number> &data;
-};
-
-
-
template <int dim, int fe_degree>
mf_data.reinit (dof, constraints, quad, data);
}
- MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+ MatrixFreeTest<dim,fe_degree,number,parallel::distributed::Vector<number> > mf (mf_data);
parallel::distributed::Vector<number> in, out, ref;
mf_data.initialize_dof_vector (in);
out.reinit (in);
data.tasks_block_size = 3;
mf_data.reinit (dof, constraints, quad, data);
- MatrixFreeTest<dim, fe_degree, number> mf (mf_data);
+ MatrixFreeTest<dim, fe_degree, number,parallel::distributed::Vector<number> > mf (mf_data);
deallog << "Norm of difference:";
// run 10 times to make a possible error more
// this tests the correctness of matrix free matrix-vector products for two
-// vectors on the same DoFHandler. Otherwise the same as matrix_free_01
+// vectors on the same DoFHandler. Otherwise the same as matrix_vector_10.cc
#include "../tests.h"
#include "../tests.h"
-#include <deal.II/matrix_free/matrix_free.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
+#include "matrix_vector_mf.h"
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
-template <int dim, int fe_degree, typename Number>
-void
-helmholtz_operator (const MatrixFree<dim,Number> &data,
- Vector<Number> &dst,
- const Vector<Number> &src,
- const std::pair<unsigned int,unsigned int> &cell_range)
-{
- FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
- const unsigned int n_q_points = fe_eval.n_q_points;
-
- for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
- {
- fe_eval.reinit (cell);
- fe_eval.read_dof_values (src);
- fe_eval.evaluate (true, true, false);
- for (unsigned int q=0; q<n_q_points; ++q)
- {
- fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
- fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
- }
- fe_eval.integrate (true,true);
- fe_eval.distribute_local_to_global (dst);
- }
-}
-
-
-
-template <int dim, int fe_degree, typename Number>
-class MatrixFreeTest
-{
-public:
- typedef VectorizedArray<Number> vector_t;
-
- MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
- data (data_in)
- {};
-
- void vmult (Vector<Number> &dst,
- const Vector<Number> &src) const
- {
- dst = 0;
- const std_cxx1x::function<void(const MatrixFree<dim,Number> &,
- Vector<Number> &,
- const Vector<Number> &,
- const std::pair<unsigned int,unsigned int> &)>
- wrap = helmholtz_operator<dim,fe_degree,Number>;
- data.cell_loop (wrap, dst, src);
- };
-
-private:
- const MatrixFree<dim,Number> &data;
-};
-
-
-
-
template <int dim, int fe_degree, typename number>
void do_test (const DoFHandler<dim> &dof,
mf.vmult (out_dist, in_dist);
- // assemble sparse matrix with (\nabla v,
- // \nabla u) + (v, 10 * u)
+ // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
SparsityPattern sparsity;
{
CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
std::pair<unsigned int,unsigned int> subrange_deg =
data.create_cell_subrange_hp (cell_range, 1);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,1,Number> (data, dst, src,
+ helmholtz_operator<dim,1,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 2);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,2,Number> (data, dst, src,
+ helmholtz_operator<dim,2,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 3);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,3,Number> (data, dst, src,
+ helmholtz_operator<dim,3,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 4);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,4,Number> (data, dst, src,
+ helmholtz_operator<dim,4,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 5);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,5,Number> (data, dst, src,
+ helmholtz_operator<dim,5,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 6);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,6,Number> (data, dst, src,
+ helmholtz_operator<dim,6,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 7);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,7,Number> (data, dst, src,
+ helmholtz_operator<dim,7,Vector<Number> > (data, dst, src,
subrange_deg);
}
--- /dev/null
+//------------------ matrix_vector_common.h ------------------------
+// $Id$
+// Version: $Name$
+//
+//------------------ matrix_vector_common.h ------------------------
+
+
+// this is a template for matrix-vector products with the Helmholtz equation
+// (zero and first derivatives) on different kinds of meshes (Cartesian,
+// general, with and without hanging nodes). It also tests the multithreading
+// in case it was enabled
+
+#include "../tests.h"
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/parallel_vector.h>
+
+
+template <int dim, int fe_degree, typename VECTOR>
+void
+helmholtz_operator (const MatrixFree<dim,typename VECTOR::value_type> &data,
+ VECTOR &dst,
+ const VECTOR &src,
+ const std::pair<unsigned int,unsigned int> &cell_range)
+{
+ typedef typename VECTOR::value_type Number;
+ FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+ const unsigned int n_q_points = fe_eval.n_q_points;
+
+ for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+ {
+ fe_eval.reinit (cell);
+ fe_eval.read_dof_values (src);
+ fe_eval.evaluate (true, true, false);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
+ fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
+ }
+ fe_eval.integrate (true,true);
+ fe_eval.distribute_local_to_global (dst);
+ }
+}
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+ typedef VectorizedArray<Number> vector_t;
+
+ MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+ data (data_in)
+ {};
+
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const
+ {
+ dst = 0;
+ const std_cxx1x::function<void(const MatrixFree<dim,typename VECTOR::value_type> &,
+ VECTOR &,
+ const VECTOR &,
+ const std::pair<unsigned int,unsigned int> &)>
+ wrap = helmholtz_operator<dim,fe_degree,VECTOR>;
+ data.cell_loop (wrap, dst, src);
+ };
+
+private:
+ const MatrixFree<dim,Number> &data;
+};
+
+
+
+
void assemble_system ();
void output_results (const unsigned int timestep_number);
+#ifdef DEAL_II_WITH_P4EST
parallel::distributed::Triangulation<dim> triangulation;
+#else
+ Triangulation<dim> triangulation;
+#endif
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
template <int dim>
SineGordonProblem<dim>::SineGordonProblem ()
:
+#ifdef DEAL_II_WITH_P4EST
triangulation (MPI_COMM_WORLD),
+#endif
fe (QGaussLobatto<1>(fe_degree+1)),
dof_handler (triangulation),
n_global_refinements (8-2*dim),
}
deallog << " Number of global active cells: "
+#ifdef DEAL_II_WITH_P4EST
<< triangulation.n_global_active_cells()
+#else
+ << triangulation.n_active_cells()
+#endif
<< std::endl;
dof_handler.distribute_dofs (fe);
std::pair<unsigned int,unsigned int> subrange_deg =
data.create_cell_subrange_hp (cell_range, 1);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,1,Number> (data, dst, src,
+ helmholtz_operator<dim,1,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 2);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,2,Number> (data, dst, src,
+ helmholtz_operator<dim,2,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 3);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,3,Number> (data, dst, src,
+ helmholtz_operator<dim,3,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 4);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,4,Number> (data, dst, src,
+ helmholtz_operator<dim,4,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 5);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,5,Number> (data, dst, src,
+ helmholtz_operator<dim,5,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 6);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,6,Number> (data, dst, src,
+ helmholtz_operator<dim,6,Vector<Number> > (data, dst, src,
subrange_deg);
subrange_deg = data.create_cell_subrange_hp (cell_range, 7);
if (subrange_deg.second > subrange_deg.first)
- helmholtz_operator<dim,7,Number> (data, dst, src,
+ helmholtz_operator<dim,7,Vector<Number> > (data, dst, src,
subrange_deg);
}
projects / dealii.git / commitdiff