from gradient to drivative

spyro/solvers/inversion.py

@@ -404,6 +404,7 @@ def __init__(self, dictionary=None, comm=None):
         self.has_gradient_mask = False
         self.gradient_mask_available = False
         self.functional_history = []
+        self._mass_diagonal = None
 
     @property
     def real_velocity_model_file(self):
@@ -811,6 +812,28 @@ def get_functional(self, c=None):
 
         return Jm
 
+    def _gradient_to_derivative(self, gradient):
+        """Convert the L2 Riesz gradient to the Euclidean derivative.
+
+        The implemented adjoint returns the gradient in the function-space
+        inner product. SciPy expects the coefficient-space derivative, so this
+        method applies the lumped mass diagonal element-wise.
+        """
+        if self._mass_diagonal is None:
+            V = self.function_space
+            u = fire.TrialFunction(V)
+            v = fire.TestFunction(V)
+            ones = fire.Function(V)
+            ones.assign(1.0)
+            mass_form = fire.inner(u, v) * fire.dx(**self.quadrature_rule)
+            self._mass_diagonal = fire.assemble(fire.action(mass_form, ones))
+
+        derivative = fire.Function(self.function_space)
+        derivative.dat.data[:] = (
+            gradient.dat.data_ro[:] * self._mass_diagonal.dat.data_ro[:]
+        )
+        return derivative
+
     def get_gradient(self, c=None, save=True, calculate_functional=True):
         """
         Calculate the gradient of the objective functional.
@@ -871,7 +894,8 @@ def return_functional_and_gradient(self, c):
             The gradient of the functional with respect to the velocity model.
         """
         self.get_gradient(c=c)
-        dJ = self.gradient.dat.data[:]
+        derivative = self._gradient_to_derivative(self.gradient)
+        dJ = derivative.dat.data_ro.copy()
         return self.functional, dJ
 
     def run_fwi(self, **kwargs):

tests/on_one_core/test_utils.py

@@ -1,7 +1,9 @@
+import firedrake as fire
 import spyro
 import scipy as sp
 import numpy as np
 import math
+from spyro.solvers.inversion import FullWaveformInversion
 
 
 def test_butter_lowpast_filter():
@@ -61,6 +63,29 @@ def test_geometry_creation():
     assert all([test0, test1, test2, test3, test4, test5, test6])
 
 
+def test_gradient_to_derivative_uses_lumped_mass_diagonal():
+    mesh = fire.UnitSquareMesh(1, 1)
+    V = fire.FunctionSpace(mesh, "CG", 1)
+
+    dummy = type("DummyFWI", (), {})()
+    dummy.function_space = V
+    dummy.quadrature_rule = {}
+    dummy._mass_diagonal = None
+
+    gradient = fire.Function(V)
+    gradient.dat.data[:] = np.arange(1, V.dim() + 1, dtype=float)
+
+    derivative = FullWaveformInversion._gradient_to_derivative(dummy, gradient)
+
+    u = fire.TrialFunction(V)
+    v = fire.TestFunction(V)
+    ones = fire.Function(V).assign(1.0)
+    mass_diagonal = fire.assemble(fire.action(fire.inner(u, v) * fire.dx, ones))
+    expected = gradient.dat.data_ro[:] * mass_diagonal.dat.data_ro[:]
+
+    assert np.allclose(derivative.dat.data_ro[:], expected)
+
+
 if __name__ == "__main__":
     test_butter_lowpast_filter()
     test_geometry_creation()