Matlab: Engine retaining and altering values from previous calcs and paths

[slock]

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.,

Code Select Expand

path=MELTSdynamic(1);
path2=path.copy;

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

EXAMPLE SCRIPT

Code Select Expand

%% 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'))



[Paula]

This is one of those "not so much a bug as a feature" situations. It's deliberate that you might want, say, to equilibrate the model and then do "Supplemental Calculator" type calculations for actual samples to compare. The allowed phase names for calcPhaseProperties etc. are a regular expression:

(phase name in MELTS) + (optional digits, as output by calcEquilibriumState, but could also be user defined) + (optional '_' followed by alphanumeric string e.g. sample name).

The code doesn't have a built in way to tell if names like "spinel1" are the result of a previous equilibration calculation or a user-defined name for an instance of that phase. It would be possible to add that check (and I will), but really the plan was that getListProperty would be used on lists built up with addNodeBefore/After calls, and/or by joining different lists together. So just make sure to use addNodeAfter before each equilibration that you want to keep the output from. Consider that a bit like using menu option 4 in standalone alphaMELTS.

If you need to do a few extra calculations to get to where you need to be at the start of the run, then you can do those in place. Then use addNodeAfter and First.clearNode to clean up the list before getListProperty. This is a bit like menu option 3 in standalone alphaMELTS, except that at the moment it's your responsibility to discard the first node before the main calculation.

Note that in the Python version you should always use clearNode() and not removeNode(), as clearNode() tidies up the list after removing the node i.e. makes sure Prev, Next, First, Last is correct for all remaining nodes. removeNode() is a private method in MATLAB so this is not an issue in that version.

Likewise the copy. If you assign path2 = path1 then path1 and path2 point to the same object. This behavior is inherited from the matlab.mixin.Copyable class that MELTSdynamic builds on. Users going from MATLAB to Python often get tripped up in a similar way. So the intended way to start a new list is:

path2 = path1.copy;
        ... or ...
path2 = path1.copyAndKeepOutput;

Both copy the system settings and input parameters to the first node in a new list. The first one resets the output, the second one does exactly what it says on the tin. For an example of copyAndKeepOutput see tutorial.m.

Hope that makes sense!
Paula