From 0fa29e7781bb876564cc382232ce2e3cb01e34b1 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 27 Oct 2012 11:31:42 +0000 Subject: [PATCH] New test. git-svn-id: https://svn.dealii.org/trunk@27226 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/petsc/petsc_mf_testmatrix.h | 159 +++++++++++++++++++++++++++ tests/petsc/solver_03_mf.cc | 88 +++++++++++++++ tests/petsc/solver_03_mf/cmp/generic | 6 + 3 files changed, 253 insertions(+) create mode 100644 tests/petsc/petsc_mf_testmatrix.h create mode 100644 tests/petsc/solver_03_mf.cc create mode 100644 tests/petsc/solver_03_mf/cmp/generic diff --git a/tests/petsc/petsc_mf_testmatrix.h b/tests/petsc/petsc_mf_testmatrix.h new file mode 100644 index 0000000000..093ea6b6d5 --- /dev/null +++ b/tests/petsc/petsc_mf_testmatrix.h @@ -0,0 +1,159 @@ +// $Id$ + +#include "../tests.h" +#include +#include +#include +#include + + +/** + * Finite difference PETSc matrix-free object on uniform grid. + * Generator for simple 5-point discretization of Laplace problem. + */ + +class PetscFDMatrix : public dealii::PETScWrappers::MatrixFree +{ + public: + /** + * Constructor specifying grid resolution. + */ + PetscFDMatrix(unsigned int size, unsigned int dim); + + /** + * Matrix-vector multiplication: + * let dst = M*src with + * M being this matrix. + * + * Source and destination must + * not be the same vector. + */ + void vmult (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const; + + /** + * Matrix-vector multiplication: let + * dst = MT*src with + * M being this matrix. This + * function does the same as @p vmult() + * but takes the transposed matrix. + * + * Source and destination must + * not be the same vector. + */ + void Tvmult (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add + * M*src on dst + * with M being this + * matrix. + * + * Source and destination must + * not be the same vector. + */ + void vmult_add (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add + * MT*src to + * dst with M being + * this matrix. This function + * does the same as @p vmult_add() + * but takes the transposed + * matrix. + * + * Source and destination must + * not be the same vector. + */ + void Tvmult_add (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const; + +private: + /** + * Number of gridpoints in x-direction. + */ + unsigned int nx; + + /** + * Number of gridpoints in y-direction. + */ + unsigned int ny; +}; + + +// --------------- inline and template functions ----------------- + +inline +PetscFDMatrix::PetscFDMatrix(unsigned int size, unsigned int dim) + : PETScWrappers::MatrixFree (dim, dim, dim, dim), + nx (size), ny (size) +{} + + + +inline +void +PetscFDMatrix::vmult_add (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const +{ + for(unsigned int i=0;i<=ny-2;i++) + { + for(unsigned int j=0;j<=nx-2; j++) + { + // Number of the row to be entered + unsigned int row = j+(nx-1)*i; + + dst(row) += 4. * src(row); // A.set(row, row, 4.); + + if (j>0) + { + dst(row-1) += -1. * src(row); // A.set(row-1, row, -1.); + dst(row) += -1. * src(row-1); // A.set(row, row-1, -1.); + } + if (i>0) + { + dst(row-(nx-1)) += -1. * src(row); // A.set(row-(nx-1), row, -1.); + dst(row) += -1. * src(row-(nx-1)); // A.set(row, row-(nx-1), -1.); + } + } + } +} + + + +inline +void +PetscFDMatrix::vmult (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const +{ + dst = 0; + vmult_add (dst, src); +} + + + +inline +void +PetscFDMatrix::Tvmult (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const +{ + dst = 0; + vmult_add (dst, src); +} + + + +inline +void +PetscFDMatrix::Tvmult_add (dealii::PETScWrappers::VectorBase &dst, + const dealii::PETScWrappers::VectorBase &src) const +{ + vmult_add (dst, src); +} + + diff --git a/tests/petsc/solver_03_mf.cc b/tests/petsc/solver_03_mf.cc new file mode 100644 index 0000000000..a914c25cd2 --- /dev/null +++ b/tests/petsc/solver_03_mf.cc @@ -0,0 +1,88 @@ +//---------------------------- petsc_solver_03_mf.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2010 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. +// +//---------------------------- petsc_solver_03_mf.cc --------------------------- + +// test the PETSc CG solver with PETSc MatrixFree class + + +#include "../tests.h" +#include "../lac/petsc_mf_testmatrix.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void +check_solve( SOLVER& solver, const MATRIX& A, + VECTOR& u, VECTOR& f, const PRECONDITION& P) +{ + deallog << "Solver type: " << typeid(solver).name() << std::endl; + + u = 0.; + f = 1.; + try + { + solver.solve(A,u,f,P); + } + catch (std::exception& e) + { + std::cout << e.what() << std::endl; + deallog << e.what() << std::endl; + abort (); + } + + deallog << "Solver stopped after " << solver.control().last_step() + << " iterations" << std::endl; +} + + +int main(int argc, char **argv) +{ + std::ofstream logfile("solver_03_mf/output"); + deallog.attach(logfile); + deallog << std::setprecision(4); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + PetscInitialize(&argc,&argv,0,0); + { + SolverControl control(100, 1.e-3); + + const unsigned int size = 32; + unsigned int dim = (size-1)*(size-1); + + deallog << "Size " << size << " Unknowns " << dim << std::endl; + + PetscFDMatrix A(size, dim); + + PETScWrappers::Vector f(dim); + PETScWrappers::Vector u(dim); + f = 1.; + A.compress (); + f.compress (); + u.compress (); + + PETScWrappers::SolverCG solver(control); + PETScWrappers::PreconditionNone preconditioner(A); + check_solve (solver, A,u,f, preconditioner); + } + PetscFinalize (); +} + diff --git a/tests/petsc/solver_03_mf/cmp/generic b/tests/petsc/solver_03_mf/cmp/generic new file mode 100644 index 0000000000..2237281800 --- /dev/null +++ b/tests/petsc/solver_03_mf/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::Size 32 Unknowns 961 +DEAL::Solver type: N6dealii13PETScWrappers8SolverCGE +DEAL::Starting value 7.750 +DEAL::Convergence step 41 value 0.0006730 +DEAL::Solver stopped after 41 iterations -- 2.39.5