diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh index 600a1bb4b3e17c14ea11dfe4baa7fc6c98f99d3f..6e0f4c9f650996a10f31d740fa580ed85b59d906 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh @@ -551,7 +551,7 @@ void FVPressure2P<TypeTag>::getStorage(EntryType& entry, const Element& element */ template<class TypeTag> void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& intersection - , const CellData& cellDataI, const bool first) + , const CellData& cellData, const bool first) { ElementPointer elementI = intersection.inside(); ElementPointer elementJ = intersection.outside(); @@ -563,17 +563,17 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters const GlobalPosition& globalPosJ = elementJ->geometry().center(); // get mobilities and fractional flow factors - Scalar lambdaWI = cellDataI.mobility(wPhaseIdx); - Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx); - Scalar fractionalWI = cellDataI.fracFlowFunc(wPhaseIdx); - Scalar fractionalNWI = cellDataI.fracFlowFunc(nPhaseIdx); + Scalar lacellDatacellData.mobility(wPhaseIdx); + Scalar lambdaNWI = cellData.mobility(nPhaseIdx); + Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx); + Scalar fractionalNWI = cellData.fracFlowFunc(nPhaseIdx); Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx); Scalar lambdaNWJ = cellDataJ.mobility(nPhaseIdx); Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx); Scalar fractionalNWJ = cellDataJ.fracFlowFunc(nPhaseIdx); // get capillary pressure - Scalar pcI = cellDataI.capillaryPressure(); + Scalar pcI = cellData.capillaryPressure(); Scalar pcJ = cellDataJ.capillaryPressure(); //get face index @@ -604,8 +604,8 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters Scalar rhoMeanNW = 0; if (compressibility_) { - rhoMeanW = 0.5 * (cellDataI.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)); - rhoMeanNW = 0.5 * (cellDataI.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)); + rhoMeanW = 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)); + rhoMeanNW = 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)); } Scalar fMeanW = 0.5 * (fractionalWI + fractionalWJ); Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWJ); @@ -617,25 +617,25 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters //if we are at the very first iteration we can't calculate phase potentials if (!first) { - potentialW = cellDataI.fluxData().potential(wPhaseIdx, isIndexI); - potentialNW = cellDataI.fluxData().potential(nPhaseIdx, isIndexI); + potentialW = cellData.fluxData().potential(wPhaseIdx, isIndexI); + potentialNW = cellData.fluxData().potential(nPhaseIdx, isIndexI); if (compressibility_) { - density_[wPhaseIdx] = (potentialW > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); - density_[nPhaseIdx] = (potentialNW > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); + density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); + density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); density_[wPhaseIdx] = (potentialW == 0.) ? rhoMeanW : density_[wPhaseIdx]; density_[nPhaseIdx] = (potentialNW == 0.) ? rhoMeanNW : density_[nPhaseIdx]; } - potentialW = cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx); - potentialNW = cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx); + potentialW = cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx); + potentialNW = cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx); if (pressureType_ == pglobal) { - potentialW = (cellDataI.globalPressure() - cellDataJ.globalPressure() - fMeanNW * (pcI - pcJ)); - potentialNW = (cellDataI.globalPressure() - cellDataJ.globalPressure() + fMeanW * (pcI - pcJ)); + potentialW = (cellData.globalPressure() - cellDataJ.globalPressure() - fMeanNW * (pcI - pcJ)); + potentialNW = (cellData.globalPressure() - cellDataJ.globalPressure() + fMeanW * (pcI - pcJ)); } potentialW += density_[wPhaseIdx] * (distVec * gravity_); @@ -650,8 +650,8 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters if (compressibility_) { - density_[wPhaseIdx] = (potentialW > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); - density_[nPhaseIdx] = (potentialNW > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); + density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); + density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); density_[wPhaseIdx] = (potentialW == 0) ? rhoMeanW : density_[wPhaseIdx]; density_[nPhaseIdx] = (potentialNW == 0) ? rhoMeanNW : density_[nPhaseIdx]; diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh index 3a7c07a4857b7f400d9a3abfbd1ca5d3d501cc84..0e9cabb59d21767e6be53fbfdf4dcc3a143160c6 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh @@ -199,19 +199,19 @@ private: * \copydetails FVPressure2P::getFlux(EntryType&,const Intersection&,const CellData&,const bool) * * Implementation of the getFlux() function for cell-cell interfaces with hanging nodes. - * In case of hanging nodes the function does not return a vector of entries but directly manipulates the global matrix! + * In case of hanging nodes the function does not return a vector of entry but directly manipulates the global matrix! * */ template<class TypeTag> -void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entries, const Intersection& intersection - , const CellData& cellDataI, const bool first) +void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entry, const Intersection& intersection + , const CellData& cellData, const bool first) { ElementPointer elementI = intersection.inside(); ElementPointer elementJ = intersection.outside(); if (elementI->level() == elementJ->level()) { - ParentType::getFlux(entries, intersection, cellDataI, first); + ParentType::getFlux(entry, intersection, cellData, first); } // hanging node situation: neighbor has higher level else if (elementJ->level() == elementI->level() + 1) @@ -226,17 +226,17 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entries, const Intersecti int globalIdxJ = problem_.variables().index(*elementJ); // get mobilities and fractional flow factors - Scalar lambdaWI = cellDataI.mobility(wPhaseIdx); - Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx); - Scalar fractionalWI = cellDataI.fracFlowFunc(wPhaseIdx); - Scalar fractionalNWI = cellDataI.fracFlowFunc(nPhaseIdx); + Scalar lambdaWI = cellData.mobility(wPhaseIdx); + Scalar lambdaNWI = cellData.mobility(nPhaseIdx); + Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx); + Scalar fractionalNWI = cellData.fracFlowFunc(nPhaseIdx); Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx); Scalar lambdaNWJ = cellDataJ.mobility(nPhaseIdx); Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx); Scalar fractionalNWJ = cellDataJ.fracFlowFunc(nPhaseIdx); // get capillary pressure - Scalar pcI = cellDataI.capillaryPressure(); + Scalar pcI = cellData.capillaryPressure(); Scalar pcJ = cellDataJ.capillaryPressure(); //get face index @@ -334,10 +334,10 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entries, const Intersecti if (compressibility_) { - rhoMeanWIJ = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l; - rhoMeanNWIJ = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l; - rhoMeanWIK = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l; - rhoMeanNWIK = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l; + rhoMeanWIJ = (lJ * cellData.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l; + rhoMeanNWIJ = (lJ * cellData.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l; + rhoMeanWIK = (lJ * cellData.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l; + rhoMeanNWIK = (lJ * cellData.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l; } Scalar densityWIJ = density_[wPhaseIdx]; @@ -373,28 +373,28 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entries, const Intersecti { Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2; - potentialWIJ = (cellDataI.globalPressure() - fMeanNWIJ * pcI - (pressJK - fMeanNWIJ * pcJK)) / l + potentialWIJ = (cellData.globalPressure() - fMeanNWIJ * pcI - (pressJK - fMeanNWIJ * pcJK)) / l + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng; - potentialNWIJ = (cellDataI.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l + potentialNWIJ = (cellData.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng; - potentialWIK = (cellDataI.globalPressure() - fMeanNWIK * pcI - (pressJK - fMeanNWIK * pcJK)) / l + potentialWIK = (cellData.globalPressure() - fMeanNWIK * pcI - (pressJK - fMeanNWIK * pcJK)) / l + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng; - potentialNWIK = (cellDataI.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l + potentialNWIK = (cellData.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng; break; } default: { - potentialWIJ = (cellDataI.pressure(wPhaseIdx) + potentialWIJ = (cellData.pressure(wPhaseIdx) - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng; - potentialNWIJ = (cellDataI.pressure(nPhaseIdx) + potentialNWIJ = (cellData.pressure(nPhaseIdx) - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng; - potentialWIK = (cellDataI.pressure(wPhaseIdx) + potentialWIK = (cellData.pressure(wPhaseIdx) - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng; - potentialNWIK = (cellDataI.pressure(nPhaseIdx) + potentialNWIK = (cellData.pressure(nPhaseIdx) - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng; break; @@ -410,58 +410,58 @@ void FVPressure2PAdaptive<TypeTag>::getFlux(EntryType& entries, const Intersecti if (compressibility_) { - densityWIJ = (potentialWIJ > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); - densityNWIJ = (potentialNWIJ > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); + densityWIJ = (potentialWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); + densityNWIJ = (potentialNWIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ; densityNWIJ = (potentialNWIJ == 0) ? rhoMeanNWIJ : densityNWIJ; - densityWIK = (potentialWIK > 0.) ? cellDataI.density(wPhaseIdx) : cellDataK.density(wPhaseIdx); - densityNWIK = (potentialNWIK > 0.) ? cellDataI.density(nPhaseIdx) : densityNWIK; + densityWIK = (potentialWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx); + densityNWIK = (potentialNWIK > 0.) ? cellData.density(nPhaseIdx) : densityNWIK; densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK; densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK; } // update diagonal entry and right hand side - entries[matrix] = (lambdaWIJ + lambdaNWIJ) * kMean / l * faceArea; - entries[rhs] = faceArea * lambdaWIJ * kMean * ng + entry[matrix] = (lambdaWIJ + lambdaNWIJ) * kMean / l * faceArea; + entry[rhs] = faceArea * lambdaWIJ * kMean * ng * ((densityWIJ) - (lJ / l) * (kI + kK) / kI * (densityWIK - densityWIJ) / 2); - entries[rhs] += faceArea * lambdaNWIJ * kMean * ng + entry[rhs] += faceArea * lambdaNWIJ * kMean * ng * ((densityNWIJ) - (lJ / l) * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2); switch (pressureType_) { case pw: { - entries[rhs] += faceArea * lambdaNWIJ * kMean * (pcJK - pcI) / l; + entry[rhs] += faceArea * lambdaNWIJ * kMean * (pcJK - pcI) / l; break; } case pn: { - entries[rhs] -= faceArea * lambdaWIJ * kMean * (pcJK - pcI) / l; + entry[rhs] -= faceArea * lambdaWIJ * kMean * (pcJK - pcI) / l; break; } } //write hanging-node-specific stuff directly into matrix and rhs! - this->f_[globalIdxI] -= entries[rhs]; - this->f_[globalIdxJ] += entries[rhs]; + this->f_[globalIdxI] -= entry[rhs]; + this->f_[globalIdxJ] += entry[rhs]; // set diagonal entry - this->A_[globalIdxI][globalIdxI] += entries[matrix]; + this->A_[globalIdxI][globalIdxI] += entry[matrix]; //set off-diagonal - this->A_[globalIdxI][globalIdxJ] -= entries[matrix]; + this->A_[globalIdxI][globalIdxJ] -= entry[matrix]; - // set entries for cell J - this->A_[globalIdxJ][globalIdxI] -= entries[matrix]; - this->A_[globalIdxJ][globalIdxJ] += entries[matrix]; + // set entry for cell J + this->A_[globalIdxJ][globalIdxI] -= entry[matrix]; + this->A_[globalIdxJ][globalIdxJ] += entry[matrix]; - //set entries to zero -> matrix already written! - entries = 0.; + //set entry to zero -> matrix already written! + entry = 0.; // std::cout<<"finished hanging node!\n"; } else { - entries = 0; + entry = 0; } return; diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh index 170a0e0265d5603269dd429e8d6579a458612305..13fa0fb82fb04535d259069d9a3f7985371dcfb7 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh @@ -276,10 +276,10 @@ private: * Calculates the velocity at a cell-cell interface from a given pressure field. * * \param intersection Intersection of two grid cells -* \param CellData Object containing all model relevant cell data +* \param cellData Object containing all model relevant cell data */ template<class TypeTag> -void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellDataI) +void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData) { ElementPointer elementI = intersection.inside(); ElementPointer elementJ = intersection.outside(); @@ -293,17 +293,17 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, const GlobalPosition& globalPosJ = elementJ->geometry().center(); // get mobilities and fractional flow factors - Scalar lambdaWI = cellDataI.mobility(wPhaseIdx); - Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx); - Scalar fractionalWI = cellDataI.fracFlowFunc(wPhaseIdx); - Scalar fractionalNWI = cellDataI.fracFlowFunc(nPhaseIdx); + Scalar lambdaWI = cellData.mobility(wPhaseIdx); + Scalar lambdaNWI = cellData.mobility(nPhaseIdx); + Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx); + Scalar fractionalNWI = cellData.fracFlowFunc(nPhaseIdx); Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx); Scalar lambdaNWJ = cellDataJ.mobility(nPhaseIdx); Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx); Scalar fractionalNWJ = cellDataJ.fracFlowFunc(nPhaseIdx); // get capillary pressure - Scalar pcI = cellDataI.capillaryPressure(); + Scalar pcI = cellData.capillaryPressure(); Scalar pcJ = cellDataJ.capillaryPressure(); //get face index @@ -332,41 +332,41 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, Scalar potentialW = 0; Scalar potentialNW = 0; - potentialW = cellDataI.fluxData().potential(wPhaseIdx, isIndexI); - potentialNW = cellDataI.fluxData().potential(nPhaseIdx, isIndexI); + potentialW = cellData.fluxData().potential(wPhaseIdx, isIndexI); + potentialNW = cellData.fluxData().potential(nPhaseIdx, isIndexI); if (compressibility_) { - density_[wPhaseIdx] = (potentialW > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); - density_[nPhaseIdx] = (potentialNW > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); + density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); + density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); density_[wPhaseIdx] = - (potentialW == 0) ? 0.5 * (cellDataI.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) : + (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) : density_[wPhaseIdx]; density_[nPhaseIdx] = - (potentialNW == 0) ? 0.5 * (cellDataI.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) : + (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) : density_[nPhaseIdx]; } if (pressureType_ == pglobal) { - potentialW = (cellDataI.globalPressure() - cellDataJ.globalPressure() + potentialW = (cellData.globalPressure() - cellDataJ.globalPressure() - 0.5 * (fractionalNWI + fractionalNWJ) * (pcI - pcJ)); - potentialNW = (cellDataI.globalPressure() - cellDataJ.globalPressure() + potentialNW = (cellData.globalPressure() - cellDataJ.globalPressure() + 0.5 * (fractionalWI + fractionalWJ) * (pcI - pcJ)); } else { - potentialW = (cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)); - potentialNW = (cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)); + potentialW = (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)); + potentialNW = (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)); } potentialW += density_[wPhaseIdx] * (distVec * gravity_); //delta z/delta x in unitOuterNormal[z] potentialNW += density_[nPhaseIdx] * (distVec * gravity_); //store potentials for further calculations (velocity, saturation, ...) - cellDataI.fluxData().setPotential(wPhaseIdx, isIndexI, potentialW); - cellDataI.fluxData().setPotential(nPhaseIdx, isIndexI, potentialNW); + cellData.fluxData().setPotential(wPhaseIdx, isIndexI, potentialW); + cellData.fluxData().setPotential(nPhaseIdx, isIndexI, potentialNW); cellDataJ.fluxData().setPotential(wPhaseIdx, isIndexJ, -potentialW); cellDataJ.fluxData().setPotential(nPhaseIdx, isIndexJ, -potentialNW); @@ -379,14 +379,14 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, if (compressibility_) { - density_[wPhaseIdx] = (potentialW > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); - density_[nPhaseIdx] = (potentialNW > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); + density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); + density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); density_[wPhaseIdx] = - (potentialW == 0) ? 0.5 * (cellDataI.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) : + (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) : density_[wPhaseIdx]; density_[nPhaseIdx] = - (potentialNW == 0) ? 0.5 * (cellDataI.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) : + (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) : density_[nPhaseIdx]; } @@ -404,8 +404,8 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, { case pw: { - velocityW *= lambdaW * (cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist; - velocityNW *= lambdaNW * (cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + velocityW *= lambdaW * (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist; + velocityNW *= lambdaNW * (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + 0.5 * (lambdaNWI + lambdaNWJ) * (pcI - pcJ) / dist; velocityW += gravityTermW; velocityNW += gravityTermNW; @@ -413,16 +413,16 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, } case pn: { - velocityW *= lambdaW * (cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + velocityW *= lambdaW * (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist - 0.5 * (lambdaWI + lambdaWJ) * (pcI - pcJ) / dist; - velocityNW *= lambdaNW * (cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist; + velocityNW *= lambdaNW * (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist; velocityW += gravityTermW; velocityNW += gravityTermNW; break; } case pglobal: { - velocityW *= (lambdaW + lambdaNW) * (cellDataI.globalPressure() - cellDataJ.globalPressure()) / dist; + velocityW *= (lambdaW + lambdaNW) * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist; velocityW += gravityTermW; velocityW += gravityTermNW; velocityNW = 0; @@ -431,15 +431,15 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, } //store velocities - cellDataI.fluxData().setVelocity(wPhaseIdx, isIndexI, velocityW); - cellDataI.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNW); - cellDataI.fluxData().setVelocityMarker(isIndexI); + cellData.fluxData().setVelocity(wPhaseIdx, isIndexI, velocityW); + cellData.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNW); + cellData.fluxData().setVelocityMarker(isIndexI); cellDataJ.fluxData().setVelocity(wPhaseIdx, isIndexJ, velocityW); cellDataJ.fluxData().setVelocity(nPhaseIdx, isIndexJ, velocityNW); cellDataJ.fluxData().setVelocityMarker(isIndexJ); -// printvector(std::cout, cellDataI.fluxData().velocity(), "velocity", "row", 4, 1, 3); +// printvector(std::cout, cellData.fluxData().velocity(), "velocity", "row", 4, 1, 3); return; } @@ -448,10 +448,10 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, * Calculates the velocity at a boundary from a given pressure field. * * \param intersection Boundary intersection -* \param CellData Object containing all model relevant cell data +* \param cellData Object containing all model relevant cell data */ template<class TypeTag> -void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellDataI) +void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData) { ElementPointer element = intersection.inside(); @@ -477,13 +477,13 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte const GlobalPosition& globalPosJ = intersection.geometry().center(); // get mobilities and fractional flow factors - Scalar lambdaWI = cellDataI.mobility(wPhaseIdx); - Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx); - Scalar fractionalWI = cellDataI.fracFlowFunc(wPhaseIdx); - Scalar fractionalNWI = cellDataI.fracFlowFunc(nPhaseIdx); + Scalar lambdaWI = cellData.mobility(wPhaseIdx); + Scalar lambdaNWI = cellData.mobility(nPhaseIdx); + Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx); + Scalar fractionalNWI = cellData.fracFlowFunc(nPhaseIdx); // get capillary pressure - Scalar pcI = cellDataI.capillaryPressure(); + Scalar pcI = cellData.capillaryPressure(); // distance vector between barycenters GlobalPosition distVec = globalPosJ - globalPosI; @@ -523,8 +523,8 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte } else { - satW = cellDataI.saturation(wPhaseIdx); - satNW = cellDataI.saturation(nPhaseIdx); + satW = cellData.saturation(wPhaseIdx); + satNW = cellData.saturation(nPhaseIdx); } Scalar pressBound = boundValues[pressureIdx]; @@ -580,35 +580,35 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte Scalar potentialW = 0; Scalar potentialNW = 0; - potentialW = cellDataI.fluxData().potential(wPhaseIdx, isIndex); - potentialNW = cellDataI.fluxData().potential(nPhaseIdx, isIndex); + potentialW = cellData.fluxData().potential(wPhaseIdx, isIndex); + potentialNW = cellData.fluxData().potential(nPhaseIdx, isIndex); if (compressibility_) { - density_[wPhaseIdx] = (potentialW > 0.) ? cellDataI.density(wPhaseIdx) : densityWBound; - density_[wPhaseIdx] = (potentialW == 0) ? 0.5 * (cellDataI.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx]; - density_[nPhaseIdx] = (potentialNW > 0.) ? cellDataI.density(nPhaseIdx) : densityNWBound; - density_[nPhaseIdx] = (potentialNW == 0) ? 0.5 * (cellDataI.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx]; + density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : densityWBound; + density_[wPhaseIdx] = (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx]; + density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : densityNWBound; + density_[nPhaseIdx] = (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx]; } //calculate potential gradient if (pressureType_ == pglobal) { - potentialW = (cellDataI.globalPressure() - pressBound - fractionalNWI * (pcI - pcBound)); - potentialNW = (cellDataI.globalPressure() - pressBound + fractionalWI * (pcI - pcBound)); + potentialW = (cellData.globalPressure() - pressBound - fractionalNWI * (pcI - pcBound)); + potentialNW = (cellData.globalPressure() - pressBound + fractionalWI * (pcI - pcBound)); } else { - potentialW = (cellDataI.pressure(wPhaseIdx) - pressWBound); - potentialNW = (cellDataI.pressure(nPhaseIdx) - pressNWBound); + potentialW = (cellData.pressure(wPhaseIdx) - pressWBound); + potentialNW = (cellData.pressure(nPhaseIdx) - pressNWBound); } potentialW += density_[wPhaseIdx] * (distVec * gravity_); potentialNW += density_[nPhaseIdx] * (distVec * gravity_); //store potential gradients for further calculations - cellDataI.fluxData().setPotential(wPhaseIdx, isIndex, potentialW); - cellDataI.fluxData().setPotential(nPhaseIdx, isIndex, potentialNW); + cellData.fluxData().setPotential(wPhaseIdx, isIndex, potentialW); + cellData.fluxData().setPotential(nPhaseIdx, isIndex, potentialNW); //do the upwinding of the mobility depending on the phase potentials Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWBound; @@ -618,10 +618,10 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte if (compressibility_) { - density_[wPhaseIdx] = (potentialW > 0.) ? cellDataI.density(wPhaseIdx) : densityWBound; - density_[wPhaseIdx] = (potentialW == 0) ? 0.5 * (cellDataI.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx]; - density_[nPhaseIdx] = (potentialNW > 0.) ? cellDataI.density(nPhaseIdx) : densityNWBound; - density_[nPhaseIdx] = (potentialNW == 0) ? 0.5 * (cellDataI.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx]; + density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : densityWBound; + density_[wPhaseIdx] = (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx]; + density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : densityNWBound; + density_[nPhaseIdx] = (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx]; } //calculate the gravity term @@ -638,8 +638,8 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte { case pw: { - velocityW *= lambdaW * (cellDataI.pressure(wPhaseIdx) - pressBound) / dist; - velocityNW *= lambdaNW * (cellDataI.pressure(wPhaseIdx) - pressBound) / dist + velocityW *= lambdaW * (cellData.pressure(wPhaseIdx) - pressBound) / dist; + velocityNW *= lambdaNW * (cellData.pressure(wPhaseIdx) - pressBound) / dist + 0.5 * (lambdaNWI + lambdaNWBound) * (pcI - pcBound) / dist; velocityW += gravityTermW; velocityNW += gravityTermNW; @@ -647,16 +647,16 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte } case pn: { - velocityW *= lambdaW * (cellDataI.pressure(nPhaseIdx) - pressBound) / dist + velocityW *= lambdaW * (cellData.pressure(nPhaseIdx) - pressBound) / dist - 0.5 * (lambdaWI + lambdaWBound) * (pcI - pcBound) / dist; - velocityNW *= lambdaNW * (cellDataI.pressure(nPhaseIdx) - pressBound) / dist; + velocityNW *= lambdaNW * (cellData.pressure(nPhaseIdx) - pressBound) / dist; velocityW += gravityTermW; velocityNW += gravityTermNW; break; } case pglobal: { - velocityW *= (lambdaW + lambdaNW) * (cellDataI.globalPressure() - pressBound) / dist; + velocityW *= (lambdaW + lambdaNW) * (cellData.globalPressure() - pressBound) / dist; velocityW += gravityTermW; velocityW += gravityTermNW; velocityNW = 0; @@ -665,9 +665,9 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte } //store velocities - cellDataI.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW); - cellDataI.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNW); - cellDataI.fluxData().setVelocityMarker(isIndex); + cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW); + cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNW); + cellData.fluxData().setVelocityMarker(isIndex); } //end dirichlet boundary @@ -688,19 +688,19 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte } //store potential gradients for further calculations - cellDataI.fluxData().setPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]); - cellDataI.fluxData().setPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]); + cellData.fluxData().setPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]); + cellData.fluxData().setPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]); - cellDataI.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW); - cellDataI.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNW); - cellDataI.fluxData().setVelocityMarker(isIndex); + cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW); + cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNW); + cellData.fluxData().setVelocityMarker(isIndex); } //end neumann boundary else { DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!"); } -// printvector(std::cout, cellDataI.fluxData().velocity(), "velocity", "row", 4, 1, 3); +// printvector(std::cout, cellData.fluxData().velocity(), "velocity", "row", 4, 1, 3); return; } } diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh index a13bc28aeb99ffb17102a3247407c0e58b1843e7..3e6d460e0e085f6f89c89f5d2aefabf83ffcb1d3 100644 --- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh +++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh @@ -133,7 +133,7 @@ public: } //Calculates the velocity at a cell-cell interface. - void calculateVelocity(const Intersection& intersection, CellData& cellDataI); + void calculateVelocity(const Intersection& intersection, CellData& cellData); /*! \brief Indicates if velocity is reconstructed in the pressure step or in the transport step * @@ -163,14 +163,14 @@ private: * Implementation of calculateVelocity() function for cell-cell interfaces with hanging nodes. */ template<class TypeTag> -void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellDataI) +void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData) { ElementPointer elementI = intersection.inside(); ElementPointer elementJ = intersection.outside(); if (elementI->level() == elementJ->level()) { - ParentType::calculateVelocity(intersection, cellDataI); + ParentType::calculateVelocity(intersection, cellData); } else if (elementJ->level() == elementI->level() + 1) { @@ -181,17 +181,17 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters const GlobalPosition& globalPosJ = elementJ->geometry().center(); // get mobilities and fractional flow factors - Scalar lambdaWI = cellDataI.mobility(wPhaseIdx); - Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx); - Scalar fractionalWI = cellDataI.fracFlowFunc(wPhaseIdx); - Scalar fractionalNWI = cellDataI.fracFlowFunc(nPhaseIdx); + Scalar lambdaWI = cellData.mobility(wPhaseIdx); + Scalar lambdaNWI = cellData.mobility(nPhaseIdx); + Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx); + Scalar fractionalNWI = cellData.fracFlowFunc(nPhaseIdx); Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx); Scalar lambdaNWJ = cellDataJ.mobility(nPhaseIdx); Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx); Scalar fractionalNWJ = cellDataJ.fracFlowFunc(nPhaseIdx); // get capillary pressure - Scalar pcI = cellDataI.capillaryPressure(); + Scalar pcI = cellData.capillaryPressure(); Scalar pcJ = cellDataJ.capillaryPressure(); //get face index @@ -278,8 +278,8 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters // calculate potential gradients // reuse potentials from fvpressure2padaptive - Scalar potentialWIJ = cellDataI.fluxData().potential(wPhaseIdx, isIndexI); - Scalar potentialNWIJ = cellDataI.fluxData().potential(nPhaseIdx, isIndexI); + Scalar potentialWIJ = cellData.fluxData().potential(wPhaseIdx, isIndexI); + Scalar potentialNWIJ = cellData.fluxData().potential(nPhaseIdx, isIndexI); Scalar potentialWIK = potentialWIJ; Scalar potentialNWIK = potentialNWIJ; // preliminary potential. The "real" ones are found below @@ -292,10 +292,10 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters if (compressibility_) { - rhoMeanWIJ = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l; - rhoMeanNWIJ = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l; - rhoMeanWIK = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l; - rhoMeanNWIK = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l; + rhoMeanWIJ = (lJ * cellData.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l; + rhoMeanNWIJ = (lJ * cellData.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l; + rhoMeanWIK = (lJ * cellData.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l; + rhoMeanNWIK = (lJ * cellData.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l; } Scalar fMeanWIJ = (lJ * fractionalWI + lI * fractionalWJ) / l; @@ -311,10 +311,10 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters if (compressibility_) { // Upwinding for finding the upwinding direction - densityWIJ = (potentialWIJ > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); - densityNWIJ = (potentialNWIJ > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); - densityWIK = (potentialWIK > 0.) ? cellDataI.density(wPhaseIdx) : cellDataK.density(wPhaseIdx); - densityNWIK = (potentialNWIK > 0.) ? cellDataI.density(nPhaseIdx) : cellDataK.density(nPhaseIdx); + densityWIJ = (potentialWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); + densityNWIJ = (potentialNWIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); + densityWIK = (potentialWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx); + densityNWIK = (potentialNWIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx); densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ; densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ; @@ -338,28 +338,28 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters { Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2; - potentialWIJ = (cellDataI.globalPressure() - fMeanNWIJ * pcI - (pressJK - fMeanNWIJ * pcJK)) / l + potentialWIJ = (cellData.globalPressure() - fMeanNWIJ * pcI - (pressJK - fMeanNWIJ * pcJK)) / l + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng; - potentialNWIJ = (cellDataI.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l + potentialNWIJ = (cellData.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng; - potentialWIK = (cellDataI.globalPressure() - fMeanNWIK * pcI - (pressJK - fMeanNWIK * pcJK)) / l + potentialWIK = (cellData.globalPressure() - fMeanNWIK * pcI - (pressJK - fMeanNWIK * pcJK)) / l + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng; - potentialNWIK = (cellDataI.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l + potentialNWIK = (cellData.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng; break; } default: { - potentialWIJ = (cellDataI.pressure(wPhaseIdx) + potentialWIJ = (cellData.pressure(wPhaseIdx) - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng; - potentialNWIJ = (cellDataI.pressure(nPhaseIdx) + potentialNWIJ = (cellData.pressure(nPhaseIdx) - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng; - potentialWIK = (cellDataI.pressure(wPhaseIdx) + potentialWIK = (cellData.pressure(wPhaseIdx) - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng; - potentialNWIK = (cellDataI.pressure(nPhaseIdx) + potentialNWIK = (cellData.pressure(nPhaseIdx) - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng; break; @@ -369,8 +369,8 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters //store potentials for further calculations (velocity, saturation, ...) // these quantities only have correct sign (needed for upwinding) // potentials are defined slightly different for adaptive scheme - cellDataI.fluxData().addPotential(wPhaseIdx, isIndexI, potentialWIJ); - cellDataI.fluxData().addPotential(nPhaseIdx, isIndexI, potentialNWIJ); + cellData.fluxData().addPotential(wPhaseIdx, isIndexI, potentialWIJ); + cellData.fluxData().addPotential(nPhaseIdx, isIndexI, potentialNWIJ); cellDataJ.fluxData().setPotential(wPhaseIdx, isIndexJ, -potentialWIJ); cellDataJ.fluxData().setPotential(nPhaseIdx, isIndexJ, -potentialNWIJ); @@ -382,12 +382,12 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters if (compressibility_) { - densityWIJ = (potentialWIJ > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); - densityNWIJ = (potentialNWIJ > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); + densityWIJ = (potentialWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx); + densityNWIJ = (potentialNWIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx); densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ; densityNWIJ = (potentialNWIJ == 0) ? rhoMeanNWIJ : densityNWIJ; - densityWIK = (potentialWIK > 0.) ? cellDataI.density(wPhaseIdx) : cellDataK.density(wPhaseIdx); - densityNWIK = (potentialNWIK > 0.) ? cellDataI.density(nPhaseIdx) : cellDataK.density(nPhaseIdx); + densityWIK = (potentialWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx); + densityNWIK = (potentialNWIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx); densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK; densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK; } @@ -407,8 +407,8 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters { case pw: { - velocityW *= lambdaWIJ * kMean * (cellDataI.pressure(wPhaseIdx) - (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l; - velocityNW *= lambdaNWIJ * kMean * (cellDataI.pressure(nPhaseIdx) - (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l; + velocityW *= lambdaWIJ * kMean * (cellData.pressure(wPhaseIdx) - (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l; + velocityNW *= lambdaNWIJ * kMean * (cellData.pressure(nPhaseIdx) - (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l; velocityW += gravityTermW; velocityNW += gravityTermNW; @@ -416,8 +416,8 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters } case pn: { - velocityNW *= lambdaNWIJ * kMean * (cellDataI.pressure(nPhaseIdx) - (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l; - velocityW *= lambdaWIJ * kMean * (cellDataI.pressure(wPhaseIdx) - (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l; + velocityNW *= lambdaNWIJ * kMean * (cellData.pressure(nPhaseIdx) - (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l; + velocityW *= lambdaWIJ * kMean * (cellData.pressure(wPhaseIdx) - (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l; velocityW += gravityTermW; velocityNW += gravityTermNW; @@ -427,7 +427,7 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters { Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2; - velocityW *= lambdaWIJ * kMean * (cellDataI.globalPressure() - pressJK) / l; + velocityW *= lambdaWIJ * kMean * (cellData.globalPressure() - pressJK) / l; velocityW += gravityTermW; velocityW += gravityTermNW; velocityNW = 0; @@ -442,8 +442,8 @@ void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& inters //times 0.5 because cell face with hanging node is called twice! Do not set marker because it should be called twice! velocityW *= 0.5; velocityNW *= 0.5; - cellDataI.fluxData().addVelocity(wPhaseIdx, isIndexI, velocityW); - cellDataI.fluxData().addVelocity(nPhaseIdx, isIndexI, velocityNW); + cellData.fluxData().addVelocity(wPhaseIdx, isIndexI, velocityW); + cellData.fluxData().addVelocity(nPhaseIdx, isIndexI, velocityNW); } return;