Incorrect output in IPhreeqc compared to Phreeqc windows

[IgorMarkelov]

Dear David:

I am using Python in OSX to calibrate and perform sensitivity analyses of phreeqc models. I performed already more than 30 sensitivity and calibration analysis using python and never had this issue before. The problem is the following:

I am running the same file in windows (in notepad with latest version of phreeqc) and os x (using python and latest version of iphreeqc) but receive an incorrect output in IPhreeqc compared to Phreeqc windows for the couple of species in my model. For instance, the output in windows for concentrations is the following vector:

m_Cr_six_aq: [5.6068e-04 4.7857e-04 4.4061e-04 4.0531e-04 3.7123e-04 3.3752e-04 3.0370e-04 2.6952e-04 2.3475e-04 1.9889e-04 1.6242e-04 1.2551e-04 8.8445e-05 5.2451e-05 1.9329e-05 5.3029e-13 5.3028e-13 5.3027e-13 5.3027e-13]

Meanwhile in OSX: m_Cr_six_aq: [5.6068e-04 0 0 0 0 0 0 0 0 0 0 ...]

This is strange behaviour only for couple of species, for others everything is fine. As you may see, the problem is somewhere in populating of the solution vector. Interestingly, for other species in the same run it works and the python method used is the same (see below).

These are methods which I am using in Python to run Phreeqc text string:

Code: [Select]

def make_selected_output(components):
    """
    Build SELECTED_OUTPUT data block
    """
    headings = "-headings    cb    H    O    "
    for i in range(len(components)):
        headings += components[i] + "\t"
    selected_output = """
    SELECTED_OUTPUT
        -reset false
    USER_PUNCH
    """
    selected_output += headings + "\n"
    #
    # charge balance, H, and O
    #
    code = '10 w = TOT("water")\n'
    code += '20 PUNCH CHARGE_BALANCE, TOTMOLE("H"), TOTMOLE("O")\n'
    #
    # All other elements
    #
    lino = 30
    for component in components:
        code += '%d PUNCH w*TOT(\"%s\")\n' % (lino, component)
        lino += 10
    selected_output += code
    return selected_output


def get_selected_output(phreeqc):
    """Return calculation result as dict.

    Header entries are the keys and the columns
    are the values as lists of numbers.
    """
    output = phreeqc.get_selected_output_array()
    print output
    header = output[0]
    conc = {}
    for head in header:
        conc[head] = []
    for row in output[1:]:
        for col, head in enumerate(header):
            conc[head].append(row[col])
    return conc


def simulate_with(string):
    phreeqc = phreeqc_mod.IPhreeqc()
    phreeqc.load_database(r"../../phreeqc/ror.dat")
    phreeqc.run_string(string)
    conc = get_selected_output(phreeqc)
    return {key: np.array(val) for (key, val) in conc.iteritems()}


I've used this combination for more than 30 different examples and it always worked before. I've attached the phreeqc file, maybe the problem is there.

P.S.:

By the way, would you be interested in creation of a cross platform Python library for calibration and sensitivity analysis? My code is very easy to use and can perform the following sens analyses:

Sobol Sensitivity Analysis (Sobol 2001, Saltelli 2002, Saltelli et al. 2010)
Method of Morris, including groups and optimal trajectories (Morris 1991, Campolongo et al. 2007)
Fourier Amplitude Sensitivity Test (FAST) (Cukier et al. 1973, Saltelli et al. 1999)
Delta Moment-Independent Measure (Borgonovo 2007, Plischke et al. 2013)
Derivative-based Global Sensitivity Measure (DGSM) (Sobol and Kucherenko 2009)
Fractional Factorial Sensitivity Analysis (Saltelli et al. 2008)

Plus:

Calibration via various optimisation methods in Scipy.

If you are interested, I could create the open source Github repository to show you my code, and invite collaborators to improve one.

Thank you in advance.

Regards,
Igor

[dlparkhurst]

This is way complicated. I really can only debug PHREEQC (given input and database file), which I think you say is working correctly. If you do the IPhreeqc equivalent of SetOutputFileOn(true), by default there should be a phreeqc.out file. Modify your Python and PHREEQC input files so that you have the minimum to create the error. If you compare Python phreeqc.out file to the Windows PHREEQC output output, you may be able to tell what is different.

If the output is identical, then I suppose it is the definition of SELECTED_OUTPUT/USER_PUNCH, or the GetSelectedOutputValue method in IPhreeqc or its equivalent that Mike Mueller wrote for Python.

[IgorMarkelov]

Thank you very much for your reply. I will.

[dlparkhurst]

In addition to the output file, also write the selected output file to see if it has the correct data or not.