From: guido Date: Tue, 4 May 1999 15:29:14 +0000 (+0000) Subject: The First MultiGrid Test works X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4d22b35cb401b0112b23feccde815eb0034ab96f;p=dealii-svn.git The First MultiGrid Test works git-svn-id: https://svn.dealii.org/trunk@1265 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index e0e83596d5..c30ef1930c 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -11,6 +11,8 @@ #include +#define FLOAT float + /** * Abstract base class for coarse grid solvers. @@ -31,8 +33,8 @@ class MGCoarseGridSolver * about the matrix is removed to * that class. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const = 0; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const = 0; }; @@ -65,9 +67,27 @@ class MGSmootherBase * things. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const = 0; + Vector &u, + const Vector &rhs) const = 0; + +}; +/** + * Smoother doing nothing. + * @author Guido Kanschat, 1999 + */ +class MGSmootherIdentity + : + public MGSmootherBase +{ + public: + /** + * Implementation of the interface in #MGSmootherBase#. + * This function does nothing. + */ + virtual void smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const; }; /** @@ -100,8 +120,8 @@ class MGTransferBase : public Subscriptor * level. */ virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; /** * Restrict a vector from level @@ -117,10 +137,44 @@ class MGTransferBase : public Subscriptor * level. */ virtual void restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; +}; + +/** + * An array with a vector for each level. + * The purpose of this class is mostly to provide range checking that is missing in + * those ultra-fashionable STL vectors. + * @author Guido Kanschat, 1999 + */ +template +class MGVector + : + public Subscriptor, + public vector +{ + public: + /** + * Constructor allowing to initialize the number of levels. + */ + MGVector(unsigned int minlevel, unsigned int maxlevel); + /** + * Safe access operator. + */ + VECTOR& operator[](unsigned int); + /** + * Safe access operator. + */ + const VECTOR& operator[](unsigned int) const; + private: + /** + * Level of first component. + */ + unsigned int minlevel; }; + + /** * An array of matrices for each level. * This class provides the functionality of #vector# combined @@ -137,7 +191,21 @@ class MGMatrix /** * Constructor allowing to initialize the number of levels. */ - MGMatrix(unsigned int n_levels); + MGMatrix(unsigned int minlevel, unsigned int maxlevel); + /** + * Safe access operator. + */ + MATRIX& operator[](unsigned int); + /** + * Safe access operator. + */ + const MATRIX& operator[](unsigned int) const; + + private: + /** + * Level of first component. + */ + unsigned int minlevel; }; @@ -168,8 +236,8 @@ class MGCoarseGridLACIteration * matrix, vectors and * preconditioner. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const; private: /** * Reference to the solver. @@ -207,23 +275,35 @@ class MGCoarseGridLACIteration * */ class MGBase { + private: MGBase(const MGBase&); const MGBase& operator=(const MGBase&); - + + protected: + /** + * Highest level of cells. + */ + unsigned int maxlevel; + + /** + * Level for coarse grid solution. + */ + unsigned int minlevel; /** * Auxiliary vector, defect. */ - vector > d; + MGVector > d; /** * Auxiliary vector, solution. */ - vector > s; + MGVector > s; + private: /** * Auxiliary vector. */ - Vector t; + Vector t; /** * Prolongation and restriction object. @@ -231,15 +311,6 @@ class MGBase SmartPointer transfer; protected: - /** - * Highest level of cells. - */ - unsigned int maxlevel; - - /** - * Level for coarse grid solution. - */ - unsigned int minlevel; /** * The actual v-cycle multigrid method. @@ -262,9 +333,9 @@ class MGBase * derived class. */ virtual void level_residual(unsigned int level, - Vector& dst, - const Vector& src, - const Vector& rhs) = 0; + Vector& dst, + const Vector& src, + const Vector& rhs) = 0; public: @@ -273,8 +344,27 @@ class MGBase */ MGBase(const MGTransferBase& transfer, unsigned int maxlevel, unsigned int minlevel); + /** + * Virtual destructor. + */ virtual ~MGBase(); - + /** + * Execute v-cycle algorithm. + * This function assumes, that + * the vector #d# is properly + * filled with the residual in + * the outer defect correction + * scheme. After execution of + * #vcycle()#, the result is in + * the vector #s#. We propose to + * write functions #copy_from_mg# + * and #copy_to_mg# in derived + * classes of #MGBase# to access + * #d# and #s#. + */ + void vcycle(const MGSmootherBase& pre_smooth, + const MGSmootherBase& post_smooth, + const MGCoarseGridSolver& cgs); }; @@ -288,16 +378,64 @@ MGCoarseGridLACIteration template void MGCoarseGridLACIteration::operator() -(unsigned int, Vector& dst, const Vector& src) const +(unsigned int, Vector& dst, const Vector& src) const { solver.solve(matrix, dst, src, precondition); } +////////////////////////////////////////////////////////////////////// + +template +MGVector::MGVector(unsigned int min, unsigned int max) + : + vector(max-min+1), + minlevel(min) +{} + + +template +VECTOR& +MGVector::operator[](unsigned int i) +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + + +template +const VECTOR& +MGVector::operator[](unsigned int i) const +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + +////////////////////////////////////////////////////////////////////// + template -MGMatrix::MGMatrix(unsigned int nlevels) +MGMatrix::MGMatrix(unsigned int min, unsigned int max) : - vector(nlevels) + vector(max-min+1), + minlevel(min) {} +template +MATRIX& +MGMatrix::operator[](unsigned int i) +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + + +template +const MATRIX& +MGMatrix::operator[](unsigned int i) const +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + + #endif diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index 1016d092c5..fb8fa27b4a 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -110,29 +110,6 @@ class MGSmoother vector > interior_boundary_dofs; }; -/** - * Smoother doing nothing. - * @author Guido Kanschat, 1999 - */ -class MGSmootherIdentity - : - public MGSmoother -{ - public: - /** - * Constructor. - */ - template - MGSmootherIdentity(const MGDoFHandler &mg_dof); - /** - * Implementation of the interface in #MGSmootherBase#. - * This function does nothing. - */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; -}; - /** * Implementation of an SOR smoother. * diff --git a/deal.II/deal.II/include/numerics/mg_smoother.h b/deal.II/deal.II/include/numerics/mg_smoother.h index 1016d092c5..fb8fa27b4a 100644 --- a/deal.II/deal.II/include/numerics/mg_smoother.h +++ b/deal.II/deal.II/include/numerics/mg_smoother.h @@ -110,29 +110,6 @@ class MGSmoother vector > interior_boundary_dofs; }; -/** - * Smoother doing nothing. - * @author Guido Kanschat, 1999 - */ -class MGSmootherIdentity - : - public MGSmoother -{ - public: - /** - * Constructor. - */ - template - MGSmootherIdentity(const MGDoFHandler &mg_dof); - /** - * Implementation of the interface in #MGSmootherBase#. - * This function does nothing. - */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; -}; - /** * Implementation of an SOR smoother. * diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index c6252fa0d8..73a80470fb 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -1,6 +1,24 @@ // $Id$ #include +#include +#include + +template +void print_vector(ostream& s, const VECTOR& v) +{ + const unsigned int n = (unsigned int)(sqrt(v.size())+.3); + unsigned int k=0; + + for (unsigned int i=0;irestrict(level, t, d[level-1]); + transfer->restrict(level, d[level-1], t); + // do recursion level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver); // do coarse grid correction - transfer->prolongate(level, s[level], s[level-1]); - + t.reinit(s[level].size()); + transfer->prolongate(level, t, s[level-1]); + s[level].add(t); + // smoothing (modify s again) post_smooth.smooth(level, s[level], d[level]); } + +////////////////////////////////////////////////////////////////////// + MGTransferBase::~MGTransferBase() {} MGSmootherBase::~MGSmootherBase() {} + + +////////////////////////////////////////////////////////////////////// + + +void +MGSmootherIdentity::smooth (const unsigned int, + Vector &, + const Vector &) const +{} + + diff --git a/deal.II/deal.II/source/multigrid/mg_smoother.cc b/deal.II/deal.II/source/multigrid/mg_smoother.cc index cc16654c94..71b80cb4b3 100644 --- a/deal.II/deal.II/source/multigrid/mg_smoother.cc +++ b/deal.II/deal.II/source/multigrid/mg_smoother.cc @@ -107,22 +107,6 @@ MGSmoother::set_zero_interior_boundary (const unsigned int level, u(*p) = 0; }; -////////////////////////////////////////////////////////////////////// - - -template -MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler &mg_dof) - : - MGSmoother(mg_dof, 0) -{} - -void -MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const -{} - - ////////////////////////////////////////////////////////////////////// template @@ -153,5 +137,3 @@ template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler &mg_ const MGMatrix >& matrix, unsigned int steps); -template MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler &mg_dof); - diff --git a/deal.II/deal.II/source/numerics/mg_smoother.cc b/deal.II/deal.II/source/numerics/mg_smoother.cc index cc16654c94..71b80cb4b3 100644 --- a/deal.II/deal.II/source/numerics/mg_smoother.cc +++ b/deal.II/deal.II/source/numerics/mg_smoother.cc @@ -107,22 +107,6 @@ MGSmoother::set_zero_interior_boundary (const unsigned int level, u(*p) = 0; }; -////////////////////////////////////////////////////////////////////// - - -template -MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler &mg_dof) - : - MGSmoother(mg_dof, 0) -{} - -void -MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const -{} - - ////////////////////////////////////////////////////////////////////// template @@ -153,5 +137,3 @@ template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler &mg_ const MGMatrix >& matrix, unsigned int steps); -template MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler &mg_dof); - diff --git a/deal.II/lac/include/lac/mgbase.h b/deal.II/lac/include/lac/mgbase.h index 3cb773a6b2..54c15ce614 100644 --- a/deal.II/lac/include/lac/mgbase.h +++ b/deal.II/lac/include/lac/mgbase.h @@ -11,6 +11,8 @@ #include +#define FLOAT float + /** * Abstract base class for coarse grid solvers. @@ -31,8 +33,8 @@ class MGCoarseGridSolver * about the matrix is removed to * that class. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const = 0; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const = 0; }; @@ -65,9 +67,27 @@ class MGSmootherBase * things. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const = 0; + Vector &u, + const Vector &rhs) const = 0; + +}; +/** + * Smoother doing nothing. + * @author Guido Kanschat, 1999 + */ +class MGSmootherIdentity + : + public MGSmootherBase +{ + public: + /** + * Implementation of the interface in #MGSmootherBase#. + * This function does nothing. + */ + virtual void smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const; }; /** @@ -100,8 +120,8 @@ class MGTransferBase : public Subscriptor * level. */ virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; /** * Restrict a vector from level @@ -117,10 +137,44 @@ class MGTransferBase : public Subscriptor * level. */ virtual void restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; +}; + +/** + * An array with a vector for each level. + * The purpose of this class is mostly to provide range checking that is missing in + * those ultra-fashionable STL vectors. + * @author Guido Kanschat, 1999 + */ +template +class MGVector + : + public Subscriptor, + public vector +{ + public: + /** + * Constructor allowing to initialize the number of levels. + */ + MGVector(unsigned int minlevel, unsigned int maxlevel); + /** + * Safe access operator. + */ + VECTOR& operator[](unsigned int); + /** + * Safe access operator. + */ + const VECTOR& operator[](unsigned int) const; + private: + /** + * Level of first component. + */ + unsigned int minlevel; }; + + /** * An array of matrices for each level. * This class provides the functionality of #vector# combined @@ -137,7 +191,21 @@ class MGMatrix /** * Constructor allowing to initialize the number of levels. */ - MGMatrix(unsigned int n_levels); + MGMatrix(unsigned int minlevel, unsigned int maxlevel); + /** + * Safe access operator. + */ + MATRIX& operator[](unsigned int); + /** + * Safe access operator. + */ + const MATRIX& operator[](unsigned int) const; + + private: + /** + * Level of first component. + */ + unsigned int minlevel; }; @@ -168,8 +236,8 @@ class MGCoarseGridLACIteration * matrix, vectors and * preconditioner. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const; private: /** * Reference to the solver. @@ -207,23 +275,35 @@ class MGCoarseGridLACIteration * */ class MGBase { + private: MGBase(const MGBase&); const MGBase& operator=(const MGBase&); - + + protected: + /** + * Highest level of cells. + */ + unsigned int maxlevel; + + /** + * Level for coarse grid solution. + */ + unsigned int minlevel; /** * Auxiliary vector, defect. */ - vector > d; + MGVector > d; /** * Auxiliary vector, solution. */ - vector > s; + MGVector > s; + private: /** * Auxiliary vector. */ - Vector t; + Vector t; /** * Prolongation and restriction object. @@ -231,15 +311,6 @@ class MGBase SmartPointer transfer; protected: - /** - * Highest level of cells. - */ - unsigned int maxlevel; - - /** - * Level for coarse grid solution. - */ - unsigned int minlevel; /** * The actual v-cycle multigrid method. @@ -262,9 +333,9 @@ class MGBase * derived class. */ virtual void level_residual(unsigned int level, - Vector& dst, - const Vector& src, - const Vector& rhs) = 0; + Vector& dst, + const Vector& src, + const Vector& rhs) = 0; public: @@ -273,8 +344,27 @@ class MGBase */ MGBase(const MGTransferBase& transfer, unsigned int maxlevel, unsigned int minlevel); + /** + * Virtual destructor. + */ virtual ~MGBase(); - + /** + * Execute v-cycle algorithm. + * This function assumes, that + * the vector #d# is properly + * filled with the residual in + * the outer defect correction + * scheme. After execution of + * #vcycle()#, the result is in + * the vector #s#. We propose to + * write functions #copy_from_mg# + * and #copy_to_mg# in derived + * classes of #MGBase# to access + * #d# and #s#. + */ + void vcycle(const MGSmootherBase& pre_smooth, + const MGSmootherBase& post_smooth, + const MGCoarseGridSolver& cgs); }; @@ -288,16 +378,64 @@ MGCoarseGridLACIteration template void MGCoarseGridLACIteration::operator() -(unsigned int, Vector& dst, const Vector& src) const +(unsigned int, Vector& dst, const Vector& src) const { solver.solve(matrix, dst, src, precondition); } +////////////////////////////////////////////////////////////////////// + +template +MGVector::MGVector(unsigned int min, unsigned int max) + : + vector(max-min+1), + minlevel(min) +{} + + +template +VECTOR& +MGVector::operator[](unsigned int i) +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + + +template +const VECTOR& +MGVector::operator[](unsigned int i) const +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + +////////////////////////////////////////////////////////////////////// + template -MGMatrix::MGMatrix(unsigned int nlevels) +MGMatrix::MGMatrix(unsigned int min, unsigned int max) : - vector(nlevels) + vector(max-min+1), + minlevel(min) {} +template +MATRIX& +MGMatrix::operator[](unsigned int i) +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + + +template +const MATRIX& +MGMatrix::operator[](unsigned int i) const +{ + Assert((i>=minlevel)&&(i::operator[](i-minlevel); +} + + #endif diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 8b3c58a7b0..a0000885a6 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -158,7 +158,7 @@ SolverCG::solve (const Matrix &A, g.add(alpha,Ad); x.add(alpha,d ); - res = sqrt(g*g); + res = g.l2_norm(); conv = control().check(it,res); if (conv) break; diff --git a/deal.II/lac/source/mgbase.cc b/deal.II/lac/source/mgbase.cc index c6252fa0d8..73a80470fb 100644 --- a/deal.II/lac/source/mgbase.cc +++ b/deal.II/lac/source/mgbase.cc @@ -1,6 +1,24 @@ // $Id$ #include +#include +#include + +template +void print_vector(ostream& s, const VECTOR& v) +{ + const unsigned int n = (unsigned int)(sqrt(v.size())+.3); + unsigned int k=0; + + for (unsigned int i=0;irestrict(level, t, d[level-1]); + transfer->restrict(level, d[level-1], t); + // do recursion level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver); // do coarse grid correction - transfer->prolongate(level, s[level], s[level-1]); - + t.reinit(s[level].size()); + transfer->prolongate(level, t, s[level-1]); + s[level].add(t); + // smoothing (modify s again) post_smooth.smooth(level, s[level], d[level]); } + +////////////////////////////////////////////////////////////////////// + MGTransferBase::~MGTransferBase() {} MGSmootherBase::~MGSmootherBase() {} + + +////////////////////////////////////////////////////////////////////// + + +void +MGSmootherIdentity::smooth (const unsigned int, + Vector &, + const Vector &) const +{} + + diff --git a/deal.II/lac/source/solver_control.cc b/deal.II/lac/source/solver_control.cc index 5fc7b74c64..85fa0e5b34 100644 --- a/deal.II/lac/source/solver_control.cc +++ b/deal.II/lac/source/solver_control.cc @@ -107,7 +107,12 @@ ReductionControl::check (const unsigned int step, }; if (check_value < reduced_tol) - return success; + { + if (log_result) + deallog << "Convergence step " << step + << " value " << check_value << endl; + return success; + } else return SolverControl::check(step, check_value); }; diff --git a/tests/lac/mg.cc b/tests/lac/mg.cc index 439d47f20d..4699550774 100644 --- a/tests/lac/mg.cc +++ b/tests/lac/mg.cc @@ -17,54 +17,73 @@ #include #include +template +void print_vector(ostream& s, const VECTOR& v) +{ + const unsigned int n = (unsigned int)(sqrt(v.size())+.3); + unsigned int k=0; + + for (unsigned int i=0;i > >matrices; + /** + * Pointer to the level matrices. + */ + SmartPointer > >matrices; public: - FDMG(unsigned int maxlevel, MGMatrix >& matrices, - FDMGTransfer& transfer) - : - MGBase(transfer, maxlevel, 0), - matrices(&matrices) - {} + FDMG(unsigned int maxlevel, MGMatrix >& matrices, + FDMGTransfer& transfer); virtual void level_residual (unsigned int level, - Vector& dst, - const Vector& src, - const Vector& rhs); + Vector& dst, + const Vector& src, + const Vector& rhs); + + void copy_to_mg(const Vector& rhs); + void copy_from_mg(Vector& lsg); }; class MGSmootherLAC : - MGSmootherBase + public MGSmootherBase { private: - SmartPointer > >matrices; + SmartPointer > >matrices; public: - MGSmootherLAC(MGMatrix >&); + MGSmootherLAC(MGMatrix >&); virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; + Vector &u, + const Vector &rhs) const; }; -typedef MGCoarseGridLACIteration , Vector >, -SparseMatrix, PreconditionRelaxation , Vector > > +typedef MGCoarseGridLACIteration , Vector >, +SparseMatrix, /*PreconditionRelaxation ,*/ + PreconditionIdentity > > Coarse; class MGPrecondition { private: FDMG& mg; - MGSmootherLAC& smooth; + MGSmootherBase& smooth; Coarse& coarse; public: - MGPrecondition(FDMG& m, MGSmootherLAC& s, Coarse& c) + MGPrecondition(FDMG& m, MGSmootherBase& s, Coarse& c) : mg(m), smooth(s), coarse(c) {} @@ -72,7 +91,11 @@ class MGPrecondition void operator () (Vector& dst, const Vector& src) const { +// print_vector(cout, src); + mg.copy_to_mg(src); + mg.vcycle(smooth,smooth,coarse); + mg.copy_from_mg(dst); } @@ -88,7 +111,7 @@ main() SolverControl control(100, 1.e-5, true); const unsigned int base = 3; - const unsigned int maxlevel = 6; + const unsigned int maxlevel = 8; const unsigned int maxsize = base * (1 << maxlevel); @@ -96,50 +119,102 @@ main() FDMGTransfer transfer(maxsize, maxsize, maxlevel); // coarse grid solver - PrimitiveVectorMemory > cgmem; - SolverControl cgcontrol(100, 1.e-12); - SolverCG , Vector > cgcg(cgcontrol,cgmem); + PrimitiveVectorMemory > cgmem; + ReductionControl cgcontrol(100, 1.e-30, 1.e-2, false, false); + SolverCG , Vector > cgcg(cgcontrol,cgmem); - vector structure(maxlevel+1); - MGMatrix > A(maxlevel+1); for (unsigned int level = 0; level <= maxlevel; ++level) { + const unsigned int minlevel = 0; + const unsigned int size = base * (1 << level); - const unsigned int dim = (size+1)*(size+1); + const unsigned int dim = (size-1)*(size-1); + deallog << "Level " << level << " size " << size << endl; // Make matrix + vector structure(maxlevel+1); + MGMatrix > A(minlevel,maxlevel); + FDMatrix testproblem(size, size); - structure[level].reinit(dim, dim, 5); - testproblem.build_structure(structure[level]); - structure[level].compress(); - A[level].reinit(structure[level]); - testproblem.laplacian(A[level]); + for(unsigned ll = minlevel; ll <= level; ++ll) + { + const unsigned int size = base * (1 << ll); + const unsigned int dim = (size-1)*(size-1); + FDMatrix testproblem(size, size); + structure[ll].reinit(dim, dim, 5); + testproblem.build_structure(structure[ll]); + structure[ll].compress(); + A[ll].reinit(structure[ll]); + testproblem.laplacian(A[ll]); + } + FDMG multigrid(level, A, transfer); - PreconditionRelaxation , Vector > - cgprec(A[0], &SparseMatrix ::template precondition_SSOR, 1.2); +// PreconditionRelaxation , Vector > + PreconditionIdentity > + cgprec;//(A[minlevel], &SparseMatrix ::template precondition_SSOR, 1.2); - Coarse coarsegrid(cgcg, A[0], cgprec); + Coarse coarsegrid(cgcg, A[minlevel], cgprec); MGSmootherLAC smoother(A); - -// MGPrecondition prec(multigrid); +// MGSmootherIdentity smoother; + + MGPrecondition precondition(multigrid, smoother, coarsegrid); + +// SolverRichardson , Vector > solver(control, mem); + SolverCG , Vector > solver(control, mem); + + Vector u(dim); + Vector f(dim); + u = 0.; + f = 1.;//(size*size); + + solver.solve(A[level], u, f, precondition); + ofstream out("T"); + testproblem.gnuplot_print(out,u); } } +FDMG::FDMG(unsigned int maxlevel, MGMatrix >& matrices, + FDMGTransfer& transfer) + : + MGBase(transfer, maxlevel, 0), + matrices(&matrices) +{ + for (unsigned int level = minlevel; level<=maxlevel ; ++level) + { + s[level].reinit(matrices[level].m()); + d[level].reinit(matrices[level].m()); + } +} + + void FDMG::level_residual (unsigned int level, - Vector& dst, - const Vector& src, - const Vector& rhs) + Vector& dst, + const Vector& src, + const Vector& rhs) { (*matrices)[level].residual(dst, src, rhs); } -MGSmootherLAC::MGSmootherLAC(MGMatrix >& matrix) +void +FDMG::copy_to_mg(const Vector& v) +{ + d[maxlevel] = v; +} + +void +FDMG::copy_from_mg(Vector& v) +{ + v = s[maxlevel]; +} + + +MGSmootherLAC::MGSmootherLAC(MGMatrix >& matrix) : matrices(&matrix) {} @@ -147,14 +222,14 @@ MGSmootherLAC::MGSmootherLAC(MGMatrix >& matrix) void MGSmootherLAC::smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const + Vector &u, + const Vector &rhs) const { - SolverControl control(0.,1); - PrimitiveVectorMemory > mem; - SolverRichardson , Vector > rich(control, mem); - PreconditionRelaxation , Vector > - prec((*matrices)[level], &SparseMatrix ::template precondition_SSOR, 1.2); + SolverControl control(1,1.e-300,false,false); + PrimitiveVectorMemory > mem; + SolverRichardson , Vector > rich(control, mem); + PreconditionRelaxation , Vector > + prec((*matrices)[level], &SparseMatrix ::template precondition_SSOR, 1.); rich.solve((*matrices)[level], u, rhs, prec); } diff --git a/tests/lac/mgbase.cc b/tests/lac/mgbase.cc index acd703a6c8..59d7287ba2 100644 --- a/tests/lac/mgbase.cc +++ b/tests/lac/mgbase.cc @@ -47,11 +47,17 @@ class MGTest MGCoarseGridSolver& solver; public: - MGTest(TransferTest& t, SmoothTest& s, MGCoarseGridTest& c) + MGTest(TransferTest& t, SmoothTest& sm, MGCoarseGridTest& cs) : MGBase(t, 9, 3), - smoother(s), - solver(c) - {} + smoother(sm), + solver(cs) + { + for (unsigned int l = minlevel; l <= maxlevel; ++l) + { + d[l].reinit(1); + s[l].reinit(1); + } + } virtual void level_residual(unsigned level, Vector& dst, diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index acaa7ed9b8..62d6bcd842 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -42,7 +42,7 @@ main() { deallog << "Size " << size << endl; - unsigned int dim = (size+1)*(size+1); + unsigned int dim = (size-1)*(size-1); // Make matrix FDMatrix testproblem(size, size); @@ -73,6 +73,7 @@ main() check_method(rich,A,u,f,prec_ssor); check_method(cg,A,u,f,prec_ssor); +// testproblem.gnuplot_print(cout, u); check_method(bicgstab,A,u,f,prec_ssor); check_method(gmres,A,u,f,prec_ssor); diff --git a/tests/lac/testmatrix.cc b/tests/lac/testmatrix.cc index 8b22d664e7..53eeee5063 100644 --- a/tests/lac/testmatrix.cc +++ b/tests/lac/testmatrix.cc @@ -11,32 +11,32 @@ FDMatrix::FDMatrix(unsigned int nx, unsigned int ny) void FDMatrix::build_structure(SparseMatrixStruct& structure) const { - for(unsigned int i=0;i<=ny;i++) + for(unsigned int i=0;i<=ny-2;i++) { - for(unsigned int j=0;j<=nx; j++) + for(unsigned int j=0;j<=nx-2; j++) { // Number of the row to be entered - unsigned int row = j+(nx+1)*i; + unsigned int row = j+(nx-1)*i; structure.add(row, row); if (j>0) { structure.add(row-1, row); structure.add(row, row-1); } - if (j0) { - structure.add(row-nx, row); - structure.add(row, row-nx); + structure.add(row-(nx-1), row); + structure.add(row, row-(nx-1)); } - if (i void FDMatrix::laplacian(SparseMatrix& A) const { - for(unsigned int i=0;i<=ny;i++) + for(unsigned int i=0;i<=ny-2;i++) { - for(unsigned int j=0;j<=nx; j++) + for(unsigned int j=0;j<=nx-2; j++) { // Number of the row to be entered - unsigned int row = j+(nx+1)*i; + unsigned int row = j+(nx-1)*i; A.set(row, row, 4.); if (j>0) @@ -59,29 +59,46 @@ FDMatrix::laplacian(SparseMatrix& A) const A.set(row-1, row, -1.); A.set(row, row-1, -1.); } - if (j0) { - A.set(row-nx, row, -1.); - A.set(row, row-nx, -1.); + A.set(row-(nx-1), row, -1.); + A.set(row, row-(nx-1), -1.); } - if (i +void +FDMatrix::gnuplot_print(ostream& s, const Vector& V) 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; + s << (j+1)/(float)nx << '\t' << (i+1)/(float)ny << '\t' << V(row) << endl; + } + s << endl; + } + s << endl; +} + FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny, unsigned int nlevels) : - structures(nlevels-1), matrices(nlevels-1) + structures(nlevels), matrices(nlevels) { unsigned int power = 1 << nlevels; @@ -91,7 +108,7 @@ FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny, nx /= power; ny /= power; - for (unsigned int level = 0; level < nlevels-1; ++ level) + for (unsigned int level = 0; level < nlevels; ++ level) { build_matrix(nx,ny,structures[level],matrices[level]); nx *= 2; @@ -101,18 +118,18 @@ FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny, void FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny, - SparseMatrixStruct& structure, SparseMatrix& matrix) + SparseMatrixStruct& structure, SparseMatrix& matrix) { - structure.reinit((nx+1)*(ny+1),(2*nx+1)*(2*ny+1),9); + structure.reinit((nx-1)*(ny-1),(2*nx-1)*(2*ny-1),9); // Traverse all points of coarse grid - for (unsigned int i=0 ; i<=ny; ++i) - for (unsigned int j=0 ; j<=nx; ++j) + for (unsigned int i=1 ; i0) - structure.add(ncoarse,nfine-(2*nx+1)-1); + structure.add(ncoarse,nfine-(2*nx-1)-1); // upper left if (i0) - structure.add(ncoarse,nfine-(2*nx+1)+1); + structure.add(ncoarse,nfine-(2*nx-1)+1); // upper right if (i0) - structure.add(ncoarse,nfine-(2*nx+1)); + structure.add(ncoarse,nfine-(2*nx-1)); // upper if (i0) - matrix.set(ncoarse,nfine-(2*nx+1)-1,.25); + matrix.set(ncoarse,nfine-(2*nx-1)-1,.25); // upper left if (i0) - matrix.set(ncoarse,nfine-(2*nx+1)+1,.25); + matrix.set(ncoarse,nfine-(2*nx-1)+1,.25); // upper right if (i0) - matrix.set(ncoarse,nfine-(2*nx+1),.5); + matrix.set(ncoarse,nfine-(2*nx-1),.5); // upper if (i &dst, - const Vector &src) const + Vector &dst, + const Vector &src) const { Assert((to_level>0) && (to_level<=matrices.size()), ExcIndexRange(to_level, 0, matrices.size()+1)); @@ -205,8 +222,8 @@ FDMGTransfer::prolongate (const unsigned int to_level, void FDMGTransfer::restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const + Vector &dst, + const Vector &src) const { Assert((from_level>0) && (from_level<=matrices.size()), ExcIndexRange(from_level, 0, matrices.size()+1)); @@ -217,3 +234,6 @@ FDMGTransfer::restrict (const unsigned int from_level, template void FDMatrix::laplacian(SparseMatrix&) const; template void FDMatrix::laplacian(SparseMatrix&) const; +template void FDMatrix::gnuplot_print(ostream& s, const Vector& V) const; +template void FDMatrix::gnuplot_print(ostream& s, const Vector& V) const; + diff --git a/tests/lac/testmatrix.h b/tests/lac/testmatrix.h index 641f949e0f..659ab39752 100644 --- a/tests/lac/testmatrix.h +++ b/tests/lac/testmatrix.h @@ -28,6 +28,9 @@ class FDMatrix template void laplacian(SparseMatrix&) const; + template + void gnuplot_print(ostream&, const Vector&) const; + private: /** * Number of gridpoints in x-direction. @@ -63,16 +66,16 @@ class FDMGTransfer * function in #MGTranferBase#. */ virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const; + Vector &dst, + const Vector &src) const; /** * Implementation of abstract * function in #MGTranferBase#. */ virtual void restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const; + Vector &dst, + const Vector &src) const; /** * Exception. @@ -88,7 +91,7 @@ class FDMGTransfer /** * Prolongation matrices. */ - vector > matrices; + vector > matrices; /** * Matrix generator. @@ -98,5 +101,5 @@ class FDMGTransfer * fine to coarse (#vmult#). */ void build_matrix(unsigned int nx, unsigned int ny, - SparseMatrixStruct&, SparseMatrix&); + SparseMatrixStruct&, SparseMatrix&); };