View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002578OpenFOAMContributionpublic2017-06-09 13:052017-07-16 21:03
Reportertniemi Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status closedResolutionsuspended 
Product Versiondev 
Summary0002578: Comments regarding chemistryModel jacobian
DescriptionNote: this is not strictly a bug report, but more of a commentary/sharing of experiences

Lately, I have been testing the EDC and TDAC with the GRI-mechanism. Generally, I have had good success in some test cases, but in others the chemistry solution tends to crash quite often.

I recently saw a presentation, where it was shown that much better stability and performance could be obtained by replacing the jacobian in the chemistryModel with an analytical version obtained from pyJac. (There were also other changes such as replacing LUSolver with Lapack) Based on this, I looked at the jacobian-function to see what might cause stability/performance problems.

If I have understood correctly, currently the jacobian does the following:
- First calculates dcdt
- Then calculates (semianalytically) d(dcdt)dc terms
- Finally calculates d(dcdt)dT numerically

However, I can see several potential issues:
1. According to ODESystem, the jacobian is supposed to return dfdx, which I think would mean d(dcdt)dt. However, here the jacobian calculates dcdt, which would correspond to dydt. (Also, dcdt contains space for dTdt which is simply ignored here.)
2. When calculating d(dcdt)/dc terms, it is assumed that the R.kf and R.kr do not have derivatives with respect to concentrations. This is not true with thirdBodyReactions or FallOffReactions.
3. d(dTdt)/dc terms and d(dTdt)/dT are assumed to be zero, which is not generally true.

I guess the jacobian doesn't need to be very accurate, so are these issues intentional optimizations or bugs?

I made some tests where I included the derivatives of R.kf and R.kr in the jacobian calculated with simple forward finite difference. I have attached a comparison of the chemFoam gri-tutorial case. As can be seen, the amount of integration steps clearly drops with high tolerances. However, the jacobian evaluation is now more costly and depending on tolerances, the overall performance may increase or decrease.

I have attached the modified chemistryModel, which I used in the tests. It is hard code optimized for GRI, but it could be easily made more general if the reaction rates would hold a flag telling whether or not they have a non-zero derivative wrt concentrations. Even better would be, if the reaction rates could calculate their own derivatives, but this would require more work.

There might also be some other optimizations possible. I tested caching the Kc value when calculating R.kr derivatives and that already shaved off 10 % of the computation time.

I also tried to address concerns 1. and 3., but they didn't seem to affect the results much.
TagsNo tags attached.

Activities

tniemi

2017-06-09 13:05

reporter  

comparison.png (7,908 bytes)   
comparison.png (7,908 bytes)   

tniemi

2017-06-09 13:06

reporter  

chemistryModel.C (20,113 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "chemistryModel.H"
#include "reactingMixture.H"
#include "UniformField.H"
#include "extrapolatedCalculatedFvPatchFields.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CompType, class ThermoType>
Foam::chemistryModel<CompType, ThermoType>::chemistryModel
(
    const fvMesh& mesh,
    const word& phaseName
)
:
    CompType(mesh, phaseName),
    ODESystem(),
    Y_(this->thermo().composition().Y()),
    reactions_
    (
        dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
    ),
    specieThermo_
    (
        dynamic_cast<const reactingMixture<ThermoType>&>
            (this->thermo()).speciesData()
    ),

    nSpecie_(Y_.size()),
    nReaction_(reactions_.size()),
    Treact_(CompType::template lookupOrDefault<scalar>("Treact", 0.0)),
    RR_(nSpecie_),
    c_(nSpecie_),
    dcdt_(nSpecie_)
{
    // Create the fields for the chemistry sources
    forAll(RR_, fieldi)
    {
        RR_.set
        (
            fieldi,
            new volScalarField::Internal
            (
                IOobject
                (
                    "RR." + Y_[fieldi].name(),
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh,
                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
            )
        );
    }

    Info<< "chemistryModel: Number of species = " << nSpecie_
        << " and reactions = " << nReaction_ << endl;
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CompType, class ThermoType>
Foam::chemistryModel<CompType, ThermoType>::~chemistryModel()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::omega
(
    const scalarField& c,
    const scalar T,
    const scalar p,
    scalarField& dcdt
) const
{
    scalar pf, cf, pr, cr;
    label lRef, rRef;

    dcdt = Zero;

    forAll(reactions_, i)
    {
        const Reaction<ThermoType>& R = reactions_[i];

        scalar omegai = omega
        (
            R, c, T, p, pf, cf, lRef, pr, cr, rRef
        );

        forAll(R.lhs(), s)
        {
            const label si = R.lhs()[s].index;
            const scalar sl = R.lhs()[s].stoichCoeff;
            dcdt[si] -= sl*omegai;
        }

        forAll(R.rhs(), s)
        {
            const label si = R.rhs()[s].index;
            const scalar sr = R.rhs()[s].stoichCoeff;
            dcdt[si] += sr*omegai;
        }
    }
}


template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omegaI
(
    const label index,
    const scalarField& c,
    const scalar T,
    const scalar p,
    scalar& pf,
    scalar& cf,
    label& lRef,
    scalar& pr,
    scalar& cr,
    label& rRef
) const
{
    const Reaction<ThermoType>& R = reactions_[index];
    scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
    return(w);
}


template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omega
(
    const Reaction<ThermoType>& R,
    const scalarField& c,
    const scalar T,
    const scalar p,
    scalar& pf,
    scalar& cf,
    label& lRef,
    scalar& pr,
    scalar& cr,
    label& rRef
) const
{
    const scalar kf = R.kf(p, T, c);
    const scalar kr = R.kr(kf, p, T, c);

    pf = 1.0;
    pr = 1.0;

    const label Nl = R.lhs().size();
    const label Nr = R.rhs().size();

    label slRef = 0;
    lRef = R.lhs()[slRef].index;

    pf = kf;
    for (label s = 1; s < Nl; s++)
    {
        const label si = R.lhs()[s].index;

        if (c[si] < c[lRef])
        {
            const scalar exp = R.lhs()[slRef].exponent;
            pf *= pow(max(0.0, c[lRef]), exp);
            lRef = si;
            slRef = s;
        }
        else
        {
            const scalar exp = R.lhs()[s].exponent;
            pf *= pow(max(0.0, c[si]), exp);
        }
    }
    cf = max(0.0, c[lRef]);

    {
        const scalar exp = R.lhs()[slRef].exponent;
        if (exp < 1.0)
        {
            if (cf > SMALL)
            {
                pf *= pow(cf, exp - 1.0);
            }
            else
            {
                pf = 0.0;
            }
        }
        else
        {
            pf *= pow(cf, exp - 1.0);
        }
    }

    label srRef = 0;
    rRef = R.rhs()[srRef].index;

    // Find the matrix element and element position for the rhs
    pr = kr;
    for (label s = 1; s < Nr; s++)
    {
        const label si = R.rhs()[s].index;
        if (c[si] < c[rRef])
        {
            const scalar exp = R.rhs()[srRef].exponent;
            pr *= pow(max(0.0, c[rRef]), exp);
            rRef = si;
            srRef = s;
        }
        else
        {
            const scalar exp = R.rhs()[s].exponent;
            pr *= pow(max(0.0, c[si]), exp);
        }
    }
    cr = max(0.0, c[rRef]);

    {
        const scalar exp = R.rhs()[srRef].exponent;
        if (exp < 1.0)
        {
            if (cr>SMALL)
            {
                pr *= pow(cr, exp - 1.0);
            }
            else
            {
                pr = 0.0;
            }
        }
        else
        {
            pr *= pow(cr, exp - 1.0);
        }
    }

    return pf*cf - pr*cr;
}


template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::derivatives
(
    const scalar time,
    const scalarField& c,
    scalarField& dcdt
) const
{
    const scalar T = c[nSpecie_];
    const scalar p = c[nSpecie_ + 1];

    for (label i = 0; i < nSpecie_; i++)
    {
        c_[i] = max(0.0, c[i]);
    }

    omega(c_, T, p, dcdt);

    // Constant pressure
    // dT/dt = ...
    scalar rho = 0.0;
    scalar cSum = 0.0;
    for (label i = 0; i < nSpecie_; i++)
    {
        const scalar W = specieThermo_[i].W();
        cSum += c_[i];
        rho += W*c_[i];
    }
    scalar cp = 0.0;
    for (label i=0; i<nSpecie_; i++)
    {
        cp += c_[i]*specieThermo_[i].cp(p, T);
    }
    cp /= rho;

    scalar dT = 0.0;
    for (label i = 0; i < nSpecie_; i++)
    {
        const scalar hi = specieThermo_[i].ha(p, T);
        dT += hi*dcdt[i];
    }
    dT /= rho*cp;

    dcdt[nSpecie_] = -dT;

    // dp/dt = ...
    dcdt[nSpecie_ + 1] = 0.0;
}


template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::jacobian
(
    const scalar t,
    const scalarField& c,
    scalarField& dcdt,
    scalarSquareMatrix& dfdc
) const
{
    const scalar T = c[nSpecie_];
    const scalar p = c[nSpecie_ + 1];

    forAll(c_, i)
    {
        c_[i] = max(c[i], 0.0);
    }

    dfdc = Zero;

    // Length of the first argument must be nSpecie_
    //omega(c_, T, p, dcdt);

    dcdt = Zero; //Does not affect much

    forAll(reactions_, ri)
    {
        const Reaction<ThermoType>& R = reactions_[ri];

        const scalar kf0 = R.kf(p, T, c_);
        const scalar kr0 = R.kr(kf0, p, T, c_);

        scalar kfp = 1;
        scalar krp = 1;

        forAll(R.lhs(), j)
        {
            const label sj = R.lhs()[j].index;
            scalar kf = kf0;
            forAll(R.lhs(), i)
            {
                const label si = R.lhs()[i].index;
                const scalar el = R.lhs()[i].exponent;
                if (i == j)
                {
                    if (el < 1.0)
                    {
                        if (c_[si] > SMALL)
                        {
                            kf *= el*pow(c_[si], el - 1.0);
                        }
                        else
                        {
                            kf = 0.0;
                        }
                    }
                    else
                    {
                        kf *= el*pow(c_[si], el - 1.0);
                    }
                }
                else
                {
                    kf *= pow(c_[si], el);
                }
            }

            forAll(R.lhs(), i)
            {
                const label si = R.lhs()[i].index;
                const scalar sl = R.lhs()[i].stoichCoeff;
                dfdc(si, sj) -= sl*kf;
            }
            forAll(R.rhs(), i)
            {
                const label si = R.rhs()[i].index;
                const scalar sr = R.rhs()[i].stoichCoeff;
                dfdc(si, sj) += sr*kf;
            }

            const scalar el = R.lhs()[j].exponent;
            kfp *= pow(c_[sj], el);
        }

        forAll(R.rhs(), j)
        {
            const label sj = R.rhs()[j].index;
            scalar kr = kr0;
            forAll(R.rhs(), i)
            {
                const label si = R.rhs()[i].index;
                const scalar er = R.rhs()[i].exponent;
                if (i == j)
                {
                    if (er < 1.0)
                    {
                        if (c_[si] > SMALL)
                        {
                            kr *= er*pow(c_[si], er - 1.0);
                        }
                        else
                        {
                            kr = 0.0;
                        }
                    }
                    else
                    {
                        kr *= er*pow(c_[si], er - 1.0);
                    }
                }
                else
                {
                    kr *= pow(c_[si], er);
                }
            }

            forAll(R.lhs(), i)
            {
                const label si = R.lhs()[i].index;
                const scalar sl = R.lhs()[i].stoichCoeff;
                dfdc(si, sj) += sl*kr;
            }
            forAll(R.rhs(), i)
            {
                const label si = R.rhs()[i].index;
                const scalar sr = R.rhs()[i].stoichCoeff;
                dfdc(si, sj) -= sr*kr;
            }

            const scalar el = R.rhs()[j].exponent;
            krp *= pow(c_[sj], el);
        }

        //Conditional check is optimization for GRI
        //The thirdBodyReactions are at the end of the reaction list
        if (ri < 284) continue;

        //Derivatives of kf and kr with respect to concentrations
        //Non-zero for eg. thirdBodyReactions
        scalarList kf0p(nSpecie_);
        scalarList kr0p(nSpecie_);

        // kf and kr derivatives using finite difference
        const scalar delta = 1e-8;
        forAll(c_,i)
        {
            const scalar temp = c_[i];
            c_[i] += delta;
            const scalar kf1 = R.kf(p, T, c_);
            kf0p[i] = (kf1 - kf0) / delta;
            kr0p[i] = (R.kr(kf1, p, T, c_) - kr0) / delta;
            c_[i] = temp;
        }

        forAll(R.lhs(), j)
        {
            const label sj = R.lhs()[j].index;
            const scalar sl = R.lhs()[j].stoichCoeff;
            forAll(c_, i)
            {
                dfdc(sj, i) -= sl*kfp*kf0p[i];
                dfdc(sj, i) += sl*krp*kr0p[i];
            }
        }

        forAll(R.rhs(), j)
        {
            const label sj = R.rhs()[j].index;
            const scalar sl = R.rhs()[j].stoichCoeff;
            forAll(c_, i)
            {
                dfdc(sj, i) -= sl*krp*kr0p[i];
                dfdc(sj, i) += sl*kfp*kf0p[i];
            }
        }
    }

    // Calculate the dcdT elements numerically
    const scalar delta = 1.0e-3;

    omega(c_, T + delta, p, dcdt_);
    for (label i=0; i<nSpecie_; i++)
    {
        dfdc(i, nSpecie_) = dcdt_[i];
    }

    omega(c_, T - delta, p, dcdt_);
    for (label i=0; i<nSpecie_; i++)
    {
        dfdc(i, nSpecie_) = 0.5*(dfdc(i, nSpecie_) - dcdt_[i])/delta;
    }

    dfdc(nSpecie_, nSpecie_) = 0;
    dfdc(nSpecie_ + 1, nSpecie_) = 0;
}


template<class CompType, class ThermoType>
Foam::tmp<Foam::volScalarField>
Foam::chemistryModel<CompType, ThermoType>::tc() const
{
    tmp<volScalarField> ttc
    (
        new volScalarField
        (
            IOobject
            (
                "tc",
                this->time().timeName(),
                this->mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            this->mesh(),
            dimensionedScalar("zero", dimTime, SMALL),
            extrapolatedCalculatedFvPatchScalarField::typeName
        )
    );

    scalarField& tc = ttc.ref();

    tmp<volScalarField> trho(this->thermo().rho());
    const scalarField& rho = trho();

    const scalarField& T = this->thermo().T();
    const scalarField& p = this->thermo().p();

    const label nReaction = reactions_.size();

    scalar pf, cf, pr, cr;
    label lRef, rRef;

    if (this->chemistry_)
    {
        forAll(rho, celli)
        {
            const scalar rhoi = rho[celli];
            const scalar Ti = T[celli];
            const scalar pi = p[celli];

            scalar cSum = 0.0;

            for (label i=0; i<nSpecie_; i++)
            {
                c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
                cSum += c_[i];
            }

            forAll(reactions_, i)
            {
                const Reaction<ThermoType>& R = reactions_[i];

                omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);

                forAll(R.rhs(), s)
                {
                    tc[celli] += R.rhs()[s].stoichCoeff*pf*cf;
                }
            }

            tc[celli] = nReaction*cSum/tc[celli];
        }
    }

    ttc.ref().correctBoundaryConditions();

    return ttc;
}


template<class CompType, class ThermoType>
Foam::tmp<Foam::volScalarField>
Foam::chemistryModel<CompType, ThermoType>::Qdot() const
{
    tmp<volScalarField> tQdot
    (
        new volScalarField
        (
            IOobject
            (
                "Qdot",
                this->mesh_.time().timeName(),
                this->mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            this->mesh_,
            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
        )
    );

    if (this->chemistry_)
    {
        scalarField& Qdot = tQdot.ref();

        forAll(Y_, i)
        {
            forAll(Qdot, celli)
            {
                const scalar hi = specieThermo_[i].Hc();
                Qdot[celli] -= hi*RR_[i][celli];
            }
        }
    }

    return tQdot;
}


template<class CompType, class ThermoType>
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
Foam::chemistryModel<CompType, ThermoType>::calculateRR
(
    const label ri,
    const label si
) const
{
    scalar pf, cf, pr, cr;
    label lRef, rRef;

    tmp<volScalarField::Internal> tRR
    (
        new volScalarField::Internal
        (
            IOobject
            (
                "RR",
                this->mesh().time().timeName(),
                this->mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            this->mesh(),
            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
        )
    );

    volScalarField::Internal& RR = tRR.ref();

    tmp<volScalarField> trho(this->thermo().rho());
    const scalarField& rho = trho();

    const scalarField& T = this->thermo().T();
    const scalarField& p = this->thermo().p();

    forAll(rho, celli)
    {
        const scalar rhoi = rho[celli];
        const scalar Ti = T[celli];
        const scalar pi = p[celli];

        for (label i=0; i<nSpecie_; i++)
        {
            const scalar Yi = Y_[i][celli];
            c_[i] = rhoi*Yi/specieThermo_[i].W();
        }

        const scalar w = omegaI
        (
            ri,
            c_,
            Ti,
            pi,
            pf,
            cf,
            lRef,
            pr,
            cr,
            rRef
        );

        RR[celli] = w*specieThermo_[si].W();
    }

    return tRR;
}


template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::calculate()
{
    if (!this->chemistry_)
    {
        return;
    }

    tmp<volScalarField> trho(this->thermo().rho());
    const scalarField& rho = trho();

    const scalarField& T = this->thermo().T();
    const scalarField& p = this->thermo().p();

    forAll(rho, celli)
    {
        const scalar rhoi = rho[celli];
        const scalar Ti = T[celli];
        const scalar pi = p[celli];

        for (label i=0; i<nSpecie_; i++)
        {
            const scalar Yi = Y_[i][celli];
            c_[i] = rhoi*Yi/specieThermo_[i].W();
        }

        omega(c_, Ti, pi, dcdt_);

        for (label i=0; i<nSpecie_; i++)
        {
            RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
        }
    }
}


template<class CompType, class ThermoType>
template<class DeltaTType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
(
    const DeltaTType& deltaT
)
{
    CompType::correct();

    scalar deltaTMin = GREAT;

    if (!this->chemistry_)
    {
        return deltaTMin;
    }

    tmp<volScalarField> trho(this->thermo().rho());
    const scalarField& rho = trho();

    const scalarField& T = this->thermo().T();
    const scalarField& p = this->thermo().p();

    scalarField c0(nSpecie_);

    forAll(rho, celli)
    {
        scalar Ti = T[celli];

        if (Ti > Treact_)
        {
            const scalar rhoi = rho[celli];
            scalar pi = p[celli];

            for (label i=0; i<nSpecie_; i++)
            {
                c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
                c0[i] = c_[i];
            }

            // Initialise time progress
            scalar timeLeft = deltaT[celli];

            // Calculate the chemical source terms
            while (timeLeft > SMALL)
            {
                scalar dt = timeLeft;
                this->solve(c_, Ti, pi, dt, this->deltaTChem_[celli]);
                timeLeft -= dt;
            }

            deltaTMin = min(this->deltaTChem_[celli], deltaTMin);

            for (label i=0; i<nSpecie_; i++)
            {
                RR_[i][celli] =
                    (c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
            }
        }
        else
        {
            for (label i=0; i<nSpecie_; i++)
            {
                RR_[i][celli] = 0;
            }
        }
    }

    return deltaTMin;
}


template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
(
    const scalar deltaT
)
{
    // Don't allow the time-step to change more than a factor of 2
    return min
    (
        this->solve<UniformField<scalar>>(UniformField<scalar>(deltaT)),
        2*deltaT
    );
}


template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
(
    const scalarField& deltaT
)
{
    return this->solve<scalarField>(deltaT);
}


// ************************************************************************* //
chemistryModel.C (20,113 bytes)   

tniemi

2017-06-13 14:21

reporter   ~0008213

Related to optimizations, I have also tried to look on how to optimize the LUDecompose method. The current implementation looked rather complicated, so I tried to create a more simplified version to see if it would be faster.

I have attached alternative versions of LUDecompose and LUBacksubstitute and an artificial benhmark testcase. In the testcase the alternative implementation is roughly 30-40 % faster than the current implementation for medium to large matrices (n>10). However, the performance depends largely on compiler optimizations, so the results may vary. I only used a single compiler (gcc 4.8.3), but several different machines. On an older machine (Intel X5650), the performance of the stock LUDecompose seemed to drastically decrease when n>500 (probably cpu cache limitation) and for these large matrices the performance difference was 300-400 %. On newer machines the performance difference stayed at 30-40 %.

I also benchmarked the alternative implementation against Eigen, which is supposedly a very optimized linear algebra library. For small matrices (n<50) there were no diffences, but with n~1000 Eigen was about 3 times faster than the alternative LUDecompose. So for large matrices there is still room to the state of the art. Of course, for such a large systems using sparse matrices would be more reasonable.

I have also attached a comparison of all the chemFoam tutorials with the stock settings and with the alternative LUSolver. For h2 and gri there is practically no difference, but on larger cases the difference is noticeable.

Note that the new implementation uses different decomposition algo, so the decomposition is not identical with the stock one.

tniemi

2017-06-13 14:21

reporter  

LUcomparison.png (3,073 bytes)   
LUcomparison.png (3,073 bytes)   

tniemi

2017-06-13 14:22

reporter  

benchmark.zip (2,029 bytes)

tniemi

2017-06-13 14:22

reporter  

alt_LUDecompose.zip (5,696 bytes)

henry

2017-06-14 10:50

manager   ~0008219

> replacing the jacobian in the chemistryModel with an analytical version

This would be a good idea and would not require interfacing to a Python library to achieve. Someone would need to spend a bit of time working out the form of the temperature derivative which would be a bit tedious but not fundamentally difficult.

> I guess the jacobian doesn't need to be very accurate, so are these issues intentional optimizations or bugs?

Simplifications rather than bugs, and not explicit optimizations. If you would like to contribute a more complete form of the evaluation of the Jacobian I would be happy to include it.

> I have attached alternative versions of LUDecompose and LUBacksubstitute

I would be happy to include these developments if you are sure that they do not adversly affect the operation of other parts of the code which rely on them. What tests have you run so far?

tniemi

2017-06-16 13:29

reporter  

alt_LUDecompose_fixed.zip (5,696 bytes)

tniemi

2017-06-16 13:39

reporter   ~0008225

> Simplifications rather than bugs, and not explicit optimizations. If you would like to contribute a more complete form of the evaluation of the Jacobian I would be happy to include it.

Unfortunately, I don't currently have enough time to do this properly. Besides chemistryModel, similar changes should also be made to TDACChemistryModel, which is slightly more complex due to the mechanism reduction.

> I would be happy to include these developments if you are sure that they do not adversly affect the operation of other parts of the code which rely on them. What tests have you run so far?

I have run the chemFoam tutorials and the Matrix-Test program. I have also checked that in the benchmark tests both the new and old decompositions give similar solutions to random linear systems. However, I did spot a slight mistake in the submitted version and I have uploaded a new one. In the previous version I was not taking the abs of the diagonal element when checking for the largest pivot. This led to unneccessary pivoting and reduced stability in cases where the diagonal was negative.

After the fix I have not observed any performance degradations or other issues. However, apart from the chemFoam tests and benchmarking tests I have not done other testing. Also, I have only tried gcc compilers, so I don't know whats the performance difference with eg. clang.

Below is a list of known differences between the implementations:

- The methods produce a slightly different decompostions and pivots. Both are valid, but they do not produce exactly the same results.
- Both methods use partial pivoting so their stability should be roughly the same.
- If the stock implementation encounters a singular matrix during pivoting, it will replace the zero with SMALL. This allows the code to always continue, although the solution of the linear system will likely be garbage. The alternative implementation will halt in these situations.

henry

2017-06-16 17:12

manager   ~0008226

While it is likely that the new LUDecompose and LUBacksubstitute will work fine for other the purposes they are used for in OpenFOAM this will need to be tested before commiting to avoid surprises for the other users of OpenFOAM.

henry

2017-07-11 09:52

manager   ~0008380

Have you completed the testing of this proposed change? Should I close this report for now and you can re-open when the testing is complete?

tniemi

2017-07-16 19:47

reporter   ~0008401

No, I haven't yet done any more work for this. You can close the report for now.

henry

2017-07-16 21:03

manager   ~0008402

Pending completion of the testing for the contribution.

Issue History

Date Modified Username Field Change
2017-06-09 13:05 tniemi New Issue
2017-06-09 13:05 tniemi File Added: comparison.png
2017-06-09 13:06 tniemi File Added: chemistryModel.C
2017-06-13 14:21 tniemi Note Added: 0008213
2017-06-13 14:21 tniemi File Added: LUcomparison.png
2017-06-13 14:22 tniemi File Added: benchmark.zip
2017-06-13 14:22 tniemi File Added: alt_LUDecompose.zip
2017-06-14 10:50 henry Note Added: 0008219
2017-06-16 13:29 tniemi File Added: alt_LUDecompose_fixed.zip
2017-06-16 13:39 tniemi Note Added: 0008225
2017-06-16 17:12 henry Note Added: 0008226
2017-07-11 09:52 henry Note Added: 0008380
2017-07-16 19:47 tniemi Note Added: 0008401
2017-07-16 21:03 henry Assigned To => henry
2017-07-16 21:03 henry Status new => closed
2017-07-16 21:03 henry Resolution open => suspended
2017-07-16 21:03 henry Note Added: 0008402