newEquilibrate.cpp

#include "cantera/equil/equil.h"
#include "cantera/equil/ChemEquil.h"
#include "cantera/equil/MultiPhaseEquil.h"
#include "cantera/equil/vcs_MultiPhaseEquil.h"
#include "cantera/base/global.h"
#include "cantera/base/stringUtils.h"

namespace Cantera
{

int new_equilibrate(thermo_t& s, const char* XY, int solver,
                    doublereal rtol, int maxsteps, int maxiter, int loglevel)
{
    // Decide which solvers to consider
    int ixy = _equilflag(XY);
    vector_int solvers;
    if (solver == -1) {
        solvers.push_back(2);
        if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) {
            solvers.push_back(1);
        }
        solvers.push_back(0);
    } else if (solver == 1) {
        if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) {
            solvers.push_back(1);
        } else {
            throw CanteraError("equilibrate",
                               std::string("MultiPhase equilibrium solver does not support property pair ") + XY);
        }
    } else if (solver == 0 || solver == 2) {
        solvers.push_back(0);
    } else {
        throw CanteraError("equilibrate",
                           "unknown solver:" + int2str(solver));
    }

    int nAttempts = 0;
    bool success = false;
    while (!solvers.empty()) {
        nAttempts++;
        solver = solvers.back();
        solvers.pop_back();

        if (solver == 0) {
            ChemEquil e;
            e.options.maxIterations = maxsteps;
            e.options.relTolerance = rtol;
            nAttempts++;
            int retnSub = e.equilibrate(s, XY, true, loglevel-1);
            if (retnSub < 0) {
                if (loglevel > 0) {
                    addLogEntry("ChemEquil solver failed.");
                }
            } else {
                if (loglevel > 0) {
                    addLogEntry("ChemEquil solver succeeded.");
                }
                return nAttempts;
            }
            s.setElementPotentials(e.elementPotentials());
        } else if (solver == 1) {
            try {
                MultiPhase m;
                m.addPhase(&s, 1.0);
                m.init();
                nAttempts++;
                std::cout << "Calling equilibrate(Multiphase)" << std::endl;
                equilibrate(m, XY, rtol, maxsteps, maxiter, loglevel-1);
                if (loglevel > 0) {
                    addLogEntry("MultiPhaseEquil solver succeeded.");
                }
                return nAttempts;
            } catch (CanteraError& err) {
                err.save();
                if (loglevel > 0) {
                    addLogEntry("MultiPhaseEquil solver failed.");
                }
                if (nAttempts < 2) {
                    if (loglevel > 0) {
                        addLogEntry("Trying single phase ChemEquil solver.");
                    }
                    solver = -1;
                } else {
                    if (loglevel > 0) {
                        endLogGroup("equilibrate");
                    }
                    throw err;
                }
            }

        }

    }

    throw CanteraError("equilibrate", "All available solvers failed");

}

}