Skip to content

Instantly share code, notes, and snippets.

@Foadsf
Created February 15, 2022 11:15
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Foadsf/7b5e5ef2b1e82f03da0f7fd104fb34d1 to your computer and use it in GitHub Desktop.
Save Foadsf/7b5e5ef2b1e82f03da0f7fd104fb34d1 to your computer and use it in GitHub Desktop.
Implementing Brent solver instead of Newton method in thermoI.H to resolve maximum number of iterations error --> https://scicomp.stackexchange.com/questions/30080/implementing-brent-solver-instead-of-newton-method-in-thermoi-h-to-resolve-maxim

Someone in this page has claimed to successfully implement Brent solver in thermoI.H instead of Newton solver to resolve the infamous

FOAM FATAL ERROR: Maximum number of iterations exceeded

issue. The part of code he is referring to is:

template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::T
(
    scalar f,
    scalar p,
    scalar T0,
    scalar (thermo<Thermo, Type>::*F)(const scalar, const scalar) const,
    scalar (thermo<Thermo, Type>::*dFdT)(const scalar, const scalar)
        const,
    scalar (thermo<Thermo, Type>::*limit)(const scalar) const
) const
{
    scalar Test = T0;
    scalar Tnew = T0;
    scalar Ttol = T0*tol_;
    int    iter = 0;

    do
    {
        Test = Tnew;
        Tnew =
            (this->*limit)
            (Test - ((this->*F)(p, Test) - f)/(this->*dFdT)(p, Test));

        if (iter++ > maxIter_)
        {
            FatalErrorIn
            (
                "thermo<Thermo, Type>::T(scalar f, scalar T0, "
                "scalar (thermo<Thermo, Type>::*F)"
                "(const scalar) const, "
                "scalar (thermo<Thermo, Type>::*dFdT)"
                "(const scalar) const, "
                "scalar (thermo<Thermo, Type>::*limit)"
                "(const scalar) const"
                ") const"
            )   << "Maximum number of iterations exceeded"
                << abort(FatalError);
        }

    } while (mag(Tnew - Test) > Ttol);

    return Tnew;
}

I was wondering if if you could help me know

  • what is the Brent solver and how it is different than Newton method?
  • how one can implement in C++ instead of the Newton solver above?
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment