api: add support for 45 degree FD approx #2326

mloubout · 2024-03-05T03:27:39Z

Add suppport for 45 degree derivatives (usually called RSFD). This only supports cases where the diagonal (45 degree axis) aligns with points at which the fields is defined.

review-notebook-app · 2024-03-05T03:27:45Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

codecov · 2024-03-05T03:34:40Z

Codecov Report

Attention: Patch coverage is 92.85714% with 16 lines in your changes are missing coverage. Please review.

Project coverage is 86.67%. Comparing base (f29cb35) to head (47be3e1).

Files	Patch %	Lines
devito/finite_differences/tools.py	69.23%	7 Missing and 1 partial ⚠️
devito/finite_differences/operators.py	71.42%	4 Missing ⚠️
devito/finite_differences/differentiable.py	81.81%	1 Missing and 1 partial ⚠️
devito/finite_differences/rsfd.py	96.22%	1 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2326      +/-   ##
==========================================
+ Coverage   86.66%   86.67%   +0.01%     
==========================================
  Files         231      232       +1     
  Lines       43260    43409     +149     
  Branches     8022     8057      +35     
==========================================
+ Hits        37490    37626     +136     
- Misses       5068     5079      +11     
- Partials      702      704       +2

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

devito/finite_differences/derivative.py

FabioLuporini · 2024-03-05T03:37:15Z

devito/finite_differences/derivative.py

@@ -210,7 +215,7 @@ def _new_from_self(self, **kwargs):
        expr = kwargs.pop('expr', self.expr)
        _kwargs = {'deriv_order': self.deriv_order, 'fd_order': self.fd_order,
                   'side': self.side, 'transpose': self.transpose, 'subs': self._ppsubs,
-                   'x0': self.x0, 'preprocessed': True}
+                   'x0': self.x0, 'preprocessed': True, 'method': self.method}
        _kwargs.update(**kwargs)
        return Derivative(expr, *self.dims, **_kwargs)


one day... we should just make use of .rebuild(...). Perhaps worth adding an issue

memo (add issue?)

devito/finite_differences/derivative.py

devito/types/dense.py

devito/finite_differences/rsfd.py

devito/finite_differences/derivative.py

EdCaunt · 2024-03-06T15:16:52Z

devito/finite_differences/differentiable.py

@@ -318,16 +329,17 @@ def div(self, shift=None, order=None):
        order: int, optional, default=None
            Discretization order for the finite differences.
            Uses `func.space_order` when not specified
-
+        method: str, optional, default='FD'
+            Discretization method. Options are 'FD' (default) and 'RSFD'


Would be worth adding a '(rotated-staggered finite-difference)' to the docstrings where this option is exposed to the user

also (nitpicking) full stop at the end (lacked by several entries)

EdCaunt · 2024-03-06T16:09:54Z

devito/finite_differences/finite_difference.py

@@ -114,8 +117,8 @@ def cross_derivative(expr, dims, fd_order, deriv_order, x0=None, **kwargs):
        finite difference. Defaults to `direct`.
    x0 : dict, optional
        Origin of the finite-difference scheme as a map dim: origin_dim.
-    symbolic : bool, optional
-        Use default or custom coefficients (weights). Defaults to False.
+    coefficients : strong, optional


Typo: string

EdCaunt · 2024-03-06T16:10:22Z

devito/finite_differences/finite_difference.py

@@ -186,8 +189,8 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
        finite difference. Defaults to `direct`.
    x0 : dict, optional
        Origin of the finite-difference scheme as a map dim: origin_dim.
-    symbolic : bool, optional
-        Use default or custom coefficients (weights). Defaults to False.
+    coefficients : strong, optional


EdCaunt · 2024-03-06T16:11:13Z

devito/finite_differences/finite_difference.py

@@ -200,30 +203,33 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
    side = None
    # First order derivative with 2nd order FD is highly non-recommended so taking


If you're changing this file, then you might as well change this line to "strongly discouraged"

"strongly discouraged. Taking the..."

EdCaunt · 2024-03-06T16:16:47Z

devito/finite_differences/operators.py

+
+def div45(func, shift=None, order=None, method='FD'):
+    """
+    Divergence of the input Function, using the 45 degrees rotated derivative.


"Using 45 degree rotated derivatives"

EdCaunt · 2024-03-06T16:23:10Z

devito/finite_differences/rsfd.py

+
+    Rotated finite differences based on:
+    https://www.sciencedirect.com/science/article/pii/S0165212599000232
+    The rotated axis (the four diagonals of a cube) are:


EdCaunt · 2024-03-06T16:40:33Z

devito/finite_differences/rsfd.py

+        Expand the expression. Defaults to True.
+    """
+    # Make sure the grid supports RSFD
+    if expr.grid.dim == 1:


if expr.grid.dim not in (2, 3): would be safer

devito/finite_differences/tools.py

EdCaunt · 2024-03-06T16:50:07Z

tests/test_derivatives.py

+    @pytest.mark.parametrize('staggered', [(True, True), (False, False),
+                                           (True, False), (False, True)])
+    @pytest.mark.parametrize('space_order', [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
+    def test_fd_space_45(self, staggered, space_order):


We should test this in 3D as well to make sure everything is implemented as intended

FabioLuporini · 2024-03-11T08:49:31Z

devito/finite_differences/derivative.py

@@ -328,6 +337,8 @@ def _eval_at(self, func):
        """
        # If an x0 already exists do not overwrite it
        x0 = self.x0 or dict(func.indices_ref._getters)
+        # expr_x0 = self.expr.indices_ref._getters


FabioLuporini · 2024-03-11T08:49:40Z

devito/finite_differences/derivative.py

@@ -210,7 +215,7 @@ def _new_from_self(self, **kwargs):
        expr = kwargs.pop('expr', self.expr)
        _kwargs = {'deriv_order': self.deriv_order, 'fd_order': self.fd_order,
                   'side': self.side, 'transpose': self.transpose, 'subs': self._ppsubs,
-                   'x0': self.x0, 'preprocessed': True}
+                   'x0': self.x0, 'preprocessed': True, 'method': self.method}
        _kwargs.update(**kwargs)
        return Derivative(expr, *self.dims, **_kwargs)


memo (add issue?)

FabioLuporini · 2024-03-11T08:51:58Z

devito/finite_differences/differentiable.py

@@ -318,16 +329,17 @@ def div(self, shift=None, order=None):
        order: int, optional, default=None
            Discretization order for the finite differences.
            Uses `func.space_order` when not specified
-
+        method: str, optional, default='FD'
+            Discretization method. Options are 'FD' (default) and 'RSFD'


also (nitpicking) full stop at the end (lacked by several entries)

FabioLuporini · 2024-03-11T08:52:37Z

devito/finite_differences/differentiable.py

@@ -659,7 +672,7 @@ def _evaluate(self, **kwargs):
        expr = self.expr._evaluate(**kwargs)

        if not kwargs.get('expand', True):
-            return self.func(expr, self.dimensions)
+            return self._rebuild(expr)


(nitpicking) func is fine, it defaults to _rebuild underneath, maybe we should just use func wherever we can

func is not fine it leads to some weird error I didn't have time to look into why I switched back to rebuild

FabioLuporini · 2024-03-11T08:53:32Z

devito/finite_differences/finite_difference.py

-    symbolic : bool, optional
-        Use default or custom coefficients (weights). Defaults to False.
+    coefficients : string, optional
+        Use taylor or custom coefficients (weights). Defaults to taylor.


(nitpicking)

As pointed out in Zoe's latest PR, the "Defaults ..." goes in the line above. We should start changing them all, obv incrementally

FabioLuporini · 2024-03-11T08:55:59Z

devito/finite_differences/operators.py

@@ -55,12 +120,12 @@ def curl(func, shift=None, order=None):
        Uses `func.space_order` when not specified
    """
    try:


you could maybe just call curl(func, shift=shift, order=order, method='RSFD')?

FabioLuporini · 2024-03-11T08:57:33Z

devito/finite_differences/tools.py

@@ -360,11 +374,20 @@ def generate_indices_staggered(expr, dim, order, side=None, x0=None):
            indices = [start, start - diff]
            indices = IndexSet(dim, indices)
        else:
-            o_min = -order//2
-            o_max = order//2
+            if x0 is None or order % 2 == 0:


could use a comment down here

FabioLuporini · 2024-03-11T08:59:28Z

devito/types/dense.py

+        if self._coefficients not in fd_weights_registry:
+            if self._coefficients == 'standard':
+                self._coefficients = 'taylor'
+                warnings.warn("The `standard` mode is deprecated and will be removed in "


so just to be sure: nobody will ever see this deprecation warning unless they were explicitly doing Function(..., coefficients='standard', ...) right?

another option is to simply just add standard to fd_weights_registry and make it point to taylor without the need to deprecate and change any code here

FabioLuporini · 2024-03-11T09:01:41Z

examples/seismic/tti/operators.py

@@ -277,31 +275,31 @@ def kernel_staggered_2d(model, u, v, **kwargs):

    if forward:
        # Stencils
-        phdx = costheta * u.dx - sintheta * u.dy


why do we need the extra 'c' now?

is this affecting use side anyhow? looks like a critical change at first glance?

If you check the 3d below you'll see the c was supposed to be there but was luckily working in 2D without it. t's mostly because this is a twisted PDE with a missing field so need some tweaks for the FD>

FabioLuporini · 2024-03-11T09:01:51Z

examples/seismic/tti/wavesolver.py

-            dims = self.model.space_dimensions
-            stagg_u = (-dims[-1])
-            stagg_v = (-dims[0], -dims[1]) if self.model.grid.dim == 3 else (-dims[0])
+            stagg_u = stagg_v = NODE


lovely clean up

EdCaunt · 2024-03-11T14:56:12Z

devito/finite_differences/finite_difference.py

@@ -198,32 +201,34 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
        ``deriv-order`` derivative of ``expr``.
    """
    side = None
-    # First order derivative with 2nd order FD is highly non-recommended so taking
+    # First order derivative with 2nd order FD is strongly discouraged so taking


Comment still a little garbled. "strongly discouraged. Taking first order FD is recommended in this case."

EdCaunt · 2024-03-11T15:03:22Z

tests/test_derivatives.py

+        Rotated derivatives require at least 2D to get access to diagonal points.
+        We create a simple 1D gradient over a 2D grid to check against a polynomial
+        """
+        # dummy axis dimension


FabioLuporini

Thanks, GTG now

…lid staggering

mloubout added the API api (symbolics, types, ...) label Mar 5, 2024

mloubout force-pushed the stagg-fd-custom branch 2 times, most recently from 9d81997 to 3ed3e19 Compare March 5, 2024 03:29

FabioLuporini self-requested a review March 5, 2024 03:35

FabioLuporini reviewed Mar 5, 2024

View reviewed changes

mloubout force-pushed the stagg-fd-custom branch 4 times, most recently from 351d79d to ae85ddf Compare March 6, 2024 03:55

EdCaunt reviewed Mar 6, 2024

View reviewed changes

mloubout force-pushed the stagg-fd-custom branch 8 times, most recently from 5f4bc51 to a7f3be3 Compare March 10, 2024 03:28

FabioLuporini reviewed Mar 11, 2024

View reviewed changes

mloubout force-pushed the stagg-fd-custom branch from d45451d to dc4eaa7 Compare March 11, 2024 14:54

EdCaunt reviewed Mar 11, 2024

View reviewed changes

FabioLuporini approved these changes Mar 12, 2024

View reviewed changes

mloubout added 6 commits March 12, 2024 10:08

api: add rsfd and shortcuts

daf9721

api: refacto coefficients for easier extension

99da6b3

api: support RSFD

660df7b

symbolics: fix default rule for symbolic coefficients

3f2b8e1

examples: add rotated fd section to userapi notebook

8ee3438

api: fix EvalDerivative and expand arithmetic

9be31c9

mloubout added 3 commits March 12, 2024 10:08

api: fix staggered FD corner case and enforce satggered case to be va…

8ca4c51

…lid staggering

api: fix symbolic indices with stagger

04d6d11

misc: cleanup docstrings and typos

47be3e1

mloubout force-pushed the stagg-fd-custom branch from dc4eaa7 to 47be3e1 Compare March 12, 2024 14:09

mloubout merged commit 8f6b6a8 into master Mar 12, 2024
31 checks passed

mloubout deleted the stagg-fd-custom branch March 12, 2024 17:21

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

api: add support for 45 degree FD approx #2326

api: add support for 45 degree FD approx #2326

mloubout commented Mar 5, 2024

review-notebook-app bot commented Mar 5, 2024

codecov bot commented Mar 5, 2024 •

edited

Loading

FabioLuporini Mar 5, 2024

FabioLuporini Mar 11, 2024

mloubout Mar 12, 2024

EdCaunt Mar 6, 2024

FabioLuporini Mar 11, 2024

EdCaunt Mar 6, 2024

EdCaunt Mar 6, 2024

EdCaunt Mar 6, 2024

EdCaunt Mar 6, 2024

EdCaunt Mar 6, 2024

EdCaunt Mar 6, 2024

EdCaunt Mar 6, 2024

EdCaunt Mar 6, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

mloubout Mar 11, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

FabioLuporini Mar 11, 2024

mloubout Mar 11, 2024

FabioLuporini Mar 11, 2024

EdCaunt Mar 11, 2024

EdCaunt Mar 11, 2024

FabioLuporini left a comment

		@@ -200,30 +203,33 @@ def generic_derivative(expr, dim, fd_order, deriv_order, matvec=direct, x0=None,
		side = None
		# First order derivative with 2nd order FD is highly non-recommended so taking

api: add support for 45 degree FD approx #2326

api: add support for 45 degree FD approx #2326

Conversation

mloubout commented Mar 5, 2024

review-notebook-app bot commented Mar 5, 2024

codecov bot commented Mar 5, 2024 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FabioLuporini left a comment

Choose a reason for hiding this comment

codecov bot commented Mar 5, 2024 •

edited

Loading