PhreeqcRM.h

Go to the documentation of this file.

 
#if !defined(PHREEQCRM_H_INCLUDED)
#define PHREEQCRM_H_INCLUDED
#ifdef USE_MPI
#include "mpi.h"
#define MP_TYPE MPI_Comm
#else
#define MP_TYPE int
#endif
 
// forward declarations
class cxxNameDouble;
class cxxSolution;
class cxxStorageBin;
class IPhreeqc;
class IPhreeqcPhast;
class PHRQ_io;
 
#if defined(SWIG) || defined(swig_python_EXPORTS)
// Forward-declare CPython's PyObject so headers can refer to it without
// including Python.h here. This must be the global PyObject type (not a
// nested class-scoped typedef) so it matches the interpreter's API and
// SWIG-generated wrappers.
typedef struct _object PyObject;
#endif
 
#include "NameDouble.h"
#include <vector>
#include <list>
#include <set>
#include <map>
#include <mutex>
#include <string>
#include <memory>
#include "RMVARS.h"
 
#include "irm_dll_export.h"
 
#if SWIG
#define IRM_DLL_EXPORT
#endif
 
//class BMI_Var;
class IRM_DLL_EXPORT PhreeqcRMStop : public std::exception
{
public:
  const char *what() const throw () {return "Failure in PhreeqcRM\n";}
};
//class LetItThrow : public std::logic_error {
//public:
//  LetItThrow(std::string error_string) : std::logic_error(error_string.c_str()) { };
//};
#include "IrmResult.h"
enum {
    METHOD_CONSTRUCT,
    METHOD_CREATEMAPPING,
    METHOD_DUMPMODULE,
    METHOD_FINDCOMPONENTS,
    METHOD_GETCONCENTRATIONS,
    METHOD_GETDENSITYCALCULATED,
    METHOD_GETERRORSTRING,
    METHOD_GETGASCOMPMOLES,
    METHOD_GETGASCOMPPRESSURES,
    METHOD_GETGASCOMPPHI,
    METHOD_GETGASPHASEVOLUME,
    METHOD_GETPRESSURE,
    METHOD_GETSATURATIONCALCULATED,
    METHOD_GETSELECTEDOUTPUT,
    METHOD_GETSOLUTIONVOLUME,
    METHOD_GETSPECIESCONCENTRATIONS,
    METHOD_GETSPECIESLOG10GAMMAS,
    METHOD_GETSPECIESLOG10MOLALITIES,
    METHOD_GETTEMPERATURE,
    METHOD_GETVISCOSITY,
    METHOD_INITIALPHREEQC2MODULE,
    METHOD_INITIALPHREEQCCELL2MODULE,
    METHOD_LOADDATABASE,
    METHOD_MPIWORKERBREAK,
    METHOD_RUNCELLS,
    METHOD_RUNFILE,
    METHOD_RUNSTRING,
    METHOD_SETCOMPONENTH2O,
    METHOD_SETCONCENTRATIONS,
    METHOD_SETDENSITYUSER,
    METHOD_SETERRORHANDLERMODE,
    METHOD_SETFILEPREFIX,
    METHOD_SETGASCOMPMOLES,
    METHOD_SETGASPHASEVOLUME,
    METHOD_SETPARTITIONUZSOLIDS,
    METHOD_SETPOROSITY,
    METHOD_SETPRESSURE,
    METHOD_SETPRINTCHEMISTRYON,
    METHOD_SETPRINTCHEMISTRYMASK,
    METHOD_SETREBALANCEBYCELL,
    METHOD_SETREBALANCEFRACTION,
    METHOD_SETREPRESENTATIVEVOLUME,
    METHOD_SETSATURATIONUSER,
    METHOD_SETSELECTEDOUTPUTON,
    METHOD_SETSPECIESSAVEON,
    METHOD_SETTEMPERATURE,
    METHOD_SETTIME,
    METHOD_SETTIMECONVERSION,
    METHOD_SETTIMESTEP,
    METHOD_SETUNITSEXCHANGE,
    METHOD_SETUNITSGASPHASE,
    METHOD_SETUNITSKINETICS,
    METHOD_SETUNITSPPASSEMBLAGE,
    METHOD_SETUNITSSOLUTION,
    METHOD_SETUNITSSSASSEMBLAGE,
    METHOD_SETUNITSSURFACE,
    METHOD_SPECIESCONCENTRATIONS2MODULE,
    METHOD_STATESAVE,
    METHOD_STATEAPPLY,
    METHOD_STATEDELETE,
    METHOD_USESOLUTIONDENSITYVOLUME
} /* MPI_METHOD */;
 
 
template<typename T>
class StaticIndexer
{
public:
    StaticIndexer(T* self)
    {
        const std::lock_guard<std::mutex> lock(StaticIndexer<T>::_InstancesLock);
        this->_Index = StaticIndexer<T>::_InstancesIndex++;
        StaticIndexer<T>::_Instances[this->_Index] = self;
    }
 
    ~StaticIndexer()
    {
        const std::lock_guard<std::mutex> lock(StaticIndexer<T>::_InstancesLock);
        auto search = StaticIndexer<T>::_Instances.find(this->_Index);
        assert(search != StaticIndexer<T>::_Instances.end());
        if (search != StaticIndexer<T>::_Instances.end())
        {
            StaticIndexer<T>::_Instances.erase(search);
        }
    }
 
    static T* GetInstance(int id)
    {
        const std::lock_guard<std::mutex> lock(_InstancesLock);
        auto search = StaticIndexer<T>::_Instances.find(size_t(id));
        if (search != StaticIndexer<T>::_Instances.end())
        {
            return search->second;
        }
        return nullptr;
    }
 
    template<typename Derived>
    static Derived* GetInstance(int id)
    {
        const std::lock_guard<std::mutex> lock(_InstancesLock);
        auto search = StaticIndexer<T>::_Instances.find(size_t(id));
        if (search != StaticIndexer<T>::_Instances.end())
        {
            if (Derived* derived = dynamic_cast<Derived*>(search->second))
            {
                return derived;
            }
        }
        return nullptr;
    }
 
    static IRM_RESULT Destroy(int id)
    {
        auto search = StaticIndexer<T>::_Instances.find(size_t(id));
        if (search != StaticIndexer<T>::_Instances.end())
        {
            assert(dynamic_cast<T*>(search->second));
            delete search->second;
            return IRM_OK;
        }
        return IRM_BADINSTANCE;
    }
 
    //template<typename Derived>
    //static IRM_RESULT Destroy(int id)
    //{
    //  auto search = StaticIndexer<T>::_Instances.find(size_t(id));
    //  if (search != StaticIndexer<T>::_Instances.end())
    //  {
    //      assert(dynamic_cast<T*>(search->second));
    //      if (Derived* derived = dynamic_cast<Derived*>(search->second))
    //      {
    //          delete derived;
    //          return IRM_OK;
    //      }
    //  }
    //  return IRM_BADINSTANCE;
    //}
 
    static void DestroyAll()
    {
        std::list<T*> items;
        for (auto pair : StaticIndexer<T>::_Instances)
        {
            assert(dynamic_cast<T*>(pair.second));
            items.push_back(pair.second);
        }
        for (auto item : items)
        {
            delete item;
        }
    }
 
    //template<typename Derived>
    //static void DestroyAll()
    //{
    //  const std::lock_guard<std::mutex> lock(_InstancesLock);
    //  std::list<Derived*> derived_items;
    //  for (auto pair : StaticIndexer<T>::_Instances)
    //  {
    //      assert(dynamic_cast<T*>(pair.second));
    //      if (Derived* derived = dynamic_cast<Derived*>(pair.second))
    //      {
    //          derived_items.push_back(derived);
    //      }
    //  }
    //  for (auto item : derived_items)
    //  {
    //      delete item;
    //  }
    //}
 
    int GetIndex()
    {
        return (int)this->_Index;
    }
 
protected:
    size_t _Index;
 
    static std::mutex           _InstancesLock;
    static std::map<size_t, T*> _Instances;
    static size_t               _InstancesIndex;
};
 
// static members of StaticIndexer
 
template<typename T>
std::mutex StaticIndexer<T>::_InstancesLock;
 
template<typename T>
std::map<size_t, T*> StaticIndexer<T>::_Instances;
 
template<typename T>
size_t StaticIndexer<T>::_InstancesIndex = 0;
 
 
 
class IRM_DLL_EXPORT PhreeqcRM : public StaticIndexer<PhreeqcRM>
{
public:
    static void             CleanupReactionModuleInstances(void);
    static int              CreateReactionModule(int nxyz, MP_TYPE nthreads);
    static IRM_RESULT       DestroyReactionModule(int n);
private:
    virtual void AddOutputVars(std::string option, std::string def);
    virtual void ClearBMISelectedOutput();
    virtual void GenerateAutoOutputVars();
    virtual void UpdateBMI(RMVARS v_enum);
public:
#ifdef USE_YAML
    static int GetGridCellCountYAML(const char* YAML_file);
#endif // #ifdef USE_YAML
 
    PhreeqcRM(int nxyz, MP_TYPE thread_count_or_communicator, PHRQ_io * io=NULL, bool delay_construct=false);
    virtual ~PhreeqcRM(void);
    IRM_RESULT                                CloseFiles(void);
    IPhreeqc * Concentrations2Utility(const std::vector< double > &c,
           const std::vector< double > &tc, const std::vector< double > &p_atm);
    IRM_RESULT                                CreateMapping(const std::vector< int > &grid2chem);
    void                                      DecodeError(int result);
    IRM_RESULT                                DumpModule(bool dump_on, bool append = false);
    void                                      ErrorHandler(int result, const std::string &e_string);
    void                                      ErrorMessage(const std::string &error_string, bool prepend = true);
    int                                       FindComponents();
    const std::vector < std::vector <int> > & GetBackwardMapping(void) {return this->backward_mapping;}
    void GetBackwardMappingSWIG(std::vector<int>& nback_output, std::vector<int>& cellnumbers_output);
    int                                       GetChemistryCellCount(void) const {return this->count_chemistry;}
    int                                       GetComponentCount(void) const {return (int) this->components.size();}
const std::vector< std::string > &          GetComponents(void) const {return this->components;}
    
    IRM_RESULT                                GetConcentrations(std::vector< double > &c_output);
    int                                       GetCurrentSelectedOutputUserNumber(void);
    std::string                               GetDatabaseFileName(void) {return this->database_file_name;}
    IRM_RESULT                                GetDensityCalculated(std::vector< double > & d_output);
    IRM_RESULT                                GetDensity(std::vector< double >& d_output);
    const std::vector < int> &                GetEndCell(void) {return this->end_cell;}
const std::vector< std::string > &          GetEquilibriumPhases(void) const { return this->EquilibriumPhasesList; }
int                                       GetEquilibriumPhasesCount(void) const { return (int) this->EquilibriumPhasesList.size(); }
 
 
    int                                       GetErrorHandlerMode(void) {return this->error_handler_mode;}
    std::string                                  GetErrorString(void);
 
const std::vector< std::string > &          GetExchangeNames(void) const { return this->ExchangeNamesList; }
const std::vector< std::string > &          GetExchangeSpecies(void) const { return this->ExchangeSpeciesNamesList; }
int                                       GetExchangeSpeciesCount(void) const { return (int) this->ExchangeSpeciesNamesList.size(); }
 
 
    std::string                               GetFilePrefix(void) {return this->file_prefix;}
const std::vector < int > &               GetForwardMapping(void) {return this->forward_mapping_root;}
 
const std::vector< std::string > &          GetGasComponents(void) const { return this->GasComponentsList; }
int                                       GetGasComponentsCount(void) const { return (int) this->GasComponentsList.size(); }
 
IRM_RESULT                                GetGasCompMoles(std::vector< double >& gas_moles_output);
 
IRM_RESULT                                GetGasCompPhi(std::vector< double >& gas_phi);
 
IRM_RESULT                                GetGasCompPressures(std::vector< double >& gas_pressure);
 
IRM_RESULT                                GetGasPhaseVolume(std::vector< double >& gas_volume_output);
 
    const std::vector < double > &            GetGfw(void) {return this->gfw;}
    int                                       GetGridCellCount(void) {return this->nxyz;}
IPhreeqc *                                GetIPhreeqcPointer(int i);
IRM_RESULT                                GetIthConcentration(int i, std::vector< double >& c_output);
IRM_RESULT                                GetIthSpeciesConcentration(int i, std::vector< double >& c_output);
 
const std::vector< std::string > &          GetKineticReactions(void) const { return this->KineticReactionsList; }
int                                       GetKineticReactionsCount(void) const { return (int) this->KineticReactionsList.size(); }
 
    int                                 GetMpiMyself(void) const {return this->mpi_myself;}
    int                                 GetMpiTasks(void) const {return this->mpi_tasks;}
    int                                       GetNthSelectedOutputUserNumber(int n);
    bool                                GetPartitionUZSolids(void) const {return this->partition_uz_solids;}
#ifdef USE_RV
    std::vector< double > &                     GetPoreVolume(void) {return this->pore_volume;}
#endif
    const std::vector< double >& GetPorosity(void);
    const std::vector< double > &                     GetPressure(void);
    const std::vector< int > &                  GetPrintChemistryMask (void) {return this->print_chem_mask_root;}
    const std::vector <bool> &                GetPrintChemistryOn(void) {return this->print_chemistry_on;}
    bool                                      GetRebalanceByCell(void) const {return this->rebalance_by_cell;}
    double                                    GetRebalanceFraction(void) const {return this->rebalance_fraction;}
IRM_RESULT               GetSaturationCalculated(std::vector< double > & sat_output);
IRM_RESULT               GetSaturation(std::vector< double >& sat_output);
    IRM_RESULT                                GetSelectedOutput(std::vector< double > &s_output);
    int                                       GetSelectedOutputColumnCount(void);
    int                                       GetSelectedOutputCount(void);
    IRM_RESULT                  GetSelectedOutputHeading(int icol, std::string &heading);
    IRM_RESULT                  GetSelectedOutputHeadings(std::vector< std::string >& headings);
    bool                                      GetSelectedOutputOn(void) {return this->selected_output_on;}
int                                       GetSelectedOutputRowCount(void);
 
int                                       GetSICount(void) const { return (int) this->SINamesList.size(); }
const std::vector< std::string > &          GetSINames(void) const { return this->SINamesList; }
 
const std::vector< std::string > &          GetSolidSolutionComponents(void) const { return this->SolidSolutionComponentsList; }
int                                       GetSolidSolutionComponentsCount(void) const { return (int) this->SolidSolutionComponentsList.size(); }
 
const std::vector< std::string > &          GetSolidSolutionNames(void) const { return this->SolidSolutionNamesList; }
 
    const std::vector< double > &               GetSolutionVolume(void);
    IRM_RESULT                                GetSpeciesConcentrations(std::vector< double > & species_conc_output);
    int                                       GetSpeciesCount(void) {return (int) this->species_names.size();}
    const std::vector< double > &               GetSpeciesD25(void) {return this->species_d_25;}
    IRM_RESULT                                GetSpeciesLog10Gammas(std::vector< double > & species_log10gammas);   
 
IRM_RESULT                                GetSpeciesLog10Molalities(std::vector< double >& species_log10molalities);    
 
    const std::vector< std::string > &          GetSpeciesNames(void) {return this->species_names;}
    bool                                      GetSpeciesSaveOn(void) {return this->species_save_on;}
 
    const std::vector<cxxNameDouble> &        GetSpeciesStoichiometry(void) {return this->species_stoichiometry;}
    void                                      GetSpeciesStoichiometrySWIG(std::vector<std::string>& species_output,
        std::vector<int>& nelt_output, std::vector<std::string>& elts_output, std::vector<double>& coef_output);
    const std::vector< double > &               GetSpeciesZ(void) {return this->species_z;}
const std::vector < int> &                GetStartCell(void) {return this->start_cell;}
 
const std::vector< std::string > &          GetSurfaceNames(void) const { return this->SurfaceNamesList; }
 
const std::vector< std::string > &          GetSurfaceSpecies(void) const { return this->SurfaceSpeciesNamesList; }
 
int                                       GetSurfaceSpeciesCount(void) const { return (int) this->SurfaceSpeciesNamesList.size(); }
 
 
const std::vector< std::string > &          GetSurfaceTypes(void) const { return this->SurfaceTypesList; }
 
    //const std::vector< double > &               GetTemperature(void) {return this->tempc;}
    const std::vector< double > &               GetTemperature(void);
    //int PhreeqcRM::GetVarItemsize(const std::string name);
    //int PhreeqcRM::GetVarNbytes(const std::string name);
 
    int                                       GetThreadCount() {return this->nthreads;}
    double                                    GetTime(void) const {return this->time;}
    double                                    GetTimeConversion(void) {return this->time_conversion;}
    double                                    GetTimeStep(void) {return this->time_step;}
    int                                       GetUnitsExchange(void) {return this->units_Exchange;}
    int                                       GetUnitsGasPhase(void) {return this->units_GasPhase;}
    int                                       GetUnitsKinetics(void) {return this->units_Kinetics;}
    int                                       GetUnitsPPassemblage(void) {return this->units_PPassemblage;}
    int                                       GetUnitsSolution(void) {return this->units_Solution;}
    int                                       GetUnitsSSassemblage(void) {return this->units_SSassemblage;}
    int                                       GetUnitsSurface(void) {return this->units_Surface;}
    const std::vector<double>& GetViscosity();
 
    
    const std::vector<IPhreeqcPhast *> &      GetWorkers() {return this->workers;}
#ifdef USE_YAML
IRM_RESULT      InitializeYAML(std::string yamlfile);
#endif
IRM_RESULT InitialPhreeqc2Concentrations(std::vector < double > & destination_c, 
    const std::vector < int >    & boundary_solution1);
    IRM_RESULT                                InitialPhreeqc2Concentrations(
                                                    std::vector < double > & destination_c,
                                                    const std::vector < int >    & boundary_solution1,
                                                    const std::vector < int >    & boundary_solution2,
                                                    const std::vector < double > & fraction1);
IRM_RESULT  InitialEquilibriumPhases2Module(const std::vector < int >& equilibrium_phases);
IRM_RESULT  InitialExchanges2Module(const std::vector < int >& exchanges);
IRM_RESULT  InitialGasPhases2Module(const std::vector < int >& gas_phases);
IRM_RESULT  InitialKinetics2Module(const std::vector < int >& kinetics);
IRM_RESULT  InitialPhreeqc2Module(const std::vector < int >& initial_conditions1);
IRM_RESULT InitialPhreeqc2Module(
    const std::vector < int >& initial_conditions1,
    const std::vector < int >& initial_conditions2,
    const std::vector < double >& fraction1);
IRM_RESULT InitialPhreeqc2SpeciesConcentrations(std::vector < double >& destination_c,
    const std::vector < int >& boundary_solution1);
IRM_RESULT InitialPhreeqc2SpeciesConcentrations(std::vector < double >& destination_c,
    const std::vector < int >& boundary_solution1,
    const std::vector < int >& boundary_solution2,
    const std::vector < double >& fraction1);
IRM_RESULT                                InitialPhreeqcCell2Module(int n,
    const std::vector< int >& cell_numbers);
 
IRM_RESULT  InitialSolidSolutions2Module(const std::vector < int >& solid_solutions);
IRM_RESULT  InitialSolutions2Module(const std::vector < int >& solutions);
IRM_RESULT  InitialSurfaces2Module(const std::vector < int >& surfaces);
 
    virtual IRM_RESULT                          LoadDatabase(const std::string &database);
    void                                      LogMessage(const std::string &str);
    int                                       MpiAbort();
    IRM_RESULT                                MpiWorker();
    IRM_RESULT                                MpiWorkerBreak();
    IRM_RESULT                                OpenFiles(void);
    void                                      OutputMessage(const std::string &str);
    IRM_RESULT                                ReturnHandler(IRM_RESULT result, const std::string& e_string);
 
    
    IRM_RESULT                                RunCells(void);
    IRM_RESULT                                RunFile(bool workers, bool initial_phreeqc, bool utility,  const std::string & chemistry_name);
    IRM_RESULT                                RunString(bool workers, bool initial_phreeqc, bool utility, const std::string & input_string);
    void                                      ScreenMessage(const std::string &str);
    IRM_RESULT                                SetComponentH2O(bool tf);
    IRM_RESULT                                SetConcentrations(const std::vector< double > &c);
    IRM_RESULT                                SetCurrentSelectedOutputUserNumber(int n_user);
    IRM_RESULT                                SetDensityUser(const std::vector< double > &density);
    IRM_RESULT                                SetDensity(const std::vector< double >& density);
    IRM_RESULT                                SetDumpFileName(const std::string & dump_name);
    IRM_RESULT                                SetErrorHandlerMode(int mode);
    IRM_RESULT                                SetErrorOn(bool tf);
    IRM_RESULT                                SetFilePrefix(const std::string & prefix);
 
IRM_RESULT                                SetGasCompMoles(const std::vector< double >& gas_moles);
IRM_RESULT                                SetGasPhaseVolume(const std::vector< double >& gas_volume);
IRM_RESULT                                SetIthConcentration(int i, std::vector< double >& c);
IRM_RESULT                                SetIthSpeciesConcentration(int i, std::vector< double >& c);
 
    IRM_RESULT                                SetMpiWorkerCallbackC(int (*fcn)(int *method, void * cookie));
    IRM_RESULT                                SetMpiWorkerCallbackCookie(void * cookie);
    IRM_RESULT                                SetMpiWorkerCallbackFortran(int (*fcn)(int *method));
    IRM_RESULT                                SetNthSelectedOutput(int n);
    IRM_RESULT                                SetPartitionUZSolids(bool tf);
    IRM_RESULT                                SetPorosity(const std::vector< double > &por);
 
    IRM_RESULT                                SetPressure(const std::vector< double > &p);
    IRM_RESULT                                SetPrintChemistryMask(const std::vector<int> & cell_mask);
    IRM_RESULT                                SetPrintChemistryOn(bool workers, bool initial_phreeqc, bool utility);
    IRM_RESULT                                SetRebalanceByCell(bool tf);
    IRM_RESULT                                SetRebalanceFraction(double f);
    IRM_RESULT                                SetRepresentativeVolume(const std::vector< double > &rv);
    IRM_RESULT                                SetSaturationUser(const std::vector< double > &sat);
    IRM_RESULT                                SetSaturation(const std::vector< double >& sat);
    IRM_RESULT                                SetScreenOn(bool tf);
    IRM_RESULT                                SetSelectedOutputOn(bool tf);
    IRM_RESULT                                SetSpeciesSaveOn(bool save_on);
    IRM_RESULT                                SetTemperature(const std::vector< double > &t);
    IRM_RESULT                                SetTime(double time);
    IRM_RESULT                                SetTimeConversion(double conv_factor);
    IRM_RESULT                                SetTimeStep(double time_step);
    IRM_RESULT                                SetUnitsExchange(int option);
    IRM_RESULT                                SetUnitsGasPhase(int option);
    IRM_RESULT                                SetUnitsKinetics(int option);
    IRM_RESULT                                SetUnitsPPassemblage(int option);
    IRM_RESULT                                SetUnitsSolution(int option);
    IRM_RESULT                                SetUnitsSSassemblage(int option);
    IRM_RESULT                                SetUnitsSurface(int option);
IRM_RESULT                                SpeciesConcentrations2Module(const std::vector< double > & species_conc);
 
IRM_RESULT StateSave(int istate);
IRM_RESULT StateApply(int istate);
IRM_RESULT StateDelete(int istate);
void UseSolutionDensityVolume(bool tf);
    void                                      WarningMessage(const std::string &warnstr);
 
#if defined(SWIG) || defined(swig_python_EXPORTS)
// Python helpers
public:
    // BasicCallback handling.
    
    void                                      set_basic_callback(PyObject* py_callable, PyObject* py_cookie = nullptr);
 
    double                                    _execute_basic_callback(double val1, double val2, const char* message);
 
    struct PythonBasicCallbackData {
        PyObject* py_callable        = nullptr;
        PyObject* py_callback_cookie = nullptr;
    };
 
private:
    PythonBasicCallbackData                   python_basic_callback_data;
#endif
 
#if defined(SWIG) || defined(swig_python_EXPORTS)
// Python helpers
public:
    // MPI Callback handling.
 
    void                                      set_mpi_worker_callback(PyObject* py_callable, PyObject* py_cookie = nullptr);
 
    int                                       _execute_mpi_worker_callback(int val);
 
    struct PythonMpiWorkerCallbackData {
        PyObject* py_callable        = nullptr;
        PyObject* py_callback_cookie = nullptr;
    };
 
private:
    PythonMpiWorkerCallbackData               python_mpi_worker_callback_data;
#endif
 
public:
// Utilities
    static std::string                        Char2TrimString(const char * str, size_t l = 0);
    static bool                               FileExists(const std::string &name);
    static void                               FileRename(const std::string &temp_name, const std::string &name,
                                                   const std::string &backup_name);
    static IRM_RESULT                         Int2IrmResult(int r, bool positive_ok);
private:
    IRM_RESULT                                CellInitialize(
                                                  int i,
                                                  int n_user_new,
                                                  const int *initial_conditions1,
                                                  const int *initial_conditions2,
                                                  double *fraction1,
                                                  std::set<std::string> &error_set);
    IRM_RESULT                                CheckCells();
    int                                       CheckSelectedOutput();
    //void                                      Collapse2Nchem(double *d_in, double *d_out);
    //void                                      Collapse2Nchem(int *i_in, int *i_out);
    IPhreeqc *                                Concentrations2UtilityH2O(const std::vector< double > &c_in,
                                                   const std::vector< double > &t_in, const std::vector< double > &p_in);
    IPhreeqc *                                Concentrations2UtilityNoH2O(const std::vector< double > &c_in,
                                                   const std::vector< double > &t_in, const std::vector< double > &p_in);
    void                                      Concentrations2Solutions(int n, std::vector< double > &c);
    void                                      Concentrations2SolutionsH2O(int n, std::vector< double > &c);
    void                                      Concentrations2SolutionsNoH2O(int n, std::vector< double > &c);
    void                                      cxxSolution2concentration(cxxSolution * cxxsoln_ptr, std::vector< double > & d, double v, double dens);
    void                                      cxxSolution2concentrationH2O(cxxSolution * cxxsoln_ptr, std::vector< double > & d, double v, double dens);
    void                                      cxxSolution2concentrationNoH2O(cxxSolution * cxxsoln_ptr, std::vector< double > & d, double v, double dens);
    void                                      GatherNchem(std::vector< double > &source, std::vector< double > &destination);
    cxxStorageBin &                           Get_phreeqc_bin(void) {return *this->phreeqc_bin;}
    IRM_RESULT                                HandleErrorsInternal(std::vector< int > & r);
    void                                      PartitionUZ(int n, int iphrq, int ihst, double new_frac);
    void                                      RebalanceLoad(void);
    void                                      RebalanceLoadPerCell(void);
    IRM_RESULT                                RunCellsThread(int i);
    IRM_RESULT                                RunFileThread(int n);
    IRM_RESULT                                RunStringThread(int n, std::string & input);
    IRM_RESULT                                RunCellsThreadNoPrint(int n);
    void                                      Scale_solids(int n, int iphrq, double frac);
    void                                      ScatterNchem(double *d_array);
    void                                      ScatterNchem(int *i_array);
    void                                      ScatterNchem(std::vector< double > &source, std::vector< double > &destination);
    void                                      ScatterNchem(std::vector< int > &source, std::vector< int > &destination);
    IRM_RESULT                                SetChemistryFileName(const char * prefix = NULL);
    IRM_RESULT                                SetDatabaseFileName(const char * db = NULL);
    void                                      SetEndCells(void);
    void                                      SetEndCellsHeterogeneous(void);
    double                                    TimeStandardTask(void);
    IRM_RESULT                                TransferCells(cxxStorageBin &t_bin, int old, int nnew);
    IRM_RESULT                                TransferCellsUZ(std::ostringstream &raw_stream, int old, int nnew);
 
private:
    //IRM_RESULT                                SetGeneric(std::vector< double > &destination, int newSize, const std::vector< double > &origin, int mpiMethod, const std::string &name, const double newValue = 0.0);
    IRM_RESULT                                SetGeneric(const std::vector< double > &source, std::vector< double > &destination_root, std::vector< double > &destination_worker, int mpiMethod, const std::string &name);
public:
 
protected:
 
#if defined(_MSC_VER)
/* disable warning C4251: 'identifier' : class 'type' needs to have dll-interface to be used by clients of class 'type2' */
#pragma warning(disable:4251)
#endif
    bool component_h2o;                      // true: use H2O, excess H, excess O, and charge;
                                             // false total H, total O, and charge
    std::string database_file_name;
    std::string chemistry_file_name;
    std::string dump_file_name;
    std::string file_prefix;
    cxxStorageBin * phreeqc_bin;
    int mpi_myself;
    int mpi_tasks;
    std::vector <std::string> components;   // list of components to be transported
    std::vector <double> gfw;               // gram formula weights converting mass to moles (1 for each component)
    double gfw_water;                       // gfw of water
    bool partition_uz_solids;
    int nxyz;                               // number of nodes
    int count_chemistry;                    // number of cells for chemistry
    double time;                            // time from transport, sec
    double time_step;                       // time step from transport, sec
    double time_conversion;                 // time conversion factor, multiply to convert to preferred time unit for output
    std::vector <double> CurrentConcentrations;
    std::vector <double> CurrentSpeciesConcentrations;
    std::vector <double> IthCurrentConcentrations;
    std::vector <double> IthCurrentSpeciesConcentrations;
    std::set<int>        IthConcentrationSet;
    std::set<int>        IthSpeciesConcentrationSet;
    std::vector <double> old_saturation_root;   // saturation fraction from previous step
    std::vector <double> old_saturation_worker;
    std::vector< double > saturation_root;      // nxyz saturation fraction
    std::vector< double > saturation_worker;        // nchem on workers saturation fraction
    std::vector< double > pressure_root;            // nxyz on root current pressure
    std::vector< double > pressure_worker;      // nchem on workers current pressure
    std::vector< double > rv_root;              // nxyz on root representative volume
    std::vector< double > rv_worker;                // nchem on workers representative volume
    std::vector< double > porosity_root;            // nxyz porosity
    std::vector< double > porosity_worker;      // nchem on workers porosity
    std::vector< double > tempc_root;             // nxyz on root temperature Celsius 
    std::vector< double > tempc_worker;         // nchem on workers temperature Celsius 
    std::vector< double > density_root;         // nxyz density
    std::vector< double > density_worker;           // nchem on workers density
    std::vector< double > viscosity_root;           // nxyz viscosity
    std::vector< double > viscosity_worker;         // nchem on workers viscosity
    std::vector< double > solution_volume_root;   // nxyz on root solution volume
    std::vector< double > solution_volume_worker;   // nchem on workers solution_volume 
    std::vector< int > print_chem_mask_root;        // nxyz print flags for output file
    std::vector< int > print_chem_mask_worker;  // nchem print flags for output file
    bool rebalance_by_cell;                 // rebalance method 0 std, 1 by_cell
    double rebalance_fraction;              // parameter for rebalancing process load for parallel
    int units_Solution;                     // 1 mg/L, 2 mol/L, 3 kg/kgs
    int units_PPassemblage;                 // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
    int units_Exchange;                     // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
    int units_Surface;                      // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
    int units_GasPhase;                     // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
    int units_SSassemblage;                 // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
    int units_Kinetics;                     // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
    std::vector <int> forward_mapping_root;                 // mapping from nxyz cells to count_chem chemistry cells
    std::vector <std::vector <int> > backward_mapping;      // mapping from count_chem chemistry cells to nxyz cells
    bool use_solution_density_volume;
    bool worker_waiting;
    // print flags
    std::vector< bool > print_chemistry_on; // print flag for chemistry output file
    bool selected_output_on;                // create selected output
 
    int error_count;
    int error_handler_mode;                 // 0, return code; 1, throw; 2 exit;
    bool need_error_check;
    std::string phreeqcrm_error_string;
 
    std::map<std::string, int> method_map;
    // threading
    int nthreads;
    std::vector<IPhreeqcPhast *> workers;
    //std::vector< int > start_cell;
    //std::vector< int > end_cell;
    PHRQ_io *phreeqcrm_io;
    bool delete_phreeqcrm_io;
 
    // mpi
#ifdef USE_MPI
    MPI_Comm phreeqcrm_comm;                                       // MPI communicator
#endif
    int (*mpi_worker_callback_fortran) (int *method);
    int (*mpi_worker_callback_c) (int *method, void *cookie);
    void *mpi_worker_callback_cookie;
 
    // mcd
    bool species_save_on;
    std::vector <std::string> species_names;
    std::vector <double> species_z;
    std::vector <double> species_d_25;
    std::vector <cxxNameDouble> species_stoichiometry;
    std::map<int, int> s_num2rm_species_num;
    std::vector< double > standard_task_vector;   // root only
 
    std::vector< int > start_cell;
    std::vector< int > end_cell;
    // reactant lists
    std::vector <std::string> ExchangeSpeciesNamesList;
    std::vector <std::string> ExchangeNamesList;
    std::vector <std::string> SurfaceSpeciesNamesList;
    std::vector <std::string> SurfaceTypesList;
    std::vector <std::string> SurfaceNamesList;
 
    std::vector <std::string> EquilibriumPhasesList;
    std::vector <std::string> GasComponentsList;
    std::vector <std::string> KineticReactionsList;
    std::vector <std::string> SolidSolutionComponentsList;
    std::vector <std::string> SolidSolutionNamesList;
    std::vector <std::string> SINamesList;
    std::set <std::string> ElementRedoxSet;
 
protected:
    static const int default_nxyz = 10;
    static const MP_TYPE default_data_for_parallel_processing;
 
    class Initializer;  // Forward declaration
    std::unique_ptr<Initializer> initializer;
 
    virtual void Construct();
 
    void set_data_for_parallel_processing(MP_TYPE value);
    void set_io(PHRQ_io* value);
    void set_nxyz(int value);
 
private:
    //friend class RM_interface;
    friend class BMIPhreeqcRM;
    friend class VarManager;
 
#if defined(_MSC_VER)
/* reset warning C4251 */
#pragma warning(default:4251)
#endif
 
};
#endif // !defined(PHREEQCRM_H_INCLUDED)