]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New test.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 27 Oct 2012 11:31:42 +0000 (11:31 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 27 Oct 2012 11:31:42 +0000 (11:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@27226 0785d39b-7218-0410-832d-ea1e28bc413d

tests/petsc/petsc_mf_testmatrix.h [new file with mode: 0644]
tests/petsc/solver_03_mf.cc [new file with mode: 0644]
tests/petsc/solver_03_mf/cmp/generic [new file with mode: 0644]

diff --git a/tests/petsc/petsc_mf_testmatrix.h b/tests/petsc/petsc_mf_testmatrix.h
new file mode 100644 (file)
index 0000000..093ea6b
--- /dev/null
@@ -0,0 +1,159 @@
+// $Id$
+
+#include "../tests.h"
+#include <deal.II/base/exceptions.h>
+#include <iostream>
+#include <deal.II/lac/petsc_matrix_free.h>
+#include <deal.II/lac/vector.h>
+
+
+/**
+ * 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 <i>dst = M*src</i> with
+                                          * <i>M</i> 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
+                                          * <i>dst = M<sup>T</sup>*src</i> with
+                                          * <i>M</i> 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
+                                          * <i>M*src</i> on <i>dst</i>
+                                          * with <i>M</i> 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
+                                          * <i>M<sup>T</sup>*src</i> to
+                                          * <i>dst</i> with <i>M</i> 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 (file)
index 0000000..a914c25
--- /dev/null
@@ -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 <cmath>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/vector_memory.h>
+#include <typeinfo>
+
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+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 (file)
index 0000000..2237281
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.