From 32939922531bf6668be15aab06dec52c410c2cbb Mon Sep 17 00:00:00 2001 From: heltai Date: Thu, 19 Oct 2006 20:32:02 +0000 Subject: [PATCH] added test for bug in umfpack wrapper. git-svn-id: https://svn.dealii.org/trunk@14027 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/umfpack/umfpack_04.cc | 218 +++++++++++++++++++++++++++ tests/umfpack/umfpack_04/cmp/generic | 37 +++++ 2 files changed, 255 insertions(+) create mode 100644 tests/umfpack/umfpack_04.cc create mode 100644 tests/umfpack/umfpack_04/cmp/generic diff --git a/tests/umfpack/umfpack_04.cc b/tests/umfpack/umfpack_04.cc new file mode 100644 index 0000000000..9e80a220b9 --- /dev/null +++ b/tests/umfpack/umfpack_04.cc @@ -0,0 +1,218 @@ +//---------------------------------------------------------------------- +// $Id: umfpack_03.cc 13855 2006-09-08 14:42:11Z heltai $ +// Version: $Name$ +// +// Copyright (C) 2002, 2003, 2004, 2005, 2006 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 the umfpack sparse direct solver on a block mass matrix + +#include "../tests.h" +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +template +void assemble_laplace (MATRIX &B, VECTOR &bb, DoFHandler &dof_handler, FiniteElement &fe) +{ + QGauss quadrature_formula(2); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + + 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) + { + fe_values.reinit (cell); + + cell_matrix = 0; + cell_rhs = 0; + + unsigned int comp_i, comp_j; + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + for (unsigned int i=0; i +void test () +{ + deallog << dim << 'd' << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube (tria,0,1); + tria.refine_global (1); + + // destroy the uniformity of the matrix by + // refining one cell + tria.begin_active()->set_refine_flag (); + tria.execute_coarsening_and_refinement (); + tria.refine_global(7-2*dim); + + FESystem fe (FE_Q(1), dim); + DoFHandler dof_handler (tria); + dof_handler.distribute_dofs (fe); + DoFRenumbering::component_wise (dof_handler); + + deallog << "Number of dofs = " << dof_handler.n_dofs() << std::endl; + + std::vector size; + DoFTools::count_dofs_per_component(dof_handler, size); + + BlockSparsityPattern b_sparsity_pattern; + SparsityPattern sparsity_pattern; + + b_sparsity_pattern.reinit(size.size(),size.size()); + for (unsigned int k=0; k Bb; + Bb.reinit (b_sparsity_pattern); + + SparseMatrix B; + B.reinit(sparsity_pattern); + + BlockVector bb(size), bx(size); + Vector b(dof_handler.n_dofs()), x(dof_handler.n_dofs()), + ub(dof_handler.n_dofs()), ubb(dof_handler.n_dofs()); + + assemble_laplace(Bb, bb, dof_handler, fe); + assemble_laplace(B, b, dof_handler, fe); + + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(size.size()), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, Bb, bx, bb); + MatrixTools::apply_boundary_values (boundary_values, B, x, b); + + ub = b; + ubb = bb; + + SolverControl control (1000, 1e-12); + + SolverCG > bcg(control, SolverCG >::AdditionalData()); + bcg.solve(Bb, bx, bb, PreconditionIdentity()); + + SolverCG > cg(control, SolverCG >::AdditionalData()); + cg.solve(B, x, b, PreconditionIdentity()); + + deallog << "Sparse Factorization" << std::endl; + SparseDirectUMFPACK umfpack, umfpackb; + umfpack.factorize(B); + umfpack.solve(ub); + deallog << "absolute norms = " + << x.l2_norm() << ' ' << ub.l2_norm() + << std::endl; + x -= ub; + + deallog << "relative norm distance = " + << x.l2_norm() / ub.l2_norm() + << std::endl; + + deallog << "Block Sparse Factorization" << std::endl; + umfpackb.factorize(Bb); + // umfpackb.factorize(B); + umfpackb.solve(ubb); + deallog << "absolute norms = " + << bx.l2_norm() << ' ' << ubb.l2_norm() + << std::endl; + b = bx; + b -= ubb; + + deallog << "relative norm distance = " + << b.l2_norm() / ubb.l2_norm() + << std::endl; +} + + +int main () +{ + std::ofstream logfile("umfpack_04/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test<1> (); + test<2> (); + test<3> (); +} diff --git a/tests/umfpack/umfpack_04/cmp/generic b/tests/umfpack/umfpack_04/cmp/generic new file mode 100644 index 0000000000..4135273bda --- /dev/null +++ b/tests/umfpack/umfpack_04/cmp/generic @@ -0,0 +1,37 @@ + +DEAL::1d +DEAL::Number of dofs = 97 +DEAL:cg::Starting value 0.215517 +DEAL:cg::Convergence step 111 value 0 +DEAL:cg::Starting value 0.215517 +DEAL:cg::Convergence step 111 value 0 +DEAL::Sparse Factorization +DEAL::absolute norms = 6.45748 6.45748 +DEAL::relative norm distance = 0 +DEAL::Block Sparse Factorization +DEAL::absolute norms = 6.45748 6.45748 +DEAL::relative norm distance = 0 +DEAL::2d +DEAL::Number of dofs = 994 +DEAL:cg::Starting value 0.595296 +DEAL:cg::Convergence step 58 value 0 +DEAL:cg::Starting value 0.595296 +DEAL:cg::Convergence step 58 value 0 +DEAL::Sparse Factorization +DEAL::absolute norms = 10.0596 10.0596 +DEAL::relative norm distance = 0 +DEAL::Block Sparse Factorization +DEAL::absolute norms = 10.0596 10.0596 +DEAL::relative norm distance = 0 +DEAL::3d +DEAL::Number of dofs = 669 +DEAL:cg::Starting value 3.12141 +DEAL:cg::Convergence step 24 value 0 +DEAL:cg::Starting value 3.12141 +DEAL:cg::Convergence step 24 value 0 +DEAL::Sparse Factorization +DEAL::absolute norms = 13.7931 13.7931 +DEAL::relative norm distance = 0 +DEAL::Block Sparse Factorization +DEAL::absolute norms = 13.7931 13.7931 +DEAL::relative norm distance = 0 -- 2.39.5