Finite difference functions

[mloubout]

@mlange05 , @navjotk , @ggorman There is a lot of new finite difference functions (not one liners) required for the adjoint TTI. You can see it in the new notebook. Do you think we want to keep the symbolic computation inside the operators, or should we create a finite_difference.py containing every schemes we could need that we just call to generate the stencils?

the notebook is TTI_3d_simplfied.ipynb

[mlange05]

Yes, we definitely want to provide finite difference related tools to generically do the expansions and a separate finite_difference.py sounds like a great idea. In fact I am currently working on exactly this in the automate-symbolic-substitution branch, were I am adding short-hand properties to the symbolic data objects that provide the right derivative expression to the user; something like:

u = TimeData(name='u', time_order=2, space_order=8, ...)
eqn = Eq(d.dt, a * (u.dxx + u.dyy))

Your new derivations look very nice, and I will work them in straight away. I assume your hand-calculated weights will provide the exact same result as the currently used as_finite_diff approach?

[mloubout]

The weights are the normal Taylor coefficient, I have a matlab function I will write in python that compute this coefficient for any order and will avoid to just define all this orders. THis will also allow to work with Non-Taylor coefficients

[mlange05]

Could sympy.calculus.finite_diff_weights not do the same trick, or is this somewhat limited?

[mlange05]

Also, I just noted that the cross-derivative functions (dxyz et al.) don't honour the order or the coefficients. Is this intended?

[mloubout]

sympy.calculus.finite_diff_weights will give the proper coefficient, but only Taylor ones, which are not necessarily the best choice. But this would be a good start. For the cross derivatives (and dxx, dyy,dzz) the order does not matter as you consider the product (commutative) as a single function and apply the finite difference to it. This is different than the "math" cross derivative that would require to develop the derivative first.

[mloubout]

Actualy even Gxx,.. and the transposes could be defined in the finite_difference par allowing to only write the discrete pde in the operator e.g for 2D TTI eqn = Eq(d.dt, a*(epsilon * u.Gxx + delta * u.Gzz) )

[mlange05]

Ok thanks, that makes sense. I just pushed this commit to my current symbolics branch, which generically derives second and cross derivatives as per your example. The implementations are slightly pythonified :snake: and modified to allow user-defined dimensions, but they still support products of functions via *args (very nice btw. :+1:). If you have some time, could you please double check that they behave as you would expect? You can get them via from devito import second_derivative, cross_derivative and use them like this:

x, y = symbols('x y')
f = Function('f')
g = Function('g')
dxx = second_derivative(f(x, y), g(x, y), order=2, dim=x)
dxy = cross_derivative(f(x, y), g(x, y), dims=(x, y))

And yes, if we can make it generic we should totally provide a utility for the Gxx terms.

[mloubout]

Looks really nice, I'll test it

[ggorman]

And add a feature test for a couple of simple analytical examples?

[mloubout]

This should be easy to include too.

[mloubout]

@ggorman I looked around, and from what people usually do, I think the simplest test we can implement is to have a trivial high order polynomial with easy analytical derivatives. This is what is currently in as_finite_diff and we could just use that. I will add that to the notebook with the function that generates the finite difference coefficient (I don't want to mess up the very pretty python class @mlange05 made)

[mloubout]

https://github.com/sympy/sympy/pull/11337/commits/372f89d03fa840ac156f410af663b863ea5f8262?diff=unified

sympy/calculus/tests/test_finite_diff.py

and further

[mloubout]

Do you think we should rewrite the Acoustic and current TTI with this finite differences to make everything consistent? If yes I can do it it shouldn't be too much work.