From 70f3c2803e29558b7749824ec50625a3ce6586ab Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 27 Oct 2016 15:41:58 +0200 Subject: [PATCH] Add MatrixFreeOperators::MassOperator --- doc/news/changes.h | 4 + include/deal.II/matrix_free/operators.h | 126 +++++++++- tests/matrix_free/mass_operator_01.cc | 228 ++++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=1.output | 13 + ...h_mpi=true.with_p4est=true.mpirun=4.output | 13 + 5 files changed, 383 insertions(+), 1 deletion(-) create mode 100644 tests/matrix_free/mass_operator_01.cc create mode 100644 tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output create mode 100644 tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output diff --git a/doc/news/changes.h b/doc/news/changes.h index f4fb4d43d7..4d4a07a74a 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -391,6 +391,10 @@ inconvenience this causes.

Specific improvements

    +
  1. New: Add MatrixFreeOperators::MassOperator representing a mass matrix. +
    + (Daniel Arndt, 2016/10/27) +
  3. Fixed: GridIn::read_vtk() accidentally only read material ids of input cells correctly if the file listed them as integers. If they were diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index c6dfc50a70..24e6e510b6 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -158,7 +158,7 @@ namespace MatrixFreeOperators virtual void compute_diagonal () = 0; /** - * Get read access to diagonal of this operator. + * Get read access to the inverse diagonal of this operator. */ const LinearAlgebra::distributed::Vector &get_matrix_diagonal_inverse() const; @@ -288,6 +288,46 @@ namespace MatrixFreeOperators + /** + * This class implements the operation of the action of a mass matrix. + * + * @author Daniel Arndt, 2016 + */ + template + class MassOperator : public Base + { + public: + + /** + * Constructor. Initializes the MatrixFree object. + */ + MassOperator (); + + /** + * For preconditioning, we store a lumped mass matrix at the diagonal entries. + */ + virtual void compute_diagonal (); + + private: + /** + * Applies the mass matrix operation on an input vector. It is + * assumed that the passed input and output vector are correctly initialized + * using initialize_dof_vector(). + */ + virtual void apply_add (LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &src) const; + + /** + * For this operator, there is just a cell contribution. + */ + void local_apply_cell (const MatrixFree &data, + LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; + }; + + + // ------------------------------------ inline functions --------------------- template @@ -717,6 +757,90 @@ namespace MatrixFreeOperators dst*= omega; } + + + //-----------------------------MassOperator---------------------------------- + + template + MassOperator:: + MassOperator () + : + Base() + {} + + + + template + void + MassOperator:: + compute_diagonal() + { + Assert((Base::data != NULL), ExcNotInitialized()); + + LinearAlgebra::distributed::Vector ones; + Base::initialize_dof_vector(Base::inverse_diagonal_entries); + Base::initialize_dof_vector(ones); + Base::inverse_diagonal_entries = 1.; + ones = 1.; + ones.update_ghost_values(); + apply_add(Base::inverse_diagonal_entries, ones); + + const IndexSet &locally_owned_elements + = Base::inverse_diagonal_entries.locally_owned_elements(); + locally_owned_elements.print(std::cout); + const std::vector constrained_dofs + = Base::data->get_constrained_dofs(); + for (unsigned int i=0; i::inverse_diagonal_entries(constrained_dofs[i]) = 1.; + + const unsigned int local_size = Base::inverse_diagonal_entries.local_size(); + for (unsigned int i=0; i::inverse_diagonal_entries.local_element(i); + Base::inverse_diagonal_entries.local_element(i) + =1./Base::inverse_diagonal_entries.local_element(i); + } + + Base::inverse_diagonal_entries.compress(VectorOperation::insert); + Base::inverse_diagonal_entries.update_ghost_values(); + } + + + + template + void + MassOperator:: + apply_add (LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &src) const + { + Base::data->cell_loop (&MassOperator::local_apply_cell, + this, dst, src); + } + + + + template + void + MassOperator:: + local_apply_cell (const MatrixFree &data, + LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const + { + FEEvaluation phi(*Base::data); + for (unsigned int cell=cell_range.first; cell +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void test () +{ + typedef double number; + + parallel::distributed::Triangulation tria (MPI_COMM_WORLD); + GridGenerator::hyper_cube (tria); + tria.refine_global(1); + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end(); + cell = tria.begin_active (); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + if (cell->center().norm()<0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + if (dim < 3 && fe_degree < 2) + tria.refine_global(2); + else + tria.refine_global(1); + if (tria.begin(tria.n_levels()-1)->is_locally_owned()) + tria.begin(tria.n_levels()-1)->set_refine_flag(); + if (tria.last()->is_locally_owned()) + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + cell = tria.begin_active (); + for (unsigned int i=0; i<10-3*dim; ++i) + { + cell = tria.begin_active (); + unsigned int counter = 0; + for (; cell!=endc; ++cell, ++counter) + if (cell->is_locally_owned()) + if (counter % (7-i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dof, relevant_set); + + ConstraintMatrix constraints (relevant_set); + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction(), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + //std::cout << "Number of cells: " << tria.n_global_active_cells() << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + data.mpi_communicator = MPI_COMM_WORLD; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + data.tasks_block_size = 7; + mf_data.reinit (dof, constraints, quad, data); + } + + MatrixFreeOperators::MassOperator mf; + mf.initialize(mf_data); + mf.compute_diagonal(); + LinearAlgebra::distributed::Vector in, out, ref; + mf_data.initialize_dof_vector (in); + out.reinit (in); + ref.reinit (in); + + for (unsigned int i=0; i quadrature_formula(fe_degree+1); + + FEValues fe_values (dof.get_fe(), quadrature_formula, + update_values | update_gradients | + update_JxW_values); + + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + sparse_matrix); + } + } + sparse_matrix.compress(VectorOperation::add); + + sparse_matrix.vmult (ref, in); + out -= ref; + const double diff_norm = out.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); + + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog << std::setprecision(4); + deallog.threshold_double(1.e-10); + + deallog.push("2d"); + test<2,1>(); + test<2,2>(); + deallog.pop(); + + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } + else + { + test<2,1>(); + test<2,2>(); + test<3,1>(); + test<3,2>(); + } +} + diff --git a/tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..25e13d5f8f --- /dev/null +++ b/tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,13 @@ + +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Norm of difference: 0 +DEAL:0:2d:: +DEAL:0:2d::Testing FE_Q<2>(2) +DEAL:0:2d::Norm of difference: 0 +DEAL:0:2d:: +DEAL:0:3d::Testing FE_Q<3>(1) +DEAL:0:3d::Norm of difference: 0 +DEAL:0:3d:: +DEAL:0:3d::Testing FE_Q<3>(2) +DEAL:0:3d::Norm of difference: 0 +DEAL:0:3d:: diff --git a/tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output b/tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output new file mode 100644 index 0000000000..25e13d5f8f --- /dev/null +++ b/tests/matrix_free/mass_operator_01.with_trilinos=true.with_mpi=true.with_p4est=true.mpirun=4.output @@ -0,0 +1,13 @@ + +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Norm of difference: 0 +DEAL:0:2d:: +DEAL:0:2d::Testing FE_Q<2>(2) +DEAL:0:2d::Norm of difference: 0 +DEAL:0:2d:: +DEAL:0:3d::Testing FE_Q<3>(1) +DEAL:0:3d::Norm of difference: 0 +DEAL:0:3d:: +DEAL:0:3d::Testing FE_Q<3>(2) +DEAL:0:3d::Norm of difference: 0 +DEAL:0:3d:: -- 2.39.5