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
class IPhreeqcPhast;

class cxxStorageBin;
//class cxxNameDouble;
#include "NameDouble.h"
class cxxSolution;
class PHRQ_io;
#include <vector>
#include <list>
#include <set>
#include <map>
#include <string>

#if defined(_WINDLL)
#define IRM_DLL_EXPORT __declspec(dllexport)
#else
#define IRM_DLL_EXPORT
#endif

class PHRQ_io;
class IPhreeqc;
class IRM_DLL_EXPORT PhreeqcRMStop : public std::exception
{
public:
  const char *what() const throw () {return "Failure in PhreeqcRM\n";}
};

#include "IrmResult.h"
enum {
    METHOD_CREATEMAPPING,
    METHOD_DUMPMODULE,
    METHOD_FINDCOMPONENTS,
    METHOD_GETCONCENTRATIONS,
    METHOD_GETDENSITY,
    METHOD_GETERRORSTRING,
    METHOD_GETPRESSURE,
    METHOD_GETSATURATION,
    METHOD_GETSELECTEDOUTPUT,
    METHOD_GETSOLUTIONVOLUME,
    METHOD_GETSPECIESCONCENTRATIONS,
    METHOD_GETTEMPERATURE,
    METHOD_INITIALPHREEQC2MODULE,
    METHOD_INITIALPHREEQCCELL2MODULE,
    METHOD_LOADDATABASE,
    METHOD_MPIWORKERBREAK,
    METHOD_RUNCELLS,
    METHOD_RUNFILE,
    METHOD_RUNSTRING,
    METHOD_SETCOMPONENTH2O,
    METHOD_SETCONCENTRATIONS,
    METHOD_SETDENSITY,
    METHOD_SETERRORHANDLERMODE,
    METHOD_SETFILEPREFIX,
    METHOD_SETPARTITIONUZSOLIDS,
    METHOD_SETPOROSITY,
    METHOD_SETPRESSURE,
    METHOD_SETPRINTCHEMISTRYON,
    METHOD_SETPRINTCHEMISTRYMASK,
    METHOD_SETREBALANCEBYCELL,
    METHOD_SETREBALANCEFRACTION,
    METHOD_SETREPRESENTATIVEVOLUME,
    METHOD_SETSATURATION,
    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_USESOLUTIONDENSITYVOLUME
} /* MPI_METHOD */;

class IRM_DLL_EXPORT PhreeqcRM
{
public:
    static void             CleanupReactionModuleInstances(void);
    static int              CreateReactionModule(int nxyz, MP_TYPE nthreads);
    static IRM_RESULT       DestroyReactionModule(int n);
    static PhreeqcRM      * GetInstance(int n);

    PhreeqcRM(int nxyz, MP_TYPE thread_count_or_communicator, PHRQ_io * io=NULL);
    ~PhreeqcRM(void);
    IRM_RESULT                                CloseFiles(void);
    IPhreeqc * Concentrations2Utility(std::vector<double> &c,
           std::vector<double> tc, std::vector<double> p_atm);
    IRM_RESULT                                CreateMapping(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;}
    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);
    std::string                               GetDatabaseFileName(void) {return this->database_file_name;}
    IRM_RESULT                                GetDensity(std::vector<double> & density);
    const std::vector < int> &                GetEndCell(void) const {return this->end_cell;}
    int                                       GetErrorHandlerMode(void) {return this->error_handler_mode;}
    std::string                                  GetErrorString(void);
    std::string                               GetFilePrefix(void) {return this->file_prefix;}
    const std::vector < int > &               GetForwardMapping(void) {return this->forward_mapping_root;}
    const std::vector < double > &            GetGfw(void) {return this->gfw;}
    int                                       GetGridCellCount(void) {return this->nxyz;}
    IPhreeqc *                                GetIPhreeqcPointer(int i);
    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> &                     GetPressure(void);
    const std::vector<int> &                  GetPrintChemistryMask (void) {return this->print_chem_mask_root;}
    const std::vector <bool> &                GetPrintChemistryOn(void) const {return this->print_chemistry_on;}
    bool                                      GetRebalanceByCell(void) const {return this->rebalance_by_cell;}
    double                                    GetRebalanceFraction(void) const {return this->rebalance_fraction;}
IRM_RESULT               GetSaturation(std::vector<double> & sat);
    IRM_RESULT                                GetSelectedOutput(std::vector<double> &so);
    int                                       GetSelectedOutputColumnCount(void);
    int                                       GetSelectedOutputCount(void);
    IRM_RESULT                                GetSelectedOutputHeading(int icol, std::string &heading);
    bool                                      GetSelectedOutputOn(void) {return this->selected_output_on;}
    int                                       GetSelectedOutputRowCount(void);
    const std::vector<double> &               GetSolutionVolume(void);
    IRM_RESULT                                GetSpeciesConcentrations(std::vector<double> & species_conc);
    int                                       GetSpeciesCount(void) {return (int) this->species_names.size();}
    const std::vector<double> &               GetSpeciesD25(void) {return this->species_d_25;}
    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;}
    const std::vector<double> &               GetSpeciesZ(void) {return this->species_z;}
    const std::vector < int> &                GetStartCell(void) const {return this->start_cell;}
    //const std::vector<double> &               GetTemperature(void) {return this->tempc;}
    const std::vector<double> &               GetTemperature(void);
    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<IPhreeqcPhast *> &      GetWorkers() {return this->workers;}
    IRM_RESULT                                InitialPhreeqc2Concentrations(
                                                    std::vector < double > & destination_c,
                                                    std::vector < int >    & boundary_solution1);
    IRM_RESULT                                InitialPhreeqc2Concentrations(
                                                    std::vector < double > & destination_c,
                                                    std::vector < int >    & boundary_solution1,
                                                    std::vector < int >    & boundary_solution2,
                                                    std::vector < double > & fraction1);
    IRM_RESULT  InitialPhreeqc2Module(std::vector < int >    & initial_conditions1);
    IRM_RESULT InitialPhreeqc2Module(
        std::vector < int >    & initial_conditions1,
        std::vector < int >    & initial_conditions2,
        std::vector < double > & fraction1);
    IRM_RESULT                                InitialPhreeqc2SpeciesConcentrations(
                                                    std::vector < double > & destination_c,
                                                    std::vector < int >    & boundary_solution1);
    IRM_RESULT                                InitialPhreeqc2SpeciesConcentrations(
                                                    std::vector < double > & destination_c,
                                                    std::vector < int >    & boundary_solution1,
                                                    std::vector < int >    & boundary_solution2,
                                                    std::vector < double > & fraction1);
    IRM_RESULT                                InitialPhreeqcCell2Module(int n, const std::vector<int> &cell_numbers);
    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                                RunCells(void);
    IRM_RESULT                                ReturnHandler(IRM_RESULT result, const std::string &e_string);
    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                                SetDensity(const std::vector<double> &density);
    IRM_RESULT                                SetDumpFileName(const std::string & dump_name);
    IRM_RESULT                                SetErrorHandlerMode(int mode);
    IRM_RESULT                                SetFilePrefix(const std::string & prefix);
    IRM_RESULT                                SetMpiWorkerCallbackC(int (*fcn)(int *method, void * cookie));
    IRM_RESULT                                SetMpiWorkerCallbackCookie(void * cookie);
    IRM_RESULT                                SetMpiWorkerCallbackFortran(int (*fcn)(int *method));
    IRM_RESULT                                SetPartitionUZSolids(bool tf);
    IRM_RESULT                                SetPorosity(const std::vector<double> &por);

    IRM_RESULT                                SetPressure(const std::vector<double> &p);
    IRM_RESULT                                SetPrintChemistryMask(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                                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(std::vector<double> & species_conc);
void UseSolutionDensityVolume(bool tf);
    void                                      WarningMessage(const std::string &warnstr);

    // 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);
protected:
    IRM_RESULT                                CellInitialize(
                                                  int i,
                                                  int n_user_new,
                                                  int *initial_conditions1,
                                                  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(std::vector<double> &c_in,
                                                   std::vector<double> &t_in, std::vector<double> &p_in);
    IPhreeqc *                                Concentrations2UtilityNoH2O(std::vector<double> &c_in,
                                                   std::vector<double> &t_in, 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);
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> 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> 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;

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

    // 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

private:
    //friend class RM_interface;
    static std::map<size_t, PhreeqcRM*> Instances;
    static size_t InstancesIndex;

#if defined(_MSC_VER)
/* reset warning C4251 */
#pragma warning(default:4251)
#endif

};
#endif // !defined(PHREEQCRM_H_INCLUDED)