diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh index 073f9a196ce5d17b5fada80d266c65cc5ed61491..c5d04d53eeaaf1c3d5035995dcc0636194572903 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh @@ -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)); diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh index deebb8c0434259835ff0ebd5ece10aa7e058c4df..7de45cacf0a19d8f4c152bfc22711ede202b4e13 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh @@ -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); diff --git a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh index 710cdaf306b7dd1c0404812305a3be2f91e0ee22..dc5120e7e3d016cd5d5f7649b0d4f5c94cf1eb31 100644 --- a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh +++ b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh @@ -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)); diff --git a/dumux/decoupled/2p/transport/fv/gravitypart.hh b/dumux/decoupled/2p/transport/fv/gravitypart.hh index ce6b65daad4e9f17d4f898fc17fef776face3896..6af0d9b17e018cdcf08382a248e8f8d41a4014f4 100644 --- a/dumux/decoupled/2p/transport/fv/gravitypart.hh +++ b/dumux/decoupled/2p/transport/fv/gravitypart.hh @@ -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)); diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh index 8082259b157c43fcd822b7fe248fe46cd390eb98..46d244df53ea936f4f6593e46b1bc464674ecd35 100644 --- a/dumux/decoupled/common/fv/fvpressure.hh +++ b/dumux/decoupled/common/fv/fvpressure.hh @@ -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);