From e00a20542f5bb3dc98c1fc56587230b1f73a7c14 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Mon, 20 Apr 2026 22:30:24 -0300
Subject: [PATCH] from gradient to drivative

---
 spyro/solvers/inversion.py      | 26 +++++++++++++++++++++++++-
 tests/on_one_core/test_utils.py | 25 +++++++++++++++++++++++++
 2 files changed, 50 insertions(+), 1 deletion(-)

diff --git a/spyro/solvers/inversion.py b/spyro/solvers/inversion.py
index 15ba8483d..63b3719b6 100644
--- a/spyro/solvers/inversion.py
+++ b/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):
diff --git a/tests/on_one_core/test_utils.py b/tests/on_one_core/test_utils.py
index f4e387443..61380f69e 100644
--- a/tests/on_one_core/test_utils.py
+++ b/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()