Total failure of Matlab when using MELTS

[slock]

I am having a weird issue where Matlab just completely breaks while doing MELTS calculations. I haven't been able to come up with a nice simple example case because it seems to be slightly path dependent but I have copied the simplest one I could find to reproduce the error which just requires a decreasing T path. If I run this script once I get the error at about i=14:


Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension.

Error in MELTSengine/calcEquilibriumState (line 425)
                properties = reshape(properties, nc, numberPhases);

Error in reproducing_Matlab_crashing_error (line 72)
      path.engine.calcEquilibriumState(1, 1);

 If I run it again Matlab dies with 'Matlab has encountered an internal error and needs to shut down'. I encountered similar errors seemingly randomly when running a number of my scripts. I am running the most recent version of alphaMELTS on Catalina. Any ideas? Let me know if this doesn't reproduce the error for you. 

Thanks a lot!
Simon


Example script:

Code Select Expand



%% SJL 9/19
%Script to reproduce MELTS failing error


%start off with a clean slate
clear all
close all
fclose('all');

%add the function directory
if isunix
    addpath('/Users/simonlock/src/MELTS_Matlab/package')
else
    addpath(strcat(test_script_dir,'\..'))
end

%turn the warning off
warning('off', 'MATLAB:loadlibrary:cppoutput');


%% PARAMS


comp=[45.5652
    0.9372
   14.0558
    1.5063
    0.1067
    8.6226
         0
   15.6014
         0
         0
   10.2931
    2.6090
         0
         0
    0.7029
         0
         0
         0
         0];

  T=949.4395;

  p=711.6284;
 
 
  %%
  path=MELTSdynamic(1);
 
  path.engine.temperature=T-1;
  path.engine.pressure=p;
 
  path.engine.setBulkComposition(comp)
 
  path.engine.calcEquilibriumState(1, 1);
 
  %%
  for i=1:1000
      disp(i)

      path.engine.temperature=path.engine.temperature-0.1;
      path.engine.pressure=path.engine.pressure;
     
     
      path.engine.calcEquilibriumState(1, 1);
  end
 
 

[Paula]

Short answer is that it's getting stuck in a local minimum, which is why it fails. Not unrelated, it's trying to exsolve multiple instances of cpx, spinel etc. that are probably not real. I think ultimately that this so much adding and dropping of phases leads to the segmentation fault, though I haven't tried running it through the debugger or memory checker to find out exactly where yet. Debugging the mixed C and MATLAB in this way is possible in the Linux version of MATLAB but it is extremely slow...

In the meantime, I'm going to outline first how to go about tracking down problems in alphaMELTS for MATLAB/Python, and then some suggested strategies to work around the specific problem that Simon describes.

So first thing is always start MATLAB from the command line on Mac or Linux so that MELTS-related output goes to the terminal. On Windows a console will open automatically. For Python it will depend on which IDE you are using and whether you are on Windows or not. If you are having trouble seeing the output for an example that does work, like the tutorial, then let me know.

Most of the time you do not need to see this output, but if something does go wrong it can help with diagnosis. It is also useful for checking whether a particular setting has been picked up. If you set a breakpoint just after the setSystemProperties() line in the script then you should be able to see a "Processed line..." message in the terminal output.

The other trick is that, by default, the latest version of alphaMELTS for MATLAB/Python will write out alphaMELTS-style output files. (There are other output options but that is for a different thread.) If the MATLAB/Python version crashes then afterwards navigate in terminal to where the output files are and do "run-alphamelts.command -x". This will make the Phase_mass_tbl.txt file without doing any extra calculations, and that file can be very useful for working out what happened.

More soon...

[Paula]

Hi Simon,

So if you run your script as is then on the last step (i=14; index=15 in the _tbl.txt files) it adds a third cpx and a third spinel. The compositions are virtually identical to the second cpx and spinel respectively, so it really shouldn't have added them. If you try to add a break point within MATLAB, e.g. with "dbstop if error", it gives a memory error just like it did when you tried restarting the script. I suspect that is all linked but I'm not sure how.

The simplest way to get the script to run to i=1000 is to add the following:

path.engine.setSystemProperties("Limit number", "spinel 2");
path.engine.setSystemProperties("Limit number", "clinopyroxene 2");

It you limit spinel but not cpx, which is limited to 3 by default anyway, then it gets to i=599 and crashes when hornblende comes in. Likewise if you limit cpx to 2 but not spinel. In both cases there is additional weirdness. I'm not sure why it is adding the extra phases in the way it is but I have some ideas, which may also suggest a more general workaround.

There are several other things to try in addition, and one or more might make your runs go more smoothly:

• Increase the temperature increment if you can. Unless you are trying to model the granite ternary minimum, something like 1 or 2^oC should work pretty well. With 0.1^oC I'm not that surprised you are hitting problems with the very small energy differences between the various pyroxenes etc.
• Possibly remove the tiny amount of Cr₂O₃ from the bulk, as it overstabilizes spinel.
• Start the run at a higher temperature. The terminal output shows that the algorithm is struggling on the first calculation, and I definitely saw one instance where it put metastable hornblende in and then dropped it on the next go.

Anyway, I will try to debug this behavior properly at some point. In the meantime, hopefully that gives you and other something to work with to make the code more reliable.

Paula

[slock]

Hi Paula,

Thank you very much for your reply! I hadn't realized that you could restrict the number of phases of a particular type, that's really useful to know. I hadn't also thought about the issues with small T (and presumably p) step sizes so thank you for that, I will bear that in mind in future. Let me know if you do track down the source of the seg fault. It would be nice if we can make it so the code fails gracefully rather than crashing matlab completely.

Thanks,

Simon

[OEM]

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