Commit 868964a0 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

made matrix types consistent

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4049 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 73af2ada
......@@ -81,7 +81,7 @@ set(DumuxIncludeDirectories
##############
# set appropriate compiler flags for debug/release compilation modes
add_definitions("-std=c++0x -Wall -Wno-sign-compare")
add_definitions("-std=c++0x -Wall -Wno-sign-compare --no-strict-aliasing")
if(CMAKE_BUILD_TYPE STREQUAL "debug")
# debug mode
......
......@@ -169,20 +169,60 @@ public:
}
reuseMatrix_ = false;
// set the entries for the ghost nodes
// calculate the global residual
residual_ = 0;
gridOperatorSpace_->residual(u, residual_);
typedef typename Matrix::block_type BlockType;
// set the entries for the ghost nodes
BlockType Id(0.0);
for (int i=0; i < BlockType::rows; ++i)
Id[i][i] = 1.0;
for (int i=0; i < ghostIndices_.size(); ++i) {
int globI = ghostIndices_[i];
(*matrix_)[globI] = 0;
(*matrix_)[globI][globI] = Id;
residual_[globI] = 0;
u[globI] = 0;
}
// typedef typename Matrix::RowIterator RowIterator;
// typedef typename Matrix::ColIterator ColIterator;
// const typename Matrix::block_type::size_type rowsInBlock = Matrix::block_type::rows;
// const typename Matrix::block_type::size_type colsInBlock = Matrix::block_type::cols;
// Scalar diagonalEntry[rowsInBlock];
// RowIterator endIBlock = matrix_->end();
// for (RowIterator iBlock = matrix_->begin(); iBlock != endIBlock; ++iBlock) {
// BlockType &diagBlock = (*iBlock)[iBlock.index()];
//
// for (int i = 0; i < rowsInBlock; ++i) {
// diagonalEntry[i] = 0;
// for (int j = 0; j < colsInBlock; ++j) {
// diagonalEntry[i] = std::max(diagonalEntry[i],
// std::abs(diagBlock[i][j]));
// }
//
// if (diagonalEntry[i] < 1e-14)
// diagonalEntry[i] = 1.0;
// }
//
// // divide right-hand side
// for (int i = 0; i < rowsInBlock; i++) {
// (residual_)[iBlock.index()][i] /= diagonalEntry[i];
// }
//
// // divide row of the jacobian
// ColIterator endJBlock = iBlock->end();
// for (ColIterator jBlock = iBlock->begin(); jBlock != endJBlock; ++jBlock) {
// for (int i = 0; i < rowsInBlock; i++) {
// for (int j = 0; j < colsInBlock; j++) {
// (*jBlock)[i][j] /= diagonalEntry[i];
// }
// }
// }
// }
// calculate the global residual
residual_ = 0;
gridOperatorSpace_->residual(u, residual_);
}
void setMatrixReuseable(bool yesno = true)
......
......@@ -344,7 +344,7 @@ public:
typedef Dune::SeqILU0<JacobianMatrix,SolVector,RhsVector> SeqPreCond;
JacobianMatrix B(A);
exchanger_.sumEntries(B);
SeqPreCond seqPreCond(B, 0.9);
SeqPreCond seqPreCond(B, 1.0);
typedef Dune::PDELab::NonoverlappingOperator<GridFunctionSpace,JacobianMatrix,SolVector,RhsVector> POP;
POP pop(gfs,A,phelper);
......@@ -387,6 +387,7 @@ class ISTLBackend_NoOverlap_Loop_Pardiso
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
typedef Dune::PDELab::ParallelISTLHelper<GridFunctionSpace> PHELPER;
public:
......@@ -426,8 +427,8 @@ public:
\param[in] r right hand side
\param[in] reduction to be achieved
*/
template<class JacobianMatrix, class SolVector, class RhsVector>
void apply(JacobianMatrix& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
template<class M, class SolVector, class RhsVector>
void apply(M& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
{
typedef Dune::SeqPardiso<JacobianMatrix,SolVector,RhsVector> SeqPreCond;
JacobianMatrix B(A);
......
......@@ -88,7 +88,7 @@ int startFromDGF(int argc, char **argv)
// create grid
// -> load the grid from file
GridPointer gridPtr = GridPointer(dgfFileName);
GridPointer gridPtr(dgfFileName);
if (mpiHelper.size() > 1)
gridPtr->loadBalance();
Dune::gridinfo(*gridPtr);
......
......@@ -26,10 +26,10 @@
#include <dumux/common/exceptions.hh>
#include <dune/istl/overlappingschwarz.hh>
#include <dune/istl/schwarz.hh>
//#include <dune/istl/schwarz.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include "dune/istl/owneroverlapcopy.hh"
//#include "dune/istl/owneroverlapcopy.hh"
#include <dune/istl/io.hh>
......@@ -182,6 +182,7 @@ class NewtonController
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
......@@ -333,8 +334,8 @@ public:
* Throws Dumux::NumericalProblem if the linear solver didn't
* converge.
*/
template <class Matrix, class Vector>
void newtonSolveLinear(const Matrix &A,
template <class Vector>
void newtonSolveLinear(const JacobianMatrix &A,
Vector &u,
const Vector &b)
{
......@@ -535,8 +536,8 @@ protected:
};
template <class Matrix, class Vector>
void solveLinear_(const Matrix &A,
template <class Vector>
void solveLinear_(const JacobianMatrix &A,
Vector &x,
const Vector &b,
Scalar residReduction)
......
Supports Markdown
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