Kernel death for magmasat

[PennyWieser]

I've got more magmasat compositions to test vs. app results for some more varied compositions from George Cooper

For some reason, I get Kernel death when running these samples en-mass. It gets down to printing "calculating sample 28" (the last one), then dies. Through literally running groups using trial and error,, I found its sample 5. Is there a way to get Magmasat to at least give an error of which sample is causing Kernel death as trying to find it made me loose my mind? Also, while sometimes it pops up kernel death, sometimes it just gets stuck forever (e.g., when running it alone). I've attached the whole dataset, and a separate file with the one that breaks it. Statia_MI_Cooper_et_al_5.xlsx Statia_MI_Cooper_et_al.xlsx

[kaylai]

Hi, sorry, I think I overlooked this thinking it was another issue. Are you still having this problem?

[PennyWieser]

Yus, I have uploaded the files and a notebook into the "testing" folder. the spreadsheet ending _5 contains the sample which causes Kernel death. This alone isn't really the issue, the issue is, when you run this sample within the rest of the samples, you don't get any outputs for any of them. So unless users wanted to go through literally line by line like I did to find the broken one, you can't use the magmasat model. Is there a way to get it to move on instead, or is kernel death inevitably going to kill all the work its done so far? In that case, can you see where it happened so you can delete that sample?

[kaylai]

Working on this. Unfortunately I think this will be a tough nut to crack. The reason the kernel dies is due to a segmentation fault in the thermoengine code, which means the C code passes a bad value somewhere and shuts everything down. There is no way to catch segmentation faults in python. However, this has been an ongoing issue with thermoengine, so for posterity I used a fault debugger built into python to at least get the traceback to show where the fault occurs. Here is the traceback:

$ python -Xfaulthandler fault_test.py 
Calculating sample 5
duanDriver-2: z = 2.71854, v = 2.1267, delv = 0, iter = 37
duanCO2Driver: z = 3.06283, v = 4.44978, delv = 0, iter = 34
Fatal Python error: Segmentation fault

Current thread 0x000000010d0225c0 (most recent call first):
  File "/Users/kiacovin/.local/lib/python3.7/site-packages/rubicon/objc/runtime.py", line 833 in __call__
  File "/Users/kiacovin/.local/lib/python3.7/site-packages/rubicon/objc/runtime.py", line 908 in __call__
  File "/opt/anaconda3/lib/python3.7/site-packages/thermoengine/equilibrate.py", line 4246 in equilibrate_tp
  File "/Users/kiacovin/Dropbox/Research/__Manuscripts in Progress/__ENKI/__TheCode/VESIcal/VESIcal.py", line 6055 in calculate_saturation_pressure
  File "/Users/kiacovin/Dropbox/Research/__Manuscripts in Progress/__ENKI/__TheCode/VESIcal/VESIcal.py", line 6747 in calculate
  File "/Users/kiacovin/Dropbox/Research/__Manuscripts in Progress/__ENKI/__TheCode/VESIcal/VESIcal.py", line 1613 in __init__
  File "/Users/kiacovin/Dropbox/Research/__Manuscripts in Progress/__ENKI/__TheCode/VESIcal/VESIcal.py", line 1313 in calculate_saturation_pressure
  File "fault_test.py", line 6 in <module>
Segmentation fault: 11

I'm still not exactly sure why this occurs sometimes, but I'll submit it as an issue on the thermoengine gitlab page.

[PennyWieser]

Hi @kaylai I think thats great - Just knowing which sample is fine. George (who's data this is) actually told me several of them crashed the Mac app too, so its not a fault with just us! Their work around in the cardiff group is to adjust temp a bit till you get a result.

[kaylai]

Oh that’s really interesting. Do you know roughly by how much they need to adjust the temp in order for it not to fail?

[PennyWieser]

When I spoke with Emma Bennett about this, she said only a few degrees. What is interesting is that George said the app failed on about half of his compositions, while VESIcal only fails on a single one. So its better than the App!

[kaylai]

Even more interesting... I pulled out just the calculate_saturation_pressure part of the routine and rand the sample (see below script). It ran fine, and I got an output! Need to figure out exactly where the segmentation fault is happening. Something with check_calibration???

import numpy as np

#--------------MELTS preamble---------------#
from thermoengine import equilibrate
# instantiate thermoengine equilibrate MELTS instance
melts = equilibrate.MELTSmodel('1.2.0')

# Suppress phases not required in the melts simulation
phases = melts.get_phase_names()
for phase in phases:
    melts.set_phase_inclusion_status({phase: False})
melts.set_phase_inclusion_status({'Fluid': True, 'Liquid': True})

bulk_comp = {'SiO2':50.63,
             'TiO2':0.95,
             'Al2O3':19.81,
             'FeO':9.90,
             'MgO':3.81,
             'MnO':0.12,
             'CaO':11.51,
             'Na2O':2.65,
             'K2O':0.42,
             'P2O5':0.10,
             'CO2':0.02426744,
             'H2O':1.6}

temperature = 950

bulk_comp_orig = bulk_comp

feasible = melts.set_bulk_composition(bulk_comp)
#Coarse search
fluid_mass = 0
pressureMPa = 2000 #NOTE that pressure is in MPa for MagmaSat calculations but reported in bars.
while fluid_mass <= 0:
    pressureMPa -= 100
    if pressureMPa <= 0:
        break

    output = melts.equilibrate_tp(temperature, pressureMPa, initialize=True)
    (status, temperature, pressureMPa, xmlout) = output[0]
    fluid_mass = melts.get_mass_of_phase(xmlout, phase_name='Fluid')

pressureMPa+=100

#Refined search 1
feasible = melts.set_bulk_composition(bulk_comp)
fluid_mass = 0
while fluid_mass <= 0:
    pressureMPa -= 10
    if pressureMPa <= 0:
        break

    output = melts.equilibrate_tp(temperature, pressureMPa, initialize=True)
    (status, temperature, pressureMPa, xmlout) = output[0]
    fluid_mass = melts.get_mass_of_phase(xmlout, phase_name='Fluid')

pressureMPa += 10

#Refined search 2
feasible = melts.set_bulk_composition(bulk_comp)
fluid_mass = 0
while fluid_mass <= 0:
    pressureMPa -= 1
    if pressureMPa <= 0:
        break

    output = melts.equilibrate_tp(temperature, pressureMPa, initialize=True)
    (status, temperature, pressureMPa, xmlout) = output[0]
    fluid_mass = melts.get_mass_of_phase(xmlout, phase_name='Fluid')

if pressureMPa != np.nan:
    satP = pressureMPa*10 #convert pressure to bars
    flmass = fluid_mass
    flsystem_wtper = 100 * fluid_mass / (fluid_mass + melts.get_mass_of_phase(xmlout, phase_name='Liquid'))
    flcomp = melts.get_composition_of_phase(xmlout, phase_name='Fluid', mode='component')
    try:
        flH2O = flcomp['Water']
    except:
        flH2O = 0.0
    try:
        flCO2 = flcomp['Carbon Dioxide']
    except:
        flCO2 = 0.0
else:
    flmass = np.nan
    flsystem_wtper = np.nan
    flH2O = np.nan
    flCO2 = np.nan
    warnmessage = 'Calculation failed.'

feasible = melts.set_bulk_composition(bulk_comp_orig) #this needs to be reset always!

print({"SaturationP_bars": satP, "FluidMass_grams": flmass, "FluidProportion_wt": flsystem_wtper,
                    "XH2O_fl": flH2O, "XCO2_fl": flCO2})