Recent posts

#31

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��␂


#32

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

#33

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

#34

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

#35

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

#36

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

#37

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

#38

alphaMELTS for MATLAB/Python / Re: MATLAB for looping T,P spa...

Last post by Ri Cao - June 28, 2023, 01:36:21 AM

Dear Paula,
Thank you very much for your prompt response.

Yes, I have emailed you my modified RC12frac.m.

I incorporated the Find Liquidus function from tutorial.m into RC12frac.m at the beginning of the script, but the problem still exists. Please have a look through it. Thank you very much.

Best regards,
Ri

#39

alphaMELTS for MATLAB/Python / Re: MATLAB for looping T,P spa...

Last post by Paula - June 26, 2023, 11:14:22 AM

Hi Ri,

Can you email me the MATLAB script and input that produced these errors, please? From what you say about removing the "Suppress: leucite" line fixing things, I suspect you've hit a genuine bug and the files would be very helpful to track it down.

Are you calling the Find Liquidus function at the start? If so, make sure to check the status string for errors. I think the RCfrac script just prints the status, rather than examining it. I would expect Find Liquidus to return without crashing MATLAB, even if the liquid composition is outside the range spanned by the MELTS liquid model, so you could use that to bypass the equilibration steps that are more likely to crash MATLAB for an illegal composition.

Thanks,
Paula

#40

alphaMELTS for MATLAB/Python / Re: MATLAB for looping T,P spa...

Last post by Ri Cao - June 20, 2023, 02:53:57 PM

Dear Paula,
Thank you very much for your prompt response and clarification.

Yes, I run the Matlab-MELTS from the terminal.

I have documented the error message from terminal there:
%%%%%%%
Assertion failed: (mxCreateString called with non-UTF-8 input: bulk                oxygen              olivine            orthopyroxene      clinopyroxene      clinopyroxene      biotite            plagioclase        rutile              spinel              water              ?
), function mxArray *(anonymous namespace)::MxCreateStringNCharsT_safe<char>::operator()(const char *, const SourceType *, size_t) const [SourceType = char], file array2.cpp, line 500.
Caught SIGABRT: usually caused by an abort() or assert()
MATLAB_maci64(50312,0x11249be00) malloc: Incorrect checksum for freed object 0x7faf38edf838: probably modified after being freed.
Corrupt value: 0x4500010001000100
MATLAB_maci64(50312,0x11249be00) malloc: *** set a breakpoint in malloc_error_break to debug
Caught SIGABRT: usually caused by an abort() or assert()
Caught SIGILL: illegal instruction
%%%%
Perhaps this is caused by so many minerals precipitating in the system that closes Matlab automatically?


Another scenario where the systems stop automatically and leave an error message. there is:
%%%%%
MATLAB_maci64(51628,0x700003183000) malloc: Incorrect checksum for freed object 0x7f8ba7b2f400: probably modified after being freed.
Corrupt value: 0x0
MATLAB_maci64(51628,0x700003183000) malloc: *** set a breakpoint in malloc_error_break to debug
%%%%%
I have a thought about this, and it seems that this is related to the inhibition of leucite precipitation. For instance, if I remove "Suppress" "leucite" in the settings, then that works.

One exercise I am currently trying to do is generate 1000 random parental magma datasets of Venus composition within the errorbar (1 sigma). I was wondering if there is a way we could ask Matlab to jump or ignore automatically during the run. For instance, perhaps there is a range for each oxide (i.e., SiO2) that we can set before each run? Any thoughts on this will be really helpful! Thank you very much.

Best regards,
Ri Cao