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

rpmodels #79

Merged
merged 12 commits into from
Mar 20, 2019
Merged

rpmodels #79

merged 12 commits into from
Mar 20, 2019

Conversation

aadm
Copy link
Contributor

@aadm aadm commented Mar 7, 2019

A bunch of rock physics models, all written from scratch from invididual papers and/or books (references in each function). All tested informally during daily job, will provide a jupyter notebook where I have recreated a few figures from cited references.

@kwinkunks
Copy link
Member

This looks tremendous, thank you Alessandro. I will make time for this and Matteo's PR asap.

Exciting to see bruges grow!

Copy link
Collaborator

@EvanBianco EvanBianco left a comment

Choose a reason for hiding this comment

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

I wonder if all these functions should go into rockphysics.elastic.py?

Alternately, perhaps these functions should go into a new submodule with the name isotropy? That might pair nicely with the one we already have called anisotropy

The last function concerns shale, but AFAICT the equation still assumes the rock is isotropic.

bruges/rockphysics/rpmodels.py Outdated Show resolved Hide resolved
return vp, vs, rho/1e3, K/1e9


def critpor(K0, G0, phi, phi_c=0.4):
Copy link
Collaborator

Choose a reason for hiding this comment

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

Should this function be renamed to better describe what it does? More over should it be pulled apart and written as two functions? It seems like this function does two things (returns Kdry and Gdry) and those two are independent.
e.g.
def calculate_criticalpor_dry_bulk_modulus(K0, phi, phi_c=0.4): returns Kdry

def calculate_criticalpor_dry_shear_modulus(G0, phi, phi_c=0.4): returns Gdry

Copy link
Contributor Author

Choose a reason for hiding this comment

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

All these functions behave in a similar way, i.e. they return K_dry and G_dry (bulk and shear elastic moduli). These two things always go together I suppose, because you then use them togetjher with density to reconstruct velocities. In the case of critical porosity it is true, you could split the function in two bu in other cases they are interlinked; so I would keep it like this for consistency with other functions and also to avoid the proliferation of too many functions.

bruges/rockphysics/rpmodels.py Outdated Show resolved Hide resolved
bruges/rockphysics/rpmodels.py Outdated Show resolved Hide resolved
@aadm
Copy link
Contributor Author

aadm commented Mar 13, 2019

This is a notebook that shows some test I made using the rock physics models included here. Essentially I have tried to reproduce some of the figures from the original sources' (papers and books) used to code the equations. It is not comprehensive though!

https://github.com/aadm/geophysical_notes/blob/master/bruges_rpmodels.ipynb

If I change the names of the functions as Evan suggested, I'll also modify this notebook. I also wonder where this notebook would go. Is there an official repo for extra documentation around bruges?

@aadm
Copy link
Contributor Author

aadm commented Mar 13, 2019

All the functions require in input non-SI units, i.e. elastic bulk and shear moduli (K, G) in GPa, pressure in MPa, velocities in m/s, density in g/cc etc. This may not be scientifically sound or "clean" but these units are what practitioners commonly use. Maybe we could check what units all the other functions in bruges require and try to find some sort of agreement at the library level.

@aadm
Copy link
Contributor Author

aadm commented Mar 13, 2019

The name of the various functions are now modified and less cryptic. So we have:

critical_porosity, hertz_mindlin, soft_sand, stiff_sand
contact_cement, constant_cement, increasing_cement
vernik_consol_sand, vernik_sand_diagenesis, vernik_shale
vernik_soft_sand_1, vernik_soft_sand_2, vels

I have decided not to add _model to the end of each of these functions to keep names a bit shorter.

@aadm
Copy link
Contributor Author

aadm commented Mar 13, 2019

For now, the file with all the new stuff is called rockphysicsmodels.py. I accept suggestions as to where to move these functions (e.g., elastic.py or isotropy.py, but honestly I dont' like either; perhaps a simple models.py so that we will have bruges.rockphysics.models?

@aadm
Copy link
Contributor Author

aadm commented Mar 13, 2019

Docstrings now have a brief description of each model; not 100% happy with these and they may need rewriting at some point.

@aadm
Copy link
Contributor Author

aadm commented Mar 14, 2019

Just to explain my latest commit, while moving my function vels to fluidsub.py as agreed with @EvanBianco I also tried your fluid substitution routines (I usually have my own routines which are different, but still based on Gassmann's equation obviously). I realized that the docstrings (as well as the names) for smith_gassmann and avseth_gassmann were misleading, at least to me. The first is simply put the direct form of Gassmann: you have as input dry-rock, mineral and fluid bulk modulus and porosity and you calculate the saturated rock bulk modulus. The other (avseth_fluidsub) is the indirect form where you bypass kdry altogether -- having as input the rock in one fluid state and knowing the elastic properties of initial and final fluid state you get the saturated rock bulk modulus with the final fluid. I have quickly rewritten the docstrings, but I believe the names should be changed as well.

@EvanBianco
Copy link
Collaborator

EvanBianco commented Mar 19, 2019

Thanks for all the changes here @aadm, things are looking great. I have suggested a bit of re-formatting some of the mathy bits for readability, but everything else is looking good as far as I can tell.

@kwinkunks, do you think we should merge these additions into a so-called rock-physics-model develop branch, first, and write some tests for it, before merging it to master?

@aadm , it's a bit tricky to verify the equations are correct, but ultimately we'll want to work towards writing tests for these functions.

Also, I think it might be helpful if we include an example section in each function's documentation, using the doctest format.

I'll defer to @kwinkunks to advise if he wants to have those examples in place before we pull the code, or if we can work on that collaboratively. @aadm, if you have some examples – of verified values that have been sanity checked, then we could add those in pretty quickly.

Comments welcome.

@aadm
Copy link
Contributor Author

aadm commented Mar 20, 2019

Thanks @EvanBianco I see the doctest format, and it would be easy to include these. The problem is finding these "verified values". I have this notebook where I have tested most of the functions as much as I can, and I only have ONE numerical examples, the others are all reproduction of published plots, so the accuracy is what it is:

https://github.com/aadm/geophysical_notes/blob/master/bruges_rpmodels.ipynb

All the years working on these bits of code and re-reading old papers, I still don't have this database of published values that I can compare the results to. What I used to do is to run the functions included in commercial packages and all the numbers match (up to a certain float), but this is obviously not a rigorous way of testing because who knows if the equations included in said packages are "correct"?

@kwinkunks kwinkunks changed the base branch from master to develop March 20, 2019 16:58
@kwinkunks kwinkunks merged commit eb49693 into agilescientific:develop Mar 20, 2019
@kwinkunks
Copy link
Member

Thank you Alessandro and Evan!

@EvanBianco
Copy link
Collaborator

@aadm, "All the years working on these bits of code and re-reading old papers, I still don't have this database of published values that I can compare the results to"

Sadly, for the state of our science, this is far too common in our industry than it should be. This complete lack of appreciation for reproducibility has become a big turn-off for me in the rock physics community, and I totally feel your pain.

This was exactly the motivation for doing the Repro Zoo, and why we all need to soldier on towards more useful technical publication systems.

Thanks again, for your efforts on adding these functions to bruges. Great contribution!

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.

None yet

5 participants