Skip to content
Snippets Groups Projects
Commit 65aaf74a authored by Christian Engwer's avatar Christian Engwer
Browse files

added test for bug #333 

[[Imported from SVN: r841]]
parent e41170b2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
istl/overlappingschwarz.hh
+ 10
− 8
View file @ 65aaf74a
......@@ -301,6 +301,8 @@ namespace Dune
/** @brief The vector type containing the subdomain to row index mapping. */
typedef std::vector<subdomain_type,allocator> subdomain_vector;
/** @brief The vector type containing the row index to subdomain mapping. */
typedef std::vector<SLList<size_type,TA>,TA> rowtodomain_vector;
enum {
//! \brief The category the precondtioner is part of.
......@@ -322,7 +324,7 @@ namespace Dune
* @param mat The matrix to precondition.
* @param rowToDomain The mapping of the rows onto the domains.
*/
SeqOverlappingSchwarz(const matrix_type& mat, const std::vector<SLList<size_type,TA>,TA>& rowToDomain,
SeqOverlappingSchwarz(const matrix_type& mat, const rowtodomain_vector& rowToDomain,
field_type relaxationFactor=1);
/*!
......@@ -357,7 +359,7 @@ namespace Dune
int maxlength;
void initialize(const std::vector<SLList<size_type,TA>,TA>& rowToDomain, const matrix_type& mat);
void initialize(const rowtodomain_vector& rowToDomain, const matrix_type& mat);
template<typename T>
struct Assigner
......@@ -485,11 +487,11 @@ namespace Dune
}
template<class M, class X, class TM, class TA>
SeqOverlappingSchwarz<M,X,TM,TA>::SeqOverlappingSchwarz(const matrix_type& mat_, const std::vector<SLList<size_type,TA>,TA>& rowToDomain,
SeqOverlappingSchwarz<M,X,TM,TA>::SeqOverlappingSchwarz(const matrix_type& mat_, const rowtodomain_vector& rowToDomain,
field_type relaxationFactor)
: mat(mat_), relax(relaxationFactor)
{
typedef typename std::vector<SLList<size_type,TA>,TA>::const_iterator RowDomainIterator;
typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
typedef typename SLList<size_type,TA>::const_iterator DomainIterator;
#ifdef DUNE_ISTL_WITH_CHECKING
assert(rowToDomain.size()==mat.N());
......@@ -555,7 +557,7 @@ namespace Dune
#endif
// Create a row to subdomain mapping
std::vector<SLList<size_type,TA>,TA> rowToDomain(mat.N());
rowtodomain_vector rowToDomain(mat.N());
int domainId=0;
......@@ -569,7 +571,7 @@ namespace Dune
}
template<class M, class X, class TM, class TA>
void SeqOverlappingSchwarz<M,X,TM,TA>::initialize(const std::vector<SLList<size_type,TA>,TA>& rowToDomain, const matrix_type& mat)
void SeqOverlappingSchwarz<M,X,TM,TA>::initialize(const rowtodomain_vector& rowToDomain, const matrix_type& mat)
{
typedef typename std::vector<SuperMatrixInitializer<matrix_type> >::iterator InitializerIterator;
typedef typename subdomain_vector::const_iterator DomainIterator;
......@@ -591,7 +593,7 @@ namespace Dune
// Set up the supermatrices according to the subdomains
typedef OverlappingSchwarzInitializer<std::vector<SuperMatrixInitializer<matrix_type> >,
std::vector<SLList<size_type,TA>,TA> > Initializer;
rowtodomain_vector > Initializer;
Initializer initializer(initializers, rowToDomain);
copyToSuperMatrix(initializer, mat);
......@@ -675,7 +677,7 @@ namespace Dune
// loop over all Matrix row entries and calculate defect.
typedef typename M::ConstColIterator col_iterator;
typedef typename subdomain_type::const_iterator domain_iterator;
typedef typename subdomain_type::const_iterator domain_iterator;
for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
typename X::block_type tmp;
......
istl/test/.gitignore
+ 1
− 0
View file @ 65aaf74a
......@@ -13,6 +13,7 @@ gmon.out
vectorcommtest
matrixtest
matrixiteratortest
overlappingschwarztest
bcrsbuildtest
superlutest
mv
istl/test/overlappingschwarztest.cc
+ 57
− 22
View file @ 65aaf74a
......@@ -10,6 +10,7 @@
#include <dune/common/timer.hh>
#include <dune/common/sllist.hh>
#include <dune/istl/overlappingschwarz.hh>
int main(int argc, char** argv)
{
......@@ -24,12 +25,12 @@ int main(int argc, char** argv)
typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
typedef Dune::FieldVector<double,BS> VectorBlock;
typedef Dune::BlockVector<VectorBlock> Vector;
typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
typedef Dune::BlockVector<VectorBlock> BVector;
typedef Dune::MatrixAdapter<BCRSMat,BVector,BVector> Operator;
BCRSMat mat;
Operator fop(mat);
Vector b(N*N), x(N*N);
BVector b(N*N), x(N*N);
setupLaplacian(mat,N);
b=0;
......@@ -48,11 +49,14 @@ int main(int argc, char** argv)
int domainsPerDim=(N+domainSize-1)/domainSize;
// set up the overlapping domains
typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector> Schwarz;
typedef Dune::SeqOverlappingSchwarz<BCRSMat,BVector> Schwarz;
typedef Schwarz::subdomain_vector subdomain_vector;
subdomain_vector domains(domainsPerDim*domainsPerDim);
// set up the rowToDomain vector
typedef Schwarz::rowtodomain_vector rowtodomain_vector;
rowtodomain_vector rowToDomain(N*N);
for(int j=0; j < N; ++j)
for(int i=0; i < N; ++i)
......@@ -60,27 +64,41 @@ int main(int argc, char** argv)
int xdomain = i/domainSize;
int ydomain = j/domainSize;
int mainDomain=ydomain*domainsPerDim+xdomain;
domains[mainDomain].insert(j*N+i);
int id=j*N+i;
domains[mainDomain].insert(id);
rowToDomain[id].push_back(mainDomain);
// check left domain
int domain = (i-overlap)/domainSize;
if(domain>=0 && domain<domainsPerDim)
domains[ydomain*domainsPerDim+domain].insert(j*N+i);
{
domains[ydomain*domainsPerDim+domain].insert(id);
rowToDomain[id].push_back(ydomain*domainsPerDim+domain);
}
//check right domain
domain = (i+overlap)/domainSize;
if(domain>=0 && domain<domainsPerDim)
domains[ydomain*domainsPerDim+domain].insert(j*N+i);
{
domains[ydomain*domainsPerDim+domain].insert(id);
rowToDomain[id].push_back(ydomain*domainsPerDim+domain);
}
// check lower domain
domain = (j-overlap)/domainSize;
if(domain>=0 && domain<domainsPerDim)
domains[domain*domainsPerDim+xdomain].insert(j*N+i);
{
domains[domain*domainsPerDim+xdomain].insert(id);
rowToDomain[id].push_back(domain*domainsPerDim+xdomain);
}
//check right domain
domain = (j+overlap)/domainSize;
if(domain>=0 && domain<domainsPerDim)
domains[domain*domainsPerDim+xdomain].insert(j*N+i);
{
domains[domain*domainsPerDim+xdomain].insert(id);
rowToDomain[id].push_back(domain*domainsPerDim+xdomain);
}
}
typedef subdomain_vector::const_iterator iterator;
......@@ -100,38 +118,55 @@ int main(int argc, char** argv)
Dune::printvector(std::cout, x, std::string("X"), std::string("X"));
}
Dune::SeqOverlappingSchwarz<BCRSMat,Vector> prec(mat, domains, 1);
Dune::Timer watch;
watch.reset();
Dune::LoopSolver<Vector> solver(fop, prec, 1e-2,100,2);
Dune::InverseOperatorResult res;
std::cout<<"Additive Schwarz (domains vector)"<<std::endl;
// b=0;
// x=100;
// setBoundary(x,b,N);
std::cout<<"Additive Schwarz"<<std::endl;
solver.apply(x,b, res);
Dune::SeqOverlappingSchwarz<BCRSMat,BVector> prec0(mat, domains, 1);
Dune::LoopSolver<BVector> solver0(fop, prec0, 1e-2,100,2);
solver0.apply(x,b, res);
Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> prec1(mat,domains, 1);
Dune::LoopSolver<Vector> solver1(fop, prec1, 1e-2,100,2);
std::cout << "Multiplicative Schwarz (domains vector)"<<std::endl;
b=0;
x=100;
//setBoundary(x,b,N);
std::cout << "Multiplicative Schwarz"<<std::endl;
Dune::SeqOverlappingSchwarz<BCRSMat,BVector,Dune::MultiplicativeSchwarzMode> prec1(mat, domains, 1);
Dune::LoopSolver<BVector> solver1(fop, prec1, 1e-2,100,2);
solver1.apply(x,b, res);
Dune::SeqSOR<BCRSMat,Vector,Vector> sor(mat, 1,1);
Dune::LoopSolver<Vector> solver2(fop, sor, 1e-2,100,2);
std::cout<<"Additive Schwarz (rowToDomain vector)"<<std::endl;
b=0;
x=100;
// setBoundary(x,b,N);
Dune::SeqOverlappingSchwarz<BCRSMat,BVector> prec2(mat, rowToDomain, 1);
Dune::LoopSolver<BVector> solver2(fop, prec2, 1e-2,100,2);
solver2.apply(x,b, res);
std::cout << "Multiplicative Schwarz (rowToDomain vector)"<<std::endl;
b=0;
x=100;
//setBoundary(x,b,N);
Dune::SeqOverlappingSchwarz<BCRSMat,BVector,Dune::MultiplicativeSchwarzMode> prec3(mat, rowToDomain, 1);
Dune::LoopSolver<BVector> solver3(fop, prec3, 1e-2,100,2);
solver3.apply(x,b, res);
std::cout << "SOR"<<std::endl;
solver2.apply(x,b, res);
b=0;
x=100;
//setBoundary(x,b,N);
Dune::SeqSOR<BCRSMat,BVector,BVector> sor(mat, 1,1);
Dune::LoopSolver<BVector> solver4(fop, sor, 1e-2,100,2);
solver4.apply(x,b, res);
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment