search

kaylai / VESIcal

A generalized python library for calculating and plotting various things related to mixed volatile (H2O-CO2) solubility in silicate melts.
MIT License
27 stars 9 forks source link

Bug in Zhang and Duan EOS #82

Open simonwmatthews opened 4 years ago

simonwmatthews commented 4 years ago

From discussion in issue #46.

The fugacity function seems to be being called with pressure = np.nan. Recreated using a moderately high CO2 concentration.

/Users/simonmatthews/repos/VESIcal/VESIcal.py:2602: RuntimeWarning: overflow encountered in exp return Pnp.exp(lnfc) /Users/simonmatthews/repos/VESIcal/VESIcal.py:4157: RuntimeWarning: invalid value encountered in double_scalars CO2_CO2 = ((44.01XCO2)/(44.01XCO2+(1-(XCO2+XCO3))FW_one))100 /Users/simonmatthews/repos/VESIcal/VESIcal.py:4158: RuntimeWarning: invalid value encountered in double_scalars CO2_CO3 = ((44.01XCO3)/(44.01XCO3+(1-(XCO2+XCO3))FW_one))*100 /Users/simonmatthews/repos/VESIcal/VESIcal.py:1603: RuntimeWarning: The pressure is outside the calibration range of the Eguchi & Dasgupta (2018) carbon model, as nan bar is not between 500.0 and 50000.0 bar. The pressure is outside the calibration range of the Kerrick and Jacobs (1981) EOS model, as nan bar is not between 1.0 and 100000.0 bar. warnings.warn(self.calib_check,RuntimeWarning)

simonwmatthews commented 4 years ago

This appears to be a bug in the Zhang and Duan fugacity model code, where an infinite fugacity is returned above a certain pressure.

simonwmatthews commented 4 years ago

(Added you Penny so you'll know when the problem is fixed)

PennyWieser commented 4 years ago

Thanks @simonwmatthews - what other models will this affect? Wondering if I should stop testing till you've fixed it. I'm also finding offsets with Dixon vs. VolatileCalc and the Dixon macro

simonwmatthews commented 4 years ago

This issue won't affect others, but I think some of the models are converting degC to K multiple times, i.e. adding too many factors of 273.15, which I will fix later. This might explain the offset.

simonwmatthews commented 4 years ago

The initial problem has been fixed, i.e. the attempts to calculate negative/nan fugacities (mostly due to pressure unit issues). However, this has highlighted that the ZD EOS is returning the wrong results. I have no idea way. The Eguchi model returns different results when used with other EOS, and so one solution is to remove the Eguchi model for the time being?

PennyWieser commented 4 years ago

I have a version of the Z+D EOS coded up in matlab from Jonathon tucker? Not sure if it would help if I share that? How much difference is there between the Eguchi results with the 2 EOS? How do I change the one its calling so that I can test it aganst the calibration dataset?

simonwmatthews commented 4 years ago

That could be useful. It might also be possible to compare the code with the ThermoEngine SWIM model, as this is based on the same form of EOS. I'm sure it must be factors of ten going wrong somewhere, but I haven't been able to figure out where. I'm moving on from this now in order to try to get other stuff done. We can always add the model back post-submission if we can't get it to work before then.

You can change the fugacity model by: eguchi = v.default_models['EguchiCarbon'] fm = fugacity_KJ81_co2() eguchi.set_fugacity_model(fm)

Then do calculations by passing the eguchi model object instead of specifying model='EguchiCarbon', i.e. v.calculate_saturation_pressure(model=eguchi,...)

simonwmatthews commented 4 years ago

I've temporarily removed EguchiCarbon from the default_models list, which means it won't be accessible without specifically initialising a model instance. Once its working it can be reinstated easily.