diff --git a/dumux/material/constraintsolvers/compositionfromfugacities.hh b/dumux/material/constraintsolvers/compositionfromfugacities.hh index 8c256fafbbabe44ffc865058def2fd4b28a622d8..560352e46b76f2511153c3147b030dfac05bff07 100644 --- a/dumux/material/constraintsolvers/compositionfromfugacities.hh +++ b/dumux/material/constraintsolvers/compositionfromfugacities.hh @@ -60,15 +60,15 @@ public: int phaseIdx, const ComponentVector &fugVec) { - if (FluidSystem::isIdealMixture(phaseIdx)) - return; - - // Pure component fugacities - for (int i = 0; i < numComponents; ++ i) { - //std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n"; - fluidState.setMoleFraction(phaseIdx, - i, - 1.0/numComponents); + if (!FluidSystem::isIdealMixture(phaseIdx)) + { + // Pure component fugacities + for (unsigned int i = 0; i < numComponents; ++ i) + { + fluidState.setMoleFraction(phaseIdx, + i, + 1.0/numComponents); + } } } @@ -207,7 +207,6 @@ protected: Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx); fluidState.setDensity(phaseIdx, rho); - return; } template <class FluidState> @@ -291,46 +290,43 @@ protected: Dune::FieldVector<Scalar, numComponents> origComp; Scalar relError = 0; Scalar sumDelta = 0; - Scalar sumx = 0; - for (int i = 0; i < numComponents; ++i) { + for (unsigned int i = 0; i < numComponents; ++i) + { origComp[i] = fluidState.moleFraction(phaseIdx, i); relError = std::max(relError, std::abs(x[i])); - - sumx += std::abs(fluidState.moleFraction(phaseIdx, i)); sumDelta += std::abs(x[i]); } -#if 1 // chop update to at most 20% change in composition const Scalar maxDelta = 0.2; if (sumDelta > maxDelta) x /= (sumDelta/maxDelta); -#endif //Scalar curDefect = calculateDefect_(fluidState, phaseIdx, targetFug); //Scalar nextDefect; //Scalar sumMoleFrac = 0.0; //ComponentVector newB(1e100); //for (int numTries = 0; numTries < 1; ++numTries) { - // change composition - for (int i = 0; i < numComponents; ++i) { - Scalar newx = origComp[i] - x[i]; -#if 1 - // only allow negative mole fractions if the target fugacity is negative - if (targetFug[i] > 0) - newx = std::max(0.0, newx); - // only allow positive mole fractions if the target fugacity is positive - else if (targetFug[i] < 0) - newx = std::min(0.0, newx); - // if the target fugacity is zero, the mole fraction must also be zero - else - newx = 0; -#endif - fluidState.setMoleFraction(phaseIdx, i, newx); - //sumMoleFrac += std::abs(newx); - } + // change composition + for (unsigned int i = 0; i < numComponents; ++i) + { + Scalar newx = origComp[i] - x[i]; + + // only allow negative mole fractions if the target fugacity is negative + if (targetFug[i] > 0) + newx = std::max(0.0, newx); + // only allow positive mole fractions if the target fugacity is positive + else if (targetFug[i] < 0) + newx = std::min(0.0, newx); + // if the target fugacity is zero, the mole fraction must also be zero + else + newx = 0; + + fluidState.setMoleFraction(phaseIdx, i, newx); + //sumMoleFrac += std::abs(newx); + } - paramCache.updateComposition(fluidState, phaseIdx); + paramCache.updateComposition(fluidState, phaseIdx); /* // if the sum of the mole fractions gets 0, we take the diff --git a/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh b/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh index d417607e8dd7a81d74d4e4b827f0af19d6843a12..4ee240cbc6bdff391112fa9a9082292890f72821 100644 --- a/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh +++ b/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh @@ -1328,44 +1328,44 @@ void FVTransport2P2C<TypeTag>::innerUpdate(TransportSolutionType& updateVec) if (verbosity_ > 0) std::cout<<" Sub-time-step size: "<<subDt<< std::endl; - bool stopTimeStep = false; - int size = problem_.gridView().size(0); - for (int i = 0; i < size; i++) + bool stopTimeStep = false; + int size = problem_.gridView().size(0); + for (int i = 0; i < size; i++) + { + EntryType newVal(0); + int transportedQuantities = GET_PROP_VALUE(TypeTag, NumEq) - 1; // NumEq - 1 pressure Eq + for (int eqNumber = 0; eqNumber < transportedQuantities; eqNumber++) { - EntryType newVal(0); - int transportedQuantities = GET_PROP_VALUE(TypeTag, NumEq) - 1; // NumEq - 1 pressure Eq - for (int eqNumber = 0; eqNumber < transportedQuantities; eqNumber++) - { - newVal[eqNumber] = totalConcentration_[eqNumber][i]; - newVal[eqNumber] += updateVec[eqNumber][i] * subDt; - } - if (!asImp_().inPhysicalRange(newVal)) - { - stopTimeStep = true; + newVal[eqNumber] = totalConcentration_[eqNumber][i]; + newVal[eqNumber] += updateVec[eqNumber][i] * subDt; + } + if (!asImp_().inPhysicalRange(newVal)) + { + stopTimeStep = true; - break; - } + break; } + } #if HAVE_MPI - int rank = 0; - if (stopTimeStep) - rank = problem_.gridView().comm().rank(); + int rank = 0; + if (stopTimeStep) + rank = problem_.gridView().comm().rank(); - rank = problem_.gridView().comm().max(rank); - problem_.gridView().comm().broadcast(&stopTimeStep,1,rank); + rank = problem_.gridView().comm().max(rank); + problem_.gridView().comm().broadcast(&stopTimeStep,1,rank); #endif - if (stopTimeStep && accumulatedDtOld > dtThreshold_) - { - problem_.timeManager().setTimeStepSize(accumulatedDtOld); - break; - } - else - { - asImp_().updateTransportedQuantity(updateVec, subDt); - } + if (stopTimeStep && accumulatedDtOld > dtThreshold_) + { + problem_.timeManager().setTimeStepSize(accumulatedDtOld); + break; + } + else + { + asImp_().updateTransportedQuantity(updateVec, subDt); + } if (accumulatedDt_ >= realDt) diff --git a/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh b/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh index 8b29523389f09523bd5f709e6ca3bd73404bc56e..f9b1b4420fb732b96c5e3a4907a1ff466605e2c1 100644 --- a/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh @@ -137,39 +137,42 @@ public: typename FluidSystem::ParameterCache paramCache; paramCache.updateAll(fluidState_); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - {// relative permeabilities - Scalar kr; - if (phaseIdx == wPhaseIdx) - kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx)); - else // ATTENTION: krn requires the liquid saturation - // as parameter! - kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx)); - mobility_[phaseIdx] = kr / fluidState_.viscosity(phaseIdx); - Valgrind::CheckDefined(mobility_[phaseIdx]); - int compIIdx = phaseIdx; - for(int compIdx = 0; compIdx < numComponents; ++compIdx) - { - int compJIdx = compIdx; - // binary diffusion coefficents - diffCoeff_[phaseIdx][compIdx] = 0.0; - if(compIIdx!= compJIdx) - diffCoeff_[phaseIdx][compIdx] = FluidSystem::binaryDiffusionCoefficient(fluidState_, - paramCache, - phaseIdx, - compIIdx, - compJIdx); - Valgrind::CheckDefined(diffCoeff_[phaseIdx][compIdx]); + {// relative permeabilities + Scalar kr; + if (phaseIdx == wPhaseIdx) + kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx)); + else // ATTENTION: krn requires the liquid saturation + // as parameter! + kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx)); + + mobility_[phaseIdx] = kr / fluidState_.viscosity(phaseIdx); + Valgrind::CheckDefined(mobility_[phaseIdx]); + int compIIdx = phaseIdx; + for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx) + { + // binary diffusion coefficents + diffCoeff_[phaseIdx][compJIdx] = 0.0; + if(compIIdx!= compJIdx) + { + diffCoeff_[phaseIdx][compJIdx] = + FluidSystem::binaryDiffusionCoefficient(fluidState_, + paramCache, + phaseIdx, + compIIdx, + compJIdx); } + Valgrind::CheckDefined(diffCoeff_[phaseIdx][compJIdx]); } + } - // porosity - porosity_ = problem.spatialParams().porosity(element, - fvGeometry, - scvIdx); - Valgrind::CheckDefined(porosity_); - // energy related quantities not contained in the fluid state + // porosity + porosity_ = problem.spatialParams().porosity(element, + fvGeometry, + scvIdx); + Valgrind::CheckDefined(porosity_); + // energy related quantities not contained in the fluid state - asImp_().updateEnergy_(priVars, problem,element, fvGeometry, scvIdx, isOldSol); + asImp_().updateEnergy_(priVars, problem,element, fvGeometry, scvIdx, isOldSol); } /*! diff --git a/dumux/porousmediumflow/sequential/cellcentered/transport.hh b/dumux/porousmediumflow/sequential/cellcentered/transport.hh index 054f0f1c9f2f329230eaa703a9fe3f2075c8d893..f8d548f695f857a2c93f2cefd806a6890652f9b1 100644 --- a/dumux/porousmediumflow/sequential/cellcentered/transport.hh +++ b/dumux/porousmediumflow/sequential/cellcentered/transport.hh @@ -590,41 +590,41 @@ void FVTransport<TypeTag>::innerUpdate(TransportSolutionType& updateVec) if (verbosity_ > 0) std::cout<<" Sub-time-step size: "<<subDt<<"\n"; - TransportSolutionType transportedQuantity; - asImp_().getTransportedQuantity(transportedQuantity); + TransportSolutionType transportedQuantity; + asImp_().getTransportedQuantity(transportedQuantity); - bool stopTimeStep = false; - int size = transportedQuantity.size(); - for (int i = 0; i < size; i++) + bool stopTimeStep = false; + int size = transportedQuantity.size(); + for (int i = 0; i < size; i++) + { + Scalar newVal = transportedQuantity[i] += updateVec[i] * subDt; + if (!asImp_().inPhysicalRange(newVal)) { - Scalar newVal = transportedQuantity[i] += updateVec[i] * subDt; - if (!asImp_().inPhysicalRange(newVal)) - { - stopTimeStep = true; + stopTimeStep = true; - break; - } + break; } + } #if HAVE_MPI - int rank = 0; - if (stopTimeStep) - rank = problem_.gridView().comm().rank(); + int rank = 0; + if (stopTimeStep) + rank = problem_.gridView().comm().rank(); - rank = problem_.gridView().comm().max(rank); - problem_.gridView().comm().broadcast(&stopTimeStep,1,rank); + rank = problem_.gridView().comm().max(rank); + problem_.gridView().comm().broadcast(&stopTimeStep,1,rank); #endif - if (stopTimeStep && accumulatedDtOld > dtThreshold_) - { - problem_.timeManager().setTimeStepSize(accumulatedDtOld); - break; - } - else - { - asImp_().setTransportedQuantity(transportedQuantity); - } + if (stopTimeStep && accumulatedDtOld > dtThreshold_) + { + problem_.timeManager().setTimeStepSize(accumulatedDtOld); + break; + } + else + { + asImp_().setTransportedQuantity(transportedQuantity); + } if (accumulatedDt_ >= realDt) diff --git a/test/material/immiscibleflash/test_immiscibleflash.cc b/test/material/immiscibleflash/test_immiscibleflash.cc index 74d353b841df44729d1a11965c73e49182f3a59f..768da13fcc3be8bc8da6eeeced1907a10de8a50f 100644 --- a/test/material/immiscibleflash/test_immiscibleflash.cc +++ b/test/material/immiscibleflash/test_immiscibleflash.cc @@ -32,9 +32,9 @@ #include <dumux/material/constraintsolvers/computefromreferencephase.hh> #include <dumux/material/constraintsolvers/immiscibleflash.hh> -#include <dumux/material/fluidstates/immisciblefluidstate.hh> +#include <dumux/material/fluidstates/immiscible.hh> -#include <dumux/material/fluidsystems/h2on2fluidsystem.hh> +#include <dumux/material/fluidsystems/h2on2.hh> #include <dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/mp/2padapter.hh> diff --git a/test/material/ncpflash/test_ncpflash.cc b/test/material/ncpflash/test_ncpflash.cc index 1060d4c72846208ccc5e041c48399334d0b70eb1..a0d73818ed4e8f528205048e648249197b8f66fa 100644 --- a/test/material/ncpflash/test_ncpflash.cc +++ b/test/material/ncpflash/test_ncpflash.cc @@ -32,9 +32,9 @@ #include <dumux/material/constraintsolvers/computefromreferencephase.hh> #include <dumux/material/constraintsolvers/ncpflash.hh> -#include <dumux/material/fluidstates/compositionalfluidstate.hh> +#include <dumux/material/fluidstates/compositional.hh> -#include <dumux/material/fluidsystems/h2on2fluidsystem.hh> +#include <dumux/material/fluidsystems/h2on2.hh> #include <dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/mp/2padapter.hh> diff --git a/test/porousmediumflow/2p/implicit/pointsources/test_ccadaptive2ppointsource.cc b/test/porousmediumflow/2p/implicit/pointsources/test_ccadaptive2ppointsource.cc index 43fdde93d0e8d5cedfc2564f10a12d632b90750f..a39ae64aa66cd251eda03c87f46c7e5a82cb7652 100644 --- a/test/porousmediumflow/2p/implicit/pointsources/test_ccadaptive2ppointsource.cc +++ b/test/porousmediumflow/2p/implicit/pointsources/test_ccadaptive2ppointsource.cc @@ -22,7 +22,9 @@ * \brief Test for the two-phase CC model with point source */ #include <config.h> +#if HAVE_UG #include "lensproblem.hh" +#endif #include <dumux/common/start.hh> /*!