Unexpected behavior during P-T grid simulation

[OEM]

Dear Colleagues!
I am trying to simulate phase diagram for Pinatubo magma composition in a wide range of pressures and temperatures using Matlab.
When I restrict the minimum temperature to 780 C results looks reasonable. For the minimum temperature of 740 C most of the plagioclase table is filled with a constant value.
740: https://disk.yandex.ru/i/o3ZM38bk-rvTJQ
780:  https://disk.yandex.ru/i/YtnCTTwMyu18pg
What am I doing wrong? The code is below:

Code Select Expand

clear all
close all
P = 2000:250:7750;
T=1000:-20:780;
bulk = zeros(19, 1);
% SiO2, TiO2, Al2O3, Fe2O3, Cr2O3, FeO,        MgO
%                         1     2     3      4        5     6          8
bulk([1:6 8 11:16])  = [64.6  0.53   16.5  0.15*4.37 0.00 0.85*4.37  2.39 ...
    5.23 4.49 1.54 0.00 6.24 0.4];
%, CaO, Na2O, K2O, P2O5, H2O CO2
%    11    12  13   14  15

bcom=100*bulk/sum(bulk);
%% carry out equilibrium crystallisation modelling
ptpath=MELTSdynamic(1);
for i=1:numel(P)
    for j=1:numel(T)
        disp([P(i),T(j)]);
        % Change P-T conditions
        ptpath=ptpath.addNodeAfter;
        ptpath.engine.set('BulkComposition',bcom);
        ptpath.engine.setSystemProperties("Log fO2 Path","+2FMQ")
        ptpath.engine.setSystemProperties("Mode", "Fractionate Fluids");
        ptpath.engine.pressure = P(i);
        ptpath.engine.temperature = T(j);
        % calculate the equilibrium
        ptpath.engine.calcEquilibriumState(1,1);
        disp(ptpath.engine.status.message);
        liq=ptpath.engine.getProperty('dispComposition','liquid1');
        if liq(19)~=0,  break; end
     end
end
len=length(ptpath.getListProperty('pressure'));
N={'P','T'};
Conditions=table(zeros(len,1),zeros(len,1));
Conditions.Properties.VariableNames=N;
Conditions.P=ptpath.getListProperty('pressure')';
Conditions.T=ptpath.getListProperty('temperature')';
Names=ptpath.systemNames; nminerals=numel(Names);
Mass=array2table(zeros(0,nminerals));
Mass.Properties.VariableNames=Names;

Mass{1:len,1}=ptpath.getListProperty('mass',Names(1))';
Mass{1:len,2}=ptpath.getListProperty('mass',Names(2))';
for j=3:width(Mass)
    Mass{1:len,j}=ptpath.getListProperty('mass',strcat(Names(j), '1'))';
    if Names(j)=="rhmoxide1"
        Mass{1:len,j}=ptpath.getListProperty('mass','rhm-oxide1')';
    end
    if Names(j)=="kfeldspar1"
        Mass{1:len,j}=ptpath.getListProperty('mass','k-feldspar1')';
    end
    if Names(j)=="kfeldspar2"
        Mass{1:len,j}=ptpath.getListProperty('mass','k-feldspar2')';
    end
end
arr=table2array(Mass);
arr(isnan(arr))=0;
plag=arr(2:end,18);
bulk=arr(2:end,1);
Plag=reshape(plag./bulk,numel(T),numel(P));
surf(P,T,Plag); shading interp; hold on
contour3(P,T,Plag,0:0.05:0.5,'Color','w');
xlabel('Pressure,bar')
ylabel('Temperature, ^oC')
zlabel('plag. content')
tabb=[Conditions,Mass];
Best regards 
Oleg Melnik

[Paula]

Try updating to the latest version here and see if it helps. The newest version has k-feldspar renamed to alkali-feldspar (see the alphaMELTS 2 README).

There is definitely a bug fix for the other problem you mentioned in this topic reply. When I added oxygen as a "phase" (in a similar sense to how "bulk" is treated as a phase) I forgot to add extra storage space. It should be better now.

Paula

[OEM]

Thank you, Paola! The code is more stable now!
Is there a manual how to set  ptpath.engine.setSystemProperties with all available options? For example, I would like to change the maximum number of iterations, because for some parameters the code does not converge.
Best
Oleg

[Paula]

Hi! Sorry for delayed response. Please could you post a code example with parameters for which is doesn't converge?

Paula

[OEM]

Here is the code, Paola.

Code Select Expand

clear all
close all
%% load in appropriate packages and functions
addpath('D:\YandexDisk\Work\_Projects\Melts_for_Matlab\package')
warning('off','all')
P=[4000 4000 4000 4000 4000 4000 4000 4000 4000 8300 8300 8300 9400 9400 9400 9600 9600 9600 9600 9600 9600 9600 9600 9700 9700 9700 9800 9800 9800 2212 2212 2238 2238 2238 2238 2245 2245 2300 2300 2300 2300 2080 2080 2080 2080 2250 2250 2250 2250 2250 2250 2237 2237 2237 2237 3890 3890 3890 3890 2174 2174 2174 2174];
T=[900 900 900 900 900 950 950 950 950 892 892 892 750 750 750 841 841 841 841 841 943 943 943 892 892 892 995 995 995 785 785 776 776 776 776 776 776 834 834 834 834 866 866 866 866 899 899 899 899 899 899 781 781 781 781 780 780 780 780 747 747 778 778];
H2Ogl=[9 8.30000000000000 7.50000000000000 6.70000000000000 4.30000000000000 8.90000000000000 6.60000000000000 5.50000000000000 4.40000000000000 9.30000000000000 8 6.40000000000000 10.5000000000000 9.70000000000000 9 11.7000000000000 10.3000000000000 9.50000000000000 8 7.20000000000000 8.50000000000000 6.80000000000000 6 10.4000000000000 12.9000000000000 8.20000000000000 7.60000000000000 5.30000000000000 4.80000000000000 7 6.69000000000000 7 7 7 7 7 7 7 6.96000000000000 5.76000000000000 6.19000000000000 7 6.21000000000000 5.04000000000000 6.07000000000000 6.96000000000000 6.24000000000000 6.04000000000000 5.66000000000000 3.78000000000000 3.09000000000000 7 7 6.50000000000000 5.83000000000000 8.50000000000000 8.08000000000000 7.27000000000000 8.15000000000000 6.38000000000000 6.35000000000000 7.01000000000000 6.90000000000000];
bulk = zeros(19, 1);
%                       SiO2, TiO2, Al2O3, Fe2O3, Cr2O3, FeO,        MgO
%                         1     2     3      4        5     6          8
bulk([1:6 8 11:16])  = [64.6  0.53   16.5 0.15*4.37 0.000 0.85*4.37  2.39 ...
    5.23 4.49 1.54 0.00   0   0 ];

%, CaO, Na2O, K2O, P2O5, H2O CO2
%    11    12  13   14  15   16
%% carry out equilibrium crystallisation modelling
ptpath=MELTSdynamic(4);
for i=1:numel(P)
        disp([P(i),T(i),H2Ogl(i)]);
        bulk(15)=H2Ogl(i);
        bulk(16)=0;
        bcom=100*bulk/sum(bulk);
        ptpath.engine.set('BulkComposition',bcom);
        ptpath.engine.setSystemProperties("Log fO2 Path","+2FMQ")
        ptpath.engine.pressure = P(i);
        ptpath.engine.temperature = T(i);
        ptpath.engine.calcEquilibriumState(1,1);
        ptpath=ptpath.addNodeAfter;
end

After several successful calculations the system stops responding and matlab closes without report generation. Windows, Matlab R2023a, update 5.
Best
Oleg

[Paula]

Hi! Sorry for the delay again. 

I wasn't entirely able to reproduce the problem (on Windows I got a crash but right at the end when it was pretty much done, on Mac no crash, and on Linux crash early on more similar to what you describe). The variability across platforms is a good indication that it's a memory problem.

I found a couple of issues that I've fixed but can't be sure they are what you were hitting. I'm going to keep testing, combine some other bug fixes and new features, and aim to post an update to GitList some time next week. I'll let you know when it's there.

Thanks for your patience!
Paula