From: Timo Heister Date: Tue, 21 Apr 2015 16:59:05 +0000 (-0400) Subject: add tests X-Git-Tag: v8.3.0-rc1~222^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8472326d7b06f55ade89893eb03d8290ea5a8698;p=dealii.git add tests --- diff --git a/tests/mpi/interpolate_to_different_mesh_01.cc b/tests/mpi/interpolate_to_different_mesh_01.cc new file mode 100644 index 0000000000..8d82629fdb --- /dev/null +++ b/tests/mpi/interpolate_to_different_mesh_01.cc @@ -0,0 +1,197 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +// test VectorTools::interpolate_to_different_mesh in parallel + +using namespace dealii; + +namespace LA +{ + using namespace dealii::LinearAlgebraTrilinos; +} + +template +class SomeFunction : public Function +{ +public: + double value (const Point &p, + const unsigned int) const + { + return 1+p(0)*p(0); + } +}; + +template +void setup(DoFHandler & dh, + FE_Q & fe, + LA::MPI::Vector & vec, + LA::MPI::Vector & lr_vec) +{ + dh.distribute_dofs (fe); + vec.reinit(dh.locally_owned_dofs(), MPI_COMM_WORLD); + IndexSet locally_relevant; + DoFTools::extract_locally_relevant_dofs (dh, locally_relevant); + lr_vec.reinit(locally_relevant, MPI_COMM_WORLD); +} + +template +void output(DoFHandler & dh, LA::MPI::Vector & v, unsigned int loop, const std::string filename_) +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + DataOut data_out; + data_out.add_data_vector (dh, v, "1"); + data_out.build_patches (1); + std::ostringstream filename; + filename << filename_ + << Utilities::int_to_string (loop, 2) + << "." + << Utilities::int_to_string (myid,2) + << ".vtu"; + + std::ofstream output (filename.str().c_str()); + data_out.write_vtu (output); + if (myid) + { + std::vector filenames; + for (unsigned int i=0; i +void test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + parallel::distributed::Triangulation tr(MPI_COMM_WORLD, + dealii::Triangulation::none, + parallel::distributed::Triangulation::no_automatic_repartitioning); + + GridGenerator::hyper_cube(tr); + tr.refine_global(2); + parallel::distributed::Triangulation tr2(MPI_COMM_WORLD, + dealii::Triangulation::none, + parallel::distributed::Triangulation::no_automatic_repartitioning); + + GridGenerator::hyper_cube(tr2); + tr2.refine_global(2); + + FE_Q fe(1); + DoFHandler dh(tr); + DoFHandler dh2(tr2); + + SomeFunction func; + + + for (unsigned int loop = 0; loop < 5; ++loop) + { + // randomly refine: + std::vector r_flags(tr.n_active_cells()*dim, false); + std::vector c_flags(tr.n_active_cells(), false); + + for (unsigned int i=0;i= 2) + { + for (unsigned int j=0;j local_errors (tr.n_active_cells()); + VectorTools::integrate_difference (dh, lr_vec1, func, local_errors, + QGauss(3), VectorTools::L2_norm); + double total_local_error = local_errors.l2_norm(); + const double total_global_error + = std::sqrt (Utilities::MPI::sum (total_local_error * total_local_error, MPI_COMM_WORLD)); + if (myid == 0) + deallog << "err: " << total_global_error << std::endl; + } + + //output(dh, lr_vec1, loop, "solutionA-"); + //output(dh2, lr_vec2, loop, "solutionB-"); + + // also update tr2 to be the same as tr1 + tr2.load_coarsen_flags(c_flags); + tr2.load_refine_flags(r_flags); + tr2.execute_coarsening_and_refinement (); + + // repartition both in the same way + tr.repartition(); + tr2.repartition(); + + // checks: + const unsigned int checksum = tr.get_checksum (); + const unsigned int checksum2 = tr2.get_checksum (); + if (myid == 0) + deallog << "Checksum: " + << checksum << " " << checksum2 + << std::endl; + } +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + test<2>(); +} diff --git a/tests/mpi/interpolate_to_different_mesh_01.with_trilinos=true.mpirun=3.output b/tests/mpi/interpolate_to_different_mesh_01.with_trilinos=true.mpirun=3.output new file mode 100644 index 0000000000..60a0396ee6 --- /dev/null +++ b/tests/mpi/interpolate_to_different_mesh_01.with_trilinos=true.mpirun=3.output @@ -0,0 +1,30 @@ + +DEAL:0::locally owned cells: 0 / 49 +DEAL:0::err: 0.0114109 +DEAL:0::Checksum: 946342831 946342831 +DEAL:0::locally owned cells: 57 / 148 +DEAL:0::err: 0.00598894 +DEAL:0::Checksum: 585311429 585311429 +DEAL:0::locally owned cells: 148 / 415 +DEAL:0::err: 0.00203731 +DEAL:0::Checksum: 2102618372 2102618372 +DEAL:0::locally owned cells: 378 / 1162 +DEAL:0::err: 0.000842838 +DEAL:0::Checksum: 1959145913 1959145913 +DEAL:0::locally owned cells: 1146 / 3505 +DEAL:0::err: 0.000338503 +DEAL:0::Checksum: 1046055222 1046055222 + +DEAL:1::locally owned cells: 0 / 49 +DEAL:1::locally owned cells: 51 / 148 +DEAL:1::locally owned cells: 140 / 415 +DEAL:1::locally owned cells: 395 / 1162 +DEAL:1::locally owned cells: 1172 / 3505 + + +DEAL:2::locally owned cells: 49 / 49 +DEAL:2::locally owned cells: 40 / 148 +DEAL:2::locally owned cells: 127 / 415 +DEAL:2::locally owned cells: 389 / 1162 +DEAL:2::locally owned cells: 1187 / 3505 + diff --git a/tests/mpi/interpolate_to_different_mesh_02.cc b/tests/mpi/interpolate_to_different_mesh_02.cc new file mode 100644 index 0000000000..e91d5efc04 --- /dev/null +++ b/tests/mpi/interpolate_to_different_mesh_02.cc @@ -0,0 +1,301 @@ +#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 +#include +#include + +#include "../tests.h" + +// test VectorTools::interpolate_to_different_mesh in parallel +// this is a slightly modified version from the example by Sam Cox from the mailing + +using namespace dealii; + +namespace LA +{ + using namespace dealii::LinearAlgebraTrilinos; +} + +template +class SeventhProblem +{ +public: + SeventhProblem (unsigned int prob_number); + ~SeventhProblem (); + void run (unsigned int cycle); +private: + void setup_system (); + void setup_second_system (); + void assemble_system (); + void solve (); + MPI_Comm mpi_communicator; + parallel::distributed::Triangulation<2>::Settings settings; + parallel::distributed::Triangulation triangulation; + DoFHandler dof_handler; + FE_Q fe; + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + ConstraintMatrix constraints; + TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::MPI::Vector locally_relevant_solution; + TrilinosWrappers::MPI::Vector interpolated_locally_relevant_solution; + TrilinosWrappers::MPI::Vector system_rhs; + parallel::distributed::Triangulation second_triangulation; + DoFHandler second_dof_handler; + FE_Q second_fe; + IndexSet second_locally_owned_dofs; + IndexSet second_locally_relevant_dofs; + TrilinosWrappers::MPI::Vector second_locally_relevant_solution; + ConditionalOStream pcout; + unsigned int prob_number; +}; +template +SeventhProblem::SeventhProblem (unsigned int prob_number) +: +mpi_communicator (MPI_COMM_WORLD), +settings (parallel::distributed::Triangulation<2,2>::no_automatic_repartitioning), +triangulation (mpi_communicator, + typename Triangulation::MeshSmoothing + (Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening), + settings), +dof_handler (triangulation), +fe (2), +second_triangulation (mpi_communicator, + typename Triangulation::MeshSmoothing + (Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening), + settings), +second_dof_handler (second_triangulation), +second_fe (2), +pcout (deallog.get_file_stream(), + (Utilities::MPI::this_mpi_process(mpi_communicator) + == 0)), +prob_number(prob_number) +{} + +template +SeventhProblem::~SeventhProblem () +{ + dof_handler.clear (); + second_dof_handler.clear (); +} + +template +void SeventhProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + locally_owned_dofs = dof_handler.locally_owned_dofs (); + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + locally_relevant_solution.reinit (locally_owned_dofs, + locally_relevant_dofs, mpi_communicator); + system_rhs.reinit (locally_owned_dofs, mpi_communicator); + system_rhs = 0; + constraints.clear (); + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + constraints); + constraints.close (); + CompressedSimpleSparsityPattern csp (locally_relevant_dofs); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + SparsityTools::distribute_sparsity_pattern (csp, + dof_handler.n_locally_owned_dofs_per_processor(), + mpi_communicator, + locally_relevant_dofs); + system_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); +} + +template +void SeventhProblem::setup_second_system () +{ + second_dof_handler.distribute_dofs (fe); + second_locally_owned_dofs = second_dof_handler.locally_owned_dofs (); + DoFTools::extract_locally_relevant_dofs (second_dof_handler, + second_locally_relevant_dofs); + second_locally_relevant_solution.reinit (second_locally_owned_dofs, + second_locally_relevant_dofs, mpi_communicator); + interpolated_locally_relevant_solution.reinit(second_locally_owned_dofs, + second_locally_relevant_dofs, mpi_communicator); +} + +template +void SeventhProblem::assemble_system () +{ + const QGauss quadrature_formula(3); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + 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(); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + cell_rhs = 0; + fe_values.reinit (cell); + for (unsigned int q_point=0; q_point + 0.5+0.25*std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) + ? 1 : -1); + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + system_matrix.compress (VectorOperation::add); + system_rhs.compress (VectorOperation::add); +} +template +void SeventhProblem::solve () +{ + LA::MPI::Vector + completely_distributed_solution (locally_owned_dofs, mpi_communicator); + SolverControl solver_control (dof_handler.n_dofs(), 1e-12); + TrilinosWrappers::SolverCG solver(solver_control, mpi_communicator); + TrilinosWrappers::PreconditionAMG preconditioner; + TrilinosWrappers::PreconditionAMG::AdditionalData data; + preconditioner.initialize(system_matrix, data); + solver.solve (system_matrix, completely_distributed_solution, system_rhs, + preconditioner); + pcout << " Solved in " << solver_control.last_step() + << " iterations." << std::endl; + constraints.distribute (completely_distributed_solution); + locally_relevant_solution = completely_distributed_solution; +} + +template +void SeventhProblem::run (unsigned int cycle) +{ + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (1); + GridGenerator::hyper_cube (second_triangulation); + second_triangulation.refine_global (1); + setup_system(); + } + else + { + + Vector estimated_error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + QGauss(3), + typename FunctionMap::type(), + locally_relevant_solution, + estimated_error_per_cell); + + parallel::distributed::GridRefinement:: + refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.5, 0.3); + std::vector r_flags; + std::vector c_flags; + triangulation.prepare_coarsening_and_refinement(); + triangulation.save_refine_flags(r_flags); + triangulation.save_coarsen_flags(c_flags); + + triangulation.execute_coarsening_and_refinement (); + + setup_system(); + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() + << std::endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + assemble_system (); + solve (); + + setup_second_system(); + second_locally_relevant_solution = locally_relevant_solution; + + VectorTools::interpolate_to_different_mesh(dof_handler, locally_relevant_solution, + second_dof_handler, interpolated_locally_relevant_solution); + second_triangulation.load_coarsen_flags(c_flags); + second_triangulation.load_refine_flags(r_flags); + second_triangulation.execute_coarsening_and_refinement(); + } +} + +void seventh_grid() +{ + + ConditionalOStream pcout(deallog.get_file_stream(), + (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) + == 0)); + + pcout << "7th Starting" << std::endl; + SeventhProblem<2> lap(1); + const unsigned int n_cycles = 5; + for (unsigned int cycle=0; cycle