From d21a049ff46ba55bec8bba7315ae967d159764d9 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 30 Apr 1999 14:09:21 +0000 Subject: [PATCH] more multigrid git-svn-id: https://svn.dealii.org/trunk@1236 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/multigrid/mg_base.h | 10 +- deal.II/deal.II/source/multigrid/mg_base.cc | 5 + deal.II/doc/parser/dealparse.y | 4 +- deal.II/lac/include/lac/mgbase.h | 10 +- deal.II/lac/source/mgbase.cc | 5 + tests/deal.II/mg2.cc | 22 +++ tests/lac/Makefile | 37 ++++- tests/lac/mg.cc | 160 ++++++++++++++++++++ tests/{deal.II/mg1.cc => lac/mgbase.cc} | 43 ++++-- tests/lac/solver.cc | 76 +++++----- tests/lac/testmatrix.cc | 137 +++++++++++++++++ tests/lac/testmatrix.h | 63 ++++++++ 12 files changed, 510 insertions(+), 62 deletions(-) create mode 100644 tests/deal.II/mg2.cc create mode 100644 tests/lac/mg.cc rename tests/{deal.II/mg1.cc => lac/mgbase.cc} (66%) diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index e864500b26..e0e83596d5 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -131,7 +131,7 @@ template class MGMatrix : public Subscriptor, - public Vector + public vector { public: /** @@ -159,7 +159,7 @@ class MGCoarseGridLACIteration * preconditioning method for later * use. */ - MGCoarseGridLACIteration(const SOLVER&, const MATRIX&, const PRECOND&); + MGCoarseGridLACIteration( SOLVER&, const MATRIX&, const PRECOND&); /** * Implementation of the abstract @@ -174,7 +174,7 @@ class MGCoarseGridLACIteration /** * Reference to the solver. */ - const SOLVER& solver; + SOLVER& solver; /** * Reference to the matrix. */ @@ -280,7 +280,7 @@ class MGBase template MGCoarseGridLACIteration -::MGCoarseGridLACIteration(const SOLVER& s, const MATRIX& m, const PRECOND& p) +::MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p) : solver(s), matrix(m), precondition(p) {} @@ -288,7 +288,7 @@ MGCoarseGridLACIteration template void MGCoarseGridLACIteration::operator() -(unsigned int level, Vector& dst, const Vector& src) const +(unsigned int, Vector& dst, const Vector& src) const { solver.solve(matrix, dst, src, precondition); } diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index 1ccc5258a6..c6252fa0d8 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -45,3 +45,8 @@ MGBase::level_mgstep(unsigned level, post_smooth.smooth(level, s[level], d[level]); } +MGTransferBase::~MGTransferBase() +{} + +MGSmootherBase::~MGSmootherBase() +{} diff --git a/deal.II/doc/parser/dealparse.y b/deal.II/doc/parser/dealparse.y index c8d0bafd2d..e84390d1d5 100644 --- a/deal.II/doc/parser/dealparse.y +++ b/deal.II/doc/parser/dealparse.y @@ -197,6 +197,8 @@ enumerator: IDENTIFIER class_declaration: class_head | class_head { + cout << "@@@" << endl << setw(OUTW) << "@Class-Definition:" << $1 << endl + << setw(OUTW) << "@In-Class:" << class_stack.top() << endl; class_name = class_stack.top() + string("::") + class_name; class_stack.push(class_name); access_stack.push(string("private")); @@ -204,8 +206,6 @@ class_declaration: class_head { class_stack.pop(); access_stack.pop(); - cout << "@@@" << endl << setw(OUTW) << "@Class-Definition:" << $1 << endl - << setw(OUTW) << "@In-Class:" << class_stack.top() << endl; } ; diff --git a/deal.II/lac/include/lac/mgbase.h b/deal.II/lac/include/lac/mgbase.h index f614a308dc..3cb773a6b2 100644 --- a/deal.II/lac/include/lac/mgbase.h +++ b/deal.II/lac/include/lac/mgbase.h @@ -131,7 +131,7 @@ template class MGMatrix : public Subscriptor, - public Vector + public vector { public: /** @@ -159,7 +159,7 @@ class MGCoarseGridLACIteration * preconditioning method for later * use. */ - MGCoarseGridLACIteration(const SOLVER&, const MATRIX&, const PRECOND&); + MGCoarseGridLACIteration( SOLVER&, const MATRIX&, const PRECOND&); /** * Implementation of the abstract @@ -174,7 +174,7 @@ class MGCoarseGridLACIteration /** * Reference to the solver. */ - const SOLVER& solver; + SOLVER& solver; /** * Reference to the matrix. */ @@ -280,7 +280,7 @@ class MGBase template MGCoarseGridLACIteration -::MGCoarseGridLACIteration(const SOLVER& s, const MATRIX& m, const PRECOND& p) +::MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p) : solver(s), matrix(m), precondition(p) {} @@ -288,7 +288,7 @@ MGCoarseGridLACIteration template void MGCoarseGridLACIteration::operator() -(unsigned int level, Vector& dst, const Vector& src) const +(unsigned int, Vector& dst, const Vector& src) const { solver.solve(matrix, dst, src, precondition); } diff --git a/deal.II/lac/source/mgbase.cc b/deal.II/lac/source/mgbase.cc index 1ccc5258a6..c6252fa0d8 100644 --- a/deal.II/lac/source/mgbase.cc +++ b/deal.II/lac/source/mgbase.cc @@ -45,3 +45,8 @@ MGBase::level_mgstep(unsigned level, post_smooth.smooth(level, s[level], d[level]); } +MGTransferBase::~MGTransferBase() +{} + +MGSmootherBase::~MGSmootherBase() +{} diff --git a/tests/deal.II/mg2.cc b/tests/deal.II/mg2.cc new file mode 100644 index 0000000000..6e3ace4dac --- /dev/null +++ b/tests/deal.II/mg2.cc @@ -0,0 +1,22 @@ +// $Id$ + +// deal_II_libraries.g=-ldeal_II_2d.g +// deal_II_libraries=-ldeal_II_2d + +#include +#include +#include +#include +#include + +main() +{ + Triangulation<2> tr; + MGDoFHandler<2> dof(&tr); + FELinear<2> fe; + + tr.create_hypercube(-1.,1.); + tr.refine_global(3); + dof.distribute_dofs(fe); +} + diff --git a/tests/lac/Makefile b/tests/lac/Makefile index 3d1a1cf7ab..1b7ce6c299 100644 --- a/tests/lac/Makefile +++ b/tests/lac/Makefile @@ -41,7 +41,7 @@ libraries = $(libs) flags = $(CXXFLAGS.o) $(CXXFLAGS) endif -all: solver +all: solver mgbase mg ############################################################ # Typical block for building a running program @@ -67,6 +67,39 @@ endif solver: $(solver-o-files) $(libraries) $(CXX) $(flags) -o $@ $^ + +############################################################ + + +mgbase-cc-files = mgbase.cc + +ifeq ($(debug-mode),on) +mgbase-o-files = $(mgbase-cc-files:.cc=.go) +else +mgbase-o-files = $(mgbase-cc-files:.cc=.o) +endif + +mgbase: $(mgbase-o-files) $(libraries) + $(CXX) $(flags) -o $@ $^ + + + +############################################################ + + +mg-cc-files = mg.cc testmatrix.cc + +ifeq ($(debug-mode),on) +mg-o-files = $(mg-cc-files:.cc=.go) +else +mg-o-files = $(mg-cc-files:.cc=.o) +endif + +mg: $(mg-o-files) $(libraries) + $(CXX) $(flags) -o $@ $^ + + + ############################################################ # Continue with other targets if needed ############################################################ @@ -99,7 +132,7 @@ veryclean: clean # Automatic generation of dependencies ############################################################ -all-cc-files = $(solver-cc-files) # $(target2-cc-files) ... +all-cc-files = $(solver-cc-files) $(mgbase-cc-files) $(mg-cc-files) # $(target2-cc-files) ... Make.depend: $(all-cc-files) @echo =====Dependecies=== Make.depend diff --git a/tests/lac/mg.cc b/tests/lac/mg.cc new file mode 100644 index 0000000000..439d47f20d --- /dev/null +++ b/tests/lac/mg.cc @@ -0,0 +1,160 @@ +// $Id$ + +// Test program for multigrid with test matrices (fd) + +#include +#include +#include "testmatrix.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class FDMG + : + public MGBase +{ + private: + SmartPointer > >matrices; + public: + FDMG(unsigned int maxlevel, MGMatrix >& matrices, + FDMGTransfer& transfer) + : + MGBase(transfer, maxlevel, 0), + matrices(&matrices) + {} + + + virtual void level_residual (unsigned int level, + Vector& dst, + const Vector& src, + const Vector& rhs); +}; + +class MGSmootherLAC + : + MGSmootherBase +{ + private: + SmartPointer > >matrices; + public: + MGSmootherLAC(MGMatrix >&); + + virtual void smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const; + +}; + +typedef MGCoarseGridLACIteration , Vector >, +SparseMatrix, PreconditionRelaxation , Vector > > +Coarse; + +class MGPrecondition +{ + private: + FDMG& mg; + MGSmootherLAC& smooth; + Coarse& coarse; + public: + MGPrecondition(FDMG& m, MGSmootherLAC& s, Coarse& c) + : + mg(m), smooth(s), coarse(c) + {} + + + void operator () (Vector& dst, const Vector& src) const + { + + } + + +}; + + +main() +{ + ofstream logfile("mg.output"); + deallog.attach(logfile); + + PrimitiveVectorMemory > mem; + SolverControl control(100, 1.e-5, true); + + const unsigned int base = 3; + const unsigned int maxlevel = 6; + + const unsigned int maxsize = base * (1 << maxlevel); + + // grid transfer + FDMGTransfer transfer(maxsize, maxsize, maxlevel); + + // coarse grid solver + PrimitiveVectorMemory > cgmem; + SolverControl cgcontrol(100, 1.e-12); + SolverCG , Vector > cgcg(cgcontrol,cgmem); + + vector structure(maxlevel+1); + MGMatrix > A(maxlevel+1); + + for (unsigned int level = 0; level <= maxlevel; ++level) + { + const unsigned int size = base * (1 << level); + + const unsigned int dim = (size+1)*(size+1); + + // Make matrix + 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]); + + FDMG multigrid(level, A, transfer); + + PreconditionRelaxation , Vector > + cgprec(A[0], &SparseMatrix ::template precondition_SSOR, 1.2); + + Coarse coarsegrid(cgcg, A[0], cgprec); + + MGSmootherLAC smoother(A); + +// MGPrecondition prec(multigrid); + } +} + +void +FDMG::level_residual (unsigned int level, + Vector& dst, + const Vector& src, + const Vector& rhs) +{ + (*matrices)[level].residual(dst, src, rhs); +} + +MGSmootherLAC::MGSmootherLAC(MGMatrix >& matrix) + : + matrices(&matrix) +{} + + +void +MGSmootherLAC::smooth (const unsigned int level, + 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); + + rich.solve((*matrices)[level], u, rhs, prec); +} diff --git a/tests/deal.II/mg1.cc b/tests/lac/mgbase.cc similarity index 66% rename from tests/deal.II/mg1.cc rename to tests/lac/mgbase.cc index 8c2efd4f32..acd703a6c8 100644 --- a/tests/deal.II/mg1.cc +++ b/tests/lac/mgbase.cc @@ -3,8 +3,7 @@ // deal_II_libraries.g=-ldeal_II_2d.g // deal_II_libraries=-ldeal_II_2d -#include -#include +#include class TransferTest : @@ -16,7 +15,7 @@ class TransferTest virtual void restrict(unsigned l, Vector& dst, const Vector& src) const; - friend class MultiGridBase; + friend class MGBase; }; class SmoothTest @@ -29,17 +28,29 @@ class SmoothTest const Vector& rhs) const; }; +class MGCoarseGridTest + : + public MGCoarseGridSolver +{ + public: + virtual void operator () (unsigned int level, + Vector& dst, + const Vector& src) const ; +}; + class MGTest : - public MultiGridBase + public MGBase { MGSmootherBase& smoother; - + MGCoarseGridSolver& solver; + public: - MGTest(TransferTest& t, SmoothTest& s) - : MultiGridBase(t, 9, 3), - smoother(s) + MGTest(TransferTest& t, SmoothTest& s, MGCoarseGridTest& c) + : MGBase(t, 9, 3), + smoother(s), + solver(c) {} virtual void level_residual(unsigned level, @@ -49,7 +60,7 @@ class MGTest void doit() { - level_mgstep(9, smoother, smoother); + level_mgstep(9, smoother, smoother, solver); } }; @@ -57,7 +68,9 @@ main() { TransferTest tr; SmoothTest s; - MGTest mg(tr, s); + MGCoarseGridTest c; + + MGTest mg(tr, s, c); mg.doit(); } @@ -83,7 +96,7 @@ SmoothTest::smooth (const unsigned int level, Vector &, const Vector&) const { - cout << "Smoothing on" << level << endl; + cout << "Smoothing on " << level << endl; } void @@ -94,3 +107,11 @@ MGTest::level_residual(unsigned l, { cout << "Residual on " << l << endl; } + +void +MGCoarseGridTest::operator() (unsigned int level, + Vector&, + const Vector&) const +{ + cout << "Solving on " << level << endl; +} diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index a70a6c4787..acaa7ed9b8 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -3,6 +3,7 @@ // Test program for linear solvers. #include +#include #include "testmatrix.h" #include #include @@ -15,17 +16,27 @@ #include #include -#define MATRIX SparseMatrix -#define VECTOR Vector +template +void +check_method( SOLVER& solver, const MATRIX& A, + VECTOR& u, VECTOR& f, const PRECONDITION& P) +{ + u = 0.; + f = 1.; + solver.solve(A,u,f,P); +} main() { - PrimitiveVectorMemory mem; + ofstream logfile("solver.output"); + deallog.attach(logfile); + + PrimitiveVectorMemory > mem; SolverControl control(100, 1.e-5, true); - SolverCG cg(control, mem); - SolverGMRES gmres(control, mem,20); - SolverBicgstab bicgstab(control, mem); - SolverRichardson rich(control, mem); + SolverCG , Vector > cg(control, mem); + SolverGMRES , Vector > gmres(control, mem,20); + SolverBicgstab , Vector > bicgstab(control, mem); + SolverRichardson , Vector > rich(control, mem); for (unsigned int size=10; size <= 10; size *= 10) { @@ -38,49 +49,40 @@ main() SparseMatrixStruct structure(dim, dim, 5); testproblem.build_structure(structure); structure.compress(); - MATRIX A(structure); + SparseMatrix A(structure); testproblem.laplacian(A); - PreconditionIdentity + PreconditionIdentity > prec_no; - PreconditionRelaxation - prec_ssor(A, &MATRIX::template precondition_SSOR, 1.2); + PreconditionRelaxation , Vector > + prec_sor(A, &SparseMatrix ::template precondition_SOR, 1.2); + PreconditionRelaxation , Vector > + prec_ssor(A, &SparseMatrix ::template precondition_SSOR, 1.2); - VECTOR f(dim); - VECTOR u(dim); + Vector f(dim); + Vector u(dim); deallog.push("no"); - f = 1.; - u = 0.; - cg.solve(A,u,f,prec_no); - - f = 1.; - u = 0.; - bicgstab.solve(A,u,f,prec_no); - - f = 1.; - u = 0.; - gmres.solve(A,u,f,prec_no); + check_method(cg,A,u,f,prec_no); + check_method(bicgstab,A,u,f,prec_no); + check_method(gmres,A,u,f,prec_no); deallog.pop(); deallog.push("ssor"); - f = 1.; - u = 0.; - rich.solve(A,u,f,prec_ssor); - - f = 1.; - u = 0.; - cg.solve(A,u,f,prec_ssor); + check_method(rich,A,u,f,prec_ssor); + check_method(cg,A,u,f,prec_ssor); + check_method(bicgstab,A,u,f,prec_ssor); + check_method(gmres,A,u,f,prec_ssor); - f = 1.; - u = 0.; - bicgstab.solve(A,u,f,prec_ssor); + deallog.pop(); + deallog.push("sor"); - f = 1.; - u = 0.; - gmres.solve(A,u,f,prec_ssor); + check_method(rich,A,u,f,prec_sor); + check_method(cg,A,u,f,prec_sor); + check_method(bicgstab,A,u,f,prec_sor); + check_method(gmres,A,u,f,prec_sor); deallog.pop(); } diff --git a/tests/lac/testmatrix.cc b/tests/lac/testmatrix.cc index 656cadef33..8b22d664e7 100644 --- a/tests/lac/testmatrix.cc +++ b/tests/lac/testmatrix.cc @@ -78,5 +78,142 @@ FDMatrix::laplacian(SparseMatrix& A) const } } +FDMGTransfer::FDMGTransfer(unsigned int nx, unsigned int ny, + unsigned int nlevels) + : + structures(nlevels-1), matrices(nlevels-1) +{ + unsigned int power = 1 << nlevels; + + Assert ((nx%power)==0, ExcDivide(nx,power)); + Assert ((ny%power)==0, ExcDivide(ny,power)); + + nx /= power; + ny /= power; + + for (unsigned int level = 0; level < nlevels-1; ++ level) + { + build_matrix(nx,ny,structures[level],matrices[level]); + nx *= 2; + ny *= 2; + } +} + +void +FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny, + SparseMatrixStruct& structure, SparseMatrix& matrix) +{ + 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) + { + // coarse grid point number + unsigned int ncoarse =j+(nx+1)*i; + // same point on fine grid + unsigned int nfine = 2*j+(4*nx+2)*i; + + structure.add(ncoarse,nfine); + // left + if (j>0) + { + structure.add(ncoarse,nfine-1); + // lower left + if (i>0) + structure.add(ncoarse,nfine-(2*nx+1)-1); + // upper left + if (i0) + structure.add(ncoarse,nfine-(2*nx+1)+1); + // upper right + if (i0) + structure.add(ncoarse,nfine-(2*nx+1)); + // upper + if (i0) + { + matrix.set(ncoarse,nfine-1,.5); + // lower left + if (i>0) + matrix.set(ncoarse,nfine-(2*nx+1)-1,.25); + // upper left + if (i0) + matrix.set(ncoarse,nfine-(2*nx+1)+1,.25); + // upper right + if (i0) + matrix.set(ncoarse,nfine-(2*nx+1),.5); + // upper + if (i &dst, + const Vector &src) const +{ + Assert((to_level>0) && (to_level<=matrices.size()), + ExcIndexRange(to_level, 0, matrices.size()+1)); + + matrices[to_level-1].Tvmult(dst,src); +} + + +void +FDMGTransfer::restrict (const unsigned int from_level, + Vector &dst, + const Vector &src) const +{ + Assert((from_level>0) && (from_level<=matrices.size()), + ExcIndexRange(from_level, 0, matrices.size()+1)); + + matrices[from_level-1].vmult(dst,src); +} + + template void FDMatrix::laplacian(SparseMatrix&) const; template void FDMatrix::laplacian(SparseMatrix&) const; diff --git a/tests/lac/testmatrix.h b/tests/lac/testmatrix.h index ffb32d8e8a..641f949e0f 100644 --- a/tests/lac/testmatrix.h +++ b/tests/lac/testmatrix.h @@ -1,6 +1,8 @@ // $Id$ #include +#include +#include /** * Finite difference matrix on uniform grid. @@ -37,3 +39,64 @@ class FDMatrix */ unsigned int ny; }; + +/** + * Grid transfer for finite differences on uniform grid. + */ + +class FDMGTransfer + : + public MGTransferBase +{ + public: + /** + * Constructor. Prepares grid + * transfer matrices for #nlevels# + * levels on an #nx# times #ny# + * grid. + */ + FDMGTransfer(unsigned int nx, unsigned int ny, + unsigned int nlevels); + + /** + * Implementation of abstract + * function in #MGTranferBase#. + */ + virtual void prolongate (const unsigned int to_level, + 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; + + /** + * Exception. + */ + DeclException2(ExcDivide, unsigned int, unsigned int, + << "Cannot divide " << arg1 << " by " << arg2); + + private: + /** + * Prolongation matrix structures. + */ + vector structures; + /** + * Prolongation matrices. + */ + vector > matrices; + + /** + * Matrix generator. + * The arguments #nx# and #ny# + * are the numbers on the COARSE level. + * Builds a transfer matrix from + * fine to coarse (#vmult#). + */ + void build_matrix(unsigned int nx, unsigned int ny, + SparseMatrixStruct&, SparseMatrix&); +}; -- 2.39.5