From 1f3ea2959f8c8f243eb2f8b8f8cf929bae815b8b Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 16 Apr 2014 19:30:38 +0000 Subject: [PATCH] Add additional tests (simplied version of matrix_free/step-48) to see what is going on git-svn-id: https://svn.dealii.org/trunk@32783 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-48/step-48.cc | 2 - tests/matrix_free/step-48.cc | 4 +- tests/matrix_free/step-48b.cc | 412 ++++++++++++++++++ .../step-48b.with_mpi=true.mpirun=4.output | 19 + tests/matrix_free/step-48c.cc | 404 +++++++++++++++++ .../step-48c.with_mpi=true.mpirun=4.output | 19 + 6 files changed, 855 insertions(+), 5 deletions(-) create mode 100644 tests/matrix_free/step-48b.cc create mode 100644 tests/matrix_free/step-48b.with_mpi=true.mpirun=4.output create mode 100644 tests/matrix_free/step-48c.cc create mode 100644 tests/matrix_free/step-48c.with_mpi=true.mpirun=4.output diff --git a/deal.II/examples/step-48/step-48.cc b/deal.II/examples/step-48/step-48.cc index 95d9d67d9d..77e1f71970 100644 --- a/deal.II/examples/step-48/step-48.cc +++ b/deal.II/examples/step-48/step-48.cc @@ -286,8 +286,6 @@ namespace Step48 ConditionalOStream pcout; void make_grid_and_dofs (); - void oldstyle_operation (); - void assemble_system (); void output_results (const unsigned int timestep_number); #ifdef DEAL_II_WITH_P4EST diff --git a/tests/matrix_free/step-48.cc b/tests/matrix_free/step-48.cc index 3add3fe5d9..b920f8b52c 100644 --- a/tests/matrix_free/step-48.cc +++ b/tests/matrix_free/step-48.cc @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2012 - 2013 by the deal.II authors +// Copyright (C) 2012 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -207,8 +207,6 @@ namespace Step48 private: void make_grid_and_dofs (); - void oldstyle_operation (); - void assemble_system (); void output_results (const unsigned int timestep_number); #ifdef DEAL_II_WITH_P4EST diff --git a/tests/matrix_free/step-48b.cc b/tests/matrix_free/step-48b.cc new file mode 100644 index 0000000000..11cb3c2edb --- /dev/null +++ b/tests/matrix_free/step-48b.cc @@ -0,0 +1,412 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// simplified version of step-48, no sine term (i.e., usual wave equation) + + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + + +namespace Step48 +{ + using namespace dealii; + + const unsigned int dimension = 2; + const unsigned int fe_degree = 4; + + + + template + class SineGordonOperation + { + public: + SineGordonOperation(const MatrixFree &data_in, + const double time_step); + + void apply (parallel::distributed::Vector &dst, + const std::vector*> &src) const; + + private: + const MatrixFree &data; + const VectorizedArray delta_t_sqr; + parallel::distributed::Vector inv_mass_matrix; + + void local_apply (const MatrixFree &data, + parallel::distributed::Vector &dst, + const std::vector*> &src, + const std::pair &cell_range) const; + }; + + + + + template + SineGordonOperation:: + SineGordonOperation(const MatrixFree &data_in, + const double time_step) + : + data(data_in), + delta_t_sqr(make_vectorized_array(time_step *time_step)) + { + VectorizedArray one = make_vectorized_array (1.); + + data.initialize_dof_vector (inv_mass_matrix); + + FEEvaluationGL fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=0; cell1e-15) + inv_mass_matrix.local_element(k) = 1./inv_mass_matrix.local_element(k); + else + inv_mass_matrix.local_element(k) = 0; + } + + + + + + template + void SineGordonOperation:: + local_apply (const MatrixFree &data, + parallel::distributed::Vector &dst, + const std::vector*> &src, + const std::pair &cell_range) const + { + AssertDimension (src.size(), 2); + FEEvaluationGL current (data), old (data); + for (unsigned int cell=cell_range.first; cell current_value = current.get_value(q); + const VectorizedArray old_value = old.get_value(q); + + current.submit_value (2.*current_value - old_value - + delta_t_sqr * std::sin(current_value),q); + current.submit_gradient (- delta_t_sqr * + current.get_gradient(q), q); + } + + current.integrate (true,true); + current.distribute_local_to_global (dst); + } + } + + + + + template + void SineGordonOperation:: + apply (parallel::distributed::Vector &dst, + const std::vector*> &src) const + { + dst = 0; + data.cell_loop (&SineGordonOperation::local_apply, + this, dst, src); + dst.scale(inv_mass_matrix); + } + + + + template + class ExactSolution : public Function + { + public: + ExactSolution (const unsigned int n_components = 1, + const double time = 0.) : Function(n_components, time) {} + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + template + double ExactSolution::value (const Point &p, + const unsigned int /* component */) const + { + double t = this->get_time (); + + const double m = 0.5; + const double c1 = 0.; + const double c2 = 0.; + const double factor = (m / std::sqrt(1.-m*m) * + std::sin(std::sqrt(1.-m*m)*t+c2)); + double result = 1.; + for (unsigned int d=0; d + class SineGordonProblem + { + public: + SineGordonProblem (); + void run (); + + private: + void make_grid_and_dofs (); + void output_results (const unsigned int timestep_number); + +#ifdef DEAL_II_WITH_P4EST + parallel::distributed::Triangulation triangulation; +#else + Triangulation triangulation; +#endif + FE_Q fe; + DoFHandler dof_handler; + ConstraintMatrix constraints; + IndexSet locally_relevant_dofs; + + MatrixFree matrix_free_data; + + parallel::distributed::Vector solution, old_solution, old_old_solution; + + const unsigned int n_global_refinements; + double time, time_step; + const double final_time; + const double cfl_number; + const unsigned int output_timestep_skip; + }; + + + + template + SineGordonProblem::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), + time (-10), + final_time (-9), + cfl_number (.1/fe_degree), + output_timestep_skip (200) + {} + + + template + void SineGordonProblem::make_grid_and_dofs () + { + GridGenerator::hyper_cube (triangulation, -15, 15); + triangulation.refine_global (n_global_refinements); + + 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); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + constraints.close(); + + QGaussLobatto<1> quadrature (fe_degree+1); + typename MatrixFree::AdditionalData additional_data; + additional_data.mpi_communicator = MPI_COMM_WORLD; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_partition; + + matrix_free_data.reinit (dof_handler, constraints, + quadrature, additional_data); + + matrix_free_data.initialize_dof_vector (solution); + old_solution.reinit (solution); + old_old_solution.reinit (solution); + } + + + + + template + void + SineGordonProblem::output_results (const unsigned int timestep_number) + { + constraints.distribute (solution); + solution.update_ghost_values(); + + Vector norm_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (dof_handler, + solution, + ZeroFunction(), + norm_per_cell, + QGauss(fe_degree+1), + VectorTools::L2_norm); + const double solution_norm = + std::sqrt(Utilities::MPI::sum (norm_per_cell.norm_sqr(), MPI_COMM_WORLD)); + + deallog << " Time:" + << std::setw(8) << std::setprecision(3) << time + << ", solution norm: " + << std::setprecision(5) << std::setw(7) << solution_norm + << std::endl; + } + + + + template + void + SineGordonProblem::run () + { + make_grid_and_dofs(); + + const double local_min_cell_diameter = + triangulation.last()->diameter()/std::sqrt(dim); + const double global_min_cell_diameter + = -Utilities::MPI::max(-local_min_cell_diameter, MPI_COMM_WORLD); + time_step = cfl_number * global_min_cell_diameter; + time_step = (final_time-time)/(int((final_time-time)/time_step)); + deallog << " Time step size: " << time_step << ", finest cell: " + << global_min_cell_diameter << std::endl << std::endl; + + + VectorTools::interpolate (dof_handler, + ExactSolution (1, time), + solution); + VectorTools::interpolate (dof_handler, + ExactSolution (1, time-time_step), + old_solution); + output_results (0); + + std::vector*> previous_solutions; + previous_solutions.push_back(&old_solution); + previous_solutions.push_back(&old_old_solution); + + SineGordonOperation sine_gordon_op (matrix_free_data, + time_step); + + unsigned int timestep_number = 1; + + for (time+=time_step; time<=final_time; time+=time_step, ++timestep_number) + { + old_old_solution.swap (old_solution); + old_solution.swap (solution); + sine_gordon_op.apply (solution, previous_solutions); + + if (timestep_number % output_timestep_skip == 0) + output_results(timestep_number / output_timestep_skip); + } + output_results(timestep_number / output_timestep_skip + 1); + + deallog << std::endl + << " Performed " << timestep_number << " time steps." + << std::endl << std::endl; + } +} + + + +int main (int argc, char **argv) +{ + Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); + + 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.depth_console(0); + deallog.threshold_double(1.e-10); + + { + deallog.push("2d"); + Step48::SineGordonProblem<2> sg_problem; + sg_problem.run (); + deallog.pop(); + } + { + deallog.push("3d"); + Step48::SineGordonProblem<3> sg_problem; + sg_problem.run (); + deallog.pop(); + } + } + else + { + deallog.depth_console(0); + { + Step48::SineGordonProblem<2> sg_problem; + sg_problem.run (); + } + { + Step48::SineGordonProblem<3> sg_problem; + sg_problem.run (); + } + } +} + diff --git a/tests/matrix_free/step-48b.with_mpi=true.mpirun=4.output b/tests/matrix_free/step-48b.with_mpi=true.mpirun=4.output new file mode 100644 index 0000000000..5cf259b9be --- /dev/null +++ b/tests/matrix_free/step-48b.with_mpi=true.mpirun=4.output @@ -0,0 +1,19 @@ + +DEAL:0:2d:: Number of global active cells: 256 +DEAL:0:2d:: Number of degrees of freedom: 4225 +DEAL:0:2d:: Time step size: 0.04762, finest cell: 1.875 +DEAL:0:2d:: +DEAL:0:2d:: Time:-10.0 , solution norm: 9.5599 +DEAL:0:2d:: Time:-8.95 , solution norm: 21.848 +DEAL:0:2d:: +DEAL:0:2d:: Performed 22 time steps. +DEAL:0:2d:: +DEAL:0:3d:: Number of global active cells: 64 +DEAL:0:3d:: Number of degrees of freedom: 4913 +DEAL:0:3d:: Time step size: 0.20000, finest cell: 7.5000 +DEAL:0:3d:: +DEAL:0:3d:: Time:-10.0 , solution norm: 29.779 +DEAL:0:3d:: Time:-8.80 , solution norm: 84.725 +DEAL:0:3d:: +DEAL:0:3d:: Performed 6 time steps. +DEAL:0:3d:: diff --git a/tests/matrix_free/step-48c.cc b/tests/matrix_free/step-48c.cc new file mode 100644 index 0000000000..c14b9ebce5 --- /dev/null +++ b/tests/matrix_free/step-48c.cc @@ -0,0 +1,404 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// simplified form for step-48 test + + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + + +namespace Step48 +{ + using namespace dealii; + + const unsigned int dimension = 2; + const unsigned int fe_degree = 4; + + + + template + class SineGordonOperation + { + public: + SineGordonOperation(const MatrixFree &data_in, + const double time_step); + + void apply (parallel::distributed::Vector &dst, + const std::vector*> &src) const; + + private: + const MatrixFree &data; + const VectorizedArray delta_t_sqr; + parallel::distributed::Vector inv_mass_matrix; + + void local_apply (const MatrixFree &data, + parallel::distributed::Vector &dst, + const std::vector*> &src, + const std::pair &cell_range) const; + }; + + + + + template + SineGordonOperation:: + SineGordonOperation(const MatrixFree &data_in, + const double time_step) + : + data(data_in), + delta_t_sqr(make_vectorized_array(time_step *time_step)) + { + VectorizedArray one = make_vectorized_array (1.); + + data.initialize_dof_vector (inv_mass_matrix); + + FEEvaluationGL fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=0; cell1e-15) + inv_mass_matrix.local_element(k) = 1./inv_mass_matrix.local_element(k); + else + inv_mass_matrix.local_element(k) = 0; + } + + + + + + template + void SineGordonOperation:: + local_apply (const MatrixFree &data, + parallel::distributed::Vector &dst, + const std::vector*> &src, + const std::pair &cell_range) const + { + AssertDimension (src.size(), 2); + FEEvaluationGL current (data), old (data); + for (unsigned int cell=cell_range.first; cell current_value = current.get_value(q); + const VectorizedArray old_value = old.get_value(q); + + current.submit_value (2.*current_value - old_value,q); + current.submit_gradient (- delta_t_sqr * + current.get_gradient(q), q); + } + + current.integrate (true,true); + current.distribute_local_to_global (dst); + } + } + + + + + template + void SineGordonOperation:: + apply (parallel::distributed::Vector &dst, + const std::vector*> &src) const + { + dst = 0; + data.cell_loop (&SineGordonOperation::local_apply, + this, dst, src); + dst.scale(inv_mass_matrix); + } + + + + template + class InitialSolution : public Function + { + public: + InitialSolution (const unsigned int n_components = 1, + const double time = 0.) : Function(n_components, time) {} + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + template + double InitialSolution::value (const Point &p, + const unsigned int /* component */) const + { + return 4.*std::exp(-p.square()*10); + } + + + + + template + class SineGordonProblem + { + public: + SineGordonProblem (); + void run (); + + private: + void make_grid_and_dofs (); + void output_results (const unsigned int timestep_number); + +#ifdef DEAL_II_WITH_P4EST + parallel::distributed::Triangulation triangulation; +#else + Triangulation triangulation; +#endif + FE_Q fe; + DoFHandler dof_handler; + ConstraintMatrix constraints; + IndexSet locally_relevant_dofs; + + MatrixFree matrix_free_data; + + parallel::distributed::Vector solution, old_solution, old_old_solution; + + const unsigned int n_global_refinements; + double time, time_step; + const double final_time; + const double cfl_number; + const unsigned int output_timestep_skip; + }; + + + + template + SineGordonProblem::SineGordonProblem () + : +#ifdef DEAL_II_WITH_P4EST + triangulation (MPI_COMM_WORLD), +#endif + fe (QGaussLobatto<1>(fe_degree+1)), + dof_handler (triangulation), + n_global_refinements (9-2*dim), + time (-10), + final_time (-9), + cfl_number (.1/fe_degree), + output_timestep_skip (200) + {} + + + template + void SineGordonProblem::make_grid_and_dofs () + { + GridGenerator::hyper_cube (triangulation, -15, 15); + triangulation.refine_global (n_global_refinements); + + 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); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + constraints.clear(); + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close(); + + QGaussLobatto<1> quadrature (fe_degree+1); + typename MatrixFree::AdditionalData additional_data; + additional_data.mpi_communicator = MPI_COMM_WORLD; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_partition; + + matrix_free_data.reinit (dof_handler, constraints, + quadrature, additional_data); + + matrix_free_data.initialize_dof_vector (solution); + old_solution.reinit (solution); + old_old_solution.reinit (solution); + } + + + + + template + void + SineGordonProblem::output_results (const unsigned int timestep_number) + { + constraints.distribute (solution); + solution.update_ghost_values(); + + Vector norm_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (dof_handler, + solution, + ZeroFunction(), + norm_per_cell, + QGauss(fe_degree+1), + VectorTools::L2_norm); + const double solution_norm = + std::sqrt(Utilities::MPI::sum (norm_per_cell.norm_sqr(), MPI_COMM_WORLD)); + + deallog << " Time:" + << std::setw(8) << std::setprecision(3) << time + << ", solution norm: " + << std::setprecision(5) << std::setw(7) << solution_norm + << std::endl; + } + + + + template + void + SineGordonProblem::run () + { + make_grid_and_dofs(); + + const double local_min_cell_diameter = + triangulation.last()->diameter()/std::sqrt(dim); + const double global_min_cell_diameter + = -Utilities::MPI::max(-local_min_cell_diameter, MPI_COMM_WORLD); + time_step = cfl_number * global_min_cell_diameter; + time_step = (final_time-time)/(int((final_time-time)/time_step)); + deallog << " Time step size: " << time_step << ", finest cell: " + << global_min_cell_diameter << std::endl << std::endl; + + + VectorTools::interpolate (dof_handler, + InitialSolution (1, time), + solution); + VectorTools::interpolate (dof_handler, + InitialSolution (1, time-time_step), + old_solution); + output_results (0); + + std::vector*> previous_solutions; + previous_solutions.push_back(&old_solution); + previous_solutions.push_back(&old_old_solution); + + SineGordonOperation sine_gordon_op (matrix_free_data, + time_step); + + unsigned int timestep_number = 1; + + for (time+=time_step; time<=final_time; time+=time_step, ++timestep_number) + { + old_old_solution.swap (old_solution); + old_solution.swap (solution); + sine_gordon_op.apply (solution, previous_solutions); + + if (timestep_number % output_timestep_skip == 0) + output_results(timestep_number / output_timestep_skip); + } + output_results(timestep_number / output_timestep_skip + 1); + + deallog << std::endl + << " Performed " << timestep_number << " time steps." + << std::endl << std::endl; + } +} + + + +int main (int argc, char **argv) +{ + Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv); + + 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.depth_console(0); + deallog.threshold_double(1.e-10); + + { + deallog.push("2d"); + Step48::SineGordonProblem<2> sg_problem; + sg_problem.run (); + deallog.pop(); + } + { + deallog.push("3d"); + Step48::SineGordonProblem<3> sg_problem; + sg_problem.run (); + deallog.pop(); + } + } + else + { + deallog.depth_console(0); + { + Step48::SineGordonProblem<2> sg_problem; + sg_problem.run (); + } + { + Step48::SineGordonProblem<3> sg_problem; + sg_problem.run (); + } + } +} + diff --git a/tests/matrix_free/step-48c.with_mpi=true.mpirun=4.output b/tests/matrix_free/step-48c.with_mpi=true.mpirun=4.output new file mode 100644 index 0000000000..9639684813 --- /dev/null +++ b/tests/matrix_free/step-48c.with_mpi=true.mpirun=4.output @@ -0,0 +1,19 @@ + +DEAL:0:2d:: Number of global active cells: 1024 +DEAL:0:2d:: Number of degrees of freedom: 16641 +DEAL:0:2d:: Time step size: 0.02381, finest cell: 0.9375 +DEAL:0:2d:: +DEAL:0:2d:: Time:-10.0 , solution norm: 1.6062 +DEAL:0:2d:: Time:-8.98 , solution norm: 1.1332 +DEAL:0:2d:: +DEAL:0:2d:: Performed 43 time steps. +DEAL:0:2d:: +DEAL:0:3d:: Number of global active cells: 512 +DEAL:0:3d:: Number of degrees of freedom: 35937 +DEAL:0:3d:: Time step size: 0.10000, finest cell: 3.7500 +DEAL:0:3d:: +DEAL:0:3d:: Time:-10.0 , solution norm: 0.78144 +DEAL:0:3d:: Time:-8.90 , solution norm: 0.57949 +DEAL:0:3d:: +DEAL:0:3d:: Performed 11 time steps. +DEAL:0:3d:: -- 2.39.5