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

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

    static void solve(FluidState &fluidState,
                      ParameterCache &paramCache,
                      int phasePresence,
                      const MMPCAuxConstraint<Evaluation> *auxConstraints,
                      unsigned numAuxConstraints,
                      bool setViscosity,
                      bool setInternalEnergy)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
        static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
                      "The scalar type of the fluid state must be 'Evaluation'");

#ifndef NDEBUG
        // currently this solver can only handle fluid systems which
        // assume ideal mixtures of all fluids. TODO: relax this
        // (requires solving a non-linear system of equations, i.e. using
        // newton method.)
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            assert(FluidSystem::isIdealMixture(phaseIdx));
        }
#endif

        // compute all fugacity coefficients
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            paramCache.updatePhase(fluidState, phaseIdx);

            // since we assume ideal mixtures, the fugacity
            // coefficients of the components cannot depend on
            // composition, i.e. the parameters in the cache are valid
            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                Evaluation fugCoeff = FsToolbox::template toLhs<Evaluation>(
                    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) {
//.........这里部分代码省略.........

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