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

void createSurfaceGasFluidSystem(FluidState& gasFluidState)
{
    static const int gasPhaseIdx = FluidSystem::gasPhaseIdx;

    // temperature
    gasFluidState.setTemperature(273.15 + 20);

    // gas pressure
    gasFluidState.setPressure(gasPhaseIdx, 1e5);

    // gas saturation
    gasFluidState.setSaturation(gasPhaseIdx, 1.0);

    //  gas composition: mostly methane, a bit of propane
    gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::H2OIdx, 0.0);
    gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C1Idx, 0.94);
    gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C3Idx, 0.06);
    gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C6Idx, 0.00);
    gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C10Idx, 0.00);
    gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C15Idx, 0.00);
    gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C20Idx, 0.00);

    // gas density
    typename FluidSystem::template ParameterCache<typename FluidState::Scalar> paramCache;
    paramCache.updatePhase(gasFluidState, gasPhaseIdx);
    gasFluidState.setDensity(gasPhaseIdx,
                             FluidSystem::density(gasFluidState, paramCache, gasPhaseIdx));
}

    static void guessInitial(FluidState &fluidState,
                             unsigned phaseIdx,
                             const ComponentVector &/*fugVec*/)
    {
        if (FluidSystem::isIdealMixture(phaseIdx))
            return;

        // Pure component fugacities
        for (unsigned i = 0; i < numComponents; ++ i) {
            //std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
            fluidState.setMoleFraction(phaseIdx,
                                       i,
                                       1.0/numComponents);
        }
    }

void guessInitial(FluidState& fluidState, unsigned phaseIdx)
{
    if (phaseIdx == FluidSystem::gasPhaseIdx) {
        fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.74785);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.0121364);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.00606028);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.00268136);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.000204256);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 8.78291e-06);
    }
    else if (phaseIdx == FluidSystem::oilPhaseIdx) {
        fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.50);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.03);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.07);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.20);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.15);
        fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 0.05);
    }
    else {
        assert(phaseIdx == FluidSystem::waterPhaseIdx);
    }
}

    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 void solveIdealMix_(FluidState &fluidState,
                               ParameterCache &paramCache,
                               int phaseIdx,
                               const ComponentVector &fugacities)
    {
        for (int i = 0; i < numComponents; ++ i) {
            const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
                                                                     paramCache,
                                                                     phaseIdx,
                                                                     i);
            const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
            Valgrind::CheckDefined(phi);
            Valgrind::CheckDefined(gamma);
            Valgrind::CheckDefined(fugacities[i]);
            fluidState.setFugacityCoefficient(phaseIdx, i, phi);
            fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
        };

        paramCache.updatePhase(fluidState, phaseIdx);

        const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
        fluidState.setDensity(phaseIdx, rho);
        return;
    }

//.........这里部分代码省略.........
                    FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
                fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
            }
        }

        // create the linear system of equations which defines the
        // mole fractions
        static const int numEq = numComponents*numPhases;
        Dune::FieldMatrix<Evaluation, numEq, numEq> M(Toolbox::createConstant(0.0));
        Dune::FieldVector<Evaluation, numEq> x(Toolbox::createConstant(0.0));
        Dune::FieldVector<Evaluation, numEq> b(Toolbox::createConstant(0.0));

        // assemble the equations expressing the fact that the
        // fugacities of each component are equal in all phases
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
            const Evaluation& entryCol1 =
                fluidState.fugacityCoefficient(/*phaseIdx=*/0, compIdx)
                *fluidState.pressure(/*phaseIdx=*/0);
            unsigned col1Idx = compIdx;

            for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
                unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
                unsigned col2Idx = phaseIdx*numComponents + compIdx;

                const Evaluation& entryCol2 =
                    fluidState.fugacityCoefficient(phaseIdx, compIdx)
                    *fluidState.pressure(phaseIdx);

                M[rowIdx][col1Idx] = entryCol1;
                M[rowIdx][col2Idx] = -entryCol2;
            }
        }

        // assemble the equations expressing the assumption that the
        // sum of all mole fractions in each phase must be 1 for the
        // phases present.
        unsigned presentPhases = 0;
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!(phasePresence & (1 << phaseIdx)))
                continue;

            unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
            presentPhases += 1;

            b[rowIdx] = Toolbox::createConstant(1.0);
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                unsigned colIdx = phaseIdx*numComponents + compIdx;

                M[rowIdx][colIdx] = Toolbox::createConstant(1.0);
            }
        }

        assert(presentPhases + numAuxConstraints == numComponents);

        // incorperate the auxiliary equations, i.e., the explicitly given mole fractions
        for (unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
            unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
            b[rowIdx] = auxConstraints[auxEqIdx].value();

            unsigned colIdx = auxConstraints[auxEqIdx].phaseIdx()*numComponents + auxConstraints[auxEqIdx].compIdx();
            M[rowIdx][colIdx] = 1.0;
        }

        // solve for all mole fractions
        try {
            Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-50);
            M.solve(x, b);
        }
        catch (const Dune::FMatrixError &e) {
            OPM_THROW(NumericalProblem,
                      "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() << "; M="<<M);
        }
        catch (...) {
            throw;
        }


        // set all mole fractions and the additional quantities in
        // the fluid state
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                unsigned rowIdx = phaseIdx*numComponents + compIdx;
                fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
            }
            paramCache.updateComposition(fluidState, phaseIdx);

            const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
            fluidState.setDensity(phaseIdx, rho);

            if (setViscosity) {
                const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
                fluidState.setViscosity(phaseIdx, mu);
            }

            if (setInternalEnergy) {
                const Evaluation& h =  FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
                fluidState.setEnthalpy(phaseIdx, h);
            }
        }
    }

    static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents> &J,
                             Dune::FieldVector<Evaluation, numComponents> &defect,
                             FluidState &fluidState,
                             ParameterCache &paramCache,
                             int phaseIdx,
                             const ComponentVector &targetFug)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        // reset jacobian
        J = 0;

        Scalar absError = 0;
        // calculate the defect (deviation of the current fugacities
        // from the target fugacities)
        for (int i = 0; i < numComponents; ++ i) {
            const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
                                                          paramCache,
                                                          phaseIdx,
                                                          i);
            const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
            fluidState.setFugacityCoefficient(phaseIdx, i, phi);

            defect[i] = targetFug[i] - f;
            absError = std::max(absError, std::abs(Toolbox::value(defect[i])));
        }

        // assemble jacobian matrix of the constraints for the composition
        static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
        for (int i = 0; i < numComponents; ++ i) {
            ////////
            // approximately calculate partial derivatives of the
            // fugacity defect of all components in regard to the mole
            // fraction of the i-th component. This is done via
            // forward differences

            // deviate the mole fraction of the i-th component
            Evaluation xI = fluidState.moleFraction(phaseIdx, i);
            fluidState.setMoleFraction(phaseIdx, i, xI + eps);
            paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);

            // compute new defect and derivative for all component
            // fugacities
            for (int j = 0; j < numComponents; ++j) {
                // compute the j-th component's fugacity coefficient ...
                const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
                                                                         paramCache,
                                                                         phaseIdx,
                                                                         j);
                // ... and its fugacity ...
                const Evaluation& f =
                    phi *
                    fluidState.pressure(phaseIdx) *
                    fluidState.moleFraction(phaseIdx, j);
                // as well as the defect for this fugacity
                const Evaluation& defJPlusEps = targetFug[j] - f;

                // use forward differences to calculate the defect's
                // derivative
                J[j][i] = (defJPlusEps - defect[j])/eps;
            }

            // reset composition to original value
            fluidState.setMoleFraction(phaseIdx, i, xI);
            paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);

            // end forward differences
            ////////
        }

        return absError;
    }

inline void testAll()
{
    typedef Opm::FluidSystems::Spe5<Scalar> FluidSystem;

    enum {
        numPhases = FluidSystem::numPhases,
        waterPhaseIdx = FluidSystem::waterPhaseIdx,
        gasPhaseIdx = FluidSystem::gasPhaseIdx,
        oilPhaseIdx = FluidSystem::oilPhaseIdx,

        numComponents = FluidSystem::numComponents,
        H2OIdx = FluidSystem::H2OIdx,
        C1Idx = FluidSystem::C1Idx,
        C3Idx = FluidSystem::C3Idx,
        C6Idx = FluidSystem::C6Idx,
        C10Idx = FluidSystem::C10Idx,
        C15Idx = FluidSystem::C15Idx,
        C20Idx = FluidSystem::C20Idx
    };

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

    typedef Opm::ThreePhaseMaterialTraits<Scalar, waterPhaseIdx, oilPhaseIdx, gasPhaseIdx> MaterialTraits;
    typedef Opm::LinearMaterial<MaterialTraits> MaterialLaw;
    typedef typename MaterialLaw::Params MaterialLawParams;

    typedef typename FluidSystem::template ParameterCache<Scalar> ParameterCache;

    ////////////
    // Initialize the fluid system and create the capillary pressure
    // parameters
    ////////////
    Scalar T = 273.15 + 20; // 20 deg Celsius
    FluidSystem::init(/*minTemperature=*/T - 1,
                      /*maxTemperature=*/T + 1,
                      /*minPressure=*/1.0e4,
                      /*maxTemperature=*/40.0e6);

    // 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();

    ////////////
    // Create a fluid state
    ////////////
    FluidState gasFluidState;
    createSurfaceGasFluidSystem<FluidSystem>(gasFluidState);

    FluidState fluidState;
    ParameterCache paramCache;

    // temperature
    fluidState.setTemperature(T);

    // oil pressure
    fluidState.setPressure(oilPhaseIdx, 4000 * 6894.7573); // 4000 PSI

    // oil saturation
    fluidState.setSaturation(oilPhaseIdx, 1.0);
    fluidState.setSaturation(gasPhaseIdx, 1.0 - fluidState.saturation(oilPhaseIdx));

    // oil composition: SPE-5 reservoir oil
    fluidState.setMoleFraction(oilPhaseIdx, H2OIdx, 0.0);
    fluidState.setMoleFraction(oilPhaseIdx, C1Idx, 0.50);
    fluidState.setMoleFraction(oilPhaseIdx, C3Idx, 0.03);
    fluidState.setMoleFraction(oilPhaseIdx, C6Idx, 0.07);
    fluidState.setMoleFraction(oilPhaseIdx, C10Idx, 0.20);
    fluidState.setMoleFraction(oilPhaseIdx, C15Idx, 0.15);
    fluidState.setMoleFraction(oilPhaseIdx, C20Idx, 0.05);

    //makeOilSaturated<Scalar, FluidSystem>(fluidState, gasFluidState);

    // set the saturations and pressures of the other phases
    for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
        if (phaseIdx != oilPhaseIdx) {
            fluidState.setSaturation(phaseIdx, 0.0);
            fluidState.setPressure(phaseIdx, fluidState.pressure(oilPhaseIdx));
        }

        // initial guess for the composition (needed by the ComputeFromReferencePhase
        // constraint solver. TODO: bug in ComputeFromReferencePhase?)
        guessInitial<FluidSystem>(fluidState, phaseIdx);
    }

    typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
    CFRP::solve(fluidState,
                paramCache,
                /*refPhaseIdx=*/oilPhaseIdx,
                /*setViscosity=*/false,
                /*setEnthalpy=*/false);

    ////////////
    // Calculate the total molarities of the components
    ////////////
//.........这里部分代码省略.........

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;
}