geodynamics / burnman Goto Github PK

all of it.

It occurred to me that we are missing 1D seismic models for ak135 and iasp91 in seismic.py. Should be a simple enough thing to include them...

Dear Burnman developers,
I have been trying to fit various datasets of mineral physics data (EoS and cp) with a script based on example_fit_eos.py and frequently encounter problems of this type:
E_th = 3. * n * constants.gas_constant * T * debye_fn_cheb(debye_T / T) File "/Users/ruedas/chemrock/burnman/burnman/eos/debye.py", line 94, in debye_fn_cheb assert(x > 0.0) # check for invalid x AssertionError
I have determined that at least for one my test examples, the cause is that the Debye temperature is allowed to assume the value 0 in a call of debye_fn_cheb via thermal_energy. I think it should be allowed to impose constraints on the fitting parameters in order to prevent the routine from crashing. As far as I can tell, the use of scipy.optimize.least_squares, which seems to offer that feature in its newest version, is not implemented in the current master. Is it planned/worth considering in the future?

Generally, it would be desirable to have functions return gracefully even in the case of failure in order to permit adjusting parameters or skipping results while continuing the executation of a more comprehensive script.

And the use of 'names'

It would be good to add both the internal/helmholtz energy and gibbs/enthalpy to static equations of state, as these are both useful in their own right and also key components of thermal equations of state.

P = -dU/dV | S = -dA/dV | T
V = dG/dP | T = dH/dP | S

At 0 K, S=0 J/K/mol.
The integration only needs to be done for one of the four potentials; Legendre transformations will take care of the others. Which integration is easiest (VdP or PdV) depends on whether the EoS is explicit in V or P.

This is a really old file which has only a single limited use. It should be replaced with more general functions utilising linear algebra. Specifically, from a bulk composition and list of potential minerals, output:

• A potential set of mineral compositions and proportions
• A set of reaction vectors which do not change the bulk composition

We'll eventually need the same functions in burnman.equilibrate anyway.

Several minerals require a magnetic contribution to the Gibbs free energy. A parameterisation has been provided by Sundman, which should be easy enough to implement...

1. Many of the computations in planetary builder loop over the layers (for gravity, pressure, mass, moment of inertia). The computations for the separate layers could be within the layer class (including a value for the starting value of the integration). This would mean a layer can also be used separately, e.g. currently example_geodynamic_adiabat.py (which produces the input to aspect) is using a layer-type-thing (the mantle), but it uses its own gravity-calculating-function (with a starting gravity value). I think it would be nicer to generalize this more, so layer is an independently usable thing.

2. Currently layer can handle constant or linear temperature. An adiabat would be nice, defined with a T0 at the top of the layer (or T_pot at 0 GPa). Within planet, T0 could be computed from the temperature at the bottom of the previous layer, and potentially a predefined jump across the layer boundary (or the thermal boundary layers). Another potential model would be a geothermal read in from a script that is then interpolated to the depths at which layer is defined.

3. Currently temperature is implemented as a separate layer class. I don't know why we don't have functions within layer that would return or set the temperature of a layer. This would make it easier for people to implement their own function for the temperature of a layer, and potentially include thermal boundary layers or an internal heat source.

Currently, solution models have many nested loops, where presumably a matrix operation would be much more efficient...

At the moment, the solid solution init takes a series of order-dependent inputs. This is ok, but would it be better to have a param dictionary to avoid confusion? If so, we should make the change soon before many people use the new functionality...

May I ask the reference of liquid water EOS in minerals/other.py?

We should make names consistent for the SolidSolution constructor. Specifically:

Change SolidSolution.type to SolidSolution.solution_type
Add name to SolidSolution.init

Then, we should add a use of the constructor to example_solid_solution.py. Something like:
garnet = burnman.SolidSolution(name = 'garnet', endmembers = [[py(), '[Mg]3[Al]2Si3O12'], [alm(), '[Fe]3[Al]2Si3O12'], [gr(), '[Ca]3[Al]2Si3O12'], [andr(), '[Ca]3[Fe]2Si3O12']], solution_type = 'asymmetric', alphas = [1.0, 1.0, 2.7, 2.7], energy_interaction = [[2.5e3, 31.e3, 53.2e3], [5.e3, 37.24e3], [2.e3]])

Not sure which is right... I don't think it gets used.

Now that we have the start of a working gibbs minimizer, we should start to think about setting up a high-level function to automatically calculate stable phases along geotherms, or within a 2D P-T grid (or 3D P-T-X, or 4D P-T-X-X etc etc).

To this end, we need the following:

• Change the composite function so that it doesn't renormalise amounts to one.

• Make sure that seismic velocity functions etc. still produce the same values.

• Some way to intelligently pick different starting guesses for compositions (at the moment, we just use whatever springs out of linalg.lstsq)

• On a related note, linalg.lstsq sometimes fails. Dunno why. Would be good also to add dependent endmembers and force each to be >0 (the dependent endmember amounts can then be folded back in to the composition)

• An algorithm to add/subtract phases at a point until the equilibrium assemblage has been found. (0D)

• An algorithm to find the next point along a P-T path where a phase becomes stable/unstable (1D)

• Algorithms to construct a higher level phase diagram (using Schreinemakers Rules; 2D)

  Inevitably, we'll find places where the gibbs minimization functions fall over, so we'll make them more robust as and when we find those problems.

Seems odd to have Grueneisen, eta, alpha etc. with EOS's that don't require them. Can throw errors...

After running tests, the directory structure is full of annoying output files. They can be removed with the line:
rm -f */*.error */*/*.error */*-e */*/*-e contrib/tutorial/uncertainty.dat

Can we write this into the test.sh script somehow? Maybe appending a flag if you want to keep them?
test.sh --keep_output would be nice.

Here is our current code coverage: http://burnman.org/coverage/

Click on the "coverage" column to sort by it. Of course we can ignore all the minerals, but things like eos/* should be tested better than it currently is.

The most recent release of burnman is now very old. We should make another minor release. I've kicked the more major changes down the road (see Issue #388), but these things should be done before we release 0.10.0:

• Merge #383
• Merge #384
• Merge #386
• Merge #387
• Merge #389
• Merge #394

When running composite minerals of composites (see example_composition, example 2), produces error:

File "example_composition.py", line 79, in
[minerals.SLB_2011.mg_fe_perovskite(0.2),
TypeError: init() takes exactly 1 argument (2 given)

Thoughts?

The Debye model for heat capacity developed by S+LB is great for a lot of simple materials. For others, landau-type transitions lead to spikes in heat capacity (particularly at low temperatures), and other contributions (e.g. magnetic) can also be important.

I think it might be quite nice to add an excess entropy/heat capacity contribution with a form that makes sense in the framework of magnetism/low T landau transitions. As S = \int Cp/T dT, it would make sense to make this
S_excess = a + bT
and thus
Cp_excess = bT (which does a reasonable job in modelling magnetic contributions)

What do we think? Both about the idea, and where it should go? HP2011 already has such terms in the polynomial for Cp and as S_0.

• updated flowchart with workflows including input/output
• Ipython pvt fitter
• planet builder example
• seismic example (earth->ASPECT->seismic data)
• phase diagram example (pyrolite?)

mineral.py includes:
@material_property
@copy_documentation(Material.grueneisen_parameter)
def grueneisen_parameter(self):
if self.temperature < 1.e-12:
return 0.
else:
return self.thermal_expansivity * self.isothermal_bulk_modulus
* self.molar_volume / self.heat_capacity_v

This actually means the many (differently) implemented grueneisen parameters in the eos's are not being called. I think we should just change this to self.method.grueneisen_parameter(..). Agreed?

At the moment, we have hardwired the isothermal method into the thermal equations of state. It would be good to separate methods, so that the user can interchange isothermal and thermal parts. This flexibility is standard in other codes. We can always add a warning if we want a user to use a specific isothermal and thermal combination.

Is there a reason we supply EoS methods as strings? This is not very flexible and doesn't allow the user to add their own. It turns out we already support passing in a classname, so the following works:

rock.set_method(burnman.slb.SLB3)

I think passing in an instance instead of a classname would be better (so that the user can have EoS with parameters. Then you would need to write

rock.set_method(burnman.slb.SLB3())

but we could also change the class hierarchy to allow things like

rock.set_method(burnman.slb.SLB(3))

Thoughts?

create_burnman_readable_perplex_table.py needs a few more details for users:

• Make it clear that the PerpleX output script requires the full path and executable name for werami, and the project name without .dat.
• Make a note of which files need to be moved to the current directory/which would work with global path names

I have no clue what changed, but all of a sudden I'm having trouble running test.sh as the tests will fail to import burman. The only way I can fix this is to force a switch in directory:
abspath = os.path.abspath(file)
dname = os.path.dirname(abspath)
os.chdir(dname)
Does anybody have a clue or a better solution?

Can we decide on a new standard blurb at the top of each file? Something to include CIG, and add 2015? Also, egocentric, I know, but can we add me? :)

We should also change the Github link from burnman.org to the CIG page.

• Document and test the ability for setting fractions in different ways. Clean it up.
• Is phase_names used somewhere? Why do we need it? Should it be public?
• Does nesting composites work?
• Deprecate or remove main.velocities_from_rocks()
• Remove volume fraction
• Remove set_composition() for now.
• Move check_pairs() to material.py
• Do we need set_averaging_scheme(), or can that just be provided on initialization?
• Should set_method() reset the state? (Yes, for everything)
• Remove composition() for now
• Don't use try/except in unroll()
• Remove evaluate_phase_properties()
• Remove calculate_elastic_properties()
• Remove fugacity and chemical_potentials for now
• TODO: add the thermodynamic potentials in Composite
• Don't set internal variables for material properties
• Document aliases and derived functions for material (Timo)
• Document aliases and derived functions for mineral (Timo)
• Document aliases and derived functions for composite (Timo)
• Copy documentation from base class.
• Remove _calculate_moduli()
• Don't use try/except to for mass_to_molar_fractions
• Test material.material_property() (Timo)
• Make aliases into properties, document them
• Material should set self.name == __tostring or whatever
• Remove mineral.composition()
• Rip out solid solution endmember names and compositions
• remove SolidSolution.composition
• Call new Reuss function in solid solution shear modulus averaging
• Look over example composite

• 3d seismic model, import to axisem (Sanne)
• mineos output inside BurnMan
• planet.layer temperature as function (Timo)
• planet builder with self-consistent conductive/convective profiles (+heat generation)
• liquid iron bug (Bob)
• liquid database (Bob)
• EOS statistics
• Gibbs minimizer
• ASPECT -> seismic data (1d averaged, 2d/3d data)
• pvt statistics (Bob)

Currently, the code is written for BM 3rd order expansion, with a second or third order expansion for the shear modulus. This is not obvious in the user code (although it is in the documentation)

Would it be better to define two parameters, and then have two separate flags: self.shear_order and self.bulk_order, to clarify what is actually being done in the calculations?

sympy is a really valuable package for a lot of thermodynamics. It would be really good to have this installed on the tester

Calling set_method should reset the internal state of the Material completely and the user should be forced to call set_state before querying any information. If not the user should get a helpful error message. Make sure this works (also with the caching mechanisms in place!).

*** End of file "/home/bob/projects/burnman/misc/ref/step_3.py.out" reached while trying to read line 1001.
*** File "step_3.py.tmp" has more lines than file "/home/bob/projects/burnman/misc/ref/step_3.py.out",
*** line 1001 is the last one read from file "step_3.py.tmp"

! step_3.py ... FAIL

Check: /home/bob/projects/burnman/tutorial/step_3.py.tmp /home/bob/projects/burnman/misc/ref/step_3.py.out

1001 <== /usr/lib/python2.7/dist-packages/matplotlib/axes/_base.py: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal

      ==>

*** End of file "/home/bob/projects/burnman/misc/ref/step_3.py.out" reached while trying to read line 1001.
*** File "step_3.py.tmp" has more lines than file "/home/bob/projects/burnman/misc/ref/step_3.py.out",
*** line 1001 is the last one read from file "step_3.py.tmp"

+++ File "step_3.py.tmp" differs from file "/home/bob/projects/burnman/misc/ref/step_3.py.out"
done

At the moment, we have tried to optimise set_state by not resetting state if P,T is the same as stored in memory.

When we're trying to fit experimental data, sometimes we want to change a parameter and look at the same P, T. I think it would be useful to have either a separate function, or an optional argument to set state to allow us to calculate at the same P, T if an important parameter has changed. Doing a full param check is probably a bit excessive, so we might also want to add a warning the first time we trip the P, T = P, T old statement.

• New tutorial
• Move old tutorial
• Add existing examples (or add a comment to the documentation if the examples are not needed)
• Improve autogen documentation in the examples, especially for things not covered in the tutorial
• Double check that the autogenerated documentation for the code base is complete
• Mention/describe the contributions and ipython directories

• update list of equations of state (and include references)

Separate question. I'd prefer to keep the code overview out of the manual and just online. I'd never use it myself in this document, and I'm afraid someone one day might print this.

Or did anybody have somethings against this? I forgot.

Not an issue, per-se, but definitely something we should add.

Currently, in the literature, most academics report the uncertainty in parameter i of a least squares fit as sqrt(cov[i][i]). Given the large covariances in equation of state fitting, this is a severe misuse of statistics. I would like to see something better, for example, for EoS fitting:

• Visualisation for the covariances
• Confidence intervals in volume at a given P, T (using the delta method?)

I am getting

$ python test_solidsolution.py -v
test_1_gibbs (__main__.test_solidsolution) ... /home/heister/Dropbox/work/burnman-git/burnman/averaging_schemes.py:576: UserWarning: Oops, called reuss_average with Xi<=0!
  warnings.warn("Oops, called reuss_average with Xi<=0!")

@bobmyhill did you make the test? Is it correct to ignore this warning here?

Currently, mineral.pressure, mineral.temperature and mineral.params can be written by the user.
This is horrible, because it allows the user to reset apparent temperature (or params, or anything else) without recalculating everything - leading to inconsistencies.

Instead, we should make pressure and temperature settable only through set_state, and params only settable through set_params. When set_params gets called, mineral.pressure and mineral.temperature gets deleted.

Same thing for solid solutions.

Fun fact: The values in the Stixrude and Lithgow-Bertelloni paper are not those used in hefesto or perplex. They are rounded severely enough that there is a noticeable discrepancy in transition pressure for some reactions (like wad-rw at 1673 K, for example, which is about 0.5 GPa off). I guess we missed this in benchmarking, as the minerals we chose to benchmark against were less affected by the rounding.

We should get the original dataset from somewhere. The PerpleX database has all the full values, so it would be a relatively easy task to write a mini program to convert stx11ver.dat to burnman-readable classes.

Move all the different equations of state to one directory, so they are easily scanned by ocular inspection.

(from a discussion from hangouts today)

The functionality for averaging in main is somewhat restrictive (fixed order of certain parameters). In the light of adding other parameters to be averaged (thermal expansivity, heat capacities, ...), this becomes a problem.

Suggestion:
Add a new class that gets in the constructor:

• a rock
• an averaging scheme
• list of p, T

The class would have individual functions that return vectors of data for v_s, v_p, moduli, density, etc.. I guess we can change the interface of averaging scheme (maybe to operate on a list of minerals and a single given p and T?), which will be used inside this new class.

Internally the class could compute the required things in a lazy fashion (compute v_s from the moduli only when needed, compute the moduli only once when needed, etc.).

I am not sure how to name this class, maybe "Model"?

Usage would look something like this:

rock = ... # type Material as before
p, T = ... # lists as before
model = Model(rock, p, T, VoigtReussHill)

plt.plot(p, model.v_s())
# and model.K(), model.G(), model.density(), ...

[vs_err,rho_err]=burnman.compare_chifactor \
            ([model.v_s(), model.density()],[seis_vs,seis_rho])

This is much cleaner than what we have right now!

Questions:

• What would a good name for this be? Model? Earth? Realization?
• We would lose the ability to look at moduli for individual phases. Is that a problem? I guess one could offer the ability to look at the individual elastic properties before they are averaged (like calculate_moduli does).
• Other thoughts?

@ian-r-rose is the plot supposed to look broken?

As well as the other calculate_moduli functions scattered around.

• ipython notebooks: planet builder (cayman)
• online clicky thingy:
  • pvt fitter with uncertainties
  • planet builder (simplified)
  • make profiles example
• put examples on website
• exoplanet core size and effects on 3D dynamics (cayman)
• ASPECT high resolution, cluster stuff
• Mg/Si ratio for real
• thermal structure for planet builder
• perplex input -> rock (Cayman, Bob)
• quasi-chemical solution models (Bob)
• Add stoichiometry calculators to Planet builder [using Bob's chemistry tools) (Cayman)

Readjust burnman to set the EoS within the minerals... and just have that automated if the EoS is set in params.
I guess some sort of warning should pop up if minerals are mixed with different EoSs.