Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added the Vinet EOS #158

Closed
wants to merge 1 commit into from

Conversation

CaymanUnterborn
Copy link
Contributor

Added the Vinet EOS from a copied birch_murnaghan.py. May need to be updated to be compatible, just let me know what needs to be done. Also included the 4th order birch murghnahan EOS and Fe minerals to use both new EOS's

@@ -20,6 +21,8 @@ def create(method):
if isinstance(method, basestring):
if method == "slb2":
return slb.SLB2()
elif method == "V":
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"V" to "Vinet"

@tjhei
Copy link
Member

tjhei commented Dec 11, 2015

please rebase and add unittests. /run-tests

'Kprime_0': 6.08,
'molar_mass': 0.055845,
'n': 1,} #number of atoms per formula unit
Mineral.__init__(self)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Github is complaining about no newline here?

@tjhei
Copy link
Member

tjhei commented Dec 13, 2015

can you take over tjhei@8f76d18 ?

This currently does not include shear modulus
"""
def __init__(self):
self.order=4
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please add a newline after this line

@sannecottaar
Copy link
Contributor

Do we want to finish this for the release?

@CaymanUnterborn
Copy link
Contributor Author

I think there may be a rebase issue although I'm not sure exactly what the "Merge remote-tracking branch 'geodynamics/master' into Planets"

@tjhei
Copy link
Member

tjhei commented Jan 5, 2016

We'll need to send you through another git bootcamp, I guess. :-) Don't do any merge/pulls unless you really want to.

@tjhei
Copy link
Member

tjhei commented Jan 5, 2016

What is the status of this, Cayman?

@CaymanUnterborn
Copy link
Contributor Author

I need to split BM into two separate files and then there's nothing else to do really. Can you explain the errors the tester is spitting out Timo so I can fix them? I still can't get anything but test.py to run on my machine.

@tjhei
Copy link
Member

tjhei commented Jan 5, 2016

please take over tjhei@fdcdeb2 then

@tjhei
Copy link
Member

tjhei commented Jan 5, 2016

(or I'll do that when it is good to merge). Anybody else have some comments?

@CaymanUnterborn
Copy link
Contributor Author

Alright, BM is split into 4th. Everything should be okay to go, tests.py passes.

raise ValueError('Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?')
return opt.brentq(func, sol[0], sol[1])

def volume_fourth_order(pressure,params):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you experimented with how fussy brentq is for this EoS? My intuition is that it would still blow up fairly easily (as with BM3).

I wonder if using the new bracket function, as above, would be more robust.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My experiments show it's also pretty fussy, particularly with brentq. I'll switch to bracket though.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using bracket should be less fussy (unless I've totally screwed that up). It would be neat to see an example -- perhaps a test?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm imagining the example will be Planet builder. I can write a test too, once I think of one...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test would not need to be complicated, just something that instantiates
a mineral and evaluates the EoS. Extra points if it is near the edge of the
range of validity for the EoS.

On Tue, Jan 5, 2016 at 11:08 AM, Cayman Unterborn notifications@github.com
wrote:

In burnman/eos/birch_murnaghan_4th.py
#158 (comment):

  • return K

+def volume(pressure, params):

  • """
  • Get the birch-murnaghan volume at a reference temperature for a given
  • pressure :math:[Pa]. Returns molar volume in :math:[m^3]
  • """
  • func = lambda x: birch_murnaghan(params['V_0']/x, params) - pressure
  • try:
  •    sol = bracket(func, params['V_0'], 1.e-2*params['V_0'])
    
  • except:
  •    raise ValueError('Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?')
    
  • return opt.brentq(func, sol[0], sol[1])

+def volume_fourth_order(pressure,params):

I'm imagining the example will be Planet builder. I can write a test too,
once I think of one...


Reply to this email directly or view it on GitHub
https://github.com/geodynamics/burnman/pull/158/files#r48882813.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't we have a test that does this already? I can write a "high_P" test that will do this although I'm not sure what that would tell us

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose test_minerals.py does to that, though it does not call all the
relevant functions in the EoS.

On Tue, Jan 5, 2016 at 11:28 AM, Cayman Unterborn notifications@github.com
wrote:

In burnman/eos/birch_murnaghan_4th.py
#158 (comment):

  • return K

+def volume(pressure, params):

  • """
  • Get the birch-murnaghan volume at a reference temperature for a given
  • pressure :math:[Pa]. Returns molar volume in :math:[m^3]
  • """
  • func = lambda x: birch_murnaghan(params['V_0']/x, params) - pressure
  • try:
  •    sol = bracket(func, params['V_0'], 1.e-2*params['V_0'])
    
  • except:
  •    raise ValueError('Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?')
    
  • return opt.brentq(func, sol[0], sol[1])

+def volume_fourth_order(pressure,params):

Don't we have a test that does this already? I can write a "high_P" test
that will do this although I'm not sure what that would tell us


Reply to this email directly or view it on GitHub
https://github.com/geodynamics/burnman/pull/158/files#r48885256.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've tried adding things to test_minerals to make it work, however it seems the fluids don't have molar masses and thus won't work if I call say their V_phi's (to test for K) or rho's. Do we not have a way to call just bulk modulus?

@CaymanUnterborn
Copy link
Contributor Author

Any word on this? My paper got accepted today so I'd like to include the planet builder in the release if possible.

Returns shear modulus :math:`G` of the mineral. :math:`[Pa]`
Currently not included in the Vinet EOS, so omitted.
"""
return 0
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you make this a floating point number?

@ian-r-rose
Copy link
Contributor

I think it is close. Can you figure out what is going on with the failing test?

@tjhei
Copy link
Member

tjhei commented Jan 7, 2016

This is just the table.py output, I can fix that if you accept the PR otherwise.

'molar_mass': 0.055845,
}
Mineral.__init__(self)
class Fe_Dewaele(Mineral):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please add two empty lines between classes

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is all ready to go. Might need a rebase though

@CaymanUnterborn
Copy link
Contributor Author

What's the word? Timo said he can fix the tests but is there anything else?

@bobmyhill
Copy link
Member

There still aren't any unit tests, which would be very useful. Also, a couple of P, T points to endmember_benchmarks.py to check for implementation changes in the future (as I've done for HP, HHPH and SLB). Finally, I think that a benchmark for each of BM4 and Vinet would be good (some figure from a paper that we can use to show the EoS has been implemented correctly). We've caught problems in both HP and SLB from plotting up diagrams from the original papers.

@CaymanUnterborn
Copy link
Contributor Author

What tests I could figure out have been added, feel free to suggest more. I've also added the benchmarking for both BM4 and Vinet.

(your choice in geotherm will not matter in this case)
or 'vinet' (vinet equation of state, if you choose to ignore temperature
(your choice in geotherm will not matter in this case)))"""
method = 'vinet'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason you switch to "vinet" here? I assume slb3 is a better default?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with Timo. slb3 should be the default - but of course vinet should be added to the list.
vinet should also be added to example_compare_all_methods.py.
I don't think there's any need to add another example... this is fine for me.

@tjhei
Copy link
Member

tjhei commented Feb 24, 2016

@bobmyhill can you take another look?

@tjhei
Copy link
Member

tjhei commented Feb 24, 2016

(and we need a rebase or merge from master before we can merge this, Cayman)

plt.title("Comparing with Figure 1 of Dewaele et al., (2006)")

plt.show()

def check_mgd_shim_duffy_kenichi():
"""
Attemmpts to recreate Shim Duffy Kenichi (2002)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(sp.)

@bobmyhill
Copy link
Member

  • I didn't think about it before, but birch_murnaghan_4th.py is not a very elegant way to deal with the fourth order EoS. Can you merge it into birch_murnaghan.py, making a new subclass (BM4)? The new functions look fine, so it wouldn't take very long.
  • Vinet looks fine.
    ... comments on examples/tests to follow ...

@@ -71,6 +97,41 @@ def test_reference_values(self):
Grun_test = i.grueneisen_parameter(pressure, temperature, rock.params['V_0'], rock.params)
self.assertFloatEqual(Grun_test, rock.params['grueneisen_0'])

def test_reference_values_noG(self):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why the different temperatures (300., 0.) in this test? The temperature shouldn't matter, surely.
If you have a reference pressure as 0. (rather than 1.e5) Pa, then this should be in the params as P_0, and you should test that the Vinet EoS handles different reference pressures in the same way that MT and BM do.

IMO, this should be made into at least two tests; one for Vinet, and one for BM4.

Another sensible (and extremely useful) test for BM4 would be to set K'' to the value inferred by BM3, and make sure that the results are the same at some high pressure. Same goes for BM3 wrt. BM2, if we don't have that test already.

@CaymanUnterborn
Copy link
Contributor Author

They were originally entirely in BM, however because BM2 and BM3, however the naming convention there doesn't match what BM4 really is. Ian suggested to just split it into two different EOS's. Is that right Ian?

@bobmyhill
Copy link
Member

What do you mean? BM4 is just an extension of BM3, just as BM3 is to BM2.

If you wanted to be more consistent in terms of expressions, you could rewrite the fourth order expansions as functions of f (as they are written here: http://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-99-3356).

@CaymanUnterborn
Copy link
Contributor Author

BM adopts a 3rd order P-V-K and a second or third order G, BM4 uses a 4th order P-V-K. I believe it was that reason that we split them.

@CaymanUnterborn
Copy link
Contributor Author

And thanks for those equations, I'll rewrite. I like that formulation much better.

@bobmyhill
Copy link
Member

Oh, yes, I remember. Hmmm, frustrating. I've opened a new issue (#234), because I still think BM should have a single home.

@tjhei
Copy link
Member

tjhei commented Feb 26, 2016

Can we maybe get this merged as is and defer the cleanup to after the release?

@CaymanUnterborn
Copy link
Contributor Author

Can you run me through exactly what I need to enter to rebase/squash correctly Timo?

@tjhei
Copy link
Member

tjhei commented Feb 26, 2016

Looking through the history I can see that you merged from master at least once. Rebasing/Squashing is really awkward in this case. Luckily, there is an easy way to squash everything into a single commit:

  1. I assume you call the official remote repo "upstream":

    $ git remote -v
    mygithub    git@github.com:tjhei/burnman.git (fetch)
    mygithub    git@github.com:tjhei/burnman.git (push)
    upstream    git@github.com:geodynamics/burnman.git (fetch)
    upstream    git@github.com:geodynamics/burnman.git (push)
    
  2. Make a backup:

    git checkout Planets
    git branch vinet_backup_before_merge
    
  3. Get up to date:

    git remote update
    git merge upstream/master
    # resolve conflict in misc/ref/table.py.out by throwing away your version:
    git checkout misc/ref/table.py.out --theirs
    git commit -a
    
  4. Create another backup:

    git branch vinet_backup_merged_before_squash
    
  5. Rebase and squash in one step:

    git reset --hard upstream/master
    git merge --squash vinet_backup_merged_before_squash
    git commit     # (make a nice commit message)
    
  6. Test, make sure it looks okay.

  7. Done, upload:

    git push -f mygithub Planets
    

why the hell is your branch called "Planets"?

@CaymanUnterborn
Copy link
Contributor Author

Because Planets are what I do Timo, the Earth is so passé these days :P

This is also where Planet builder was built so that's why it has that name.

I'll do this once I get back home this evening.

commit 2454445d3bfb8f08d4e97a02217b2dce82436c9d
Merge: 74be206 cb7e671
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Feb 27 00:04:40 2016 -0500

    Merge remote-tracking branch 'upstream/master' into Planets

commit 74be206
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Thu Feb 25 16:13:43 2016 -0500

    Cleaned up tests, more to do

commit b47c22f
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Wed Feb 24 14:29:51 2016 -0500

    Fixed default method in example. Reformatted BM4 to f notation

commit cb8e698
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Wed Feb 24 11:08:25 2016 -0500

    Fixed comments

commit a83a921
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Mon Feb 22 19:43:16 2016 -0500

    Added tests and benchmarks for Vinet and BM4

commit e151d90
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Fri Feb 12 12:21:40 2016 -0500

    Added vinet to example_user_input

commit 8e0a7df
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Thu Jan 7 15:04:30 2016 -0500

    Cosmetic fixes

commit 167e261
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 14:14:54 2016 -0500

    Added K back in.

commit 5c515f7
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 14:11:09 2016 -0500

    Redundant Code removal

commit dded24a
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 14:06:41 2016 -0500

    Removed K, bracketed V

commit be7f28b
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 13:59:56 2016 -0500

    Cleanup

commit cfcd96c
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 13:46:29 2016 -0500

    Split BM into 4th order

commit 2c96d8a
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 05:22:09 2016 -0500

    Fixed decmials

commit 3359fe0
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 05:18:47 2016 -0500

    Removed unnecessary Fe mineral

commit ec8c5d5
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 05:16:57 2016 -0500

    Removing duplicate mineral

commit ca412fa
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 04:48:48 2016 -0500

    Cleaned up failed test

commit bea8270
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 04:44:34 2016 -0500

    Added a liquid Fe EOS for BM4 Calculations

commit 7ba35ed
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 04:41:41 2016 -0500

    Fixed decimals

commit d579344
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 04:35:31 2016 -0500

    Added New line Timo Requested

commit b0af75f
Merge: a8325da d5d3d03
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Tue Jan 5 04:24:34 2016 -0500

    Merge remote-tracking branch 'geodynamics/master' into Planets

commit a8325da
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sun Dec 13 14:09:58 2015 -0800

    Fixed uppercase/lowercase issue

commit 9de0593
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Dec 12 23:02:23 2015 -0800

    Fixed typo causing test failure

commit cf492ab
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Dec 12 22:27:12 2015 -0800

    Lowercase V in vines and updated test table

commit 200f3b6
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Dec 12 22:14:58 2015 -0800

    Revert "Split Bm4 into it’s own EOS and reverted BM.py"

    This reverts commit 29a8b20.

commit 018d5cb
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Dec 12 19:18:30 2015 -0800

    Split Bm4 into it’s own EOS and reverted BM.py

commit bb7f6d8
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Dec 12 18:15:27 2015 -0800

    Fixed no new line issue

commit fa7efd5
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Dec 12 18:06:08 2015 -0800

    Fixed an import

commit 3e830f3
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Sat Dec 12 18:02:53 2015 -0800

    Included the 4th order Birch murnaghan EOS. Added 2 Epsilon Fe minerals
    in minerals.other and a liquid Fe EOS to use the 4th order BM EOS

commit d274159
Author: Cayman Unterborn <kmanunterborn@gmail.com>
Date:   Thu Dec 10 18:54:27 2015 -0800

    Added the Vinet EOS

    Fixed if to elif
@CaymanUnterborn
Copy link
Contributor Author

Should be squashed/rebased.

@tjhei
Copy link
Member

tjhei commented Feb 27, 2016

@bobmyhill @sannecottaar is this okay for the release as is? We can fix tests/design at some point.

I will fix the test before merging of course.

@tjhei
Copy link
Member

tjhei commented Mar 1, 2016

@bobmyhill @sannecottaar ?

@sannecottaar
Copy link
Contributor

I agree with Bob that BM should be organized differently, but otherwise I
think this looks good for now.
Thanks, Cayman.

On Tue, Mar 1, 2016 at 5:36 PM, Timo Heister notifications@github.com
wrote:

@bobmyhill https://github.com/bobmyhill @sannecottaar
https://github.com/sannecottaar ?


Reply to this email directly or view it on GitHub
#158 (comment).

@bobmyhill
Copy link
Member

Oops, sorry, meant to comment and didn't. Yes, this is okay for the release. There's nothing wrong with the code; it just needs some tidying. I'm not a big fan of including the Anderson and Ahrens as a test case, as their EoS is an adiabatic one starting from the melting point at 100 kPa, not an isothermal one starting at 0 or 300 K. But if it recreates the figure in the paper (Figure B1), then it's ok.

tjhei added a commit that referenced this pull request Mar 2, 2016
Add Vinet and BM4
@tjhei
Copy link
Member

tjhei commented Mar 15, 2016

I guess this can be closed. @lordkman are you going to reorganize things for after the release?

@tjhei tjhei closed this Mar 15, 2016
@CaymanUnterborn
Copy link
Contributor Author

I will indeed. Already preparing some new things to add and I'll just combine that into it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants