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

 void assign(const FluidState& fs)
 {
     typedef typename FluidState::Scalar FsScalar;
     typedef Opm::MathToolbox<FsScalar> FsToolbox;
     for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
         saturation_[phaseIdx] = FsToolbox::template decay<Scalar>(fs.saturation(phaseIdx));
     }
 }

    static Evaluation krg(const Params &params,
                          const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& Sw = 1 - FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
        return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
    }

    static Evaluation pcnw(const Params &params,
                           const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(waterPhaseIdx));
        return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
    }

    static Evaluation krw(const Params &params,
                          const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& Sw = FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
        return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
    }

    static Evaluation krn(const Params &params, const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& Sw =
            1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
        return twoPhaseSatKrn(params, Sw);
    }

    static Evaluation pcgn(const Params &params,
                           const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const auto& Sw = 1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(gasPhaseIdx));
        return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
    }

    static Evaluation pcgn(const Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
        typedef MathToolbox<Evaluation> Toolbox;

        Scalar PC_VG_REG = 0.01;

        // sum of liquid saturations
        const auto& St =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx))
            + FsToolbox::template toLhs<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));

        Evaluation Se = (St - params.Swrx())/(1. - params.Swrx());

        // regularization
        if (Se < 0.0)
            Se=0.0;
        if (Se > 1.0)
            Se=1.0;

        if (Se>PC_VG_REG && Se<1-PC_VG_REG)
        {
            const Evaluation& x = Toolbox::pow(Se,-1/params.vgM()) - 1;
            return Toolbox::pow(x, 1.0 - params.vgM())/params.vgAlpha();
        }

        // value and derivative at regularization point
        Scalar Se_regu;
        if (Se<=PC_VG_REG)
            Se_regu = PC_VG_REG;
        else
            Se_regu = 1-PC_VG_REG;
        const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
        const Evaluation& pc = Toolbox::pow(x, 1.0/params.vgN())/params.vgAlpha();
        const Evaluation& pc_prime =
            Toolbox::pow(x, 1/params.vgN()-1)
            * std::pow(Se_regu,-1/params.vgM()-1)
            / (-params.vgM())
            / params.vgAlpha()
            / (1 - params.Sgr() - params.Swrx())
            / params.vgN();

        // evaluate tangential
        return ((Se-Se_regu)*pc_prime + pc)/params.betaGN();
    }

    static Scalar pcnw(const Params &params, const FluidState &fs)
    {
        Scalar S = fs.saturation(Traits::wPhaseIdx);
        Valgrind::CheckDefined(S);

        Scalar wPhasePressure =
            S*params.pcMaxSat(Traits::wPhaseIdx) +
            (1.0 - S)*params.pcMinSat(Traits::wPhaseIdx);

        S = fs.saturation(Traits::nPhaseIdx);
        Valgrind::CheckDefined(S);

        Scalar nPhasePressure =
            S*params.pcMaxSat(Traits::nPhaseIdx) +
            (1.0 - S)*params.pcMinSat(Traits::nPhaseIdx);

        return nPhasePressure - wPhasePressure;
    }

    static Evaluation pcnw(const Params& params, const FluidState& fs)
    {
        const Evaluation& Sw =
            Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));

        assert(0.0 <= Sw && Sw <= 1.0);

        return twoPhaseSatPcnw(params, Sw);
    }

    static Evaluation pcnw(const Params &params, const FluidState &fluidState)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& Sw =
            FsToolbox::template toLhs<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));

        return twoPhaseSatPcnw(params, Sw);
    }

    static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
    pcgn(const Params &params, const FluidState &fs)
    {
        Scalar S = fs.saturation(Traits::nPhaseIdx);
        Valgrind::CheckDefined(S);

        Scalar nPhasePressure =
            S*params.pcMaxSat(Traits::nPhaseIdx) +
            (1.0 - S)*params.pcMinSat(Traits::nPhaseIdx);

        S = fs.saturation(Traits::gPhaseIdx);
        Valgrind::CheckDefined(S);

        Scalar gPhasePressure =
            S*params.pcMaxSat(Traits::gPhaseIdx) +
            (1.0 - S)*params.pcMinSat(Traits::gPhaseIdx);

        return gPhasePressure - nPhasePressure;
    }

    static void relativePermeabilities(ContainerT &values,
                                       const Params &params,
                                       const FluidState &state)
    {
        for (int phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
            Scalar S = state.saturation(phaseIdx);
            Valgrind::CheckDefined(S);

            values[phaseIdx] = std::max(std::min(S,1.0),0.0);
        }
    }

    static Evaluation krn(const Params &params,
                          const FluidState &fluidState)
    {
        typedef MathToolbox<Evaluation> Toolbox;
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        // the Eclipse docu is inconsistent here: In some places the connate water
        // saturation is represented by "Swl", in others "Swco" is used.
        Scalar Swco = params.Swl();

        Scalar Sowcr = params.Sowcr();
        Scalar Sogcr = params.Sogcr();
        Scalar Som = std::min(Sowcr, Sogcr); // minimum residual oil saturation

        Scalar eta = params.eta(); // exponent of the beta term

        const Evaluation& Sw = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(waterPhaseIdx));
        const Evaluation& So = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(oilPhaseIdx));
        const Evaluation& Sg = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));

        Evaluation SSw;
        if (Sw > Swco)
            SSw = (Sw - Swco)/(1 - Swco - Som);
        else
            SSw = 0.0;

        Evaluation SSo;
        if (So > Som)
            SSo = (So - Som)/(1 - Swco - Som);
        else
            SSo = 0.0;

        Evaluation SSg = Sg/(1 - Swco - Som);

        Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco);
        Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
        Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg);

        Evaluation beta = Toolbox::pow(SSo/((1 - SSw)*(1 - SSg)), eta);
        return beta*krow*krog/krocw;
    }

    static void capillaryPressures(ContainerT &values,
                                   const Params &params,
                                   const FluidState &fluidState)
    {
        Scalar p_cgw = ThreePLaw::pCGW(params, fluidState.saturation(wPhaseIdx));
        Scalar p_cnw = ThreePLaw::pCNW(params, fluidState.saturation(wPhaseIdx));
        Scalar p_cgn = ThreePLaw::pCGN(params,
                                       fluidState.saturation(wPhaseIdx)
                                       + fluidState.saturation(nPhaseIdx));
        Scalar p_cAlpha = ThreePLaw::pCAlpha(params,
                                             fluidState.saturation(nPhaseIdx));
        Scalar p_cnw1 = 0.0;
        Valgrind::CheckDefined(p_cgw);
        Valgrind::CheckDefined(p_cnw);
        Valgrind::CheckDefined(p_cgn);
        Valgrind::CheckDefined(p_cAlpha);

        values[gPhaseIdx] = 0;
        values[nPhaseIdx] = - (p_cAlpha*p_cgn + (1 - p_cAlpha)*(p_cgw - p_cnw1));
        values[wPhaseIdx] = values[nPhaseIdx] - (p_cAlpha*p_cnw + (1 - p_cAlpha)*p_cnw1);
    }

    static void capillaryPressures(ContainerT &values,
                                   const Params &params,
                                   const FluidState &fluidState)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        switch (params.approach()) {
        case EclTwoPhaseGasOil: {
            const Evaluation& So =
                FsToolbox::template decay<Evaluation>(fluidState.saturation(oilPhaseIdx));

            values[oilPhaseIdx] = 0.0;
            values[gasPhaseIdx] = GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), So);
            break;
        }

        case EclTwoPhaseOilWater: {
            const Evaluation& Sw =
                FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));

            values[waterPhaseIdx] = 0.0;
            values[oilPhaseIdx] = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
            break;
        }

        case EclTwoPhaseGasWater: {
            const Evaluation& Sw =
                FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));

            values[waterPhaseIdx] = 0.0;
            values[gasPhaseIdx] =
                OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw)
                + GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), 0.0);
            break;
        }

        }
    }