C++ FluidState::moleFraction方法代码示例

Scalar computeSumxg(FluidState& resultFluidState,
                    const FluidState& prestineFluidState,
                    const FluidState& gasFluidState,
                    Scalar additionalGas)
{
    static const int oilPhaseIdx = FluidSystem::oilPhaseIdx;
    static const int gasPhaseIdx = FluidSystem::gasPhaseIdx;
    static const int numComponents = FluidSystem::numComponents;

    typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
    typedef Opm::NcpFlash<Scalar, FluidSystem> Flash;

    resultFluidState.assign(prestineFluidState);

    // add a bit of additional gas components
    ComponentVector totalMolarities;
    for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++ compIdx)
        totalMolarities =
            prestineFluidState.molarity(oilPhaseIdx, compIdx)
            + additionalGas*gasFluidState.moleFraction(gasPhaseIdx, compIdx);

    // "flash" the modified fluid state
    typename FluidSystem::ParameterCache paramCache;
    Flash::solve(resultFluidState, totalMolarities);

    Scalar sumxg = 0;
    for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx)
        sumxg += resultFluidState.moleFraction(gasPhaseIdx, compIdx);

    return sumxg;
}

    static Scalar thermalConductivity(const FluidState &fluidState,
                                      const ParameterCache &paramCache,
                                      int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);

        Scalar temperature  = fluidState.temperature(phaseIdx) ;
        Scalar pressure = fluidState.pressure(phaseIdx);

        if (phaseIdx == lPhaseIdx){// liquid phase
            return H2O::liquidThermalConductivity(temperature, pressure);
        }
        else{// gas phase
            Scalar lambdaDryAir = Air::gasThermalConductivity(temperature, pressure);

            if (useComplexRelations){
                Scalar xAir = fluidState.moleFraction(phaseIdx, AirIdx);
                Scalar xH2O = fluidState.moleFraction(phaseIdx, H2OIdx);
                Scalar lambdaAir = xAir * lambdaDryAir;

                // Assuming Raoult's, Daltons law and ideal gas
                // in order to obtain the partial density of water in the air phase
                Scalar partialPressure  = pressure * xH2O;

                Scalar lambdaH2O =
                    xH2O
                    * H2O::gasThermalConductivity(temperature, partialPressure);
                return lambdaAir + lambdaH2O;
            }
            else
                return lambdaDryAir; // conductivity of Nitrogen [W / (m K ) ]
        }
    }

void makeOilSaturated(FluidState& fluidState, const FluidState& gasFluidState)
{
    static const int gasPhaseIdx = FluidSystem::gasPhaseIdx;

    FluidState prestineFluidState;
    prestineFluidState.assign(fluidState);

    Scalar sumxg = 0;
    for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx)
        sumxg += fluidState.moleFraction(gasPhaseIdx, compIdx);

    // Newton method
    Scalar tol = 1e-8;
    Scalar additionalGas = 0; // [mol]
    for (int i = 0; std::abs(sumxg - 1) > tol; ++i) {
        if (i > 50)
            throw std::runtime_error("Newton method did not converge after 50 iterations");

        Scalar eps = std::max(1e-8, additionalGas*1e-8);

        Scalar f = 1 - computeSumxg<Scalar, FluidSystem>(prestineFluidState,
                                                         fluidState,
                                                         gasFluidState,
                                                         additionalGas);
        Scalar fStar = 1 - computeSumxg<Scalar, FluidSystem>(prestineFluidState,
                                                             fluidState,
                                                             gasFluidState,
                                                             additionalGas + eps);
        Scalar fPrime = (fStar - f)/eps;

        additionalGas -= f/fPrime;
    };
}

    static Scalar update_(FluidState &fluidState,
                          ParameterCache &paramCache,
                          Dune::FieldVector<Evaluation, numComponents> &x,
                          Dune::FieldVector<Evaluation, numComponents> &b,
                          int phaseIdx,
                          const Dune::FieldVector<Evaluation, numComponents> &targetFug)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        // store original composition and calculate relative error
        Dune::FieldVector<Evaluation, numComponents> origComp;
        Scalar relError = 0;
        Evaluation sumDelta = Toolbox::createConstant(0.0);
        Evaluation sumx = Toolbox::createConstant(0.0);
        for (int i = 0; i < numComponents; ++i) {
            origComp[i] = fluidState.moleFraction(phaseIdx, i);
            relError = std::max(relError, std::abs(Toolbox::value(x[i])));

            sumx += Toolbox::abs(fluidState.moleFraction(phaseIdx, i));
            sumDelta += Toolbox::abs(x[i]);
        }

        // chop update to at most 20% change in composition
        const Scalar maxDelta = 0.2;
        if (sumDelta > maxDelta)
            x /= (sumDelta/maxDelta);

        // change composition
        for (int i = 0; i < numComponents; ++i) {
            Evaluation newx = origComp[i] - x[i];
            // only allow negative mole fractions if the target fugacity is negative
            if (targetFug[i] > 0)
                newx = Toolbox::max(0.0, newx);
            // only allow positive mole fractions if the target fugacity is positive
            else if (targetFug[i] < 0)
                newx = Toolbox::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);
        }

        paramCache.updateComposition(fluidState, phaseIdx);

        return relError;
    }

    static Scalar density(const FluidState &fluidState,
                          const ParameterCache &paramCache,
                          int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        if (phaseIdx == lPhaseIdx) {
            // use normalized composition for to calculate the density
            // (the relations don't seem to take non-normalized
            // compositions too well...)
            Scalar xlBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(lPhaseIdx, BrineIdx)));
            Scalar xlCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(lPhaseIdx, CO2Idx)));
            Scalar sumx = xlBrine + xlCO2;
            xlBrine /= sumx;
            xlCO2 /= sumx;

            Scalar result = liquidDensity_(temperature,
                                           pressure,
                                           xlBrine,
                                           xlCO2);

            Valgrind::CheckDefined(result);
            return result;
        }

        assert(phaseIdx == gPhaseIdx);

        // use normalized composition for to calculate the density
        // (the relations don't seem to take non-normalized
        // compositions too well...)
        Scalar xgBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(gPhaseIdx, BrineIdx)));
        Scalar xgCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(gPhaseIdx, CO2Idx)));
        Scalar sumx = xgBrine + xgCO2;
        xgBrine /= sumx;
        xgCO2 /= sumx;

        Scalar result = gasDensity_(temperature,
                                    pressure,
                                    xgBrine,
                                    xgCO2);
        Valgrind::CheckDefined(result);
        return result;
    }

    void updateMix(const FluidState &fs)
    {
        Scalar sumx = 0.0;
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
            sumx += fs.moleFraction(phaseIdx, compIdx);
        sumx = std::max(Scalar(1e-10), sumx);

        // Calculate the Peng-Robinson parameters of the mixture
        //
        // See: R. Reid, et al.: The Properties of Gases and Liquids,
        // 4th edition, McGraw-Hill, 1987, p. 82
        Scalar newA = 0;
        Scalar newB = 0;
        for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
            const Scalar moleFracI = fs.moleFraction(phaseIdx, compIIdx);
            Scalar xi = std::max(Scalar(0), std::min(Scalar(1.0), moleFracI));
            Valgrind::CheckDefined(xi);

            for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
                const Scalar moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
                Scalar xj = std::max(Scalar(0), std::min(Scalar(1), moleFracJ));
                Valgrind::CheckDefined(xj);

                // mixing rule from Reid, page 82
                newA +=  xi * xj * aCache_[compIIdx][compJIdx];

                assert(std::isfinite(newA));
            }

            // mixing rule from Reid, page 82
            newB += std::max(Scalar(0), xi) * this->pureParams_[compIIdx].b();
            assert(std::isfinite(newB));
        }

        // assert(newB > 0);
        this->setA(newA);
        this->setB(newB);

        Valgrind::CheckDefined(this->a());
        Valgrind::CheckDefined(this->b());

    }

 void assign(const FluidState& fs)
 {
     for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
         averageMolarMass_[phaseIdx] = 0;
         sumMoleFractions_[phaseIdx] = 0;
         for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
             moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
             averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
             sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
         }
     }
 }

    void assign(const FluidState& fs)
    {
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            averageMolarMass_[phaseIdx] = 0;
            sumMoleFractions_[phaseIdx] = 0;
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                moleFraction_[phaseIdx][compIdx] =
                    Opm::decay<Scalar>(fs.moleFraction(phaseIdx, compIdx));

                averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
                sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
            }
        }
    }

    static LhsEval thermalConductivity(const FluidState &fluidState,
                                       const ParameterCache &/*paramCache*/,
                                       unsigned phaseIdx)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        assert(0 <= phaseIdx  && phaseIdx < numPhases);

        const LhsEval& temperature =
            FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& pressure =
            FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));

        if (phaseIdx == liquidPhaseIdx)
            return H2O::liquidThermalConductivity(temperature, pressure);
        else { // gas phase
            const LhsEval& lambdaDryAir = Air::gasThermalConductivity(temperature, pressure);

            if (useComplexRelations){
                const LhsEval& xAir =
                    FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, AirIdx));
                const LhsEval& xH2O =
                    FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
                LhsEval lambdaAir = xAir*lambdaDryAir;

                // Assuming Raoult's, Daltons law and ideal gas
                // in order to obtain the partial density of water in the air phase
                LhsEval partialPressure  = pressure*xH2O;

                LhsEval lambdaH2O =
                    xH2O*H2O::gasThermalConductivity(temperature, partialPressure);
                return lambdaAir + lambdaH2O;
            }
            else
                return lambdaDryAir; // conductivity of Nitrogen [W / (m K ) ]
        }
    }

void checkSame(const FluidState &fsRef, const FluidState &fsFlash)
{
    enum { numPhases = FluidState::numPhases };
    enum { numComponents = FluidState::numComponents };

    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
        Scalar error;

        // check the pressures
        error = 1 - fsRef.pressure(phaseIdx)/fsFlash.pressure(phaseIdx);
        if (std::abs(error) > 1e-6) {
            std::cout << "pressure error phase " << phaseIdx << ": "
                      << fsFlash.pressure(phaseIdx)  << " flash vs "
                      << fsRef.pressure(phaseIdx) << " reference"
                      << " error=" << error << "\n";
        }

        // check the saturations
        error = fsRef.saturation(phaseIdx) - fsFlash.saturation(phaseIdx);
        if (std::abs(error) > 1e-6)
            std::cout << "saturation error phase " << phaseIdx << ": "
                      << fsFlash.saturation(phaseIdx) << " flash vs "
                      << fsRef.saturation(phaseIdx) << " reference"
                      << " error=" << error << "\n";

        // check the compositions
        for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
            error = fsRef.moleFraction(phaseIdx, compIdx) - fsFlash.moleFraction(phaseIdx, compIdx);
            if (std::abs(error) > 1e-6)
                std::cout << "composition error phase " << phaseIdx << ", component " << compIdx << ": "
                          << fsFlash.moleFraction(phaseIdx, compIdx) << " flash vs "
                          << fsRef.moleFraction(phaseIdx, compIdx) << " reference"
                          << " error=" << error << "\n";
        }
    }
}

    static LhsEval viscosity(const FluidState &fluidState,
                             const ParameterCache &/*paramCache*/,
                             unsigned phaseIdx)
    {
        typedef Opm::MathToolbox<LhsEval> LhsToolbox;
        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;

        assert(0 <= phaseIdx  && phaseIdx < numPhases);

        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));

        if (phaseIdx == liquidPhaseIdx)
        {
            // assume pure water for the liquid phase
            // TODO: viscosity of mixture
            // couldn't find a way to solve the mixture problem
            return H2O::liquidViscosity(T, p);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            if(!useComplexRelations){
                return Air::gasViscosity(T, p);
            }
            else //using a complicated version of this fluid system
            {
                /* Wilke method. See:
                 *
                 * See: R. Reid, et al.: The Properties of Gases and Liquids,
                 * 4th edition, McGraw-Hill, 1987, 407-410 or
                 * 5th edition, McGraw-Hill, 2000, p. 9.21/22
                 *
                 */

                LhsEval muResult = 0;
                const LhsEval mu[numComponents] = {
                    H2O::gasViscosity(T, H2O::vaporPressure(T)),
                    Air::gasViscosity(T, p)
                };

                // molar masses
                const Scalar M[numComponents] =  {
                    H2O::molarMass(),
                    Air::molarMass()
                };

                for (unsigned i = 0; i < numComponents; ++i) {
                    LhsEval divisor = 0;
                    for (unsigned j = 0; j < numComponents; ++j) {
                        LhsEval phiIJ =
                            1 +
                            LhsToolbox::sqrt(mu[i]/mu[j]) * // 1 + (mu[i]/mu[j]^1/2
                            std::pow(M[j]/M[i], 1./4.0);   // (M[i]/M[j])^1/4

                        phiIJ *= phiIJ;
                        phiIJ /= std::sqrt(8*(1 + M[i]/M[j]));
                        divisor += FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, j))*phiIJ;
                    }
                    const auto& xAlphaI = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, i));
                    muResult += xAlphaI*mu[i]/divisor;
                }
                return muResult;
            }
        }
        OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
    }

    static Scalar computeFugacityCoefficient(const FluidState &fs,
                                             const Params &params,
                                             int phaseIdx,
                                             int compIdx)
    {
        // note that we normalize the component mole fractions, so
        // that their sum is 100%. This increases numerical stability
        // considerably if the fluid state is not physical.
        Scalar Vm = params.molarVolume(phaseIdx);

        // Calculate b_i / b
        Scalar bi_b = params.bPure(phaseIdx, compIdx) / params.b(phaseIdx);

        // Calculate the compressibility factor
        Scalar RT = R*fs.temperature(phaseIdx);
        Scalar p = fs.pressure(phaseIdx); // molar volume in [bar]
        Scalar Z = p*Vm/RT; // compressibility factor

        // Calculate A^* and B^* (see: Reid, p. 42)
        Scalar Astar = params.a(phaseIdx)*p/(RT*RT);
        Scalar Bstar = params.b(phaseIdx)*p/(RT);

        // calculate delta_i (see: Reid, p. 145)
        Scalar sumMoleFractions = 0.0;
        for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
            sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
        Scalar deltai = 2*std::sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
        Scalar tmp = 0;
        for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
            tmp +=
                fs.moleFraction(phaseIdx, compJIdx)
                / sumMoleFractions
                * std::sqrt(params.aPure(phaseIdx, compJIdx))
                * (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
        };
        deltai *= tmp;

        Scalar base =
            (2*Z + Bstar*(u + std::sqrt(u*u - 4*w))) /
            (2*Z + Bstar*(u - std::sqrt(u*u - 4*w)));
        Scalar expo =  Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai);

        Scalar fugCoeff =
            std::exp(bi_b*(Z - 1))/std::max(1e-9, Z - Bstar) *
            std::pow(base, expo);

        ////////
        // limit the fugacity coefficient to a reasonable range:
        //
        // on one side, we want the mole fraction to be at
        // least 10^-3 if the fugacity is at the current pressure
        //
        fugCoeff = std::min(1e10, fugCoeff);
        //
        // on the other hand, if the mole fraction of the component is 100%, we want the
        // fugacity to be at least 10^-3 Pa
        //
        fugCoeff = std::max(1e-10, fugCoeff);
        ///////////

        return fugCoeff;
    }

Scalar bringOilToSurface(FluidState& surfaceFluidState, Scalar alpha, const FluidState& reservoirFluidState, bool guessInitial)
{
    enum {
        numPhases = FluidSystem::numPhases,
        waterPhaseIdx = FluidSystem::waterPhaseIdx,
        gasPhaseIdx = FluidSystem::gasPhaseIdx,
        oilPhaseIdx = FluidSystem::oilPhaseIdx,

        numComponents = FluidSystem::numComponents
    };

    typedef Opm::NcpFlash<Scalar, FluidSystem> Flash;
    typedef Opm::ThreePhaseMaterialTraits<Scalar, waterPhaseIdx, oilPhaseIdx, gasPhaseIdx> MaterialTraits;
    typedef Opm::LinearMaterial<MaterialTraits> MaterialLaw;
    typedef typename MaterialLaw::Params MaterialLawParams;
    typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;

    const Scalar refPressure = 1.0135e5; // [Pa]

    // set the parameters for the capillary pressure law
    MaterialLawParams matParams;
    for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
        matParams.setPcMinSat(phaseIdx, 0.0);
        matParams.setPcMaxSat(phaseIdx, 0.0);
    }
    matParams.finalize();

    // retieve the global volumetric component molarities
    surfaceFluidState.setTemperature(273.15 + 20);

    ComponentVector molarities;
    for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
        molarities[compIdx] = reservoirFluidState.molarity(oilPhaseIdx, compIdx);

    if (guessInitial) {
        // we start at a fluid state with reservoir oil.
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
            for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
                surfaceFluidState.setMoleFraction(phaseIdx,
                                                  compIdx,
                                                  reservoirFluidState.moleFraction(phaseIdx, compIdx));
            }
            surfaceFluidState.setDensity(phaseIdx, reservoirFluidState.density(phaseIdx));
            surfaceFluidState.setPressure(phaseIdx, reservoirFluidState.pressure(phaseIdx));
            surfaceFluidState.setSaturation(phaseIdx, 0.0);
        }
        surfaceFluidState.setSaturation(oilPhaseIdx, 1.0);
        surfaceFluidState.setSaturation(gasPhaseIdx, 1.0 - surfaceFluidState.saturation(oilPhaseIdx));
    }

    typename FluidSystem::template ParameterCache<Scalar> paramCache;
    paramCache.updateAll(surfaceFluidState);

    // increase volume until we are at surface pressure. use the
    // newton method for this
    ComponentVector tmpMolarities;
    for (int i = 0;; ++i) {
        if (i >= 20)
            throw Opm::NumericalIssue("Newton method did not converge after 20 iterations");

        // calculate the deviation from the standard pressure
        tmpMolarities = molarities;
        tmpMolarities /= alpha;
        Flash::template solve<MaterialLaw>(surfaceFluidState, matParams, paramCache, tmpMolarities);
        Scalar f = surfaceFluidState.pressure(gasPhaseIdx) - refPressure;

        // calculate the derivative of the deviation from the standard
        // pressure
        Scalar eps = alpha*1e-10;
        tmpMolarities = molarities;
        tmpMolarities /= alpha + eps;
        Flash::template solve<MaterialLaw>(surfaceFluidState, matParams, paramCache, tmpMolarities);
        Scalar fStar = surfaceFluidState.pressure(gasPhaseIdx) - refPressure;
        Scalar fPrime = (fStar - f)/eps;

        // newton update
        Scalar delta = f/fPrime;
        alpha -= delta;
        if (std::abs(delta) < std::abs(alpha)*1e-9) {
            break;
        }
    }

    // calculate the final result
    tmpMolarities = molarities;
    tmpMolarities /= alpha;
    Flash::template solve<MaterialLaw>(surfaceFluidState, matParams, paramCache, tmpMolarities);
    return alpha;
}

    static Scalar viscosity(const FluidState &fluidState,
                            const ParameterCache &paramCache,
                            int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);

        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p = fluidState.pressure(phaseIdx);

        if (phaseIdx == lPhaseIdx)
        {
            // assume pure water for the liquid phase
            // TODO: viscosity of mixture
            // couldn't find a way to solve the mixture problem
            return H2O::liquidViscosity(T, p);
        }
        else if (phaseIdx == gPhaseIdx)
        {
            if(!useComplexRelations){
                return Air::gasViscosity(T, p);
            }
            else //using a complicated version of this fluid system
            {
                /* Wilke method. See:
                 *
                 * See: R. Reid, et al.: The Properties of Gases and Liquids,
                 * 4th edition, McGraw-Hill, 1987, 407-410 or
                 * 5th edition, McGraw-Hill, 2000, p. 9.21/22
                 *
                 */

                Scalar muResult = 0;
                const Scalar mu[numComponents] = {
                    H2O::gasViscosity(T,
                                      H2O::vaporPressure(T)),
                    Air::gasViscosity(T, p)
                };

                // molar masses
                const Scalar M[numComponents] =  {
                    H2O::molarMass(),
                    Air::molarMass()
                };

                for (int i = 0; i < numComponents; ++i)
                {
                    Scalar divisor = 0;
                    for (int j = 0; j < numComponents; ++j)
                    {
                        Scalar phiIJ = 1 + std::sqrt(mu[i]/mu[j]) * // 1 + (mu[i]/mu[j]^1/2
                            std::pow(M[j]/M[i], 1./4.0);   // (M[i]/M[j])^1/4

                        phiIJ *= phiIJ;
                        phiIJ /= std::sqrt(8*(1 + M[i]/M[j]));
                        divisor += fluidState.moleFraction(phaseIdx, j)*phiIJ;
                    }
                    muResult += fluidState.moleFraction(phaseIdx, i)*mu[i] / divisor;
                }
                return muResult;
            }
        }
        OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
    }

    static Scalar density(const FluidState &fluidState,
                          const ParameterCache &paramCache,
                          int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);

        Scalar T = fluidState.temperature(phaseIdx);
        Scalar p;
        if (isCompressible(phaseIdx))
            p = fluidState.pressure(phaseIdx);
        else {
            // random value which will hopefully cause things to blow
            // up if it is used in a calculation!
            p = - 1e100;
            Valgrind::SetUndefined(p);
        }


        Scalar sumMoleFrac = 0;
        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
            sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);

        if (phaseIdx == lPhaseIdx)
        {
            if (!useComplexRelations)
                // assume pure water
                return H2O::liquidDensity(T, p);
            else
            {
                // See: Ochs 2008 (2.6)
                Scalar rholH2O = H2O::liquidDensity(T, p);
                Scalar clH2O = rholH2O/H2O::molarMass();

                return
                    clH2O
                    * (H2O::molarMass()*fluidState.moleFraction(lPhaseIdx, H2OIdx)
                           +
                           Air::molarMass()*fluidState.moleFraction(lPhaseIdx, AirIdx))
                   / sumMoleFrac;
            }
        }
        else if (phaseIdx == gPhaseIdx)
        {
            if (!useComplexRelations)
                // for the gas phase assume an ideal gas
                return
                    IdealGas::molarDensity(T, p)
                    * fluidState.averageMolarMass(gPhaseIdx)
                    / std::max(1e-5, sumMoleFrac);

            Scalar partialPressureH2O =
                fluidState.moleFraction(gPhaseIdx, H2OIdx)  *
                fluidState.pressure(gPhaseIdx);

            Scalar partialPressureAir =
                fluidState.moleFraction(gPhaseIdx, AirIdx)  *
                fluidState.pressure(gPhaseIdx);

            return
                H2O::gasDensity(T, partialPressureH2O) +
                Air::gasDensity(T, partialPressureAir);
        }
        OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
    }