Recent posts

#81

alphaMELTS for MATLAB/Python / Re: Unexpected behavior during...

Last post by Paula - September 09, 2023, 04:46:26 PM

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

Paula

#82

alphaMELTS for MATLAB/Python / Re: Unexpected behavior during...

Last post by OEM - August 30, 2023, 01:30:25 AM

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

#83

alphaMELTS for MATLAB/Python / Re: Unexpected behavior during...

Last post by Paula - August 28, 2023, 10:45:29 AM

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

#84

alphaMELTS for MATLAB/Python / Re: Total failure of Matlab wh...

Last post by OEM - August 27, 2023, 07:45:50 AM

My matlab code crashes also. Here is an interesting observation from the log file. It tries to add a phase to the table.

Code Select Expand

Assertion in struct mxArray_tag *__cdecl `anonymous-namespace'::MxCreateStringNCharsT_safe<char>::operator ()(const char *,const char *,unsigned __int64) const at B:\matlab\foundation\matrix\src\matrix\array2.cpp line 500:
mxCreateString called with non-UTF-8 input: bulk                oxygen              liquid              sphene              clinopyroxene       clinopyroxene       clinopyroxene       clinoamphibole      clinoamphibole      clinoamphibole      biotite             alkali-feldspar     plagioclase         quartz              tridymite           cristobalite        leucite             rutile              spinel              rhm-oxide           ortho-oxide         water               ��Q��␂


#85

alphaMELTS for MATLAB/Python / Unexpected behavior during P-T...

Last post by OEM - August 27, 2023, 02:19:23 AM

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

#86

alphaMELTS for MATLAB/Python / Re: Fractional Crystallization...

Last post by agion - August 15, 2023, 05:20:27 AM

Paula,

Thank you for the advice. Your first suggestion for the output flags is working well. I have not noticed any problems with the run speed and all the calculations have been stable.

Thanks,
Austin

#87

alphaMELTS for MATLAB/Python / Re: Fractional Crystallization...

Last post by Paula - July 30, 2023, 02:44:57 AM

From MELTSengine.m (these refer to runMode and outputFlag, which are the first and second codes that you input):

% Calculation mode is set in parent MELTSdyanmic list: 1. rhyolite-MELTS 1.0.2; 2. pMELTS; 3. rhyolite-MELTS 1.1.0; 4. rhyolite-MELTS 1.2.0.
% Output is post-equilibration / pre-fractionation (1) by default, but calculation can also be batch, no fractionation (0), or output can be after fractionation (2).

The value of outputFlag affects what ends up in quantities labelled as OUTPUT only in the MELTSengine.m comments. Quantities that are labelled "INPUT and OUTPUT" will always be updated based on the last thing to happen. So by default (outputFlag=1) the bulkComposition variable will be the values after fractionation, whereas dispComposition("bulk") will be the values after equilibration and before fractionation (similar to what goes into the text output).

When you add a node it will inherit the (post-fractionation) bulkComposition values and the dispComposition("bulk") values will be reset (unless you use copyAndKeepOutput). If this isn't the order of events that you need then you can try calling calcEquilibriumState twice and add the new node in between e.g:

• 1st time outputFlag = 0: only equilibrium state is calculated and everything is pre-fractionation
• 2nd time outputFlag = 1: bulkComposition is post-fractionation, dispComposition("bulk") is post-equilibration / pre-fractionation

Or:

• 1st time outputFlag = 1: bulkComposition is post-fractionation, dispComposition("bulk") is post-equilibration / pre-fractionation
• 2nd time outputFlag = 2: everything is post-fractionation (so you can see the small amounts of fractionating phases that are retained by the system)

Either way, the equilibration in the second call to calcEquilibriumState is trivial and shouldn't affect the run speed (but if it does, please let me know).

Hope that helps,
Paula

#88

alphaMELTS for MATLAB/Python / Fractional Crystallization in ...

Last post by agion - July 17, 2023, 07:36:50 AM

When running a fractional crystallization model in MELTS for MATLAB is it possible to determine when during the process the solids are removed from the system? At the moment it seems like the solids are removed after each equilibration, but is it possible to change this to remove the solids after adding the next node? Or should I simply perform the calculation with fractionation turned off and manually remove the solids from the bulk composition at each node?

Thanks,
Austin

#89

Software tools, support and development / Re: C++ API for alphaMELTS 2 o...

Last post by Paula - July 04, 2023, 09:01:39 AM

I believe ENKI has some C++ code but I don't know the current status.

The alphaMELTS for MATLAB/Python package includes a C library of functions that you could call from C++. Check the MELTSdynamic.h header file for the function prototypes.

Note: some of the variable names in MELTSdynamic.h could be a little deceptive e.g. *entropy is actually pointer to whatever reference quantity is relevant for the current calculation mode, so it could just as easily be reference enthalpy or reference volume, or it may be ignored for isobaric/isothermal calculations.

If you are comfortable reading either Python or MATLAB code then look at meltsengine.py or MELTSengine.m to see exactly what is passed and what is 'returned' (everything is passed by reference, though MATLAB kind of tries to hide this). Search for 'self.status.libalphamelts' in the Python file or 'calllib' in the MATLAB file to locate where the C library functions are invoked.

At one point I was planning to wrap the C library functions in C++, as calls to MATLAB's C++ Library interface resemble calls to Python's ctypes more closely. That had the potential of making it easier to develop and maintain both versions simultaneously, among other advantages.

At the time the MATLAB C++ Library interface was quite new and I couldn't find a way to force it to reload the C++ library, which is a useful feature of the current alphaMELTS for MATLAB version. It looks from the MATLAB documentation like it is possible to reload the library now, so I may revisit writing simple C++ wrappers soon and would welcome any input you have.

Thanks,
Paula

#90

Software tools, support and development / C++ API for alphaMELTS 2 or ME...

Last post by Gabriel Girela - June 28, 2023, 06:26:27 AM

Hello,

I wanted to know if there is any C++ API so I can use rhyolite-MELTS inside another C++ code.

If this is not the case, is there a way of using alphaMELTS2 but passing the parameters like composition, pressure, temperature, etc. as terminal arguments, instead of through the input melts file?

Thank you very much,

Gabriel Girela