From 0fc3ab3cf01e5c808fa91a9937ef2bc2b3a0e7fa Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Mon, 18 Apr 2016 17:35:56 +0200 Subject: [PATCH] Test for periodic boundary conditions using DG and MeshWorker, update other tests --- tests/bits/step-4_dg_periodic.cc | 322 ++++++++++++++++++ tests/bits/step-4_dg_periodic.output | 57 ++++ tests/mpi/periodicity_01.cc | 22 +- ...odicity_01.with_petsc=true.mpirun=3.output | 276 +++++++-------- ...odicity_01.with_petsc=true.mpirun=5.output | 276 +++++++-------- ...odicity_01.with_petsc=true.mpirun=7.output | 276 +++++++-------- 6 files changed, 781 insertions(+), 448 deletions(-) create mode 100644 tests/bits/step-4_dg_periodic.cc create mode 100644 tests/bits/step-4_dg_periodic.output diff --git a/tests/bits/step-4_dg_periodic.cc b/tests/bits/step-4_dg_periodic.cc new file mode 100644 index 0000000000..a252698fd8 --- /dev/null +++ b/tests/bits/step-4_dg_periodic.cc @@ -0,0 +1,322 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 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. +// +// --------------------------------------------------------------------- + + + +// solves a 2D Poisson equation with FE_DGQ elements and periodic boundary conditions + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include + + +template +class MatrixIntegrator : public MeshWorker::LocalIntegrator +{ +public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; +}; + +template +void MatrixIntegrator +::cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, + info.fe_values()); +} + +template +void MatrixIntegrator +::boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const +{ + const unsigned int deg = info.fe_values(0).get_fe().degree; + LocalIntegrators::Laplace + ::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0), + LocalIntegrators::Laplace:: + compute_penalty(dinfo, dinfo, deg, deg)); +} + +template +void MatrixIntegrator +::face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const +{ + const unsigned int deg = info1.fe_values(0).get_fe().degree; + LocalIntegrators::Laplace + ::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, + dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix, + info1.fe_values(0), info2.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); +} + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void solve (); + void output_results (const unsigned int cycle) const; + void check_periodicity (const unsigned int cycle) const; + + + Triangulation triangulation; + FE_DGQ fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + + + +template +Step4::Step4 () + : + fe (1), + dof_handler (triangulation) +{} + + +template +void Step4::make_grid () +{ + GridGenerator::hyper_cube (triangulation, -1, 1, true); + typedef typename dealii::Triangulation::cell_iterator CellIteratorTria; + std::vector > periodic_faces; + const unsigned int b_id1 = 2; + const unsigned int b_id2 = 3; + const unsigned int direction = 1; + + dealii::GridTools::collect_periodic_faces (triangulation, b_id1, b_id2, + direction, periodic_faces, dealii::Tensor<1,dim>()); + triangulation.add_periodicity(periodic_faces); + triangulation.refine_global (1); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_flux_sparsity_pattern (dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + MappingQGeneric mapping(1); + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + MeshWorker::Assembler::MatrixSimple > assembler; + assembler.initialize(system_matrix); + MatrixIntegrator integrator; + MeshWorker::integration_loop(dof_handler.begin_active(), + dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + for (unsigned int i=0; i +void Step4::solve () +{ + solution = 0; + SolverControl solver_control (1000, 1e-10, false, false); + SolverCG<> solver (solver_control); + + PreconditionJacobi > preconditioner; + preconditioner.initialize(system_matrix); + + solver.solve (system_matrix, solution, system_rhs, + preconditioner); +} + + + +template +void Step4::output_results (const unsigned int cycle) const +{ + std::string filename = "solution-"+dealii::Utilities::int_to_string(cycle,2); + + dealii::DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "u"); + + data_out.build_patches (fe.degree); + std::ofstream output ((filename + ".vtk").c_str()); + data_out.write_vtk (output); +} + + + +template +void Step4::check_periodicity (const unsigned int cycle) const +{} + +template <> +void Step4<2>::check_periodicity (const unsigned int cycle) const +{ + unsigned int n_points = 4; + for (unsigned int i = 0; i value1(1); + Vector value2(1); + + Point<2> point1; + point1(0)=2*(1.*i/n_points+eps)-1; + point1(1)=-1.; + Point<2> point2; + point2(0)=2*(1.*i/n_points+eps)-1; + point2(1)=1.; + + VectorTools::point_value (dof_handler, solution, point1, value1); + VectorTools::point_value (dof_handler, solution, point2, value2); + + const double rel_error = std::abs((value2[0]-value1[0])/value1[0]); + const double rel_tol = 1./std::pow(2., cycle); + + if (rel_error < rel_tol) + deallog << point1 << "\t pass" << std::endl; + else + { + deallog << point1 << "\t fail" << std::endl; + deallog << point1 << "\t" << value1[0] << "\t" + << value2[0] << "\t" + << rel_error << std::endl; + all_passed = false; + } + } + AssertThrow(all_passed, ExcInternalError()); +} + + + +template +void Step4::run() +{ + for (unsigned int cycle = 0; cycle < 4; ++cycle) + { + if (cycle == 0) + make_grid(); + else + triangulation.refine_global(1); + + setup_system(); + solve(); + //output_results(cycle); + deallog.push(Utilities::int_to_string(dof_handler.n_dofs(),5)); + check_periodicity(cycle); + deallog.pop(); + } +} + + +int main (int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + try + { + Step4<2> test; + test.run(); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/step-4_dg_periodic.output b/tests/bits/step-4_dg_periodic.output new file mode 100644 index 0000000000..00ac41fe93 --- /dev/null +++ b/tests/bits/step-4_dg_periodic.output @@ -0,0 +1,57 @@ + +DEAL:00016::-0.468750 -1.00000 pass +DEAL:00016::0.0312500 -1.00000 pass +DEAL:00016::0.531250 -1.00000 pass +DEAL:00064::-0.734375 -1.00000 pass +DEAL:00064::-0.484375 -1.00000 pass +DEAL:00064::-0.234375 -1.00000 pass +DEAL:00064::0.0156250 -1.00000 pass +DEAL:00064::0.265625 -1.00000 pass +DEAL:00064::0.515625 -1.00000 pass +DEAL:00064::0.765625 -1.00000 pass +DEAL:00256::-0.867188 -1.00000 pass +DEAL:00256::-0.742188 -1.00000 pass +DEAL:00256::-0.617188 -1.00000 pass +DEAL:00256::-0.492188 -1.00000 pass +DEAL:00256::-0.367188 -1.00000 pass +DEAL:00256::-0.242188 -1.00000 pass +DEAL:00256::-0.117188 -1.00000 pass +DEAL:00256::0.00781250 -1.00000 pass +DEAL:00256::0.132812 -1.00000 pass +DEAL:00256::0.257812 -1.00000 pass +DEAL:00256::0.382812 -1.00000 pass +DEAL:00256::0.507812 -1.00000 pass +DEAL:00256::0.632812 -1.00000 pass +DEAL:00256::0.757812 -1.00000 pass +DEAL:00256::0.882812 -1.00000 pass +DEAL:01024::-0.933594 -1.00000 pass +DEAL:01024::-0.871094 -1.00000 pass +DEAL:01024::-0.808594 -1.00000 pass +DEAL:01024::-0.746094 -1.00000 pass +DEAL:01024::-0.683594 -1.00000 pass +DEAL:01024::-0.621094 -1.00000 pass +DEAL:01024::-0.558594 -1.00000 pass +DEAL:01024::-0.496094 -1.00000 pass +DEAL:01024::-0.433594 -1.00000 pass +DEAL:01024::-0.371094 -1.00000 pass +DEAL:01024::-0.308594 -1.00000 pass +DEAL:01024::-0.246094 -1.00000 pass +DEAL:01024::-0.183594 -1.00000 pass +DEAL:01024::-0.121094 -1.00000 pass +DEAL:01024::-0.0585938 -1.00000 pass +DEAL:01024::0.00390625 -1.00000 pass +DEAL:01024::0.0664062 -1.00000 pass +DEAL:01024::0.128906 -1.00000 pass +DEAL:01024::0.191406 -1.00000 pass +DEAL:01024::0.253906 -1.00000 pass +DEAL:01024::0.316406 -1.00000 pass +DEAL:01024::0.378906 -1.00000 pass +DEAL:01024::0.441406 -1.00000 pass +DEAL:01024::0.503906 -1.00000 pass +DEAL:01024::0.566406 -1.00000 pass +DEAL:01024::0.628906 -1.00000 pass +DEAL:01024::0.691406 -1.00000 pass +DEAL:01024::0.753906 -1.00000 pass +DEAL:01024::0.816406 -1.00000 pass +DEAL:01024::0.878906 -1.00000 pass +DEAL:01024::0.941406 -1.00000 pass diff --git a/tests/mpi/periodicity_01.cc b/tests/mpi/periodicity_01.cc index aef92fd4c8..8a2d3c70ea 100644 --- a/tests/mpi/periodicity_01.cc +++ b/tests/mpi/periodicity_01.cc @@ -255,7 +255,7 @@ namespace Step40 dof_handler.n_dofs(), dof_handler.n_locally_owned_dofs()); - SolverControl solver_control (dof_handler.n_dofs(), 1e-12); + SolverControl solver_control (dof_handler.n_dofs(), 1e-12, false, false); PETScWrappers::SolverCG solver(solver_control, mpi_communicator); @@ -273,9 +273,6 @@ namespace Step40 PETScWrappers::PreconditionJacobi(system_matrix)); #endif - pcout << " Solved in " << solver_control.last_step() - << " iterations." << std::endl; - constraints.distribute (completely_distributed_solution); locally_relevant_solution = completely_distributed_solution; @@ -353,13 +350,17 @@ namespace Step40 if (Utilities::MPI::this_mpi_process(mpi_communicator)==0) { - pcout << point1 << "\t" << get_real_assert_zero_imag(value1[0]) << std::endl; if (std::abs(value2[0]-value1[0])>1e-8) { + pcout << point1 << "\t" << "fail" << std::endl; std::cout<1e-8) { + pcout << point1 << "\t fail check 0" << std::endl; std::cout<1e-8) { + pcout << point3 << "\t fail check 1" << std::endl; std::cout<