Commit 1edee536 authored by Markus Wolff's avatar Markus Wolff
Browse files

same clean-up and small modifications

   - clean-up in 2p models
   - added possibility to fix a pressure somewhere in the domain



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7430 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 45fa65ee
......@@ -163,6 +163,9 @@ public:
}
// std::cout<<"Pressure defect = "<<pressureNorm<<"; "<<numIter<<" Iterations needed for initial pressure field"<<std::endl;
}
storePressureSolution();
return;
}
......@@ -334,7 +337,7 @@ public:
if (!compressibility_)
{
const Element& element = *(problem_.gridView().template end<0> ());
const Element& element = *(problem_.gridView().template begin<0> ());
FluidState fluidState;
fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
......
......@@ -125,7 +125,7 @@ public:
if (!compressibility_)
{
const Element& element = *(problem_.gridView().template end<0> ());
const Element& element = *(problem_.gridView().template begin<0> ());
FluidState fluidState;
fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
......@@ -654,6 +654,10 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
velocityNW /= density_[nPhaseIdx];
}
//store potential gradients for further calculations
cellDataI.fluxData().setPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]);
cellDataI.fluxData().setPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]);
cellDataI.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
cellDataI.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNW);
cellDataI.fluxData().setVelocityMarker(isIndex);
......
......@@ -340,7 +340,7 @@ public:
if (!compressibility_)
{
const Element& element = *(problem_.gridView().template end<0> ());
const Element& element = *(problem_.gridView().template begin<0> ());
FluidState fluidState;
fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
......
......@@ -179,7 +179,7 @@ public:
GravityPart (Problem& problem)
: ConvectivePart<TypeTag>(problem), problem_(problem), preComput_(GET_PROP_VALUE(TypeTag, PrecomputedConstRels))
{
const Element& element = *(problem_.gridView().template end<0> ());
const Element& element = *(problem_.gridView().template begin<0> ());
FluidState fluidState;
fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
......
......@@ -142,13 +142,32 @@ public:
}
//@}
/*! set a pressure to be fixed at a certain cell */
void setFixPressureAtIndex(Scalar pressure, int globalIdx)
{
hasFixPressureAtIndex_ = true;
fixPressureAtIndex_ = pressure;
idxFixPressureAtIndex_ = globalIdx;
}
/*! unset a fixed pressure at a certain cell */
void unsetFixPressureAtIndex(int globalIdx)
{
hasFixPressureAtIndex_ = false;
fixPressureAtIndex_ = 0.0;
idxFixPressureAtIndex_ = 0.0;
}
//! Constructs a FVPressure object
/**
* \param problem a problem class object
*/
FVPressure(Problem& problem) :
problem_(problem), size_(problem.gridView().size(0)),
pressure_(size_), A_(size_, size_, (2 * dim + 1) * size_, Matrix::random), f_(size_)
pressure_(size_), A_(size_, size_, (2 * dim + 1) * size_, Matrix::random), f_(size_),
fixPressureAtIndex_(0),
idxFixPressureAtIndex_(0),
hasFixPressureAtIndex_(false)
{}
private:
......@@ -170,6 +189,9 @@ private:
protected:
Matrix A_;
RHSVector f_;
Scalar fixPressureAtIndex_;
Scalar idxFixPressureAtIndex_;
bool hasFixPressureAtIndex_;
};
//! initializes the matrix to store the system of equations
......@@ -325,6 +347,14 @@ void FVPressure<TypeTag>::solve()
if (verboseLevelSolver)
std::cout << __FILE__ <<": solve for pressure" << std::endl;
//set a fixed pressure for a certain cell
if (hasFixPressureAtIndex_)
{
A_[idxFixPressureAtIndex_] = 0;
A_[idxFixPressureAtIndex_][idxFixPressureAtIndex_] = 1;
f_[idxFixPressureAtIndex_] = fixPressureAtIndex_;
}
Solver solver(problem_);
solver.solve(A_, pressure_, f_);
// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
......
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