#include <vector>
+#define FLOAT float
+
/**
* Abstract base class for coarse grid solvers.
* about the matrix is removed to
* that class.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const = 0;
+ virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src) const = 0;
};
* things.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const = 0;
+ Vector<FLOAT> &u,
+ const Vector<FLOAT> &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<FLOAT> &u,
+ const Vector<FLOAT> &rhs) const;
};
/**
* level.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<FLOAT> &dst,
+ const Vector<FLOAT> &src) const = 0;
/**
* Restrict a vector from level
* level.
*/
virtual void restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<FLOAT> &dst,
+ const Vector<FLOAT> &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 VECTOR>
+class MGVector
+ :
+ public Subscriptor,
+ public vector<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<MATRIX># combined
/**
* 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;
};
* matrix, vectors and
* preconditioner.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const;
+ virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src) const;
private:
/**
* Reference to the solver.
* */
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<Vector<float> > d;
+ MGVector<Vector<FLOAT> > d;
/**
* Auxiliary vector, solution.
*/
- vector<Vector<float> > s;
+ MGVector<Vector<FLOAT> > s;
+ private:
/**
* Auxiliary vector.
*/
- Vector<float> t;
+ Vector<FLOAT> t;
/**
* Prolongation and restriction object.
SmartPointer<const MGTransferBase> transfer;
protected:
- /**
- * Highest level of cells.
- */
- unsigned int maxlevel;
-
- /**
- * Level for coarse grid solution.
- */
- unsigned int minlevel;
/**
* The actual v-cycle multigrid method.
* derived class.
*/
virtual void level_residual(unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs) = 0;
+ Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src,
+ const Vector<FLOAT>& rhs) = 0;
public:
*/
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);
};
template<class SOLVER, class MATRIX, class PRECOND>
void
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<FLOAT>& dst, const Vector<FLOAT>& src) const
{
solver.solve(matrix, dst, src, precondition);
}
+//////////////////////////////////////////////////////////////////////
+
+template<class VECTOR>
+MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+ :
+ vector<VECTOR>(max-min+1),
+ minlevel(min)
+{}
+
+
+template<class VECTOR>
+VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i)
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<VECTOR>::operator[](i-minlevel);
+}
+
+
+template<class VECTOR>
+const VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i) const
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<VECTOR>::operator[](i-minlevel);
+}
+
+//////////////////////////////////////////////////////////////////////
+
template<class MATRIX>
-MGMatrix<MATRIX>::MGMatrix(unsigned int nlevels)
+MGMatrix<MATRIX>::MGMatrix(unsigned int min, unsigned int max)
:
- vector<MATRIX>(nlevels)
+ vector<MATRIX>(max-min+1),
+ minlevel(min)
{}
+template<class MATRIX>
+MATRIX&
+MGMatrix<MATRIX>::operator[](unsigned int i)
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<MATRIX>::operator[](i-minlevel);
+}
+
+
+template<class MATRIX>
+const MATRIX&
+MGMatrix<MATRIX>::operator[](unsigned int i) const
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<MATRIX>::operator[](i-minlevel);
+}
+
+
#endif
vector<vector<int> > interior_boundary_dofs;
};
-/**
- * Smoother doing nothing.
- * @author Guido Kanschat, 1999
- */
-class MGSmootherIdentity
- :
- public MGSmoother
-{
- public:
- /**
- * Constructor.
- */
- template<int dim>
- MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof);
- /**
- * Implementation of the interface in #MGSmootherBase#.
- * This function does nothing.
- */
- virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
-};
-
/**
* Implementation of an SOR smoother.
*
vector<vector<int> > interior_boundary_dofs;
};
-/**
- * Smoother doing nothing.
- * @author Guido Kanschat, 1999
- */
-class MGSmootherIdentity
- :
- public MGSmoother
-{
- public:
- /**
- * Constructor.
- */
- template<int dim>
- MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof);
- /**
- * Implementation of the interface in #MGSmootherBase#.
- * This function does nothing.
- */
- virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
-};
-
/**
* Implementation of an SOR smoother.
*
// $Id$
#include <lac/mgbase.h>
+#include <iostream>
+#include <cmath>
+
+template<class VECTOR>
+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<n;++i)
+ {
+ for (unsigned int j=0;j<n;++j)
+// s << ((FLOAT)i)/n << '\t' << ((FLOAT)j)/n << '\t' << v(k++) << endl;
+ s << '\t' << v(k++);
+ s << endl;
+ }
+ s << endl;
+}
MGBase::~MGBase ()
{};
MGBase::MGBase(const MGTransferBase& transfer,
unsigned maxlevel, unsigned minlevel)
:
- d(maxlevel-minlevel),
- s(maxlevel-minlevel),
- transfer(&transfer),
- maxlevel(maxlevel), minlevel(minlevel)
-{}
+ maxlevel(maxlevel), minlevel(minlevel),
+ d(minlevel,maxlevel),
+ s(minlevel,maxlevel),
+ transfer(&transfer)
+{
+}
+
+void
+MGBase::vcycle(const MGSmootherBase& pre_smooth,
+ const MGSmootherBase& post_smooth,
+ const MGCoarseGridSolver& coarse_grid_solver)
+{
+ level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
+}
+
void
MGBase::level_mgstep(unsigned level,
const MGSmootherBase& post_smooth,
const MGCoarseGridSolver& coarse_grid_solver)
{
+ s[level] = 0.;
+
if (level == minlevel)
{
coarse_grid_solver(level, s[level], d[level]);
return;
}
-
- // smoothing of the residual by modifying s
- pre_smooth.smooth(level, s[level], d[level]);
+ // smoothing of the residual by modifying s
+ pre_smooth.smooth(level, s[level], d[level]);
// t = d-As
+ t.reinit(s[level].size());
level_residual(level, t, s[level], d[level]);
-
+
// make t rhs of lower level
- transfer->restrict(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<FLOAT> &,
+ const Vector<FLOAT> &) const
+{}
+
+
u(*p) = 0;
};
-//////////////////////////////////////////////////////////////////////
-
-
-template<int dim>
-MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof)
- :
- MGSmoother(mg_dof, 0)
-{}
-
-void
-MGSmootherIdentity::smooth (const unsigned int,
- Vector<float> &,
- const Vector<float> &) const
-{}
-
-
//////////////////////////////////////////////////////////////////////
template<int dim>
const MGMatrix<SparseMatrix<float> >& matrix,
unsigned int steps);
-template MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<deal_II_dimension> &mg_dof);
-
u(*p) = 0;
};
-//////////////////////////////////////////////////////////////////////
-
-
-template<int dim>
-MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<dim> &mg_dof)
- :
- MGSmoother(mg_dof, 0)
-{}
-
-void
-MGSmootherIdentity::smooth (const unsigned int,
- Vector<float> &,
- const Vector<float> &) const
-{}
-
-
//////////////////////////////////////////////////////////////////////
template<int dim>
const MGMatrix<SparseMatrix<float> >& matrix,
unsigned int steps);
-template MGSmootherIdentity::MGSmootherIdentity(const MGDoFHandler<deal_II_dimension> &mg_dof);
-
#include <vector>
+#define FLOAT float
+
/**
* Abstract base class for coarse grid solvers.
* about the matrix is removed to
* that class.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const = 0;
+ virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src) const = 0;
};
* things.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const = 0;
+ Vector<FLOAT> &u,
+ const Vector<FLOAT> &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<FLOAT> &u,
+ const Vector<FLOAT> &rhs) const;
};
/**
* level.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<FLOAT> &dst,
+ const Vector<FLOAT> &src) const = 0;
/**
* Restrict a vector from level
* level.
*/
virtual void restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<FLOAT> &dst,
+ const Vector<FLOAT> &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 VECTOR>
+class MGVector
+ :
+ public Subscriptor,
+ public vector<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<MATRIX># combined
/**
* 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;
};
* matrix, vectors and
* preconditioner.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const;
+ virtual void operator() (unsigned int level, Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src) const;
private:
/**
* Reference to the solver.
* */
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<Vector<float> > d;
+ MGVector<Vector<FLOAT> > d;
/**
* Auxiliary vector, solution.
*/
- vector<Vector<float> > s;
+ MGVector<Vector<FLOAT> > s;
+ private:
/**
* Auxiliary vector.
*/
- Vector<float> t;
+ Vector<FLOAT> t;
/**
* Prolongation and restriction object.
SmartPointer<const MGTransferBase> transfer;
protected:
- /**
- * Highest level of cells.
- */
- unsigned int maxlevel;
-
- /**
- * Level for coarse grid solution.
- */
- unsigned int minlevel;
/**
* The actual v-cycle multigrid method.
* derived class.
*/
virtual void level_residual(unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs) = 0;
+ Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src,
+ const Vector<FLOAT>& rhs) = 0;
public:
*/
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);
};
template<class SOLVER, class MATRIX, class PRECOND>
void
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<FLOAT>& dst, const Vector<FLOAT>& src) const
{
solver.solve(matrix, dst, src, precondition);
}
+//////////////////////////////////////////////////////////////////////
+
+template<class VECTOR>
+MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+ :
+ vector<VECTOR>(max-min+1),
+ minlevel(min)
+{}
+
+
+template<class VECTOR>
+VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i)
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<VECTOR>::operator[](i-minlevel);
+}
+
+
+template<class VECTOR>
+const VECTOR&
+MGVector<VECTOR>::operator[](unsigned int i) const
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<VECTOR>::operator[](i-minlevel);
+}
+
+//////////////////////////////////////////////////////////////////////
+
template<class MATRIX>
-MGMatrix<MATRIX>::MGMatrix(unsigned int nlevels)
+MGMatrix<MATRIX>::MGMatrix(unsigned int min, unsigned int max)
:
- vector<MATRIX>(nlevels)
+ vector<MATRIX>(max-min+1),
+ minlevel(min)
{}
+template<class MATRIX>
+MATRIX&
+MGMatrix<MATRIX>::operator[](unsigned int i)
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<MATRIX>::operator[](i-minlevel);
+}
+
+
+template<class MATRIX>
+const MATRIX&
+MGMatrix<MATRIX>::operator[](unsigned int i) const
+{
+ Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
+ return vector<MATRIX>::operator[](i-minlevel);
+}
+
+
#endif
g.add(alpha,Ad);
x.add(alpha,d );
- res = sqrt(g*g);
+ res = g.l2_norm();
conv = control().check(it,res);
if (conv) break;
// $Id$
#include <lac/mgbase.h>
+#include <iostream>
+#include <cmath>
+
+template<class VECTOR>
+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<n;++i)
+ {
+ for (unsigned int j=0;j<n;++j)
+// s << ((FLOAT)i)/n << '\t' << ((FLOAT)j)/n << '\t' << v(k++) << endl;
+ s << '\t' << v(k++);
+ s << endl;
+ }
+ s << endl;
+}
MGBase::~MGBase ()
{};
MGBase::MGBase(const MGTransferBase& transfer,
unsigned maxlevel, unsigned minlevel)
:
- d(maxlevel-minlevel),
- s(maxlevel-minlevel),
- transfer(&transfer),
- maxlevel(maxlevel), minlevel(minlevel)
-{}
+ maxlevel(maxlevel), minlevel(minlevel),
+ d(minlevel,maxlevel),
+ s(minlevel,maxlevel),
+ transfer(&transfer)
+{
+}
+
+void
+MGBase::vcycle(const MGSmootherBase& pre_smooth,
+ const MGSmootherBase& post_smooth,
+ const MGCoarseGridSolver& coarse_grid_solver)
+{
+ level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
+}
+
void
MGBase::level_mgstep(unsigned level,
const MGSmootherBase& post_smooth,
const MGCoarseGridSolver& coarse_grid_solver)
{
+ s[level] = 0.;
+
if (level == minlevel)
{
coarse_grid_solver(level, s[level], d[level]);
return;
}
-
- // smoothing of the residual by modifying s
- pre_smooth.smooth(level, s[level], d[level]);
+ // smoothing of the residual by modifying s
+ pre_smooth.smooth(level, s[level], d[level]);
// t = d-As
+ t.reinit(s[level].size());
level_residual(level, t, s[level], d[level]);
-
+
// make t rhs of lower level
- transfer->restrict(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<FLOAT> &,
+ const Vector<FLOAT> &) const
+{}
+
+
};
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);
};
#include <lac/precondition.h>
#include <lac/mgbase.h>
+template<class VECTOR>
+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<n;++i)
+ {
+ for (unsigned int j=0;j<n;++j)
+// s << ((FLOAT)i)/n << '\t' << ((FLOAT)j)/n << '\t' << v(k++) << endl;
+ s << '\t' << v(k++);
+ s << endl;
+ }
+ s << endl;
+}
+
class FDMG
:
public MGBase
{
private:
- SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ /**
+ * Pointer to the level matrices.
+ */
+ SmartPointer<MGMatrix<SparseMatrix<FLOAT> > >matrices;
public:
- FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
- FDMGTransfer& transfer)
- :
- MGBase(transfer, maxlevel, 0),
- matrices(&matrices)
- {}
+ FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<FLOAT> >& matrices,
+ FDMGTransfer& transfer);
virtual void level_residual (unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs);
+ Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src,
+ const Vector<FLOAT>& rhs);
+
+ void copy_to_mg(const Vector<double>& rhs);
+ void copy_from_mg(Vector<double>& lsg);
};
class MGSmootherLAC
:
- MGSmootherBase
+ public MGSmootherBase
{
private:
- SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ SmartPointer<MGMatrix<SparseMatrix<FLOAT> > >matrices;
public:
- MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+ MGSmootherLAC(MGMatrix<SparseMatrix<FLOAT> >&);
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
+ Vector<FLOAT> &u,
+ const Vector<FLOAT> &rhs) const;
};
-typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<float> , Vector<float> >,
-SparseMatrix<float>, PreconditionRelaxation<SparseMatrix<float> , Vector<float> > >
+typedef MGCoarseGridLACIteration<SolverCG<SparseMatrix<FLOAT> , Vector<FLOAT> >,
+SparseMatrix<FLOAT>, /*PreconditionRelaxation<SparseMatrix<FLOAT> ,*/
+ PreconditionIdentity<Vector<FLOAT> > >
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)
{}
void operator () (Vector<double>& dst, const Vector<double>& src) const
{
+// print_vector(cout, src);
+ mg.copy_to_mg(src);
+ mg.vcycle(smooth,smooth,coarse);
+ mg.copy_from_mg(dst);
}
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);
FDMGTransfer transfer(maxsize, maxsize, maxlevel);
// coarse grid solver
- PrimitiveVectorMemory<Vector<float> > cgmem;
- SolverControl cgcontrol(100, 1.e-12);
- SolverCG<SparseMatrix<float> , Vector<float> > cgcg(cgcontrol,cgmem);
+ PrimitiveVectorMemory<Vector<FLOAT> > cgmem;
+ ReductionControl cgcontrol(100, 1.e-30, 1.e-2, false, false);
+ SolverCG<SparseMatrix<FLOAT> , Vector<FLOAT> > cgcg(cgcontrol,cgmem);
- vector<SparseMatrixStruct > structure(maxlevel+1);
- MGMatrix<SparseMatrix<float> > 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<SparseMatrixStruct > structure(maxlevel+1);
+ MGMatrix<SparseMatrix<FLOAT> > 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<SparseMatrix<float> , Vector<float> >
- cgprec(A[0], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.2);
+// PreconditionRelaxation<SparseMatrix<FLOAT> , Vector<FLOAT> >
+ PreconditionIdentity<Vector<FLOAT> >
+ cgprec;//(A[minlevel], &SparseMatrix<FLOAT> ::template precondition_SSOR<FLOAT>, 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<SparseMatrix<FLOAT> , Vector<double> > solver(control, mem);
+ SolverCG<SparseMatrix<FLOAT> , Vector<double> > solver(control, mem);
+
+ Vector<double> u(dim);
+ Vector<double> 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<SparseMatrix<FLOAT> >& 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<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs)
+ Vector<FLOAT>& dst,
+ const Vector<FLOAT>& src,
+ const Vector<FLOAT>& rhs)
{
(*matrices)[level].residual(dst, src, rhs);
}
-MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+void
+FDMG::copy_to_mg(const Vector<double>& v)
+{
+ d[maxlevel] = v;
+}
+
+void
+FDMG::copy_from_mg(Vector<double>& v)
+{
+ v = s[maxlevel];
+}
+
+
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<FLOAT> >& matrix)
:
matrices(&matrix)
{}
void
MGSmootherLAC::smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const
+ Vector<FLOAT> &u,
+ const Vector<FLOAT> &rhs) const
{
- SolverControl control(0.,1);
- 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.2);
+ 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.);
rich.solve((*matrices)[level], u, rhs, prec);
}
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<float>& dst,
{
deallog << "Size " << size << endl;
- unsigned int dim = (size+1)*(size+1);
+ unsigned int dim = (size-1)*(size-1);
// Make matrix
FDMatrix testproblem(size, size);
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);
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 (j<nx)
+ if (j<nx-2)
{
structure.add(row+1, row);
structure.add(row, row+1);
}
if (i>0)
{
- structure.add(row-nx, row);
- structure.add(row, row-nx);
+ structure.add(row-(nx-1), row);
+ structure.add(row, row-(nx-1));
}
- if (i<ny)
+ if (i<ny-2)
{
- structure.add(row+nx, row);
- structure.add(row, row+nx);
+ structure.add(row+(nx-1), row);
+ structure.add(row, row+(nx-1));
}
}
}
void
FDMatrix::laplacian(SparseMatrix<number>& 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)
A.set(row-1, row, -1.);
A.set(row, row-1, -1.);
}
- if (j<nx)
+ if (j<nx-2)
{
A.set(row+1, row, -1.);
A.set(row, row+1, -1.);
}
if (i>0)
{
- 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<ny)
+ if (i<ny-2)
{
- 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.);
}
}
}
}
+template<typename number>
+void
+FDMatrix::gnuplot_print(ostream& s, const Vector<number>& 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;
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;
void
FDMGTransfer::build_matrix(unsigned int nx, unsigned int ny,
- SparseMatrixStruct& structure, SparseMatrix<float>& matrix)
+ SparseMatrixStruct& structure, SparseMatrix<FLOAT>& 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 ; i<ny; ++i)
+ for (unsigned int j=1 ; j<nx; ++j)
{
// coarse grid point number
- unsigned int ncoarse =j+(nx+1)*i;
+ unsigned int ncoarse =j+(nx-1)*i -(nx-1) -1;
// same point on fine grid
- unsigned int nfine = 2*j+(4*nx+2)*i;
+ unsigned int nfine = 2*j+(4*nx-2)*i -(2*nx-1) -1;
structure.add(ncoarse,nfine);
// left
structure.add(ncoarse,nfine-1);
// lower left
if (i>0)
- structure.add(ncoarse,nfine-(2*nx+1)-1);
+ structure.add(ncoarse,nfine-(2*nx-1)-1);
// upper left
if (i<ny)
- structure.add(ncoarse,nfine+(2*nx+1)-1);
+ structure.add(ncoarse,nfine+(2*nx-1)-1);
}
// right
if (j<nx)
structure.add(ncoarse, nfine+1);
// lower right
if (i>0)
- structure.add(ncoarse,nfine-(2*nx+1)+1);
+ structure.add(ncoarse,nfine-(2*nx-1)+1);
// upper right
if (i<ny)
- structure.add(ncoarse,nfine+(2*nx+1)+1);
+ structure.add(ncoarse,nfine+(2*nx-1)+1);
}
// lower
if (i>0)
- structure.add(ncoarse,nfine-(2*nx+1));
+ structure.add(ncoarse,nfine-(2*nx-1));
// upper
if (i<ny)
- structure.add(ncoarse,nfine+(2*nx+1));
+ structure.add(ncoarse,nfine+(2*nx-1));
}
structure.compress();
matrix.reinit(structure);
- for (unsigned int i=0 ; i<=ny; ++i)
- for (unsigned int j=0 ; j<=nx; ++j)
+ for (unsigned int i=1 ; i<ny; ++i)
+ for (unsigned int j=1 ; j<nx; ++j)
{
// coarse grid point number
- unsigned int ncoarse =j+(nx+1)*i;
+ unsigned int ncoarse =j+(nx-1)*i -(nx-1) -1;
// same point on fine grid
- unsigned int nfine = 2*j+(4*nx+2)*i;
+ unsigned int nfine =2*j+(4*nx-2)*i -(2*nx-1) -1;
matrix.set(ncoarse,nfine,1.);
// left
matrix.set(ncoarse,nfine-1,.5);
// lower left
if (i>0)
- matrix.set(ncoarse,nfine-(2*nx+1)-1,.25);
+ matrix.set(ncoarse,nfine-(2*nx-1)-1,.25);
// upper left
if (i<ny)
- matrix.set(ncoarse,nfine+(2*nx+1)-1,.25);
+ matrix.set(ncoarse,nfine+(2*nx-1)-1,.25);
}
// right
if (j<nx)
matrix.set(ncoarse,nfine+1,.5);
// lower right
if (i>0)
- matrix.set(ncoarse,nfine-(2*nx+1)+1,.25);
+ matrix.set(ncoarse,nfine-(2*nx-1)+1,.25);
// upper right
if (i<ny)
- matrix.set(ncoarse,nfine+(2*nx+1)+1,.25);
+ matrix.set(ncoarse,nfine+(2*nx-1)+1,.25);
}
// lower
if (i>0)
- matrix.set(ncoarse,nfine-(2*nx+1),.5);
+ matrix.set(ncoarse,nfine-(2*nx-1),.5);
// upper
if (i<ny)
- matrix.set(ncoarse,nfine+(2*nx+1),.5);
+ matrix.set(ncoarse,nfine+(2*nx-1),.5);
}
}
void
FDMGTransfer::prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ Vector<FLOAT> &dst,
+ const Vector<FLOAT> &src) const
{
Assert((to_level>0) && (to_level<=matrices.size()),
ExcIndexRange(to_level, 0, matrices.size()+1));
void
FDMGTransfer::restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ Vector<FLOAT> &dst,
+ const Vector<FLOAT> &src) const
{
Assert((from_level>0) && (from_level<=matrices.size()),
ExcIndexRange(from_level, 0, matrices.size()+1));
template void FDMatrix::laplacian(SparseMatrix<double>&) const;
template void FDMatrix::laplacian(SparseMatrix<float>&) const;
+template void FDMatrix::gnuplot_print(ostream& s, const Vector<double>& V) const;
+template void FDMatrix::gnuplot_print(ostream& s, const Vector<float>& V) const;
+
template <typename number>
void laplacian(SparseMatrix<number>&) const;
+ template <typename number>
+ void gnuplot_print(ostream&, const Vector<number>&) const;
+
private:
/**
* Number of gridpoints in x-direction.
* function in #MGTranferBase#.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const;
+ Vector<FLOAT> &dst,
+ const Vector<FLOAT> &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<FLOAT> &dst,
+ const Vector<FLOAT> &src) const;
/**
* Exception.
/**
* Prolongation matrices.
*/
- vector<SparseMatrix<float> > matrices;
+ vector<SparseMatrix<FLOAT> > matrices;
/**
* Matrix generator.
* fine to coarse (#vmult#).
*/
void build_matrix(unsigned int nx, unsigned int ny,
- SparseMatrixStruct&, SparseMatrix<float>&);
+ SparseMatrixStruct&, SparseMatrix<FLOAT>&);
};