Rework main/averaging

[tjhei]

(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?

[sannecottaar]

Can 'averaging_scheme' also become an input to the Composite? Do we need a new class for this?
Composite could just have functions for the moduli, velocities, and depth (giving P and T), and pressures (giving depths and T). These only work if you have fractions and averaging_scheme defined. I like the individual functions, this way, if you change the fractions or averaging_scheme, it would also get recomputed.

Also the misfit functions in main.py could go into tools.py or a new misfitfunctions.py. These two things would empty out main.py.

[ian-r-rose]

I generally agree, Sanne. It would certainly simplify things to be able to ask a Composite for its averaged properties.

[bobmyhill]

I also agree. I made a start on this in #138 yesterday, but the elasticity/thermodynamics averaging isn't finished yet.

One question - before now, there has always been the option to input a bunch of P,T points, but with composite (and mineral, and solid solution) only one P,T point is input during set_state. Is this a problem?

I think depths should not be an input to composite - perhaps some higher level function could be added to burnman.tools to calculate pressure-depth curves for a sphere/infinite half plane of constant composition/layered composition? Eventually, of course, this is something for which equilibrium calculations will be necessary...

[bobmyhill]

I like the idea of being able to compute things individually or as a group. If we do this for Composite, maybe it would be good to also do this for Mineral and SolidSolution? The only issue is that it will need a rewrite, removing attributes and replacing with functions, such that we have things like

pyrope.set_state(P,T) # essentially does nothing apart from updating P, T
pyrope.gibbs() # actually does the calculation

This is a pretty hefty change, but it would vastly increase efficiency...

[tjhei]

The reason for the current design is how we designed Material: For a given P,T you can unroll() a material into an array of phases and fractions, but the phases/minerals can change depending on P and T. This basic design allows for things like switching between two mineral phases based on P for example.

I don't see how we can make this work for a list of P,T values without rewriting pretty much everything.

[bobmyhill]

Was that an answer to my first comment? If so, then I understand. I don't think it needs to be changed - I always calculate things for individual P,T points anyway.

For my second comment, any thoughts? Or alternative ways to make things faster? It seems like this is our last chance to make any big changes...

[tjhei]

I was commenting on

One question - before now, there has always been the option to input a bunch of P,T points, but with composite (and mineral, and solid solution) only one P,T point is input during set_state. Is this a problem?

And the answer is that I don't see an easy way to change this.

For my second comment, any thoughts?

What are you referring to?

[bobmyhill]

Ok, we're on the same page w.r.t the first comment :)

My second comment was this:
"I like the idea of being able to compute things individually or as a group. If we do this for Composite, maybe it would be good to also do the same for Mineral and SolidSolution? The only issue is that it will force a rewrite, removing attributes and replacing with functions, such that we have things like

pyrope.set_state(P,T) # essentially does nothing apart from updating P, T
pyrope.gibbs() # actually does the calculation"

Essentially, my question is how to optimise efficiency (calculating only the required variables) whilst ensuring that the user doesn't end up taking values for variables which haven't been updated since before the last set_state. We discussed this last December, but I guess it got shelved.

Sanne's point about having functions to compute elastic moduli etc. might be a nice way to solve the efficiency problem. Another alternative might be to delete attributes each time set_state is called and reinitialise them with some optional functions. The user input for each method would look something like this:

Current method

rock.set_state(P, T)
rock.gibbs
rock.G

Possibility 1

rock.set_state(P, T)
rock.gibbs()
rock.G()

Possibility 2

rock.set_state(P, T, "thermodynamic")
print rock.gibbs
rock.set_state(P, T, "all")
print rock.G

There are issues with all three. The current method is elegant, but does force the user to compute everything for each set_state (slow even for simple calculations, which is why I implemented calcgibbs). The second method is also elegant, and will be much faster for most purposes, but does force a rewrite, and might be slower if the user needs to compute a lot of related variables. The third method will be backwards compatible, and will potentially be very efficient, but will need careful deletion of attributes for each set_state to make sure the user doesn't use old values by accident.

I'm happy to stick with the old method, but thought I'd bring this up now. Any thoughts?

[tjhei]

I played with computing properties on demand inside EOS here 1. This way things are not computed when you don't need them and we avoida ton of recomputation of things too. I guess we never continued with this.

There is also python properties (see 2) that allows you to write a function get_G() that is called behind the scenes whenever you access rock.G.

[bobmyhill]

Oooh, python properties are cool! I didn't realise those existed, thanks Timo! What do you think, might that be a good way to go? From the user point of view, there would be no change, which is definitely a plus!

[bobmyhill]

An additional point - I don't think we should have setters and deleters, so in this case, attributes would be as simple as:

@property
def x(self):
    self._x = 1.0
    return self._x

Short, simple, easy to maintain - I like it!

[tjhei]

What do you think, might that be a good way to go?

I think we should keep the set_state interface. Down the road we can use properties, lazy evaluation, and caching in Material and EOS. The good news is that we don't need to worry about optimization now and don't have to change the interface later.

For now we should get the changes to the interfaces of Composite that you need right.

[tjhei]

Before we do optimizations a benchmark needs to be set up to measure the improvements and to be able to check accuracy of the changes.