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

examples: Add an example notebook for ADER-FD schemes #2338

Merged
merged 10 commits into from
May 9, 2024
Merged

examples: Add an example notebook for ADER-FD schemes #2338

merged 10 commits into from
May 9, 2024

Conversation

EdCaunt
Copy link
Contributor

@EdCaunt EdCaunt commented Mar 27, 2024

Currently work in progress. Adds a notebook showing how to solve the first-order formulation of the acoustic wave equation with 4th-order ADER timestepping.

Lmk your thoughts.

Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

Copy link

codecov bot commented Mar 27, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 86.76%. Comparing base (fcfdfa3) to head (50be3f8).

Additional details and impacted files
@@           Coverage Diff           @@
##           master    #2338   +/-   ##
=======================================
  Coverage   86.76%   86.76%           
=======================================
  Files         233      233           
  Lines       43740    43740           
  Branches     8075     8075           
=======================================
  Hits        37952    37952           
  Misses       5079     5079           
  Partials      709      709           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@@ -0,0 +1,247 @@
{
Copy link
Contributor

@mloubout mloubout Mar 27, 2024

Choose a reason for hiding this comment

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

Intro and refs?


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah, still need to add these

@@ -0,0 +1,247 @@
{
Copy link
Contributor

@mloubout mloubout Mar 27, 2024

Choose a reason for hiding this comment

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

Line #10.    # dv.grad(dv.div(v)) does not get flattened and will need to be explicitly written out

Not sure what that means but does not make sense to have to go to those sympy matrices.


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Maybe I'm overthinking? I was under the impression that dv.grad(dv.div(v)) is going to result in bigger stencils than doing the grad of the divergence algebraically and implementing that

@@ -0,0 +1,247 @@
{
Copy link
Contributor

@mloubout mloubout Mar 27, 2024

Choose a reason for hiding this comment

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

Same as above.


Reply via ReviewNB

@@ -0,0 +1,247 @@
{
Copy link
Contributor

@mloubout mloubout Mar 27, 2024

Choose a reason for hiding this comment

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

Isn't that second order just the OT4 acoustic kernel but for the first order equations?


Reply via ReviewNB

Copy link
Contributor Author

@EdCaunt EdCaunt Mar 27, 2024

Choose a reason for hiding this comment

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

OT4?

Copy link
Contributor

Choose a reason for hiding this comment

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

biharmonic = field.biharmonic(1/model.m) if kernel == 'OT4' else 0

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Same idea I think

@mloubout mloubout added the examples examples label Mar 27, 2024
"\n",
"# First time derivatives\n",
"pdt = rho*c2*dv.div(v)\n",
"vdt = b*dv.grad(p)\n",
Copy link
Contributor

Choose a reason for hiding this comment

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

This is quite surprising, it's missing the time leap-frogging

p.forward = f(v, p))
v.forward=f(p.forward, v)

Copy link
Contributor Author

@EdCaunt EdCaunt Mar 27, 2024

Choose a reason for hiding this comment

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

Yeah, you don't need time leapfrogging for ADER. You also don't need any staggering to solve first-order systems of equations. This is part of the appeal imo on top of high-order timestepping

@rhodrin
Copy link
Contributor

rhodrin commented Apr 23, 2024

What's the reason this is currently located in misc vs the seismic tutorials?

Locating in the later and referring previous tutorials where appropriate would seem more suitable to me?

@@ -0,0 +1,327 @@
{
Copy link
Contributor

@rhodrin rhodrin Apr 23, 2024

Choose a reason for hiding this comment

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

I would recommend either defining the state variables or referencing previous tutorials appropriately. Unless this tutorial is only intended for readers who are familiar with seismic modelling, make it as simple for target audience as is reasonably possible - that is, make sure everything is clearly defined.

I'm not a big fan of the phrasing here:

"one can derive expressions for these higher-order time derivatives in terms of spatial derivatives. High-order..."

Just say what you're going to do in the tutorial and reference the paper sited above as appropriate.


Reply via ReviewNB

Copy link
Contributor

Choose a reason for hiding this comment

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

Also, U expansion expression is missing $\Delta t$ terms.

@@ -0,0 +1,327 @@
{
Copy link
Contributor

@rhodrin rhodrin Apr 23, 2024

Choose a reason for hiding this comment

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

Nitpicking: required imports comment wouldn't hurt.


Reply via ReviewNB

@@ -0,0 +1,327 @@
{
Copy link
Contributor

@rhodrin rhodrin Apr 23, 2024

Choose a reason for hiding this comment

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

No need for "as usual". Ref. the difference of this setup vs the one in the existing first order notebook.


Reply via ReviewNB

@@ -0,0 +1,327 @@
{
Copy link
Contributor

@rhodrin rhodrin Apr 23, 2024

Choose a reason for hiding this comment

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

State units or dimensionless etc.


Reply via ReviewNB

@@ -0,0 +1,327 @@
{
Copy link
Contributor

@rhodrin rhodrin Apr 23, 2024

Choose a reason for hiding this comment

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

Conventional case? Conventional to whom?

"Note that if one assumes non-constant material parameters when deriving the higher time derivatives, it will be necessary..." - reference or brief explanation.


Reply via ReviewNB

@@ -0,0 +1,327 @@
{
Copy link
Contributor

@rhodrin rhodrin Apr 23, 2024

Choose a reason for hiding this comment

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

Don't we define the biharmonic method in differentiable?

Not necessarily for this PR, but with a couple of appropriate sympy functions deriving expressions of arbitrary order could be 'automated'.


Reply via ReviewNB

@@ -0,0 +1,327 @@
{
Copy link
Contributor

@rhodrin rhodrin Apr 23, 2024

Choose a reason for hiding this comment

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

Add labels etc.

Would be nice to add a comparison vs staggered 1st order acoustic (with a difference plot).


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This was useful, since it turns out this highlights the higher stability limit of ADER vs staggered

@rhodrin
Copy link
Contributor

rhodrin commented Apr 23, 2024

Still need to check the mathematical details in some places.

But in general, throughout the tutorial recommend expanding and clarifying statements a little more and adding a comparison vs staggered acoustic.

@@ -0,0 +1,565 @@
{
Copy link
Contributor

@mloubout mloubout May 7, 2024

Choose a reason for hiding this comment

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

Could you add the math for say 2nd order ADER?


Reply via ReviewNB

@@ -0,0 +1,565 @@
{
Copy link
Contributor

@mloubout mloubout May 7, 2024

Choose a reason for hiding this comment

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

Line #9.    

I still don't quite get why using grad(div) and such doesn't work.


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Because Devito derivatives are discrete, so grad(div()) is the gradient stencil applied to a divergence stencil. This will contain f[0].dx.dx which is not the same stencil as f[0].dx2

@@ -0,0 +1,565 @@
{
Copy link
Contributor

@mloubout mloubout May 7, 2024

Choose a reason for hiding this comment

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

Line #7.    # Update equations (3rd-order ADER timestepping)

I would prefer writing the PDE and use solve rather than by hand. It would also make the method more explicit


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That would obfuscate how the timestepping works imo

@@ -0,0 +1,565 @@
{
Copy link
Contributor

@mloubout mloubout May 7, 2024

Choose a reason for hiding this comment

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

Isn't taht worrying that the dirichlet boundary isn't reflecting?


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

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

No. If you imposed all the higher-order conditions to get a free surface, you would expect it to reflect, but the Dirichlet boundary condition when you hit the zero padding isn't really a free surface or anything physically meaningful. There are reflections back into the domain in the ADER case, they're just much lower amplitude. Due to the properties of the discretisation and the fact that the padding is held at zero, energy is diffused at the edges of the domain, hence the lack of reflections

@mloubout
Copy link
Contributor

mloubout commented May 7, 2024

Mostly minor comments. Neess reabase as well.

@mloubout mloubout merged commit 19d4e33 into master May 9, 2024
31 checks passed
@mloubout mloubout deleted the ader branch May 9, 2024 14:38
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
examples examples
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants