From a4e032381920a49042b8c9f4be6bb0b9eac3a49b Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 23 Nov 2012 10:04:01 +0000 Subject: [PATCH] VectorTools::interpolate should not call compress on parallel::distributed::Vector as this zeros out those entries that are not locally active but included via constraints. Test step-48 and step-37. git-svn-id: https://svn.dealii.org/trunk@27672 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/lac/parallel_vector.templates.h | 1 + .../deal.II/numerics/vector_tools.templates.h | 26 +- tests/matrix_free/step-37.cc | 660 ++++++++++++++++++ tests/matrix_free/step-37/cmp/generic | 31 + tests/mpi/step-48.cc | 426 +++++++++++ tests/mpi/step-48/ncpu_1/cmp/generic | 19 + tests/mpi/step-48/ncpu_10/cmp/generic | 19 + tests/mpi/step-48/ncpu_4/cmp/generic | 19 + 8 files changed, 1200 insertions(+), 1 deletion(-) create mode 100644 tests/matrix_free/step-37.cc create mode 100644 tests/matrix_free/step-37/cmp/generic create mode 100644 tests/mpi/step-48.cc create mode 100644 tests/mpi/step-48/ncpu_1/cmp/generic create mode 100644 tests/mpi/step-48/ncpu_10/cmp/generic create mode 100644 tests/mpi/step-48/ncpu_4/cmp/generic diff --git a/deal.II/include/deal.II/lac/parallel_vector.templates.h b/deal.II/include/deal.II/lac/parallel_vector.templates.h index cadadfcaca..d841e0a1f5 100644 --- a/deal.II/include/deal.II/lac/parallel_vector.templates.h +++ b/deal.II/include/deal.II/lac/parallel_vector.templates.h @@ -190,6 +190,7 @@ namespace parallel { AssertDimension (local_range().first, c.local_range().first); AssertDimension (local_range().second, c.local_range().second); + AssertDimension (vector_view.size(), c.vector_view.size()); vector_view = c.vector_view; if (call_update_ghost_values == true) update_ghost_values(); diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 8546036d53..6936970500 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -66,6 +66,30 @@ DEAL_II_NAMESPACE_OPEN namespace VectorTools { + namespace + { + template + void perform_compress_insert (VEC &vec) + { + vec.compress(::dealii::VectorOperation::insert); + } + + template + void perform_compress_insert (::dealii::parallel::distributed::Vector &vec) + { + // should not do compress when inserting + // elements as the ghosts are fixed and some + // will not be set at all. Instead, zero the + // ghost data afterwards to get consistent + // data. Otherwise, at least with the layout + // in r27671 the tests/mpi/step-48/ncpu_10 + // will signal incorrect results because of + // incorrect interpolation. + vec.zero_out_ghosts(); + } + + } + template void interpolate (const Mapping &mapping, @@ -259,7 +283,7 @@ namespace VectorTools = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]]; } } - vec.compress(::dealii::VectorOperation::insert); + perform_compress_insert(vec); } diff --git a/tests/matrix_free/step-37.cc b/tests/matrix_free/step-37.cc new file mode 100644 index 0000000000..005c227d17 --- /dev/null +++ b/tests/matrix_free/step-37.cc @@ -0,0 +1,660 @@ +//----------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//----------------------------------------------------------------------------- + +// test for correctness of step-37 (without output and only small sizes) + + +#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 + +#include +#include + + +namespace Step37 +{ + using namespace dealii; + + const unsigned int degree_finite_element = 2; + + + + template + class Coefficient : public Function + { + public: + Coefficient () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + template + number value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; + }; + + + + template + template + number Coefficient::value (const Point &p, + const unsigned int /*component*/) const + { + return 1. / (0.05 + 2.*p.square()); + } + + + + template + double Coefficient::value (const Point &p, + const unsigned int component) const + { + return value(p,component); + } + + + + template + void Coefficient::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const + { + Assert (values.size() == points.size(), + ExcDimensionMismatch (values.size(), points.size())); + Assert (component == 0, + ExcIndexRange (component, 0, 1)); + + const unsigned int n_points = points.size(); + for (unsigned int i=0; i(points[i],component); + } + + + + + + template + class LaplaceOperator : public Subscriptor + { + public: + LaplaceOperator (); + + void clear(); + + void reinit (const MGDoFHandler &dof_handler, + const ConstraintMatrix &constraints, + const unsigned int level = numbers::invalid_unsigned_int); + + unsigned int m () const; + unsigned int n () const; + + void vmult (Vector &dst, + const Vector &src) const; + void Tvmult (Vector &dst, + const Vector &src) const; + void vmult_add (Vector &dst, + const Vector &src) const; + void Tvmult_add (Vector &dst, + const Vector &src) const; + + number el (const unsigned int row, + const unsigned int col) const; + void set_diagonal (const Vector &diagonal); + + std::size_t memory_consumption () const; + + private: + void local_apply (const MatrixFree &data, + Vector &dst, + const Vector &src, + const std::pair &cell_range) const; + + void evaluate_coefficient(const Coefficient &function); + + MatrixFree data; + AlignedVector > coefficient; + + Vector diagonal_values; + bool diagonal_is_available; + }; + + + + template + LaplaceOperator::LaplaceOperator () + : + Subscriptor() + {} + + + + template + unsigned int + LaplaceOperator::m () const + { + return data.get_vector_partitioner()->size(); + } + + + + template + unsigned int + LaplaceOperator::n () const + { + return data.get_vector_partitioner()->size(); + } + + + + template + void + LaplaceOperator::clear () + { + data.clear(); + diagonal_is_available = false; + diagonal_values.reinit(0); + } + + + + template + void + LaplaceOperator::reinit (const MGDoFHandler &dof_handler, + const ConstraintMatrix &constraints, + const unsigned int level) + { + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + additional_data.level_mg_handler = level; + additional_data.mapping_update_flags = (update_gradients | update_JxW_values | + update_quadrature_points); + data.reinit (dof_handler, constraints, QGauss<1>(fe_degree+1), + additional_data); + evaluate_coefficient(Coefficient()); + } + + + + template + void + LaplaceOperator:: + evaluate_coefficient (const Coefficient &coefficient_function) + { + const unsigned int n_cells = data.get_size_info().n_macro_cells; + FEEvaluation phi (data); + coefficient.resize (n_cells * phi.n_q_points); + for (unsigned int cell=0; cell + void + LaplaceOperator:: + local_apply (const MatrixFree &data, + Vector &dst, + const Vector &src, + const std::pair &cell_range) const + { + FEEvaluation phi (data); + AssertDimension (coefficient.size(), + data.get_size_info().n_macro_cells * phi.n_q_points); + + for (unsigned int cell=cell_range.first; cell + void + LaplaceOperator::vmult (Vector &dst, + const Vector &src) const + { + dst = 0; + vmult_add (dst, src); + } + + + + template + void + LaplaceOperator::Tvmult (Vector &dst, + const Vector &src) const + { + dst = 0; + vmult_add (dst,src); + } + + + + template + void + LaplaceOperator::Tvmult_add (Vector &dst, + const Vector &src) const + { + vmult_add (dst,src); + } + + + + template + void + LaplaceOperator::vmult_add (Vector &dst, + const Vector &src) const + { + data.cell_loop (&LaplaceOperator::local_apply, this, dst, src); + + const std::vector & + constrained_dofs = data.get_constrained_dofs(); + for (unsigned int i=0; i + number + LaplaceOperator::el (const unsigned int row, + const unsigned int col) const + { + Assert (row == col, ExcNotImplemented()); + Assert (diagonal_is_available == true, ExcNotInitialized()); + return diagonal_values(row); + } + + + + template + void + LaplaceOperator::set_diagonal(const Vector &diagonal) + { + AssertDimension (m(), diagonal.size()); + + diagonal_values = diagonal; + + const std::vector & + constrained_dofs = data.get_constrained_dofs(); + for (unsigned int i=0; i + std::size_t + LaplaceOperator::memory_consumption () const + { + return (data.memory_consumption () + + coefficient.memory_consumption() + + diagonal_values.memory_consumption() + + MemoryConsumption::memory_consumption(diagonal_is_available)); + } + + + + + template + class LaplaceProblem + { + public: + LaplaceProblem (); + void run (); + + private: + void setup_system (); + void assemble_system (); + void assemble_multigrid (); + void solve (); + void output_results (const unsigned int cycle) const; + + typedef LaplaceOperator SystemMatrixType; + typedef LaplaceOperator LevelMatrixType; + + Triangulation triangulation; + FE_Q fe; + MGDoFHandler mg_dof_handler; + ConstraintMatrix constraints; + + SystemMatrixType system_matrix; + MGLevelObject mg_matrices; + FullMatrix coarse_matrix; + MGLevelObject mg_constraints; + + Vector solution; + Vector system_rhs; + }; + + + + template + LaplaceProblem::LaplaceProblem () + : + fe (degree_finite_element), + mg_dof_handler (triangulation) + {} + + + + + template + void LaplaceProblem::setup_system () + { + system_matrix.clear(); + mg_matrices.clear(); + mg_constraints.clear(); + + mg_dof_handler.distribute_dofs (fe); + + deallog << "Number of degrees of freedom: " + << mg_dof_handler.n_dofs() + << std::endl; + + constraints.clear(); + VectorTools::interpolate_boundary_values (mg_dof_handler, + 0, + ZeroFunction(), + constraints); + constraints.close(); + + system_matrix.reinit (mg_dof_handler, constraints); + solution.reinit (mg_dof_handler.n_dofs()); + system_rhs.reinit (mg_dof_handler.n_dofs()); + + const unsigned int nlevels = triangulation.n_levels(); + mg_matrices.resize(0, nlevels-1); + mg_constraints.resize (0, nlevels-1); + + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (1); + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + std::vector > boundary_indices(triangulation.n_levels()); + MGTools::make_boundary_list (mg_dof_handler, + dirichlet_boundary, + boundary_indices); + for (unsigned int level=0; level::iterator bc_it = boundary_indices[level].begin(); + for ( ; bc_it != boundary_indices[level].end(); ++bc_it) + mg_constraints[level].add_line(*bc_it); + + mg_constraints[level].close(); + mg_matrices[level].reinit(mg_dof_handler, + mg_constraints[level], + level); + } + coarse_matrix.reinit (mg_dof_handler.n_dofs(0), + mg_dof_handler.n_dofs(0)); + } + + + + + template + void LaplaceProblem::assemble_system () + { + QGauss quadrature_formula(fe.degree+1); + FEValues fe_values (fe, quadrature_formula, + update_values | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + std::vector local_dof_indices (dofs_per_cell); + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename DoFHandler::active_cell_iterator cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell->get_dof_indices (local_dof_indices); + fe_values.reinit (cell); + for (unsigned int i=0; i + void LaplaceProblem::assemble_multigrid () + { + coarse_matrix = 0; + QGauss quadrature_formula(fe.degree+1); + FEValues fe_values (fe, quadrature_formula, + update_gradients | update_inverse_jacobians | + 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(); + + std::vector local_dof_indices (dofs_per_cell); + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_diagonal (dofs_per_cell); + + const unsigned int n_levels = triangulation.n_levels(); + std::vector > diagonals (n_levels); + for (unsigned int level=0; level cell_no(triangulation.n_levels()); + typename MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + const unsigned int level = cell->level(); + cell->get_mg_dof_indices (local_dof_indices); + fe_values.reinit (cell); + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int i=0; i + void LaplaceProblem::solve () + { + GrowingVectorMemory<> vector_memory; + + MGTransferPrebuilt > mg_transfer; + mg_transfer.build_matrices(mg_dof_handler); + + MGCoarseGridHouseholder > mg_coarse; + mg_coarse.initialize(coarse_matrix); + + typedef PreconditionChebyshev > SMOOTHER; + MGSmootherPrecondition > + mg_smoother(vector_memory); + + typename SMOOTHER::AdditionalData smoother_data; + smoother_data.smoothing_range = 10.; + smoother_data.degree = 6; + smoother_data.eig_cg_n_iterations = 10; + mg_smoother.initialize(mg_matrices, smoother_data); + + MGMatrix > + mg_matrix(&mg_matrices); + + Multigrid > mg(mg_dof_handler, + mg_matrix, + mg_coarse, + mg_transfer, + mg_smoother, + mg_smoother); + PreconditionMG, + MGTransferPrebuilt > > + preconditioner(mg_dof_handler, mg, mg_transfer); + + const std::size_t multigrid_memory + = (mg_matrices.memory_consumption() + + mg_transfer.memory_consumption() + + coarse_matrix.memory_consumption()); + + SolverControl solver_control (1000, 1e-12*system_rhs.l2_norm()); + SolverCG<> cg (solver_control); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + } + + + + template + void LaplaceProblem::run () + { + for (unsigned int cycle=0; cycle<3; ++cycle) + { + deallog << "Cycle " << cycle << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, 0., 1.); + triangulation.refine_global (3-dim); + } + triangulation.refine_global (1); + setup_system (); + assemble_system (); + assemble_multigrid (); + solve (); + deallog << std::endl; + }; + } +} + + + +int main () +{ + std::ofstream logfile("step-37/output"); + deallog.attach(logfile); + deallog << std::setprecision (3); + deallog.threshold_double(1e-10); + deallog.depth_console(0); + + { + deallog.push("2d"); + Step37::LaplaceProblem<2> laplace_problem; + laplace_problem.run(); + deallog.pop(); + } + { + deallog.push("3d"); + Step37::LaplaceProblem<3> laplace_problem; + laplace_problem.run(); + deallog.pop(); + } + + + return 0; +} + diff --git a/tests/matrix_free/step-37/cmp/generic b/tests/matrix_free/step-37/cmp/generic new file mode 100644 index 0000000000..4ac0a4971c --- /dev/null +++ b/tests/matrix_free/step-37/cmp/generic @@ -0,0 +1,31 @@ + +DEAL:2d::Cycle 0 +DEAL:2d::Number of degrees of freedom: 81 +DEAL:2d:cg::Starting value 0.132 +DEAL:2d:cg::Convergence step 5 value 0 +DEAL:2d:: +DEAL:2d::Cycle 1 +DEAL:2d::Number of degrees of freedom: 289 +DEAL:2d:cg::Starting value 0.0677 +DEAL:2d:cg::Convergence step 5 value 0 +DEAL:2d:: +DEAL:2d::Cycle 2 +DEAL:2d::Number of degrees of freedom: 1089 +DEAL:2d:cg::Starting value 0.0343 +DEAL:2d:cg::Convergence step 6 value 0 +DEAL:2d:: +DEAL:3d::Cycle 0 +DEAL:3d::Number of degrees of freedom: 125 +DEAL:3d:cg::Starting value 0.125 +DEAL:3d:cg::Convergence step 5 value 0 +DEAL:3d:: +DEAL:3d::Cycle 1 +DEAL:3d::Number of degrees of freedom: 729 +DEAL:3d:cg::Starting value 0.0479 +DEAL:3d:cg::Convergence step 5 value 0 +DEAL:3d:: +DEAL:3d::Cycle 2 +DEAL:3d::Number of degrees of freedom: 4913 +DEAL:3d:cg::Starting value 0.0176 +DEAL:3d:cg::Convergence step 4 value 0 +DEAL:3d:: diff --git a/tests/mpi/step-48.cc b/tests/mpi/step-48.cc new file mode 100644 index 0000000000..d8092d49d8 --- /dev/null +++ b/tests/mpi/step-48.cc @@ -0,0 +1,426 @@ +//----------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//----------------------------------------------------------------------------- + +// test for correctness of step-48 (without output and only small sizes) + + +#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 oldstyle_operation (); + void assemble_system (); + void output_results (const unsigned int timestep_number) const; + + parallel::distributed::Triangulation triangulation; + 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 () + : + triangulation (MPI_COMM_WORLD), + 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); + { + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + end_cell = triangulation.end(); + for ( ; cell != end_cell; ++cell) + if (cell->is_locally_owned()) + if (cell->center().norm() < 11) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + + cell = triangulation.begin_active(); + end_cell = triangulation.end(); + for ( ; cell != end_cell; ++cell) + if (cell->is_locally_owned()) + if (cell->center().norm() < 6) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + } + + deallog << " Number of global active cells: " + << triangulation.n_global_active_cells() + << 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) const + { + parallel::distributed::Vector locally_relevant_solution; + locally_relevant_solution.reinit (dof_handler.locally_owned_dofs(), + locally_relevant_dofs, + MPI_COMM_WORLD); + locally_relevant_solution.copy_from (solution); + locally_relevant_solution.update_ghost_values (); + constraints.distribute (locally_relevant_solution); + + Vector norm_per_cell (triangulation.n_active_cells()); + VectorTools::integrate_difference (dof_handler, + locally_relevant_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::System::get_this_mpi_process (MPI_COMM_WORLD); + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + std::ofstream logfile(output_file_for_mpi("step-48").c_str()); + 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/mpi/step-48/ncpu_1/cmp/generic b/tests/mpi/step-48/ncpu_1/cmp/generic new file mode 100644 index 0000000000..f0b2c39aba --- /dev/null +++ b/tests/mpi/step-48/ncpu_1/cmp/generic @@ -0,0 +1,19 @@ + +DEAL:0:2d:: Number of global active cells: 964 +DEAL:0:2d:: Number of degrees of freedom: 16033 +DEAL:0:2d:: Time step size: 0.01176, finest cell: 0.4688 +DEAL:0:2d:: +DEAL:0:2d:: Time:-10.0 , solution norm: 9.5599 +DEAL:0:2d:: Time:-8.99 , solution norm: 21.965 +DEAL:0:2d:: +DEAL:0:2d:: Performed 86 time steps. +DEAL:0:2d:: +DEAL:0:3d:: Number of global active cells: 176 +DEAL:0:3d:: Number of degrees of freedom: 14001 +DEAL:0:3d:: Time step size: 0.047619, finest cell: 1.8750 +DEAL:0:3d:: +DEAL:0:3d:: Time:-10.0 , solution norm: 29.559 +DEAL:0:3d:: Time:-8.95 , solution norm: 90.729 +DEAL:0:3d:: +DEAL:0:3d:: Performed 22 time steps. +DEAL:0:3d:: diff --git a/tests/mpi/step-48/ncpu_10/cmp/generic b/tests/mpi/step-48/ncpu_10/cmp/generic new file mode 100644 index 0000000000..f0b2c39aba --- /dev/null +++ b/tests/mpi/step-48/ncpu_10/cmp/generic @@ -0,0 +1,19 @@ + +DEAL:0:2d:: Number of global active cells: 964 +DEAL:0:2d:: Number of degrees of freedom: 16033 +DEAL:0:2d:: Time step size: 0.01176, finest cell: 0.4688 +DEAL:0:2d:: +DEAL:0:2d:: Time:-10.0 , solution norm: 9.5599 +DEAL:0:2d:: Time:-8.99 , solution norm: 21.965 +DEAL:0:2d:: +DEAL:0:2d:: Performed 86 time steps. +DEAL:0:2d:: +DEAL:0:3d:: Number of global active cells: 176 +DEAL:0:3d:: Number of degrees of freedom: 14001 +DEAL:0:3d:: Time step size: 0.047619, finest cell: 1.8750 +DEAL:0:3d:: +DEAL:0:3d:: Time:-10.0 , solution norm: 29.559 +DEAL:0:3d:: Time:-8.95 , solution norm: 90.729 +DEAL:0:3d:: +DEAL:0:3d:: Performed 22 time steps. +DEAL:0:3d:: diff --git a/tests/mpi/step-48/ncpu_4/cmp/generic b/tests/mpi/step-48/ncpu_4/cmp/generic new file mode 100644 index 0000000000..f0b2c39aba --- /dev/null +++ b/tests/mpi/step-48/ncpu_4/cmp/generic @@ -0,0 +1,19 @@ + +DEAL:0:2d:: Number of global active cells: 964 +DEAL:0:2d:: Number of degrees of freedom: 16033 +DEAL:0:2d:: Time step size: 0.01176, finest cell: 0.4688 +DEAL:0:2d:: +DEAL:0:2d:: Time:-10.0 , solution norm: 9.5599 +DEAL:0:2d:: Time:-8.99 , solution norm: 21.965 +DEAL:0:2d:: +DEAL:0:2d:: Performed 86 time steps. +DEAL:0:2d:: +DEAL:0:3d:: Number of global active cells: 176 +DEAL:0:3d:: Number of degrees of freedom: 14001 +DEAL:0:3d:: Time step size: 0.047619, finest cell: 1.8750 +DEAL:0:3d:: +DEAL:0:3d:: Time:-10.0 , solution norm: 29.559 +DEAL:0:3d:: Time:-8.95 , solution norm: 90.729 +DEAL:0:3d:: +DEAL:0:3d:: Performed 22 time steps. +DEAL:0:3d:: -- 2.39.5