diff --git a/CMakeLists.txt b/CMakeLists.txt
index 73db5f224332a593b49df066e80325e9d84deb73..1243b78e2e59956aaf773cce92cc85e4c202b85f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -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
diff --git a/dumux/boxmodels/common/pdelabboxassembler.hh b/dumux/boxmodels/common/pdelabboxassembler.hh
index 79671d5f5f3547c146f86d0752aae40e09a28e0c..28e1356eb19135928d7db4ca6b47d0c48c58115a 100644
--- a/dumux/boxmodels/common/pdelabboxassembler.hh
+++ b/dumux/boxmodels/common/pdelabboxassembler.hh
@@ -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)
diff --git a/dumux/common/pdelabpreconditioner.hh b/dumux/common/pdelabpreconditioner.hh
index e16846d9bfbca75f4a16ef311424520b4a428996..211d7384c505823c9e2d474e6b95cbf8687b947a 100644
--- a/dumux/common/pdelabpreconditioner.hh
+++ b/dumux/common/pdelabpreconditioner.hh
@@ -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);
diff --git a/dumux/common/start.hh b/dumux/common/start.hh
index c882783e68d8cefd01ae2a017e94deba00acf0f1..9f5df7b5a935bef94eb73b485e3bf21a405f3a75 100644
--- a/dumux/common/start.hh
+++ b/dumux/common/start.hh
@@ -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);
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index ff5ab5413990abee9273549431f1ed1a69574f1d..18eceed428d4486c3ea2a8516de83dc2ca0ecc6b 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -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)