Commit 1b072751 authored by Benjamin Faigle's avatar Benjamin Faigle
Browse files

- Vector typedefs of the transport model is now inherited (for 2p2cni transport)

permeability tensor (not only scalar) is written out in initialization output
- Upwind direction in the mpfa case is now only set on unique intersection
- Upwind direction in transport considers sign of potential and not information of the fluxdata, as in the adaptive case the latter might not be correct
- Increment in pressure (for derivatives in fluid volume) is now private variable of the pressure model. Change of its standard value might alter test results
reviewed by Philipp

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12400 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 6b0078ab
......@@ -126,8 +126,8 @@ public:
void get1pFluxOnBoundary(EntryType&,
const Intersection&, const CellData&);
//initialize mult-physics-specific pressure model stuff
void initialize()
//initialize multi-physics-specific pressure model constituents
void initialize(bool solveTwice = false)
{
// assign whole domain to most complex subdomain, => 2p
int size = this->problem().gridView().size(0);
......@@ -398,7 +398,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
// compressibility term: 1p domain, so no dv_dp calculated
if (true)
{
Scalar incp = 1e-2;
Scalar& incp = this->incp_;
// numerical derivative of fluid volume with respect to pressure
PhaseVector p_(incp);
......@@ -429,7 +429,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
// dV_dp > 0 is unphysical: Try inverse increment for secant
if (cellDataI.dv_dp()>0)
{
Dune::dinfo <<__FILE__<< "dv_dp still larger 0 after inverting secant"<< std::endl;
Dune::dinfo <<__FILE__<< "dv_dp still larger 0 after inverting secant. regularize"<< std::endl;
cellDataI.dv_dp() *= -1;
}
}
......
......@@ -341,10 +341,16 @@ public:
if(dim >=3)
permZ_[globalIdx] = problem_.spatialParams().intrinsicPermeability(*eIt)[2][2];
}
*permPtrY = permY_;
*permPtrZ = permZ_;
initializationOutputWriter_.attachCellData(*permPtrY, "permeability Y");
initializationOutputWriter_.attachCellData(*permPtrZ, "permeability Z");
if(dim >=2)
{
*permPtrY = permY_;
initializationOutputWriter_.attachCellData(*permPtrY, "permeability Y");
}
if(dim >=3)
{
*permPtrZ = permZ_;
initializationOutputWriter_.attachCellData(*permPtrZ, "permeability Z");
}
}
#endif
......@@ -386,6 +392,8 @@ protected:
Scalar ErrorTermLowerBound_; //!< Handling of error term: lower bound for error dampening
Scalar ErrorTermUpperBound_; //!< Handling of error term: upper bound for error dampening
Scalar incp_ = 1e1; //!< Increment for the volume derivative w.r.t pressure
static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
private:
Problem& problem()
......@@ -777,7 +785,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
if(fabs(massIncrement[compIdx]) < 1e-8 * cellData.density(compIdx))
massIncrement[compIdx] = 1e-8* cellData.density(compIdx); // as phaseIdx = compIdx
}
Scalar incp = 1e-2;
Scalar& incp = incp_;
/**********************************
* c) Secant method for derivatives
......
......@@ -108,10 +108,12 @@ class FVTransport2P2C
typedef Dune::FieldVector<Scalar, NumComponents> ComponentVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
protected:
public:
//! @copydoc FVPressure::EntryType
typedef Dune::FieldVector<Scalar, 2> EntryType;
typedef Dune::FieldVector<Scalar, 2> TimeStepFluxType;
protected:
struct LocalTimesteppingData
{
Dune::FieldVector<EntryType, 2*dim> faceFluxes;
......@@ -167,7 +169,7 @@ public:
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
if(problem().vtkOutputLevel()>=3)
if(problem().vtkOutputLevel()>3)
{
typedef typename GET_PROP(TypeTag, SolutionTypes)::ScalarSolution ScalarSolutionType;
int size = problem_.gridView().size(0);
......@@ -300,14 +302,6 @@ protected:
Scalar accumulatedDt_; //! Current time-interval in sub-time-stepping routine
const Scalar dtThreshold_; //! Threshold for sub-time-stepping routine
std::vector<LocalTimesteppingData> timeStepData_; //! Stores data for sub-time-stepping
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
//! \copydoc Dumux::IMPETProblem::asImp_()
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
void updatedTargetDt_(Scalar &dt);
......@@ -320,6 +314,16 @@ private:
Scalar subCFLFactor_;
bool localTimeStepping_;
int verbosity_;
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
//! \copydoc Dumux::IMPETProblem::asImp_()
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
};
//! \brief Calculate the update vector and determine timestep size
......@@ -705,11 +709,11 @@ void FVTransport2P2C<TypeTag>::getFlux(ComponentVector& fluxEntries,
doUpwinding[phaseIdx] = false;
else // i.e. restrictFluxInTransport == 1
{
//check if harmonic weithing is necessary
if (!wasUpwindCell && (cellDataJ.mobility(phaseIdx) != 0. // check if outflow induce neglected (i.e. mob=0) phase flux
//check if harmonic weighting is necessary
if (potential[phaseIdx] > 0. && (cellDataJ.mobility(phaseIdx) != 0. // check if outflow induce neglected (i.e. mob=0) phase flux
or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementPtrI->father() == neighborPtr->father())))
lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
else if (wasUpwindCell && (cellDataI.mobility(phaseIdx) != 0. // check if inflow induce neglected phase flux
else if (potential[phaseIdx] < 0. && (cellDataI.mobility(phaseIdx) != 0. // check if inflow induce neglected phase flux
or (cellDataI.wasRefined() && cellDataJ.wasRefined() && elementPtrI->father() == neighborPtr->father())))
lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
else
......
......@@ -493,17 +493,21 @@ void FVTransport2P2CAdaptive<TypeTag>::getMpfaFlux(Dune::FieldVector<Scalar, 2>&
if(potential[phaseIdx] > 0.)
{
lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, true);
if(elementI->level()>neighborPointer->level())
cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, true);
}
else if(potential[phaseIdx] < 0.)
{
lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, false);
if(elementI->level()>neighborPointer->level())
cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, false);
}
else
{
doUpwinding[phaseIdx] = false;
if(elementI->level()>neighborPointer->level())
cellDataI.setUpwindCell(intersectionIterator->indexInInside(), contiEqIdx, false);
else
cellDataJ.setUpwindCell(intersectionIterator->indexInOutside(), contiEqIdx, false);
}
}
......@@ -525,10 +529,10 @@ void FVTransport2P2CAdaptive<TypeTag>::getMpfaFlux(Dune::FieldVector<Scalar, 2>&
doUpwinding[phaseIdx] = false;
else // i.e. restrictFluxInTransport == 1
{
//check if harmonic weithing is necessary
if (!cellIwasUpwindCell && cellDataJ.mobility(phaseIdx) != 0.) // check if outflow induce neglected (i.e. mob=0) phase flux
//check if harmonic weighting is necessary
if (potential[phaseIdx] > 0. && cellDataJ.mobility(phaseIdx) != 0.) // check if outflow induce neglected (i.e. mob=0) phase flux
lambda[phaseIdx] = cellDataI.mobility(phaseIdx);
else if (cellIwasUpwindCell && cellDataI.mobility(phaseIdx) != 0.) // check if inflow induce neglected phase flux
else if (potential[phaseIdx] < 0. && cellDataI.mobility(phaseIdx) != 0.) // check if inflow induce neglected phase flux
lambda[phaseIdx] = cellDataJ.mobility(phaseIdx);
else
doUpwinding[phaseIdx] = false;
......
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