]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
verify consistency of block relaxation and block preconditiioning methods
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 20 May 2011 02:07:54 +0000 (02:07 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 20 May 2011 02:07:54 +0000 (02:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@23728 0785d39b-7218-0410-832d-ea1e28bc413d

tests/lac/solver_relaxation_01.cc
tests/lac/solver_relaxation_02.cc [new file with mode: 0644]
tests/lac/solver_relaxation_02/cmp/generic [new file with mode: 0644]

index 36aae8a74da3aeccebc4a2bd26e5e280c2a16b4f..b5a680f11c567fde846fe866db1d7006925b8ac6 100644 (file)
@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 1998 - 2005, 2010 by the deal.II authors
+//    Copyright (C) 1998 - 2005, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -48,7 +48,8 @@ check_solve( SOLVER& solver, const MATRIX& A,
 
 int main()
 {
-  std::ofstream logfile("solver_relaxation_01/output");
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
 //  logfile.setf(std::ios::fixed);
   deallog << std::setprecision(4);
   deallog.attach(logfile);
diff --git a/tests/lac/solver_relaxation_02.cc b/tests/lac/solver_relaxation_02.cc
new file mode 100644 (file)
index 0000000..2f36510
--- /dev/null
@@ -0,0 +1,153 @@
+//----------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2011 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.
+//
+//----------------------------------------------------------------------
+
+// Compare preconditioned Richardson with block relaxation. All output diffs
+// should be zero.
+
+#include "../tests.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/solver_relaxation.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/relaxation_block.h>
+
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+double
+check_solve( SOLVER& solver, const MATRIX& A,
+            VECTOR& u, VECTOR& f, const PRECONDITION& P)
+{
+  double result = 0.;
+  u = 0.;
+  f = 1.;
+  try 
+    {
+      solver.solve(A,u,f,P);
+    }
+  catch (SolverControl::NoConvergence& e)
+    {
+      result = e.last_residual;
+    }
+  return result;
+}
+
+int main()
+{
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
+//  logfile.setf(std::ios::fixed);
+  deallog << std::setprecision(4);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  
+  SolverControl control(100, 1.e-3);
+  SolverRichardson<> rich(control);
+  SolverRelaxation<> relax(control);
+  
+  for (unsigned int size=33; size <= 33; size *= 3)
+    {
+      unsigned int dim = (size-1)*(size-1);
+
+      deallog << "Size " << size << " Unknowns " << dim << std::endl;
+      
+                                      // Make matrix
+      FDMatrix testproblem(size, size);
+      SparsityPattern structure(dim, dim, 5);
+      testproblem.five_point_structure(structure);
+      structure.compress();
+      SparseMatrix<double>  A(structure);
+      testproblem.five_point(A);
+      
+      for (unsigned int blocksize = 2; blocksize < 32; blocksize <<= 1)
+       {
+         deallog << "Block size " << blocksize << std::endl;
+
+         const unsigned int n_blocks = dim/blocksize;
+         BlockList blist;
+         blist.initialize(n_blocks);
+         std::vector<unsigned int> indices(blocksize);
+         
+         for (unsigned int block=0;block<n_blocks;++block)
+           {
+             for (unsigned int i=0;i<blocksize;++i)
+               indices[i] = i+block*blocksize;
+             blist.add(block, indices);
+           }
+         
+         RelaxationBlock<SparseMatrix<double>,double>::AdditionalData relax_data(blist, 0.8);
+         PreconditionBlock<SparseMatrix<double>,double>::AdditionalData prec_data(blocksize, 0.8);
+         
+         PreconditionBlockJacobi<SparseMatrix<double>,double> prec_jacobi;
+         prec_jacobi.initialize(A, prec_data);
+         PreconditionBlockSOR<SparseMatrix<double>,double> prec_sor;
+         prec_sor.initialize(A, prec_data);
+         PreconditionBlockSSOR<SparseMatrix<double>,double> prec_ssor;
+         prec_ssor.initialize(A, prec_data);
+
+         RelaxationBlockJacobi<SparseMatrix<double>,double> relax_jacobi;
+         relax_jacobi.initialize(A, relax_data);
+         RelaxationBlockSOR<SparseMatrix<double>,double> relax_sor;
+         relax_sor.initialize(A, relax_data);
+         RelaxationBlockSSOR<SparseMatrix<double>,double> relax_ssor;
+         relax_ssor.initialize(A, relax_data);
+         
+         Vector<double>  f(dim);
+         Vector<double>  u(dim);
+         Vector<double> res(dim);
+         
+         f = 1.;
+         u = 1.;
+         
+         try
+           {
+             double r1, r2;
+
+             deallog.push("Jacobi");
+             r1 = check_solve(rich,A,u,f,prec_jacobi);
+             r2 = check_solve(relax,A,u,f,prec_jacobi);
+             deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+             r2 = check_solve(relax,A,u,f,relax_jacobi);
+             deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+             deallog.pop();
+             
+             deallog.push("SOR   ");
+             r1 = check_solve(rich,A,u,f,prec_sor);
+             r2 = check_solve(relax,A,u,f,prec_sor);
+             deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+             r2 = check_solve(relax,A,u,f,relax_sor);
+             deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+             deallog.pop();
+             
+             deallog.push("SSOR  ");
+             r1 = check_solve(rich,A,u,f,prec_ssor);
+             r2 = check_solve(relax,A,u,f,prec_ssor);
+             deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+             r2 = check_solve(relax,A,u,f,relax_ssor);
+             deallog << "diff " << std::fabs(r1-r2)/r1 << std::endl;
+             deallog.pop();
+           }
+         catch (std::exception& e)
+           {
+             std::cerr << "Exception: " << e.what() << std::endl;
+           }
+       }
+    }  
+}
+
diff --git a/tests/lac/solver_relaxation_02/cmp/generic b/tests/lac/solver_relaxation_02/cmp/generic
new file mode 100644 (file)
index 0000000..1f0ca05
--- /dev/null
@@ -0,0 +1,102 @@
+
+DEAL::Size 33 Unknowns 1024
+DEAL::Block size 2
+DEAL:Jacobi:Richardson::Starting value 32.00
+DEAL:Jacobi:Richardson::Failure step 100 value 16.51
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 16.51
+DEAL:Jacobi::diff 0
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 16.51
+DEAL:Jacobi::diff 0
+DEAL:SOR   :Richardson::Starting value 32.00
+DEAL:SOR   :Richardson::Failure step 100 value 12.27
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 12.27
+DEAL:SOR   ::diff 0
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 12.27
+DEAL:SOR   ::diff 0
+DEAL:SSOR  :Richardson::Starting value 32.00
+DEAL:SSOR  :Richardson::Failure step 100 value 5.660
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 5.660
+DEAL:SSOR  ::diff 0
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 5.660
+DEAL:SSOR  ::diff 0
+DEAL::Block size 4
+DEAL:Jacobi:Richardson::Starting value 32.00
+DEAL:Jacobi:Richardson::Failure step 100 value 15.27
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 15.27
+DEAL:Jacobi::diff 0
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 15.27
+DEAL:Jacobi::diff 0
+DEAL:SOR   :Richardson::Starting value 32.00
+DEAL:SOR   :Richardson::Failure step 100 value 10.60
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 10.60
+DEAL:SOR   ::diff 0
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 10.60
+DEAL:SOR   ::diff 0
+DEAL:SSOR  :Richardson::Starting value 32.00
+DEAL:SSOR  :Richardson::Failure step 100 value 4.196
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 4.196
+DEAL:SSOR  ::diff 0
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 4.196
+DEAL:SSOR  ::diff 0
+DEAL::Block size 8
+DEAL:Jacobi:Richardson::Starting value 32.00
+DEAL:Jacobi:Richardson::Failure step 100 value 14.29
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 14.29
+DEAL:Jacobi::diff 0
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 14.29
+DEAL:Jacobi::diff 0
+DEAL:SOR   :Richardson::Starting value 32.00
+DEAL:SOR   :Richardson::Failure step 100 value 9.455
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 9.455
+DEAL:SOR   ::diff 0
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 9.455
+DEAL:SOR   ::diff 0
+DEAL:SSOR  :Richardson::Starting value 32.00
+DEAL:SSOR  :Richardson::Failure step 100 value 3.335
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 3.335
+DEAL:SSOR  ::diff 0
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 3.335
+DEAL:SSOR  ::diff 0
+DEAL::Block size 16
+DEAL:Jacobi:Richardson::Starting value 32.00
+DEAL:Jacobi:Richardson::Failure step 100 value 13.60
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 13.60
+DEAL:Jacobi::diff 0
+DEAL:Jacobi:Relaxation::Starting value 32.00
+DEAL:Jacobi:Relaxation::Failure step 100 value 13.60
+DEAL:Jacobi::diff 0
+DEAL:SOR   :Richardson::Starting value 32.00
+DEAL:SOR   :Richardson::Failure step 100 value 8.752
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 8.752
+DEAL:SOR   ::diff 0
+DEAL:SOR   :Relaxation::Starting value 32.00
+DEAL:SOR   :Relaxation::Failure step 100 value 8.752
+DEAL:SOR   ::diff 0
+DEAL:SSOR  :Richardson::Starting value 32.00
+DEAL:SSOR  :Richardson::Failure step 100 value 2.874
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 2.874
+DEAL:SSOR  ::diff 0
+DEAL:SSOR  :Relaxation::Starting value 32.00
+DEAL:SSOR  :Relaxation::Failure step 100 value 2.874
+DEAL:SSOR  ::diff 0

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.