]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
MG works on locally refined meshes, needs improvement and cleaning, error in tests...
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 28 Jun 1999 15:19:27 +0000 (15:19 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 28 Jun 1999 15:19:27 +0000 (15:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@1487 0785d39b-7218-0410-832d-ea1e28bc413d

tests/deal.II/helmholtz1.th
tests/deal.II/helmholtz1mg.th
tests/deal.II/mg.cc
tests/lac/Makefile
tests/lac/mg.cc
tests/lac/mgbase.cc
tests/lac/solver.cc
tests/lac/testmatrix.cc
tests/lac/testmatrix.h

index 5abed88b92692b4919da256db13725eb4ad1720e..99fc58809a3eccae11ea4bdb5dcc1b21ed6f002f 100644 (file)
@@ -50,7 +50,7 @@ Helmholtz::build_all(MATRIX& A,
                  double u = fe.shape_value(j,k);
                  
                  elmat(i,j) += fe.JxW(k) *
-                               (1.*u*v + du*dv)
+                               (0.1*u*v + du*dv)
                                ;
                }
            }
index 450e0a4b5ed1e0bcf60da9b223743e68bf4d53f7..11d952711ee9dae06753a5eee54ac02658b06c5d 100644 (file)
@@ -46,7 +46,7 @@ Helmholtz::build_mgmatrix(MGMatrix<MATRIX>& A,
                  double u = fe.shape_value(j,k);
                  
                  elmat(i,j) += fe.JxW(k) *
-                               (1.*u*v + du*dv)
+                               (0.1*u*v + du*dv)
                                ;
                }
            }
index 356fdfbfd39a95a6ef599fe4317799867c3a2dd0..d1c38f74fd62ed160151fc687669fae91cb3ac7e 100644 (file)
@@ -41,13 +41,13 @@ class MGSmootherLAC
   public MGSmootherBase
 {
   private:
-    SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+    SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
   public:
-    MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+    MGSmootherLAC(MGMatrix<SparseMatrix<double> >&);
     
     virtual void smooth (const unsigned int level,
-                        Vector<float> &u,
-                        const Vector<float> &rhs) const;
+                        Vector<double> &u,
+                        const Vector<double> &rhs) const;
     
 };
 
@@ -107,7 +107,7 @@ main()
          
          u = 0.;
          vector<SparseMatrixStruct> mgstruct(tr.n_levels());
-         MGMatrix<SparseMatrix<float> > mgA(0,tr.n_levels()-1);
+         MGMatrix<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
          for (unsigned int i=0;i<tr.n_levels();++i)
            {
              mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
@@ -120,11 +120,11 @@ main()
          equation.build_mgmatrix(mgA, mgdof, quadrature);
          
          SolverControl cgcontrol(10,0., false, false);
-         PrimitiveVectorMemory<Vector<float> > cgmem;
-         SolverCG<SparseMatrix<float>, Vector<float> > cgsolver(cgcontrol, cgmem);
-         PreconditionIdentity<Vector<float> > cgprec;
-         MGCoarseGridLACIteration<SolverCG<SparseMatrix<float>, Vector<float> >,
-           SparseMatrix<float>, PreconditionIdentity<Vector<float> > >
+         PrimitiveVectorMemory<Vector<double> > cgmem;
+         SolverCG<SparseMatrix<double>, Vector<double> > cgsolver(cgcontrol, cgmem);
+         PreconditionIdentity<Vector<double> > cgprec;
+         MGCoarseGridLACIteration<SolverCG<SparseMatrix<double>, Vector<double> >,
+           SparseMatrix<double>, PreconditionIdentity<Vector<double> > >
            coarse(cgsolver, mgA[0], cgprec);
          
          MGSmootherLAC smoother(mgA);
@@ -149,21 +149,21 @@ RHSFunction<dim>::operator() (const Point<dim>&) const
   return 1.;
 }
 
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<double> >& matrix)
                :
                matrices(&matrix)
 {}
 
 void
 MGSmootherLAC::smooth (const unsigned int level,
-                      Vector<float> &u,
-                      const Vector<float> &rhs) const
+                      Vector<double> &u,
+                      const Vector<double> &rhs) const
 {
-  SolverControl control(1,1.e-300,false,false);
-  PrimitiveVectorMemory<Vector<float> > mem;
-  SolverRichardson<SparseMatrix<float> , Vector<float>  > rich(control, mem);
-  PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
-    prec((*matrices)[level], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.);
+  SolverControl control(2,1.e-300,false,false);
+  PrimitiveVectorMemory<Vector<double> > mem;
+  SolverRichardson<SparseMatrix<double> , Vector<double>  > rich(control, mem);
+  PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+    prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
 
   rich.solve((*matrices)[level], u, rhs, prec);
 }
index ac8fcdb906deeb0e712153c0a3b7da7e6cedbf93..138b42e782f03d949db081aa3c67777a26d45baa 100644 (file)
@@ -148,7 +148,7 @@ veryclean: clean
 # Automatic generation of dependencies
 ############################################################
 
-all-cc-files = $(solver-cc-files) $(mgbase-cc-files) $(mg-cc-files) $(benchmark-cc-files)
+all-cc-files = $(wildcard *.cc)
 
 Make.depend: $(all-cc-files)
        @echo =====Dependecies=== Make.depend
index c667839d0b4aa1865a64debb7bae61f6eb4f4f07..0cba014f426d49b90c1062eb0b4f64cbf4dc5433 100644 (file)
@@ -18,7 +18,7 @@
 #include <lac/mgbase.h>
 
 #define TYPE  long double
-#define ACCURACY 1.e-18
+#define ACCURACY 1.e-20
 
 template<class VECTOR>
 void print_vector(ostream& s, const VECTOR& v)
@@ -44,16 +44,16 @@ class FDMG
                                     /**
                                      * Pointer to the level matrices.
                                      */
-    SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+    SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
   public:
-    FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
+    FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<double> >& matrices,
         FDMGTransfer& transfer);
     
     
-    virtual void level_residual (unsigned int level,
-                                Vector<float>& dst,
-                                const Vector<float>& src,
-                                const Vector<float>& rhs);
+    virtual void level_vmult (unsigned int level,
+                                Vector<double>& dst,
+                                const Vector<double>& src,
+                                const Vector<double>& rhs);
 
     void copy_to_mg(const Vector<TYPE>& rhs);
     void copy_from_mg(Vector<TYPE>& lsg);
@@ -64,19 +64,19 @@ class MGSmootherLAC
   public MGSmootherBase
 {
   private:
-    SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+    SmartPointer<MGMatrix<SparseMatrix<double> > >matrices;
   public:
-    MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+    MGSmootherLAC(MGMatrix<SparseMatrix<double> >&);
     
     virtual void smooth (const unsigned int level,
-                        Vector<float> &u,
-                        const Vector<float> &rhs) const;
+                        Vector<double> &u,
+                        const Vector<double> &rhs) const;
     
 };
 
-typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<float> , Vector<float>  >,
-SparseMatrix<float>, /*PreconditionRelaxation<SparseMatrix<float> ,*/
-  PreconditionIdentity<Vector<float> > >
+typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<double> , Vector<double>  >,
+SparseMatrix<double>, /*PreconditionRelaxation<SparseMatrix<double> ,*/
+  PreconditionIdentity<Vector<double> > >
 Coarse;
 
 
@@ -86,7 +86,7 @@ main()
   deallog.attach(logfile);
   
   PrimitiveVectorMemory<Vector<TYPE>  > mem;
-  SolverControl control(100, ACCURACY, true);
+  SolverControl control(10000, ACCURACY);
 
   const unsigned int base = 3;
   const unsigned int maxlevel = 8;
@@ -97,9 +97,9 @@ main()
   FDMGTransfer transfer(maxsize, maxsize, maxlevel);
 
                                   // coarse grid solver
-  PrimitiveVectorMemory<Vector<float> > cgmem;
+  PrimitiveVectorMemory<Vector<double> > cgmem;
   ReductionControl cgcontrol(100, 1.e-30, 1.e-2, false, false);
-  SolverCG<SparseMatrix<float> , Vector<float> > cgcg(cgcontrol,cgmem);
+  SolverCG<SparseMatrix<double> , Vector<double> > cgcg(cgcontrol,cgmem);
 
   
   for (unsigned int level = 0; level <= maxlevel; ++level)
@@ -113,7 +113,7 @@ main()
 
                                       // Make matrix
       vector<SparseMatrixStruct >  structure(maxlevel+1);
-      MGMatrix<SparseMatrix<float> > A(minlevel,maxlevel);
+      MGMatrix<SparseMatrix<double> > A(minlevel,maxlevel);
 
       FDMatrix testproblem(size, size);
 
@@ -131,9 +131,9 @@ main()
       
       FDMG multigrid(level, A, transfer);
 
-//      PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
-      PreconditionIdentity<Vector<float> >
-       cgprec;//(A[minlevel], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+//      PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+      PreconditionIdentity<Vector<double> >
+       cgprec;//(A[minlevel], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.2);
       
       Coarse coarsegrid(cgcg, A[minlevel], cgprec);
       
@@ -143,8 +143,8 @@ main()
       PreconditionMG<FDMG, Vector<TYPE> >
        precondition(multigrid, smoother, smoother, coarsegrid);
 
-//      SolverRichardson<SparseMatrix<float> , Vector<TYPE> > solver(control, mem);
-      SolverCG<SparseMatrix<float> , Vector<TYPE> > solver(control, mem);
+//      SolverRichardson<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
+      SolverCG<SparseMatrix<double> , Vector<TYPE> > solver(control, mem);
 
       Vector<TYPE> u(dim);
       Vector<TYPE> f(dim);
@@ -157,7 +157,7 @@ main()
     }
 }
 
-FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
+FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<double> >& matrices,
           FDMGTransfer& transfer)
                :
                MGBase(transfer, 0, maxlevel),
@@ -172,12 +172,12 @@ FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
 
 
 void
-FDMG::level_residual (unsigned int level,
-                     Vector<float>& dst,
-                     const Vector<float>& src,
-                     const Vector<float>& rhs)
+FDMG::level_vmult (unsigned int level,
+                     Vector<double>& dst,
+                     const Vector<double>& src,
+                     const Vector<double>&)
 {
-  (*matrices)[level].residual(dst, src, rhs);
+  (*matrices)[level].vmult(dst, src);
 }
 
 void
@@ -193,7 +193,7 @@ FDMG::copy_from_mg(Vector<TYPE>& v)
 }
 
 
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<double> >& matrix)
                :
                matrices(&matrix)
 {}
@@ -201,14 +201,14 @@ MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
 
 void
 MGSmootherLAC::smooth (const unsigned int level,
-                      Vector<float> &u,
-                      const Vector<float> &rhs) const
+                      Vector<double> &u,
+                      const Vector<double> &rhs) const
 {
   SolverControl control(1,1.e-300,false,false);
-  PrimitiveVectorMemory<Vector<float> > mem;
-  SolverRichardson<SparseMatrix<float> , Vector<float>  > rich(control, mem);
-  PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
-    prec((*matrices)[level], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.);
+  PrimitiveVectorMemory<Vector<double> > mem;
+  SolverRichardson<SparseMatrix<double> , Vector<double>  > rich(control, mem);
+  PreconditionRelaxation<SparseMatrix<double> , Vector<double> >
+    prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
 
   rich.solve((*matrices)[level], u, rhs, prec);
 }
index 59d7287ba2ccbac9fa139da07bb6b1bf683c1073..8a47cd22a02845be2f06e086f9156a14f6b86e87 100644 (file)
@@ -10,11 +10,11 @@ class TransferTest
   public MGTransferBase
 {
     virtual void prolongate(unsigned l,
-                           Vector<float>& dst,
-                           const Vector<float>& src) const;
+                           Vector<double>& dst,
+                           const Vector<double>& src) const;
     virtual void restrict(unsigned l,
-                         Vector<float>& dst,
-                         const Vector<float>& src) const;
+                         Vector<double>& dst,
+                         const Vector<double>& src) const;
     friend class MGBase;
 };
 
@@ -24,8 +24,8 @@ class SmoothTest
 {
   public:
         virtual void smooth (const unsigned int  level,
-                            Vector<float>      &u,
-                            const Vector<float>& rhs) const;
+                            Vector<double>      &u,
+                            const Vector<double>& rhs) const;
 };
 
 class MGCoarseGridTest
@@ -34,8 +34,8 @@ class MGCoarseGridTest
 {
   public:
     virtual void operator () (unsigned int level,
-                             Vector<float>& dst,
-                             const Vector<float>& src) const ;
+                             Vector<double>& dst,
+                             const Vector<double>& src) const ;
 };
 
 
@@ -48,7 +48,7 @@ class MGTest
 
   public:
     MGTest(TransferTest& t, SmoothTest& sm, MGCoarseGridTest& cs)
-                   : MGBase(t, 9, 3),
+                   : MGBase(t, 3, 9),
                      smoother(sm),
                      solver(cs)
       {
@@ -59,10 +59,10 @@ class MGTest
          }
       }
     
-    virtual void level_residual(unsigned level,
-                               Vector<float>& dst,
-                               const Vector<float>& src,
-                               const Vector<float>& rhs);
+    virtual void level_vmult(unsigned level,
+                               Vector<double>& dst,
+                               const Vector<double>& src,
+                               const Vector<double>& rhs);
     
     void doit()
       {
@@ -83,41 +83,41 @@ main()
 
 void
 TransferTest::prolongate(unsigned l,
-                        Vector<float>&,
-                        const Vector<float>&) const
+                        Vector<double>&,
+                        const Vector<double>&) const
 {
   cout << "Prolongating " << l-1 << " to " << l << endl;
 }
 
 void
 TransferTest::restrict(unsigned l,
-                      Vector<float>&,
-                      const Vector<float>&) const
+                      Vector<double>&,
+                      const Vector<double>&) const
 {
   cout << "Restricting  " << l << " to " << l-1 << endl;
 }
 
 void
 SmoothTest::smooth (const unsigned int  level,
-                   Vector<float>      &,
-                   const Vector<float>&) const
+                   Vector<double>      &,
+                   const Vector<double>&) const
 {
   cout << "Smoothing on " << level << endl;
 }
 
 void
-MGTest::level_residual(unsigned l,
-                      Vector<float>&,
-                      const Vector<float>&,
-                      const Vector<float>&)
+MGTest::level_vmult(unsigned l,
+                      Vector<double>&,
+                      const Vector<double>&,
+                      const Vector<double>&)
 {
   cout << "Residual on  " << l << endl;
 }
 
 void
 MGCoarseGridTest::operator() (unsigned int level,
-                             Vector<float>&,
-                             const Vector<float>&) const
+                             Vector<double>&,
+                             const Vector<double>&) const
 {
   cout << "Solving on   " << level << endl;
 }
index 62d6bcd8420f5c23ab886780b3b6e5a0dfc72954..b4f7d9a096e75c178fae01d1bb9c7ad315a7bee2 100644 (file)
@@ -4,6 +4,7 @@
 
 #include <cmath>
 #include <fstream>
+#include <iomanip>
 #include "testmatrix.h"
 #include <base/logstream.h>
 #include <lac/sparsematrix.h>
@@ -32,13 +33,13 @@ main()
   deallog.attach(logfile);
   
   PrimitiveVectorMemory<Vector<double>  > mem;
-  SolverControl control(100, 1.e-5, true);
+  SolverControl control(50, 1.e-5, true);
   SolverCG<SparseMatrix<float> , Vector<double>  > cg(control, mem);
   SolverGMRES<SparseMatrix<float> , Vector<double>  > gmres(control, mem,20);
   SolverBicgstab<SparseMatrix<float> , Vector<double>  > bicgstab(control, mem);
   SolverRichardson<SparseMatrix<float> , Vector<double>  > rich(control, mem);
 
-  for (unsigned int size=10; size <= 10; size *= 10)
+  for (unsigned int size=4; size <= 40; size *= 3)
     {
       deallog << "Size " << size << endl;
       
@@ -61,6 +62,18 @@ main()
       
       Vector<double>  f(dim);
       Vector<double>  u(dim);
+      Vector<double> res(dim);
+
+      f = 1.;
+      u = 1.;
+      
+      A.residual(res,u,f);
+      A.SOR(res);
+      res.add(1.,u);
+      A.SOR_step(u,f);
+      res.add(-1.,u);
+    
+      deallog << "SOR-diff:" << res*res << endl;
       
       deallog.push("no");
 
index 5100cbe8d786ed8c7840bf19891e9e7328364c59..ac7840fb1a0280014f5fccac058b7e03fd121db6 100644 (file)
@@ -118,7 +118,7 @@ FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny,
 
 void
 FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
-                          SparseMatrixStruct& structure, SparseMatrix<float>& matrix)
+                          SparseMatrixStruct& structure, SparseMatrix<double>& matrix)
 {
   structure.reinit((nx-1)*(ny-1),(2*nx-1)*(2*ny-1),9);
   
@@ -210,8 +210,8 @@ FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
 
 void
 FDMGTransfer::prolongate (const unsigned int   to_level,
-                         Vector<float>       &dst,
-                         const Vector<float> &src) const
+                         Vector<double>       &dst,
+                         const Vector<double> &src) const
 {
   Assert((to_level>0) && (to_level<=matrices.size()),
         ExcIndexRange(to_level, 0, matrices.size()+1));
@@ -222,8 +222,8 @@ FDMGTransfer::prolongate (const unsigned int   to_level,
 
 void
 FDMGTransfer::restrict (const unsigned int   from_level,
-                       Vector<float>       &dst,
-                       const Vector<float> &src) const
+                       Vector<double>       &dst,
+                       const Vector<double> &src) const
 {
   Assert((from_level>0) && (from_level<=matrices.size()),
         ExcIndexRange(from_level, 0, matrices.size()+1));
index 8099555feea2ffbcbcc21cce71b5a9f44fc5cefa..13e9ebf7beda7b833e2d889f6b85e4960e91f19a 100644 (file)
@@ -66,16 +66,16 @@ class FDMGTransfer
                                      * function in #MGTranferBase#.
                                      */
     virtual void prolongate (const unsigned int   to_level,
-                            Vector<float>       &dst,
-                            const Vector<float> &src) const;
+                            Vector<double>       &dst,
+                            const Vector<double> &src) const;
 
                                     /**
                                      * Implementation of abstract
                                      * function in #MGTranferBase#.
                                      */
     virtual void restrict (const unsigned int   from_level,
-                          Vector<float>       &dst,
-                          const Vector<float> &src) const;
+                          Vector<double>       &dst,
+                          const Vector<double> &src) const;
 
                                     /**
                                      * Exception.
@@ -91,7 +91,7 @@ class FDMGTransfer
                                     /**
                                      * Prolongation matrices.
                                      */
-    vector<SparseMatrix<float> > matrices;
+    vector<SparseMatrix<double> > matrices;
 
                                     /**
                                      * Matrix generator.
@@ -101,5 +101,5 @@ class FDMGTransfer
                                      * fine to coarse (#vmult#).
                                      */
     void build_matrix(unsigned int nx, unsigned int ny,
-                     SparseMatrixStruct&, SparseMatrix<float>&);
+                     SparseMatrixStruct&, SparseMatrix<double>&);
 };

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.