From 868964a0ed9b28bcee430e8ef06a9dd775b611c1 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 11 Aug 2010 09:02:01 +0000
Subject: [PATCH] made matrix types consistent
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4049 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
CMakeLists.txt | 2 +-
dumux/boxmodels/common/pdelabboxassembler.hh | 48 ++++++++++++++++++--
dumux/common/pdelabpreconditioner.hh | 7 +--
dumux/common/start.hh | 2 +-
dumux/nonlinear/newtoncontroller.hh | 13 +++---
5 files changed, 57 insertions(+), 15 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 73db5f2243..1243b78e2e 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 79671d5f5f..28e1356eb1 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 e16846d9bf..211d7384c5 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 c882783e68..9f5df7b5a9 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 ff5ab54139..18eceed428 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)
--
GitLab