PROBLEM:

I've come across issues with

1) the MELTS engine remembering properties (e.g., mass, g, h) of phases that are no longer stable (the list of phases are still correct) from previous calculations on the same Node of a path.

2) calculations on a newly defined path affecting the properties of phases on the previous path.

Below is an example that walks through both of these issues in a simple example.

POSSIBLE WORKAROUNDS:

I have come up with a couple of work arounds for these problems in case other people encounter them.

1) The first issue means that you cannot rely on getListProperty calls for finding the properties of individual phases. To get the property of a phase during a run you have to loop through each Node of a path and extract the properties of the phase only if it is in the solidNames or liquidNames lists. Similarly, when doing calls for the properties of phases only do so for phases in the solidNames or liquidNames lists.

2) The second issue can be solved by defining new paths in a somewhat round about route e.g.,

`path=MELTSdynamic(1);`

path2=path.copy;

The copy function seems to reset all the values and allows for starting a fresh path.

EXAMPLE SCRIPT

`%% SJL 2/21`

%Script to reproduce issues with inherited phase values from previous

%calculations and paths

%start off with a clean slate and turn off warning

clear all

close all

fclose('all');

if isunix

addpath('/Users/simonlock/src/MELTS_Matlab/package')

addpath('/Users/simonlock/Documents/alphaMELTS_for_Matlab_functions')

else

addpath(strcat(test_script_dir,'\..'))

end

warning('off', 'MATLAB:loadlibrary:cppoutput');

%% Example of the engine retaining properties from previous calculations

disp('start a path')

path1=MELTSdynamic(1);

disp('The first path has the properties you would expect. Empty solid names:')

disp(path1.engine.solidNames)

disp('and empty e.g., mass MELTSmap')

disp(path1.engine.mass)

disp(' ')

disp('Now equilibrate this path at a T-p where olivine is the only stable solid phases')

bulk = zeros(19, 1);

bulk([1:9 11:16]) = [45.0 0.201 4.45 0.0 0.384 8.05 0.0 37.8 0.0 3.55 0.36 0.0 0.0 0.0 0.0];

bulk=bulk*100/sum(bulk);

path1.engine.setBulkComposition(bulk)

path1.engine.temperature=1700;

path1.engine.pressure=1000;

path1.engine.calcEquilibriumState(1, 1);

disp('This time we have a solid names')

disp(path1.engine.solidNames)

disp('and mass map')

disp(path1.engine.mass)

disp('with the expected sizes and values')

disp(' ')

disp('Now equilibrate the system at a lower temperature where spinel is also stable')

path1.engine.temperature=1350;

path1.engine.calcEquilibriumState(1, 1);

disp('Here again the answers are what you would expect. Solid names:')

disp(path1.engine.solidNames)

disp('mass map')

disp(path1.engine.mass)

disp('and you can ask, and get, mass values for each phase')

fprintf('Olivine: %f \n', path1.engine.mass('olivine1'))

fprintf('Spinel: %f \n',path1.engine.mass('spinel1'))

disp(' ')

disp('BUT if we now equilibrate at the original temperature were spinel goes away...')

path1.engine.temperature=1700;

path1.engine.calcEquilibriumState(1, 1);

disp('we still have values for spinel in the mass map, even though it is not in the solidNames list')

disp('Solid Names')

disp(path1.engine.solidNames)

disp('mass map')

disp(path1.engine.mass)

disp('and you can still get mass values from the mass map for both spinel and olivine even though spinel is not there')

fprintf('Olivine: %f \n', path1.engine.mass('olivine1'))

fprintf('Spinel: %f \n',path1.engine.mass('spinel1'))

disp(' and spinel turns up in getListProperty')

disp(path1.getListProperty('mass','spinel1'))

disp(' ')

disp('Note that the spinel is not contrbuting to the total mass recorded under "bulk" (or other properties as far as I have found)')

mass_tot=path1.getListProperty('mass','liquid1');

for i=1:length((path1.engine.solidNames))

mass_tot=mass_tot+path1.engine.mass(path1.engine.solidNames(i));

end

fprintf('Mass from sum of phases in solidNames: %f \n', mass_tot)

fprintf('Mass from mass("bulk"): %f \n', path1.engine.mass('bulk'))

fprintf('Including the mass from spinel: %f \n',mass_tot+path1.engine.mass('spinel1'))

%% Issues persist when initializing a new path

disp(' ')

disp('Now the issues with inheritance also continue if you start a new path')

disp('start a new path')

path2=MELTSdynamic(1);

%XXXXXXXXXXXXXXXXXXXX

%initializing this way can stop these issues

%path=MELTSdynamic(1);

%path2=path.copy;

%XXXXXXXXXXXXXXXXXXXX

disp('Even when just initialized the path has inherited the mass etc. properties from the previous path(s) you have been working on')

disp('Solid names:')

disp(path2.engine.solidNames)

disp('Mass map:')

disp(path2.engine.mass)

disp('Bulk mass')

disp(path2.engine.mass('bulk'))

disp('Along with that pesky spinel mass')

fprintf('Olivine: %f \n', path2.engine.mass('olivine1'))

fprintf('Spinel: %f \n',path2.engine.mass('spinel1'))

disp(' ')

disp('Now equilibrate this path at the high (but different) T where only olivine exists')

path2.engine.setBulkComposition(bulk)

path2.engine.temperature=1710;

path2.engine.pressure=1000;

path2.engine.calcEquilibriumState(1, 1);

disp('We get what we got before')

disp('Solid names:')

disp(path2.engine.solidNames)

disp('Mass map:')

disp(path2.engine.mass)

disp('the pesky spinel is still around in the mass map even though it has never been stable in this path')

fprintf('Olivine: %f \n', path2.engine.mass('olivine1'))

fprintf('Spinel: %f \n',path2.engine.mass('spinel1'))

disp('AND doing the calc on path2 updates the values of mass etc. for both path1 and path2!')

fprintf('Olivine path 1: %f \n', path1.engine.mass('olivine1'))

fprintf('Olivine path 2: %f \n', path2.engine.mass('olivine1'))

disp(' ')

disp('Adding another node solves these problems')

path2=path2.addNodeAfter;

disp('The maps for mass etc. are now empty, although the list of solid Names is inherited')

disp('Solid names:')

disp(path2.engine.solidNames)

disp('Mass map:')

disp(path2.engine.mass)

disp(' ')

disp('After equilibration (at a slightly different temperature) we gain the correct values for the calc')

path2.engine.temperature=1700;

path2.engine.pressure=1000;

path2.engine.calcEquilibriumState(1, 1);

disp('Solid names:')

disp(path2.engine.solidNames)

disp('Mass map:')

disp(path2.engine.mass)

disp('and path1 has not been altered by this calculation. e.g., mass')

fprintf('Olivine path 1: %f \n',path1.engine.mass('olivine1'))

fprintf('Olivine path 2: %f \n',path2.engine.mass('olivine1'))

disp(' ')

disp('Spinel turns up in the list properties for the first step where we inherited some values from the previous calcs, but not in the new Node')

path2.getListProperty('mass','spinel1')

disp('and asking for spinel at the current step leads to an error (as it should)')

disp(path2.engine.mass('spinel1'))