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

void completeReferenceFluidState(FluidState &fs,
                                 typename MaterialLaw::Params &matParams,
                                 int refPhaseIdx)
{
    enum { numPhases = FluidSystem::numPhases };
    typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;

    int otherPhaseIdx = 1 - refPhaseIdx;

    // calculate the other saturation
    fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));

    // calulate the capillary pressure
    PhaseVector pC;
    MaterialLaw::capillaryPressures(pC, matParams, fs);
    fs.setPressure(otherPhaseIdx,
                   fs.pressure(refPhaseIdx)
                   + (pC[otherPhaseIdx] - pC[refPhaseIdx]));

    // set all phase densities
    typename FluidSystem::ParameterCache paramCache;
    paramCache.updateAll(fs);
    for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
        Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
        fs.setDensity(phaseIdx, rho);
    }
}

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 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 void solve(FluidState& fluidState,
                      typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
                      unsigned refPhaseIdx,
                      bool setViscosity,
                      bool setEnthalpy)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        // compute the density and enthalpy of the
        // reference phase
        paramCache.updatePhase(fluidState, refPhaseIdx);
        fluidState.setDensity(refPhaseIdx,
                              FluidSystem::density(fluidState,
                                                   paramCache,
                                                   refPhaseIdx));

        if (setEnthalpy)
            fluidState.setEnthalpy(refPhaseIdx,
                                   FluidSystem::enthalpy(fluidState,
                                                         paramCache,
                                                         refPhaseIdx));

        if (setViscosity)
            fluidState.setViscosity(refPhaseIdx,
                                    FluidSystem::viscosity(fluidState,
                                                           paramCache,
                                                           refPhaseIdx));

        // compute the fugacities of all components in the reference phase
        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
            fluidState.setFugacityCoefficient(refPhaseIdx,
                                              compIdx,
                                              FluidSystem::fugacityCoefficient(fluidState,
                                                                               paramCache,
                                                                               refPhaseIdx,
                                                                               compIdx));
        }

        // compute all quantities for the non-reference phases
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (phaseIdx == refPhaseIdx)
                continue; // reference phase is already calculated

            ComponentVector fugVec;
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                const auto& fug = fluidState.fugacity(refPhaseIdx, compIdx);
                fugVec[compIdx] = FsToolbox::template decay<Evaluation>(fug);
            }

            CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fugVec);

            if (setViscosity)
                fluidState.setViscosity(phaseIdx,
                                        FluidSystem::viscosity(fluidState,
                                                               paramCache,
                                                               phaseIdx));

            if (setEnthalpy)
                fluidState.setEnthalpy(phaseIdx,
                                       FluidSystem::enthalpy(fluidState,
                                                             paramCache,
                                                             phaseIdx));
        }
    }

    static void solve(FluidState &fluidState,
                      ParameterCache &paramCache,
                      int phaseIdx,
                      const ComponentVector &targetFug)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        // use a much more efficient method in case the phase is an
        // ideal mixture
        if (FluidSystem::isIdealMixture(phaseIdx)) {
            solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
            return;
        }

        //Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-25);

        // save initial composition in case something goes wrong
        Dune::FieldVector<Evaluation, numComponents> xInit;
        for (int i = 0; i < numComponents; ++i) {
            xInit[i] = fluidState.moleFraction(phaseIdx, i);
        }

        /////////////////////////
        // Newton method
        /////////////////////////

        // Jacobian matrix
        Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
        // solution, i.e. phase composition
        Dune::FieldVector<Evaluation, numComponents> x;
        // right hand side
        Dune::FieldVector<Evaluation, numComponents> b;

        paramCache.updatePhase(fluidState, phaseIdx);

        // maximum number of iterations
        const int nMax = 25;
        for (int nIdx = 0; nIdx < nMax; ++nIdx) {
            // calculate Jacobian matrix and right hand side
            linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
            Valgrind::CheckDefined(J);
            Valgrind::CheckDefined(b);

            /*
            std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
            for (int i = 0; i < FluidSystem::numComponents; ++i)
                std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
            std::cout << "\n";
            std::cout << FluidSystem::phaseName(phaseIdx) << "Phase phi: ";
            for (int i = 0; i < FluidSystem::numComponents; ++i)
                std::cout << fluidState.fugacityCoefficient(phaseIdx, i) << " ";
            std::cout << "\n";
            */

            // Solve J*x = b
            x = Toolbox::createConstant(0.0);
            try { J.solve(x, b); }
            catch (Dune::FMatrixError e)
            { throw Opm::NumericalIssue(e.what()); }

            //std::cout << "original delta: " << x << "\n";

            Valgrind::CheckDefined(x);

            /*
            std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
            for (int i = 0; i < FluidSystem::numComponents; ++i)
                std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
            std::cout << "\n";
            std::cout << "J: " << J << "\n";
            std::cout << "rho: " << fluidState.density(phaseIdx) << "\n";
            std::cout << "delta: " << x << "\n";
            std::cout << "defect: " << b << "\n";

            std::cout << "J: " << J << "\n";

            std::cout << "---------------------------\n";
            */

            // update the fluid composition. b is also used to store
            // the defect for the next iteration.
            Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);

            if (relError < 1e-9) {
                const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
                fluidState.setDensity(phaseIdx, rho);

                //std::cout << "num iterations: " << nIdx << "\n";
                return;
            }
        }

        OPM_THROW(Opm::NumericalIssue,
                  "Calculating the " << FluidSystem::phaseName(phaseIdx)
                  << "Phase composition failed. Initial {x} = {"
                  << xInit
                  << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx)
                  << ", T = " << fluidState.temperature(phaseIdx));
    }

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