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

An Object Oriented interface for seismic wave finite-differences simulations #137

Closed
wants to merge 26 commits into
base: master
from

Conversation

Projects
None yet
2 participants
@leouieda
Member

leouieda commented Nov 8, 2014

Implement "simulator" objects to run the finite-differences time steppers. Re-uses the Cython time-steps and wraps them in classes, instead of the old generator functions (using yield statements).

Based on this prototype.

Fixes #70

Storage

The simulation classes store each time frame in an HDF5 file. When you need the simulation results, the class digs into the cache file and returns the data requested. Animation methods load the snapshots from the file and make animations using matplotlib.animate.

The cache file stores all data needed to resume a simulation later. The from_cache method rebuilds a simulation object from the cache file and computations can be resumed right away, Source functions are also stored in the cache file as pickles.

Rich display

Simulation and source wavelet classes have rich display capabilities for use with IPython notebooks. The snapshot method loads a time step from the cache file and plots it. The _repr_png_ for the simulations is the last time step. There is an option to convert the animation to a video and embed it in the IPython notebook.

Sample usage

A sample notebook using the code from this PR is on the tutorials repo. It will be updated along with the PR.

Testing

  • We need an analytic solution against which we can test the simulation results.
  • I would also like to have Python+numpy versions of the time stepping function so that we can test them against each other (and also use as a benchmark). It wouldn't hurt having a numba implementation as well.

TODO

  • Allow the user to specify the boundary conditions. Right now they are fixed as free-surface on top and absorbing on the other sides. (Can/should be left for a future PR)
  • Get rid of redundancy in plotting functions. There is a lot of repeated code.
  • Clean-up the way things are stored in the cache file. My use of HDF5 is far from ideal.
  • Try to adapt the elastic P-SV time step so that I don't have to copy (and make C-contiguous) the arrays when initializing the simulation (in the beginning of run). (Can/should be left for a future PR)
  • Figure out what to do with old code. Should we delete it or only deprecate it (and stop using it in the cookbook). I think deleting wouldn't be too bad since not many people use it anyway. Plus, the simulation classes are just so much better.

Checklist

  • Make tests for new code
  • Create/update docstrings
  • Update requirements
  • Code follows PEP8 style conventions
  • Code and docs have been spellchecked
  • Changelog entry
  • Documentation builds properly
  • All tests pass
  • Can be merged
Copied over the simulation classes from prototypes
Still need PEP8, examples, tests, and figuring out what to do with the
old code.

leouieda added a commit to fatiando/tutorials that referenced this pull request Nov 8, 2014

Showing off the new simulation class from PR 137
Will update this as the work on the PR progresses.

fatiando/fatiando#137

@leouieda leouieda added the enhancement label Nov 8, 2014

@leouieda leouieda added this to the 0.4 milestone Nov 8, 2014

@leouieda

This comment has been minimized.

Member

leouieda commented Nov 8, 2014

@eusoubrasileiro I've finally managed to put this stuff in a PR. If you want to
try to do any of the tasks above, I'd love it. Don't worry too much though.

If you do want to work on this, create a new branch on your fork starting from
this one (wavefd-simulation-class). I can then pull the changes from your
branch and update the PR.

Please don't do much that is outside of the task list above.
Extra work (like setting a sampling rate, 3D, migration, etc)
can be made in a another PR once this is merged. This way we can get this done
as fast as possible and build upon it.

I really want to know your comments on this!

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Nov 10, 2014

@leouieda cool 👍

Sure! When I get home I will take a deep look on it. Point by point.

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Nov 10, 2014

It's just a start I will write more aftwards (lunch break)

We have the analytical solution for scalar waves in Alford 1974 Geophysics here 2nd and 4th order. I know it's just for acoustic constant density (scalar) but for making short tests It works for me wavefd.scalar. For you, I am not sure but you might use the divergence of your field of waves for P-SV, not sure have to test. For the SH case I will try to look here for some seis books to find examples, but since you are a smart fella while I write this you might have found the answer already...

Also based on this paper I think we should include some criterias to garantee that our waves have minimum numerical dispersion.

  1. minimize numerical dispersion based on cell size (dx, dy, dz). Alford makes a dispersion analysis based on phase and group velocities as function of points per wavelength. He finds that for 4th order space 2nd time scalar we should try always grid points per wavelength > 5. Well I believe that might be possible to find something like that for elastic waves too.
  2. avoid spatial alias (another paper here expression 3 for simple plane waves)

render).

I write more soon.

@leouieda

This comment has been minimized.

Member

leouieda commented Nov 10, 2014

Analytic solutions: Great! I'll have a look at the paper when I can.
You had already implemented this long ago hadn't you?
For SH waves, the equations are pretty much the same as the acoustic type.
I will check with Leandro DiBartolo (from ON) to see if he has any analytic
solutions handy.

Numerical dispersion: That would be great! Numerical dispersion is a
pain. But it can fit on another quick PR later on maybe? I don't want to
delay this one too much so that you can merge your migration etc soon as
well. Plus, I'll need this for the geophysics course next semester. That
is why I want to get the bare minimum of these classes running.

There is so much to do with this that we could get stuck working on this
PR for years. Maybe you could add the spacial aliasing check to the Issue
list as a feature request? Just so we don't forget it and have a record of
it (you can reference this PR in the issue so that it links to the
discussion).

Sorry for being so controlling. But we really need to focus and leave some
ideas for later. I've updated the TODO list to marked 2 tasks as "Should be
left for later".

Think that the sooner this gets merged, the sooner we can use to build more
cool stuff!

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Nov 10, 2014

@leouieda master you are right. Soon after I wrote it I though : Oh well all that we could do later. I agree completely and totally I am learning to be more focused.

I will be working on the subsequent tasks! great!

@leouieda leouieda modified the milestones: 0.4, 1.0 Nov 24, 2014

@leouieda leouieda referenced this pull request Nov 26, 2014

Open

Mission statement #145

@eusoubrasileiro eusoubrasileiro referenced this pull request Dec 14, 2014

Closed

Documenting and testing wavefd-simulation-class #137 #153

12 of 15 tasks complete

leouieda and others added some commits Dec 14, 2014

Fixing doc strings.
Some misunderstanding correct due ambiguous variable names.
wavefd-simulation fixing.
Docstrings umbiguous and method names.
Fixing doc strings and tests.
Need to make:

 1) Refactor scalar

 2) Refactor cookbook

 3) Fix doc strings (don't forget development guide)

4) PEP 8

5) Check docs

6) Tests coverage acceptance

New requirements installation:

h5py;
all IPython notebook dependencies;
and mencoder;
wavefd-simulation docs/refactor/cookbook PR.
Package docstring, classes docstrings, refactor Scalar.
Changelog, conf.py docs, developer guide, docstring.
 Refactoring Scalar ready (need testing and tests)
 conf.py mocking solving problems on make docs.

@leouieda leouieda referenced this pull request Jan 22, 2015

Closed

Ementa e plano de aulas #1

0 of 7 tasks complete
@leouieda

This comment has been minimized.

Member

leouieda commented Jul 29, 2015

It might be worth checking this out for the ASCII progress bar dask/dask#485

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Aug 8, 2015

Copying from here 153 and getting back to it now in master.

I don't want to be hard with python but it seams its documentation features for abstract/interfaces/derived class or methods lack functionality. For documenting methods that are going to be used by derived classes I am now copying on all classes with equal methods. Python doesn't check if your method is not complying with the mother method it s overriding.

I wrote few test cases after watching the software carpentry lessons.

I will push more stuff as they come. It takes time a long time to understand someone-else's code without or with very few comments!!! @leouieda luckily I am learning some things.

Copying from yours 137
I need to reference this most recent ipython notebook usage procedures and the old one

  • Make tests for new code (I've done some ... working one more and maybe some doctests)
  • Create/update doc-strings (working on ... very very tough not being my code! )
  • Update requirements (I think it's okay now)
  • Code follows PEP8 style conventions (working on ..)
  • Code and docs have been spell-checked (for the end...)
  • Changelog entry
  • Documentation builds properly (builds properly but need more...)
  • All tests pass (so far passed ... working on more!)
  • Can be merged

Progress updates:
Updated (20/01/2015):
Porting my scalar code and cookbook.
Soon I will try to refactor all old cookbook code to get rid of everything old.
Need to do:

  • Re-factor scalar + dumping region
  • Re-factor scalar cookbook (do not re-factor the others) implicit optimizing plotting and animation code.
  • add padding support again needed for dumping region (I'll leave it for future for sh and psv)
  • Fix doc strings (don't forget development guide)
  • Tests coverage acceptance
  • Update installation requirements : h5py; all IPython notebook dependencies (pip recommended for latest); animate requires libavcodec-extra-54 and libav-tools (In Ubuntu 14.04 LTS)

About documenting classes:

  • Base class abstract methods should not have docstrings.
@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Aug 8, 2015

@leouieda after a lot of work to understand your code I still could not understand everything. hahaha still missing documentation.

Ive written some more tests.
My old Scalar class is totally documented.
Sources are fully documented.

Well I have fixed at least (I guess) most of the docstrings.
But I haven't finished reviewing the documentation for the other classes.
I haven't fixed the cookbooks or get rid of all old code.

Your help needed

What I need from you is to take I look on the way I am going with the docs. I haven't finished but...
I suggest that the common methods of the simulation classes should be on wavefd docstring in two topics: Simulation Storage Access and Simulation (rich) display.
Doing so we would not document the methods on those topics on the body of the classes, leaving the classes cleaner. My worry is that the module wavefd is getting really huge and so its docs.

@leouieda

This comment has been minimized.

Member

leouieda commented Aug 8, 2015

@eusoubrasileiro I know, this code is very badly documented and a terrible hack. We should talk about this code and make a plan of action before anything.

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Aug 8, 2015

Wow you are fast. Have I done anything wrong? hehe I thought I was suppose to bring back things here instead of using my old branch. It's up to you master. Just trying to start it over

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Aug 8, 2015

The "terrible hack" expression is new to me hehe Ok I will stop here and wait some time so we can talk about it. In any case I can revert this commit if you want

@leouieda

This comment has been minimized.

Member

leouieda commented Aug 8, 2015

I don't even remember a lot of this code. You shouldn't worry too much about the docstring of the base classes. They are mostly not for users and so it's not a problem if they aren't on the website. They shouldn't be on the website.

I'm more worried about cleaning up the code and not breaking the functionality. This is why we need tests first. The tests should run some simple scenarios and a notebook could be a more complete test of the rich display part. That way we can change the code without worrying that we break something.

@leouieda

This comment has been minimized.

Member

leouieda commented Aug 8, 2015

No, don't need to revert. See if you can figure out the code but it would be easier to talk (can be skype, but not this weekend).

@leouieda

This comment has been minimized.

Member

leouieda commented Aug 8, 2015

Try to see if you can get some tests and a notebook running (like the one in the tutorials repo). You can bring the notebook over to this repository in a tutorials folder in the root directory.

@eusoubrasileiro eusoubrasileiro self-assigned this Aug 8, 2015

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Nov 14, 2015

BAD NEWS
somehow it seams that this branch from my side is somehow tainted. I cant standard but Travis CI doesn't like it, although its fully updated with the master.

New todo list:

  • fix old cookbook recipes that are broken (e. g. cookbook/seismic_wavefd_scalar.py)
  • issue that will not be solved here. The new 'sphinx_bootstrap_theme' really made all documentation really confusing and unclear to understand. I rather use the old one, but we need to discuss it someday...
  • add padding support again needed for dumping region (I'll leave it for future for sh, psv)
  • sim.animate(every=5, cutoff=1, blit=True, fps=20, dpi=50, embed=True) not working, bug :-( ??
    Seams to be my machine....
  • should we change the add_point_source method to use real x, z coordinates? Or should we use indexes directly? I suggest both options the first being x and z, indexes being optional parameters tuple. Including back the area parameter.
    After thinking and testing a lot: we should leave the x, z and area (x1, x2, z1, z2) out of the classes (leave it spatially unaware). That will avoid loosing one of the greatest benefits of using HDF5: index or array slicing over the simulation through the __getitem__(self, index): method

Another ideas from @leouieda :

It might be worth checking this out for the ASCII progress bar dask/dask#485

Implemented and working in python and ipython. The original implementation had a second thread to update the status bar automatically. It "theoretically" can work in ipython notebook according to this answer on stack-overflow. I've tried that in many ways but could not make it work properly . So the only way to work everywhere ipython, python and ipython notebook was to forget about the second thread and update the status bar manually. But I found another way of implementing it as a generator of any iterable. This implementation is much cleaner and simpler.

I've been thinking of replacing the hacked version of movie embedding that I made with moviepy. I'll move some of the functions that are independent of the >simulation to a utils module in fatiando.vis. I'm still investigating but just to let you know.
More here about making a movie from a numpy array.

I am starting to move the code already to vis.utils

Trying to add more tests and removing old code.
Animate func is not working, working on tutorial notebook as well.

eusoubrasileiro added some commits Dec 23, 2015

added ipython notebook tutorials and new progress bar
added vis.utils file for starting migrating plotting feature from wavefd class
@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented on fatiando/seismic/wavefd.py in c443401 Dec 27, 2015

IPython notebook is changing widgets I started to fix the references for the new version 4.0.1.

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Dec 30, 2015

Some thoughts on simulation classes refactoring (class design)

I see that all simulation classes have 3 categories of functionalities

  1. Simulation math (variables and algorithm code for the specific partial differential equation)

    e. g. methods __init__(), add_point_source, add_any_source, _time_step
    e. g. variables dx, dz, dt, it, density, shape, velocity, Mx, Mz, mu, svel, pvel, sources etc.

  2. HDF5 simulation recording (simulation serialization, persistence) and access

    e. g. methods create_tmp_cache, get_cache, init_panels, init_cache, cache_panels, expand_cache, __getitem__, run(it is totally done inside a hdf5 file)
    e. g. variables *h5py.File, hdf5 data sets, hdf5 attributes etc.

  3. Plotting, visualization and understanding results (including ipython notebook features)

    e. g. methods _repr_png_, explore, snapshot, animate
    e. g. variables ... you got the idea

I am thinking that developing following this hierarchy top 1 to 3 bottom would be much easier and also logical. Having a class full of methods and variables is not easy to maintain, understand or expand.

One initial step we are taking is to move visualization stuff to fatiando.vis. We could link our class methods to functions there and avoid duplicates making it more reusable.

Some other brain storm thoughts:

We could design every simulation class based on 3 classes. First simulation algorithm, having a subclass? of arguments package. Second simulation hdf5 persistence (derived from the first) using the arguments package to more easily write data data sets and attributes. Third visualization stuff (derived from the second). ?

Or use python 3, multiple inheritance ?
Simulation(Algorithm, Hdf5Simulation) and write in the end just the visualization part... ?

well those have been my thoughts until now...

@leouieda

This comment has been minimized.

Member

leouieda commented Jan 13, 2016

@eusoubrasileiro thanks so much for putting work into this! I haven't been able to give you the proper attention, sorry 💩

You mentioned you're having trouble with the branch and Travis CI. This PR is getting very confusing quite fast. Don't worry too much right now but we should consider starting a new branch for this before merging. We can just copy/paste the code from here to the new branch. That would make the history cleaner and eliminate all the merges from master that are causing all this trouble. But we can leave that for later.

Once I submit my thesis we should get together and discuss your design ideas!

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Jan 14, 2016

@leouieda don't worry man! just need to know you are alive heheh you are! hahah

Good work there! Focus on your thesis! cya later!

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Feb 11, 2016

@leouieda

This comment has been minimized.

Member

leouieda commented Feb 11, 2016

I remember seeing that a long time ago when first building this in Fatiando. I don't really like the API for that code. It seems very Fortranish. But we should have a better look at their source code. Maybe there is something in there to salvage. The project appears to be dead now, the last update was in 2013.

@leouieda

This comment has been minimized.

Member

leouieda commented Feb 11, 2016

I've been thinking about the volume of code here and I finally agree with you that this should be in separate files (sorry for taking so long to see it). I think the best way would be to make wavefd a package inside seismic with files to separate things neatly. But the user wouldn't need to see this. We can import the public facing classes in wavefd/__init__.py and the user can still do from fatiando.seismic.wavefd import ElasticSH for example.

However, this PR is already out of hand. I'll have to take a look at this in March because I'll be using it in my classes again. Hopefully we can organize things and finish this off then.

Thanks again for all the help with this!

@eusoubrasileiro

This comment has been minimized.

Contributor

eusoubrasileiro commented Feb 11, 2016

holaa

yes I think that project is dead as well... good to know you knew about it beforehand
I will look at it with the same spirit sometime...

hehehe about the wavefd package. No problema man! there we go!
GOOD! I like the idea! Exactly what I think we need!

we are growing and improving the package that's what matters! 👍

@leouieda

This comment has been minimized.

Member

leouieda commented Apr 28, 2016

I'm closing this Pull Request because it's been inactive for a long time. The branch is probably very out of sync with master.

If you want to continue working of this, please sync your fork, update your branch, and re-submit the PR. Let me know if you need any help doing this!

Sorry for the hassle! This is not a critique of your work, just a way to keep things fresh in the project. I hope you'll consider re-submitting.

This is a message template.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment