diff --git a/dumux/material/constraintsolvers/compositionfromfugacities.hh b/dumux/material/constraintsolvers/compositionfromfugacities.hh index 62d89f7b3ad55ce12c45c726fda6c9a25fa54a53..560352e46b76f2511153c3147b030dfac05bff07 100644 --- a/dumux/material/constraintsolvers/compositionfromfugacities.hh +++ b/dumux/material/constraintsolvers/compositionfromfugacities.hh @@ -307,25 +307,26 @@ protected: //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)