From fee2e4a134a3e3bd11ff5c2e09ca85da76e887b0 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Wed, 28 Feb 2024 07:52:01 +0000
Subject: [PATCH 01/12] fix the optmisation problem for firedrake and pyrol.

---
 .../stokes-topology-rol-firedrake.py              | 13 ++++++-------
 pyadjoint/optimization/rol_solver.py              | 15 +++++++++++----
 2 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/examples/stokes-topology/stokes-topology-rol-firedrake.py b/examples/stokes-topology/stokes-topology-rol-firedrake.py
index 6fe2d33c..bfed4fcf 100644
--- a/examples/stokes-topology/stokes-topology-rol-firedrake.py
+++ b/examples/stokes-topology/stokes-topology-rol-firedrake.py
@@ -20,14 +20,16 @@
 # lazy evaluation mode.
 parameters["pyop2_options"]["lazy_evaluation"] = False
 
-from firedrake_adjoint import *
+from firedrake.adjoint import *
+from firedrake.__future__ import interpolate
 
 try:
     import ROL
 except ImportError:
     info_red("""This example depends on ROL.""")
     raise
-
+# Starting the annotation.
+continue_annotation()
 mu = Constant(1.0)                   # viscosity
 alphaunderbar = 2.5 * mu / (100**2)  # parameter for \alpha
 alphabar = 2.5 * mu / (0.01**2)      # parameter for \alpha
@@ -79,7 +81,7 @@ def forward(rho):
     return w
 
 if __name__ == "__main__":
-    rho = interpolate(Constant(float(V)/delta), A)
+    rho = assemble(interpolate(Constant(float(V)/delta), A))
     w   = forward(rho)
     (u, p) = split(w)
 
@@ -128,11 +130,9 @@ def eval_cb(j, rho):
 
     solver = ROLSolver(problem, params, inner_product="L2")
     rho_opt = solver.solve()
-
+    get_working_tape().clear_tape()
     q.assign(0.1)
     rho.assign(rho_opt)
-    get_working_tape().clear_tape()
-
     w = forward(rho)
     (u, p) = split(w)
 
@@ -154,4 +154,3 @@ def eval_cb(j, rho):
     rho_opt = solver.solve()
     rho_viz.assign(rho_opt)
     controls.write(rho_viz)
-
diff --git a/pyadjoint/optimization/rol_solver.py b/pyadjoint/optimization/rol_solver.py
index 06222478..cd81eaad 100644
--- a/pyadjoint/optimization/rol_solver.py
+++ b/pyadjoint/optimization/rol_solver.py
@@ -71,8 +71,8 @@ def scale(self, alpha):
         def riesz_map(self, derivs):
             dat = []
             opts = {"riesz_representation": self.inner_product}
-            for deriv in Enlist(derivs):
-                dat.append(deriv._ad_convert_type(deriv, options=opts))
+            for f, deriv in zip(self.dat, derivs):
+                dat.append(f._ad_convert_type(deriv, options=opts))
             return dat
 
         def dot(self, yy):
@@ -82,6 +82,12 @@ def dot(self, yy):
                 res += x._ad_dot(y, options=opts)
             return res
 
+        def dual(self):
+            dat = []
+            for x in self.dat:
+                dat.append(x.riesz_representation())
+            return ROLVector(dat, inner_product=self.inner_product)
+
         def norm(self):
             return self.dot(self) ** 0.5
 
@@ -123,8 +129,9 @@ def applyJacobian(self, jv, v, x, tol):
             self.con.jacobian_action(x.dat, v.dat[0], jv.dat)
 
         def applyAdjointJacobian(self, jv, v, x, tol):
-            self.con.jacobian_adjoint_action(x.dat, v.dat, jv.dat[0])
-            jv.dat = jv.riesz_map(jv.dat)
+            tmp = jv.dual().clone()
+            self.con.jacobian_adjoint_action(x.dat, v.dat, tmp.dat[0])
+            jv.dat = jv.riesz_map(tmp.dat)
 
         def applyAdjointHessian(self, ahuv, u, v, x, tol):
             self.con.hessian_action(x.dat, u.dat[0], v.dat, ahuv.dat[0])

From 3e03e99a6e82309ad9844c5af553f6dc6f0096cf Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Wed, 28 Feb 2024 13:45:00 +0000
Subject: [PATCH 02/12] Create the riesz representation in overloaded type.

---
 pyadjoint/optimization/rol_solver.py | 15 +++++++++++++++
 pyadjoint/overloaded_type.py         | 14 ++++++++++++++
 2 files changed, 29 insertions(+)

diff --git a/pyadjoint/optimization/rol_solver.py b/pyadjoint/optimization/rol_solver.py
index cd81eaad..bda3c785 100644
--- a/pyadjoint/optimization/rol_solver.py
+++ b/pyadjoint/optimization/rol_solver.py
@@ -83,6 +83,13 @@ def dot(self, yy):
             return res
 
         def dual(self):
+            """Create a dual vector.
+
+            This is a `ROLVector` in the dual space of the current `ROLVector`.
+
+            Returns:
+                ROLVector: A `ROLVector` in dual space.
+            """
             dat = []
             for x in self.dat:
                 dat.append(x.riesz_representation())
@@ -129,6 +136,14 @@ def applyJacobian(self, jv, v, x, tol):
             self.con.jacobian_action(x.dat, v.dat[0], jv.dat)
 
         def applyAdjointJacobian(self, jv, v, x, tol):
+            """
+
+            Args:
+                jv (ROLVector): The result of the adjoint action in primal space.
+                v (ROLVector):
+                x (ROLVector):
+                tol (float):
+            """
             tmp = jv.dual().clone()
             self.con.jacobian_adjoint_action(x.dat, v.dat, tmp.dat[0])
             jv.dat = jv.riesz_map(tmp.dat)
diff --git a/pyadjoint/overloaded_type.py b/pyadjoint/overloaded_type.py
index 1d286aee..dd49af3b 100644
--- a/pyadjoint/overloaded_type.py
+++ b/pyadjoint/overloaded_type.py
@@ -116,6 +116,20 @@ def _ad_convert_type(self, value, options={}):
         """
         raise NotImplementedError(f"OverloadedType._ad_convert_type not defined for class {type(self)}.")
 
+    def _riesz_representation(self, options={}):
+        """This method must be overridden.
+
+        Should implement a way to return the Riesz representation of the overloaded object.
+
+        Args:
+            options (dict): A dictionary with options that may be supplied by the user. If the Riesz representation
+                functionality offers some options on how to compute it, this is the dictionary that should be used.
+
+        Returns:
+            OverloadedType: The Riesz representation of the overloaded object.
+        """
+        raise NotImplementedError(f"OverloadedType._riesz_representation not defined for class {type(self)}.")
+    
     def _ad_create_checkpoint(self):
         """This method must be overridden.
 

From 8ba78d8f97c415f2438e6c0e003cccca08e70c19 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Thu, 29 Feb 2024 17:23:46 +0000
Subject: [PATCH 03/12] Add riesz representation overload type.

---
 pyadjoint/optimization/rol_solver.py | 20 +++++---------------
 pyadjoint/overloaded_type.py         | 16 ++++++++++++++++
 2 files changed, 21 insertions(+), 15 deletions(-)

diff --git a/pyadjoint/optimization/rol_solver.py b/pyadjoint/optimization/rol_solver.py
index bda3c785..b54bba49 100644
--- a/pyadjoint/optimization/rol_solver.py
+++ b/pyadjoint/optimization/rol_solver.py
@@ -82,17 +82,15 @@ def dot(self, yy):
                 res += x._ad_dot(y, options=opts)
             return res
 
-        def dual(self):
+        def dual(self) -> "ROLVector":
             """Create a dual vector.
 
-            This is a `ROLVector` in the dual space of the current `ROLVector`.
-
-            Returns:
-                ROLVector: A `ROLVector` in dual space.
+            This is a `ROLVector` in the dual space of the current `ROLVector`, which is in the primal space.
             """
             dat = []
+            opts = {"riesz_map": self.inner_product}
             for x in self.dat:
-                dat.append(x.riesz_representation())
+                dat.append(x._riesz_representation(options=opts))
             return ROLVector(dat, inner_product=self.inner_product)
 
         def norm(self):
@@ -136,15 +134,7 @@ def applyJacobian(self, jv, v, x, tol):
             self.con.jacobian_action(x.dat, v.dat[0], jv.dat)
 
         def applyAdjointJacobian(self, jv, v, x, tol):
-            """
-
-            Args:
-                jv (ROLVector): The result of the adjoint action in primal space.
-                v (ROLVector):
-                x (ROLVector):
-                tol (float):
-            """
-            tmp = jv.dual().clone()
+            tmp = jv.dual()
             self.con.jacobian_adjoint_action(x.dat, v.dat, tmp.dat[0])
             jv.dat = jv.riesz_map(tmp.dat)
 
diff --git a/pyadjoint/overloaded_type.py b/pyadjoint/overloaded_type.py
index dd49af3b..29e260ed 100644
--- a/pyadjoint/overloaded_type.py
+++ b/pyadjoint/overloaded_type.py
@@ -64,6 +64,22 @@ def register_overloaded_type(overloaded_type, classes=None):
     return overloaded_type
 
 
+def _riesz_representation(self, riesz_map={}):
+        """This method must be overridden.
+
+        Should implement a way to convert `self` onto its dual space with respect to
+        a certain Riesz map.
+
+        Args:
+            riesz_map (dict): A dictionary with the Riesz map. This dictionary should have for instance
+            the norm used to have the riesz representation of self.
+
+        Returns:
+            OverloadedType: An `self` in its dual space.
+
+        """
+        raise NotImplementedError(f"OverloadedType._ad_convert_type not defined for class {type(self)}.")
+
 class OverloadedType(object):
     """Base class for OverloadedType types.
 

From b2750e082fc653826766a8d3a0b1221a767da291 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Thu, 29 Feb 2024 17:27:49 +0000
Subject: [PATCH 04/12] Test daiane's repository.

---
 .circleci/config.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 6b29bc45..9cc8571c 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -59,7 +59,7 @@ jobs:
       - run:
           name: Install dolfin-adjoint
           command: |
-            git clone https://github.com/dolfin-adjoint/dolfin-adjoint ~/dolfin-adjoint
+            git clone https://github.com/Ig-dolci/dolfin-adjoint/tree/main ~/dolfin-adjoint
             python3 -m pip install -e ~/dolfin-adjoint[all]
 
       - run:

From fc89445fba14439e1f50fd5799f5af3aa8ded8c8 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Thu, 29 Feb 2024 17:32:02 +0000
Subject: [PATCH 05/12] wip.

---
 .circleci/config.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 9cc8571c..65ffe108 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -59,7 +59,7 @@ jobs:
       - run:
           name: Install dolfin-adjoint
           command: |
-            git clone https://github.com/Ig-dolci/dolfin-adjoint/tree/main ~/dolfin-adjoint
+            git clone https://github.com/Ig-dolci/dolfin-adjoint.git ~/dolfin-adjoint
             python3 -m pip install -e ~/dolfin-adjoint[all]
 
       - run:

From fdc0d8aad11e9d2bf292f27c20d1b74d1af41f89 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Thu, 29 Feb 2024 17:34:24 +0000
Subject: [PATCH 06/12] wip.

---
 pyadjoint/overloaded_type.py | 17 -----------------
 1 file changed, 17 deletions(-)

diff --git a/pyadjoint/overloaded_type.py b/pyadjoint/overloaded_type.py
index 29e260ed..04d1ad27 100644
--- a/pyadjoint/overloaded_type.py
+++ b/pyadjoint/overloaded_type.py
@@ -63,23 +63,6 @@ def register_overloaded_type(overloaded_type, classes=None):
         _overloaded_types[classes] = overloaded_type
     return overloaded_type
 
-
-def _riesz_representation(self, riesz_map={}):
-        """This method must be overridden.
-
-        Should implement a way to convert `self` onto its dual space with respect to
-        a certain Riesz map.
-
-        Args:
-            riesz_map (dict): A dictionary with the Riesz map. This dictionary should have for instance
-            the norm used to have the riesz representation of self.
-
-        Returns:
-            OverloadedType: An `self` in its dual space.
-
-        """
-        raise NotImplementedError(f"OverloadedType._ad_convert_type not defined for class {type(self)}.")
-
 class OverloadedType(object):
     """Base class for OverloadedType types.
 

From db6a8305f24f5cf5f834fa9756865d87778a8689 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Thu, 29 Feb 2024 17:40:23 +0000
Subject: [PATCH 07/12] wip.

---
 examples/stokes-topology/stokes-topology-rol-firedrake.py | 8 ++++----
 pyadjoint/overloaded_type.py                              | 3 ++-
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/examples/stokes-topology/stokes-topology-rol-firedrake.py b/examples/stokes-topology/stokes-topology-rol-firedrake.py
index bfed4fcf..aa0532ba 100644
--- a/examples/stokes-topology/stokes-topology-rol-firedrake.py
+++ b/examples/stokes-topology/stokes-topology-rol-firedrake.py
@@ -28,8 +28,7 @@
 except ImportError:
     info_red("""This example depends on ROL.""")
     raise
-# Starting the annotation.
-continue_annotation()
+
 mu = Constant(1.0)                   # viscosity
 alphaunderbar = 2.5 * mu / (100**2)  # parameter for \alpha
 alphabar = 2.5 * mu / (0.01**2)      # parameter for \alpha
@@ -81,6 +80,7 @@ def forward(rho):
     return w
 
 if __name__ == "__main__":
+    continue_annotation()
     rho = assemble(interpolate(Constant(float(V)/delta), A))
     w   = forward(rho)
     (u, p) = split(w)
@@ -130,9 +130,9 @@ def eval_cb(j, rho):
 
     solver = ROLSolver(problem, params, inner_product="L2")
     rho_opt = solver.solve()
-    get_working_tape().clear_tape()
     q.assign(0.1)
-    rho.assign(rho_opt)
+    rho.interpolate(rho_opt)
+    get_working_tape().clear_tape()
     w = forward(rho)
     (u, p) = split(w)
 
diff --git a/pyadjoint/overloaded_type.py b/pyadjoint/overloaded_type.py
index 04d1ad27..96e59f71 100644
--- a/pyadjoint/overloaded_type.py
+++ b/pyadjoint/overloaded_type.py
@@ -63,6 +63,7 @@ def register_overloaded_type(overloaded_type, classes=None):
         _overloaded_types[classes] = overloaded_type
     return overloaded_type
 
+
 class OverloadedType(object):
     """Base class for OverloadedType types.
 
@@ -128,7 +129,7 @@ def _riesz_representation(self, options={}):
             OverloadedType: The Riesz representation of the overloaded object.
         """
         raise NotImplementedError(f"OverloadedType._riesz_representation not defined for class {type(self)}.")
-    
+
     def _ad_create_checkpoint(self):
         """This method must be overridden.
 

From 11696acdb65f839dbe0cd653c2a5410ebfaa49fe Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Mon, 4 Mar 2024 09:52:35 +0000
Subject: [PATCH 08/12] wip

---
 .../stokes-topology/stokes-topology-rol-firedrake.py  | 11 ++++++-----
 pyadjoint/optimization/rol_solver.py                  |  4 +---
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/examples/stokes-topology/stokes-topology-rol-firedrake.py b/examples/stokes-topology/stokes-topology-rol-firedrake.py
index aa0532ba..6fe2d33c 100644
--- a/examples/stokes-topology/stokes-topology-rol-firedrake.py
+++ b/examples/stokes-topology/stokes-topology-rol-firedrake.py
@@ -20,8 +20,7 @@
 # lazy evaluation mode.
 parameters["pyop2_options"]["lazy_evaluation"] = False
 
-from firedrake.adjoint import *
-from firedrake.__future__ import interpolate
+from firedrake_adjoint import *
 
 try:
     import ROL
@@ -80,8 +79,7 @@ def forward(rho):
     return w
 
 if __name__ == "__main__":
-    continue_annotation()
-    rho = assemble(interpolate(Constant(float(V)/delta), A))
+    rho = interpolate(Constant(float(V)/delta), A)
     w   = forward(rho)
     (u, p) = split(w)
 
@@ -130,9 +128,11 @@ def eval_cb(j, rho):
 
     solver = ROLSolver(problem, params, inner_product="L2")
     rho_opt = solver.solve()
+
     q.assign(0.1)
-    rho.interpolate(rho_opt)
+    rho.assign(rho_opt)
     get_working_tape().clear_tape()
+
     w = forward(rho)
     (u, p) = split(w)
 
@@ -154,3 +154,4 @@ def eval_cb(j, rho):
     rho_opt = solver.solve()
     rho_viz.assign(rho_opt)
     controls.write(rho_viz)
+
diff --git a/pyadjoint/optimization/rol_solver.py b/pyadjoint/optimization/rol_solver.py
index b54bba49..b6b3c614 100644
--- a/pyadjoint/optimization/rol_solver.py
+++ b/pyadjoint/optimization/rol_solver.py
@@ -83,9 +83,7 @@ def dot(self, yy):
             return res
 
         def dual(self) -> "ROLVector":
-            """Create a dual vector.
-
-            This is a `ROLVector` in the dual space of the current `ROLVector`, which is in the primal space.
+            """Create a new `ROLVector` in the dual space of the current `self`.
             """
             dat = []
             opts = {"riesz_map": self.inner_product}

From 6bee025eada473f2bdeaf98fcf76e746a7e3ce26 Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Mon, 4 Mar 2024 10:24:01 +0000
Subject: [PATCH 09/12] wip

---
 .circleci/config.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 65ffe108..6b29bc45 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -59,7 +59,7 @@ jobs:
       - run:
           name: Install dolfin-adjoint
           command: |
-            git clone https://github.com/Ig-dolci/dolfin-adjoint.git ~/dolfin-adjoint
+            git clone https://github.com/dolfin-adjoint/dolfin-adjoint ~/dolfin-adjoint
             python3 -m pip install -e ~/dolfin-adjoint[all]
 
       - run:

From 61611c15557ef6008736c5da5bbfe172025665ee Mon Sep 17 00:00:00 2001
From: Iglesia Dolci <daia.dolci@gmail.com>
Date: Mon, 4 Mar 2024 17:46:42 +0000
Subject: [PATCH 10/12] wip

---
 pyadjoint/optimization/rol_solver.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/pyadjoint/optimization/rol_solver.py b/pyadjoint/optimization/rol_solver.py
index b6b3c614..89d592b5 100644
--- a/pyadjoint/optimization/rol_solver.py
+++ b/pyadjoint/optimization/rol_solver.py
@@ -137,8 +137,9 @@ def applyAdjointJacobian(self, jv, v, x, tol):
             jv.dat = jv.riesz_map(tmp.dat)
 
         def applyAdjointHessian(self, ahuv, u, v, x, tol):
-            self.con.hessian_action(x.dat, u.dat[0], v.dat, ahuv.dat[0])
-            ahuv.dat = ahuv.riesz_map(ahuv.dat)
+            tmp = ahuv.dual()
+            self.con.hessian_action(x.dat, u.dat[0], v.dat, tmp.dat[0])
+            ahuv.dat = ahuv.riesz_map(tmp.dat)
 
     class ROLSolver(OptimizationSolver):
         """

From ce06dd45bd1efa2fe267718fa1a58a111eb67faa Mon Sep 17 00:00:00 2001
From: Daiane Iglesia Dolci <63597005+Ig-dolci@users.noreply.github.com>
Date: Wed, 13 Mar 2024 11:08:38 +0000
Subject: [PATCH 11/12] rename overloadtype function

---
 pyadjoint/overloaded_type.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyadjoint/overloaded_type.py b/pyadjoint/overloaded_type.py
index 96e59f71..61b5f47e 100644
--- a/pyadjoint/overloaded_type.py
+++ b/pyadjoint/overloaded_type.py
@@ -116,7 +116,7 @@ def _ad_convert_type(self, value, options={}):
         """
         raise NotImplementedError(f"OverloadedType._ad_convert_type not defined for class {type(self)}.")
 
-    def _riesz_representation(self, options={}):
+    def _ad_riesz_representation(self, options={}):
         """This method must be overridden.
 
         Should implement a way to return the Riesz representation of the overloaded object.

From 203d905dcd06221a1fa3f26910199f7d9bfb05dc Mon Sep 17 00:00:00 2001
From: Daiane Iglesia Dolci <63597005+Ig-dolci@users.noreply.github.com>
Date: Wed, 13 Mar 2024 11:08:55 +0000
Subject: [PATCH 12/12] wip

---
 pyadjoint/optimization/rol_solver.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyadjoint/optimization/rol_solver.py b/pyadjoint/optimization/rol_solver.py
index 89d592b5..290c080d 100644
--- a/pyadjoint/optimization/rol_solver.py
+++ b/pyadjoint/optimization/rol_solver.py
@@ -88,7 +88,7 @@ def dual(self) -> "ROLVector":
             dat = []
             opts = {"riesz_map": self.inner_product}
             for x in self.dat:
-                dat.append(x._riesz_representation(options=opts))
+                dat.append(x._ad_riesz_representation(options=opts))
             return ROLVector(dat, inner_product=self.inner_product)
 
         def norm(self):