From 70f28f45a3637f3dd6c5f4baf0eff1329a674ec6 Mon Sep 17 00:00:00 2001
From: Markus Wolff <markus.wolff@twt-gmbh.de>
Date: Fri, 17 May 2013 13:07:53 +0000
Subject: [PATCH] More than one pressure value can now be fixed inside the
domain if a finite volume pressure model is used. Additionally, the pressure
field is now initialized using the problem.initial() function. Giving
meaningful initializations can improve the convergence of iterative linear
solvers which start the iterations from the values given in the solution
vector
-reviewed by Bernd
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10693 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/decoupled/common/fv/fvpressure.hh | 76 +++++++++++++++++--------
1 file changed, 53 insertions(+), 23 deletions(-)
diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh
index 8752ea840d..5e9a809217 100644
--- a/dumux/decoupled/common/fv/fvpressure.hh
+++ b/dumux/decoupled/common/fv/fvpressure.hh
@@ -22,6 +22,7 @@
// dumux environment
#include "dumux/common/math.hh"
#include <dumux/decoupled/common/pressureproperties.hh>
+#include <map>
/**
* @file
* @brief Finite Volume Diffusion Model
@@ -65,6 +66,11 @@ template<class TypeTag> class FVPressure
typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) RHSVector;
typedef typename GET_PROP_TYPE(TypeTag, PressureSolutionVector) PressureSolution;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
protected:
/*! Type of the vector of entries
@@ -84,6 +90,11 @@ protected:
};
+ enum
+ {
+ pressEqIdx = Indices::pressEqIdx,
+ };
+
//!Initialize the global matrix of the system of equations to solve
void initializeMatrix();
@@ -108,6 +119,21 @@ protected:
const PressureSolution& pressure() const
{ return pressure_;}
+ //!Initialization of the pressure solution vector: Initialization with meaningful values may result in better convergence of the linear solver!
+ void initializePressure()
+ {
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
+ {
+ PrimaryVariables initValues;
+ problem_.initial(initValues, *eIt);
+
+ int globalIdx = problem_.variables().index(*eIt);
+
+ pressure_[globalIdx] = initValues[pressEqIdx];
+ }
+ }
+
public:
/*! \brief Function which calculates the source entry
*
@@ -181,13 +207,13 @@ public:
*/
void initialize()
{
- size_ = problem_.gridView().size(0);//resize to make sure the final grid size (after the problem was completely built) is used!
- A_.setSize(size_, size_);
+ int size = problem_.gridView().size(0);//resize to make sure the final grid size (after the problem was completely built) is used!
+ A_.setSize(size, size);
A_.setBuildMode(Matrix::random);
- f_.resize(size_);
- pressure_.resize(size_);
+ f_.resize(size);
+ pressure_.resize(size);
+ initializePressure();
asImp_().initializeMatrix();// initialize sparse matrix
- pressure_ = 0;
}
/*! \brief Pressure update
@@ -202,6 +228,11 @@ public:
return;
}
+ void calculateVelocity()
+ {
+ DUNE_THROW(Dune::NotImplemented,"Velocity calculation not implemented in pressure model!");
+ }
+
/*! \brief Function for serialization of the pressure field.
*
* Function needed for restart option. Writes the pressure of a grid element to a restart file.
@@ -238,9 +269,7 @@ public:
*/
void setFixPressureAtIndex(Scalar pressure, int globalIdx)
{
- hasFixPressureAtIndex_ = true;
- fixPressureAtIndex_ = pressure;
- idxFixPressureAtIndex_ = globalIdx;
+ fixPressure_.insert(std::make_pair(globalIdx, pressure));
}
/*! \brief Reset the fixed pressure state
@@ -251,9 +280,12 @@ public:
*/
void unsetFixPressureAtIndex(int globalIdx)
{
- hasFixPressureAtIndex_ = false;
- fixPressureAtIndex_ = 0.0;
- idxFixPressureAtIndex_ = 0.0;
+ fixPressure_.erase(globalIdx);
+ }
+
+ void resetFixPressureAtIndex()
+ {
+ fixPressure_.clear();
}
/*! \brief Constructs a FVPressure object
@@ -261,10 +293,7 @@ public:
* \param problem A problem class object
*/
FVPressure(Problem& problem) :
- problem_(problem), size_(problem.gridView().size(0)),
- fixPressureAtIndex_(0),
- idxFixPressureAtIndex_(0),
- hasFixPressureAtIndex_(false)
+ problem_(problem)
{}
private:
@@ -278,7 +307,6 @@ private:
Problem& problem_;
- int size_;
PressureSolution pressure_;
std::string solverName_;
@@ -287,9 +315,7 @@ protected:
Matrix A_;//!<Global stiffnes matrix (sparse matrix which is build by the <tt> initializeMatrix()</tt> function)
RHSVector f_;//!<Right hand side vector
private:
- Scalar fixPressureAtIndex_;
- Scalar idxFixPressureAtIndex_;
- bool hasFixPressureAtIndex_;
+ std::map<int, Scalar> fixPressure_;
};
//!Initialize the global matrix of the system of equations to solve
@@ -470,11 +496,15 @@ void FVPressure<TypeTag>::solve()
std::cout << __FILE__ << ": solve for pressure" << std::endl;
//set a fixed pressure for a certain cell
- if (hasFixPressureAtIndex_)
+ if (fixPressure_.size() > 0)
{
- A_[idxFixPressureAtIndex_] = 0;
- A_[idxFixPressureAtIndex_][idxFixPressureAtIndex_] = 1;
- f_[idxFixPressureAtIndex_] = fixPressureAtIndex_;
+ typename std::map<int, Scalar>::iterator it = fixPressure_.begin();
+ for (;it != fixPressure_.end();++it)
+ {
+ A_[it->first] = 0;
+ A_[it->first][it->first] = 1;
+ f_[it->first] = it->second;
+ }
}
// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
--
GitLab