Skip to content
Snippets Groups Projects
Commit cf51dc1c authored by Markus Blatt's avatar Markus Blatt
Browse files

Made epsilon in setAnisotropic changeable. 

Removed N as template paramter as it is of now actual use.

[[Imported from SVN: r401]]
parent 895c3355
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
istl/paamg/test/amgtest.cc
+ 61
− 30
View file @ cf51dc1c
......@@ -10,7 +10,11 @@ int main(int argc, char** argv)
{
MPI_Init(&argc, &argv);
const int BS=1, N=1000;
const int BS=1;
int N=250;
if(argc>1)
N = atoi(argv[1]);
int procs, rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
......@@ -24,11 +28,20 @@ int main(int argc, char** argv)
int n;
BCRSMat mat = setupAnisotropic2d<N,BS>(indices, &n);
BCRSMat mat = setupAnisotropic2d<BS>(N, indices, &n, 1);
/*
if(N<20)
Dune::printmatrix(std::cout, mat, "A", "row");
*/
Vector b(indices.size()), x(indices.size());
b=1;
b=0;
x=100;
Dune::MatrixAdapter<BCRSMat,Vector,Vector> fop(mat);
Dune::Timer watch;
RemoteIndices remoteIndices(indices,indices,MPI_COMM_WORLD);
remoteIndices.rebuild<false>();
......@@ -49,40 +62,58 @@ int main(int argc, char** argv)
typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
Criterion;
Dune::Timer watch;
watch.reset();
Criterion criterion(1,100);
criterion.setMaxDistance(3);
Criterion criterion(15,10);
criterion.setMaxDistance(2);
hierarchy.build(criterion);
std::cout<<"Building hierarchy took "<<watch.elapsed()<<" seconds"<<std::endl;
double buildtime = watch.elapsed();
std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
if(N<20)
Dune::printmatrix(std::cout, hierarchy.matrices().coarsest()->matrix(), "A", "row");
hierarchy.coarsenVector(vh);
typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
typedef Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
SmootherArgs smootherArgs;
smootherArgs.iterations = 2;
Dune::MatrixAdapter<BCRSMat,Vector,Vector> op(hierarchy.matrices().coarsest()->matrix());
Dune::SeqSSOR<BCRSMat,Vector,Vector> cssor(hierarchy.matrices().coarsest()->matrix(),1,1.0);
Dune::SeqSSOR<BCRSMat,Vector,Vector> ssor(hierarchy.matrices().finest()->matrix(),1,1.0);
typedef Dune::LoopSolver<Vector> CoarseSolver;
CoarseSolver csolver(op,cssor,1E-12,8000,0);
Dune::SeqScalarProduct<Vector> sp;
typedef Dune::Amg::AMG<MHierarchy,Vector,Vector,CoarseSolver,Smoother> AMG;
AMG amg(hierarchy, csolver, smootherArgs, 1, 1);
Dune::CGSolver<Vector> amgCG(fop,amg,10e-8,80,2);
/*
hierarchy.coarsenVector(vh);
typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
typedef Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
SmootherArgs smootherArgs;
Dune::MatrixAdapter<BCRSMat,Vector,Vector> op(hierarchy.matrices().coarsest()->matrix());
Dune::SeqSSOR<BCRSMat,Vector,Vector> cssor(hierarchy.matrices().coarsest()->matrix(),1,1.0);
Dune::SeqSSOR<BCRSMat,Vector,Vector> ssor(hierarchy.matrices().finest()->matrix(),1,1.0);
typedef Dune::LoopSolver<Vector> CoarseSolver;
CoarseSolver csolver(op,cssor,1E-12,8000,0);
Dune::SeqScalarProduct<Vector> sp;
typedef Dune::Amg::AMG<MHierarchy,Vector,Vector,CoarseSolver,Smoother> AMG;
AMG amg(hierarchy, csolver, smootherArgs, 1, 1);
typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
Operator fop(hierarchy.matrices().finest()->matrix());
Dune::CGSolver<Vector> amgCG(fop,amg,10e-8,8000,2);
Dune::CGSolver<Vector> cg(fop,ssor,10e-8,8000,2);
Dune::SeqSSOR<BCRSMat,Vector,Vector> fssor(mat,1,1.0);
Dune::CGSolver<Vector> cg(fop,fssor,10e-8,8,2);
*/
watch.reset();
Dune::InverseOperatorResult r;
amgCG.apply(x,b,r);
watch.reset();
Dune::InverseOperatorResult r;
amgCG.apply(x,b,r);
double solvetime = watch.elapsed();
std::cout<<"AMG solving took "<<solvetime<<" seconds"<<std::endl;
std::cout<<"AMG solving took "<<watch.elapsed()<<" seconds"<<std::endl;
std::cout<<"AMG building took "<<(buildtime/r.elapsed*r.iterations)<<" iterations"<<std::endl;
std::cout<<"AMG building together with slving took "<<buildtime+solvetime<<std::endl;
/*
watch.reset();
cg.apply(x,b,r);
......
istl/paamg/test/anisotropic.hh
+ 16
− 17
View file @ cf51dc1c
......@@ -16,20 +16,20 @@ typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
typedef Dune::Interface<ParallelIndexSet> Interface;
typedef Dune::BufferedCommunicator<ParallelIndexSet> Communicator;
template<int N, class M>
void setupPattern(M& mat, ParallelIndexSet& indices, int overlapStart, int overlapEnd,
template<class M>
void setupPattern(int N, M& mat, ParallelIndexSet& indices, int overlapStart, int overlapEnd,
int start, int end);
template<int N, class M>
void fillValues(M& mat, int overlapStart, int overlapEnd, int start, int end);
template<class M>
void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end);
template<int N, int BS>
Dune::BCRSMatrix<Dune::FieldMatrix<double,BS,BS> > setupAnisotropic2d(ParallelIndexSet& indices, int *nout);
template<int BS>
Dune::BCRSMatrix<Dune::FieldMatrix<double,BS,BS> > setupAnisotropic2d(int N, ParallelIndexSet& indices, int *nout, double eps=1.0);
template<int N, class M>
void setupPattern(M& mat, ParallelIndexSet& indices, int overlapStart, int overlapEnd,
template<class M>
void setupPattern(int N, M& mat, ParallelIndexSet& indices, int overlapStart, int overlapEnd,
int start, int end)
{
int n = overlapEnd - overlapStart;
......@@ -43,7 +43,7 @@ void setupPattern(M& mat, ParallelIndexSet& indices, int overlapStart, int overl
GridFlag flag = owner;
bool isPublic = false;
if(i<start || i>= end)
if((i<start && i > 0) || (i>= end && i < N-1))
flag=overlap;
if(i<start+1 || i>= end-1)
......@@ -76,11 +76,10 @@ void setupPattern(M& mat, ParallelIndexSet& indices, int overlapStart, int overl
indices.endResize();
}
template<int N, class M>
void fillValues(M& mat, int overlapStart, int overlapEnd, int start, int end)
template<class M>
void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end, double eps)
{
typedef typename M::block_type Block;
double eps=.01;
Block dval = 0, bone=0, bmone=0, beps=0;
for(typename Block::RowIterator b = dval.begin(); b != dval.end(); ++b)
......@@ -106,7 +105,7 @@ void fillValues(M& mat, int overlapStart, int overlapEnd, int start, int end)
ColIterator endi = (*i).end();
if(x<start || x >= end) {
if(x<start || x >= end || x==0 || x==N-1 || y==0 || y==N-1) {
// overlap node is dirichlet border
ColIterator j = (*i).begin();
......@@ -129,8 +128,8 @@ void fillValues(M& mat, int overlapStart, int overlapEnd, int start, int end)
}
}
template<int N, int BS>
Dune::BCRSMatrix<Dune::FieldMatrix<double,BS,BS> > setupAnisotropic2d(ParallelIndexSet& indices, int *nout)
template<int BS>
Dune::BCRSMatrix<Dune::FieldMatrix<double,BS,BS> > setupAnisotropic2d(int N, ParallelIndexSet& indices, int *nout, double eps)
{
int procs, rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
......@@ -172,8 +171,8 @@ Dune::BCRSMatrix<Dune::FieldMatrix<double,BS,BS> > setupAnisotropic2d(ParallelIn
BCRSMat mat(noKnown*N, noKnown*N, noKnown*N*5, BCRSMat::row_wise);
setupPattern<N>(mat, indices, overlapStart, overlapEnd, start, end);
fillValues<N>(mat, overlapStart, overlapEnd, start, end);
setupPattern(N, mat, indices, overlapStart, overlapEnd, start, end);
fillValues(N, mat, overlapStart, overlapEnd, start, end, eps);
// Dune::printmatrix(std::cout,mat,"aniso 2d","row",9,1);
......
istl/paamg/test/galerkintest.cc
+ 2
− 2
View file @ cf51dc1c
......@@ -30,7 +30,7 @@ void testCoarsenIndices()
typedef Dune::BCRSMatrix<Block> BCRSMat;
int n;
BCRSMat mat = setupAnisotropic2d<N,BS>(indices, &n);
BCRSMat mat = setupAnisotropic2d<BS>(N, indices, &n);
RemoteIndices remoteIndices(indices,indices,MPI_COMM_WORLD);
remoteIndices.rebuild<false>();
......@@ -128,6 +128,6 @@ void testCoarsenIndices()
int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
testCoarsenIndices<5,1>();
testCoarsenIndices<4,1>();
MPI_Finalize();
}
istl/paamg/test/hierarchytest.cc
+ 1
− 1
View file @ cf51dc1c
......@@ -25,7 +25,7 @@ int main(int argc, char** argv)
int n;
BCRSMat mat = setupAnisotropic2d<N,BS>(indices, &n);
BCRSMat mat = setupAnisotropic2d<BS>(N, indices, &n);
Vector b(indices.size());
RemoteIndices remoteIndices(indices,indices,MPI_COMM_WORLD);
......
istl/test/matrixutilstest.cc
+ 2
− 61
View file @ cf51dc1c
......@@ -4,66 +4,7 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/stdstreams.hh>
template<int N, class B>
void setupSparsityPattern(Dune::BCRSMatrix<B>& A)
{
for (typename Dune::BCRSMatrix<B>::CreateIterator i = A.createbegin(); i != A.createend(); ++i) {
int x = i.index()%N; // x coordinate in the 2d field
int y = i.index()/N; // y coordinate in the 2d field
if(y>0)
// insert lower neighbour
i.insert(i.index()-N);
if(x>0)
// insert left neighbour
i.insert(i.index()-1);
// insert diagonal value
i.insert(i.index());
if(x<N-1)
//insert right neighbour
i.insert(i.index()+1);
if(y<N-1)
// insert upper neighbour
i.insert(i.index()+N);
}
}
template<int N, class B>
void setupLaplacian(Dune::BCRSMatrix<B>& A)
{
setupSparsityPattern<N>(A);
B diagonal = 0, bone=0, beps=0;
for(typename B::RowIterator b = diagonal.begin(); b != diagonal.end(); ++b)
b->operator[](b.index())=4;
for(typename B::RowIterator b=bone.begin(); b != bone.end(); ++b)
b->operator[](b.index())=-1.0;
for (typename Dune::BCRSMatrix<B>::RowIterator i = A.begin(); i != A.end(); ++i) {
int x = i.index()%N; // x coordinate in the 2d field
int y = i.index()/N; // y coordinate in the 2d field
i->operator[](i.index())=diagonal;
if(y>0)
i->operator[](i.index()-N)=bone;
if(y<N-1)
i->operator[](i.index()+N)=bone;
if(x>0)
i->operator[](i.index()-1)=bone;
if(x < N-1)
i->operator[](i.index()+1)=bone;
}
}
#include "laplacian.hh"
int main(int argc, char** argv)
{
......@@ -81,7 +22,7 @@ int main(int argc, char** argv)
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > BMatrix;
BMatrix laplace(N*N,N*N, N*N*5, BMatrix::row_wise);
setupLaplacian<N>(laplace);
setupLaplacian2d<N>(laplace);
if(N*N*5-4*2-(N-2)*4!=countNonZeros(laplace)) {
++ret;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment