From c3d5e14c2cb9c64a45230bbf7cbc1ad4221fa8a7 Mon Sep 17 00:00:00 2001
From: Davor Dundovic <33790330+ddundo@users.noreply.github.com>
Date: Wed, 6 Mar 2024 19:04:27 +0100
Subject: [PATCH] Avoid interpolation deprecation warning (#129)

Closes #105.
---
 demos/burgers-hessian.py                      |  2 +-
 demos/burgers.py                              |  2 +-
 demos/burgers1.py                             |  2 +-
 demos/burgers2.py                             |  2 +-
 demos/burgers_ee.py                           |  2 +-
 demos/burgers_oo.py                           |  2 +-
 demos/burgers_time_integrated.py              |  2 +-
 demos/gray_scott_split.py                     | 17 +++++++++--------
 demos/solid_body_rotation.py                  |  2 +-
 demos/solid_body_rotation_split.py            |  2 +-
 goalie/__init__.py                            |  1 +
 goalie/utility.py                             |  4 +++-
 test/test_interpolation.py                    |  9 +++++----
 test/test_utility.py                          | 16 ++++++++--------
 test_adjoint/examples/burgers.py              |  4 +++-
 test_adjoint/examples/steady_flow_past_cyl.py |  4 +++-
 16 files changed, 41 insertions(+), 32 deletions(-)

diff --git a/demos/burgers-hessian.py b/demos/burgers-hessian.py
index b6d1ada2..98ff7311 100644
--- a/demos/burgers-hessian.py
+++ b/demos/burgers-hessian.py
@@ -74,7 +74,7 @@ def solver(index, ic):
 def get_initial_condition(mesh_seq):
     fs = mesh_seq.function_spaces["u"][0]
     x, y = SpatialCoordinate(mesh_seq[0])
-    return {"u": interpolate(as_vector([sin(pi * x), 0]), fs)}
+    return {"u": assemble(interpolate(as_vector([sin(pi * x), 0]), fs))}
 
 
 n = 32
diff --git a/demos/burgers.py b/demos/burgers.py
index 6b8c88c8..5dda6a64 100644
--- a/demos/burgers.py
+++ b/demos/burgers.py
@@ -129,7 +129,7 @@ def solver(index, ic):
 def get_initial_condition(mesh_seq):
     fs = mesh_seq.function_spaces["u"][0]
     x, y = SpatialCoordinate(mesh_seq[0])
-    return {"u": interpolate(as_vector([sin(pi * x), 0]), fs)}
+    return {"u": assemble(interpolate(as_vector([sin(pi * x), 0]), fs))}
 
 
 # Now that we have the above functions defined, we move onto the
diff --git a/demos/burgers1.py b/demos/burgers1.py
index d077bf88..170b42be 100644
--- a/demos/burgers1.py
+++ b/demos/burgers1.py
@@ -75,7 +75,7 @@ def solver(index, ic):
 def get_initial_condition(mesh_seq):
     fs = mesh_seq.function_spaces["u"][0]
     x, y = SpatialCoordinate(mesh_seq[0])
-    return {"u": interpolate(as_vector([sin(pi * x), 0]), fs)}
+    return {"u": assemble(interpolate(as_vector([sin(pi * x), 0]), fs))}
 
 
 # In line with the
diff --git a/demos/burgers2.py b/demos/burgers2.py
index fe23279e..ada6c7fc 100644
--- a/demos/burgers2.py
+++ b/demos/burgers2.py
@@ -75,7 +75,7 @@ def solver(index, ic):
 def get_initial_condition(mesh_seq):
     fs = mesh_seq.function_spaces["u"][0]
     x, y = SpatialCoordinate(mesh_seq[0])
-    return {"u": interpolate(as_vector([sin(pi * x), 0]), fs)}
+    return {"u": assemble(interpolate(as_vector([sin(pi * x), 0]), fs))}
 
 
 def get_qoi(mesh_seq, solutions, i):
diff --git a/demos/burgers_ee.py b/demos/burgers_ee.py
index c9ce4d67..bc41c8bc 100644
--- a/demos/burgers_ee.py
+++ b/demos/burgers_ee.py
@@ -90,7 +90,7 @@ def solver(index, ic):
 def get_initial_condition(mesh_seq):
     fs = mesh_seq.function_spaces["u"][0]
     x, y = SpatialCoordinate(mesh_seq[0])
-    return {"u": interpolate(as_vector([sin(pi * x), 0]), fs)}
+    return {"u": assemble(interpolate(as_vector([sin(pi * x), 0]), fs))}
 
 
 def get_qoi(mesh_seq, solutions, i):
diff --git a/demos/burgers_oo.py b/demos/burgers_oo.py
index 9b3c88fd..8e3b2722 100644
--- a/demos/burgers_oo.py
+++ b/demos/burgers_oo.py
@@ -82,7 +82,7 @@ def solver(index, ic):
     def get_initial_condition(self):
         fs = self.function_spaces["u"][0]
         x, y = SpatialCoordinate(self[0])
-        return {"u": interpolate(as_vector([sin(pi * x), 0]), fs)}
+        return {"u": assemble(interpolate(as_vector([sin(pi * x), 0]), fs))}
 
     @annotate_qoi
     def get_qoi(self, solutions, i):
diff --git a/demos/burgers_time_integrated.py b/demos/burgers_time_integrated.py
index 379edf45..f787faa9 100644
--- a/demos/burgers_time_integrated.py
+++ b/demos/burgers_time_integrated.py
@@ -44,7 +44,7 @@ def form(index, solutions):
 def get_initial_condition(mesh_seq):
     fs = mesh_seq.function_spaces["u"][0]
     x, y = SpatialCoordinate(mesh_seq[0])
-    return {"u": interpolate(as_vector([sin(pi * x), 0]), fs)}
+    return {"u": assemble(interpolate(as_vector([sin(pi * x), 0]), fs))}
 
 
 # The solver needs to be modified slightly in order to take
diff --git a/demos/gray_scott_split.py b/demos/gray_scott_split.py
index f028bfad..02115112 100644
--- a/demos/gray_scott_split.py
+++ b/demos/gray_scott_split.py
@@ -32,16 +32,17 @@ def get_initial_condition(mesh_seq):
     x, y = SpatialCoordinate(mesh_seq[0])
     fs_a = mesh_seq.function_spaces["a"][0]
     fs_b = mesh_seq.function_spaces["b"][0]
-    a_init = Function(fs_a, name="a")
-    b_init = Function(fs_b, name="b")
-    b_init.interpolate(
-        conditional(
-            And(And(1 <= x, x <= 1.5), And(1 <= y, y <= 1.5)),
-            0.25 * sin(4 * pi * x) ** 2 * sin(4 * pi * y) ** 2,
-            0,
+    b_init = assemble(
+        interpolate(
+            conditional(
+                And(And(1 <= x, x <= 1.5), And(1 <= y, y <= 1.5)),
+                0.25 * sin(4 * pi * x) ** 2 * sin(4 * pi * y) ** 2,
+                0,
+            ),
+            fs_b,
         )
     )
-    a_init.interpolate(1 - 2 * b_init)
+    a_init = assemble(interpolate(1 - 2 * b_init, fs_a))
     return {"a": a_init, "b": b_init}
 
 
diff --git a/demos/solid_body_rotation.py b/demos/solid_body_rotation.py
index 4935f9b0..6aa3657f 100644
--- a/demos/solid_body_rotation.py
+++ b/demos/solid_body_rotation.py
@@ -91,7 +91,7 @@ def get_initial_condition(mesh_seq, field="c"):
     bell = bell_initial_condition(x, y)
     cone = cone_initial_condition(x, y)
     slot_cyl = slot_cyl_initial_condition(x, y)
-    return {field: interpolate(bell + cone + slot_cyl, fs)}
+    return {field: assemble(interpolate(bell + cone + slot_cyl, fs))}
 
 
 # Now let's set up the time interval of interest. The `"GOALIE_REGRESSION_TEST"` flag
diff --git a/demos/solid_body_rotation_split.py b/demos/solid_body_rotation_split.py
index a3fc300b..35a724fe 100644
--- a/demos/solid_body_rotation_split.py
+++ b/demos/solid_body_rotation_split.py
@@ -54,7 +54,7 @@ def get_initial_condition_split(mesh_seq):
         "slot_cyl": slot_cyl_initial_condition,
     }
     return {
-        f: interpolate(init[f](x, y), fs[0])
+        f: assemble(interpolate(init[f](x, y), fs[0]))
         for f, fs in mesh_seq.function_spaces.items()
     }
 
diff --git a/goalie/__init__.py b/goalie/__init__.py
index 7652f8b5..d943c9c4 100644
--- a/goalie/__init__.py
+++ b/goalie/__init__.py
@@ -9,6 +9,7 @@
 from goalie.point_seq import *  # noqa
 from goalie.interpolation import *  # noqa
 from goalie.error_estimation import *  # noqa
+from firedrake.__future__ import interpolate  # noqa
 
 import numpy as np  # noqa
 import os  # noqa
diff --git a/goalie/utility.py b/goalie/utility.py
index 9923010c..89cc2c9a 100644
--- a/goalie/utility.py
+++ b/goalie/utility.py
@@ -1,10 +1,12 @@
 """
 Utility functions and classes for mesh adaptation.
 """
+
 from collections import OrderedDict
 import firedrake
 import firedrake.mesh as fmesh
 from firedrake.petsc import PETSc
+from firedrake.__future__ import interpolate
 import mpi4py
 import numpy as np
 import os
@@ -49,7 +51,7 @@ def Mesh(arg, **kwargs) -> firedrake.mesh.MeshGeometry:
 
     # Cell size
     if dim == 2 and mesh.coordinates.ufl_element().cell == ufl.triangle:
-        mesh.delta_x = firedrake.interpolate(ufl.CellDiameter(mesh), P0)
+        mesh.delta_x = firedrake.assemble(interpolate(ufl.CellDiameter(mesh), P0))
 
     return mesh
 
diff --git a/test/test_interpolation.py b/test/test_interpolation.py
index 53473ad7..99d449b1 100644
--- a/test/test_interpolation.py
+++ b/test/test_interpolation.py
@@ -1,6 +1,7 @@
 """
 Test interpolation schemes.
 """
+
 from firedrake import *
 from goalie import *
 from goalie.utility import function2cofunction
@@ -76,7 +77,7 @@ def test_wrong_number_sub_spaces(self):
 
     def test_project_same_space(self):
         Vs = FunctionSpace(self.source_mesh, "CG", 1)
-        source = interpolate(self.sinusoid(), Vs)
+        source = assemble(interpolate(self.sinusoid(), Vs))
         target = Function(Vs)
         project(source, target)
         expected = source
@@ -84,7 +85,7 @@ def test_project_same_space(self):
 
     def test_project_same_space_adjoint(self):
         Vs = FunctionSpace(self.source_mesh, "CG", 1)
-        source = interpolate(self.sinusoid(), Vs)
+        source = assemble(interpolate(self.sinusoid(), Vs))
         source = function2cofunction(source)
         target = Cofunction(Vs.dual())
         project(source, target)
@@ -119,7 +120,7 @@ def test_project_same_space_mixed_adjoint(self):
     def test_project_same_mesh(self):
         Vs = FunctionSpace(self.source_mesh, "CG", 1)
         Vt = FunctionSpace(self.source_mesh, "DG", 0)
-        source = interpolate(self.sinusoid(), Vs)
+        source = assemble(interpolate(self.sinusoid(), Vs))
         target = Function(Vt)
         project(source, target)
         expected = Function(Vt).project(source)
@@ -128,7 +129,7 @@ def test_project_same_mesh(self):
     def test_project_same_mesh_adjoint(self):
         Vs = FunctionSpace(self.source_mesh, "CG", 1)
         Vt = FunctionSpace(self.source_mesh, "DG", 0)
-        source = interpolate(self.sinusoid(), Vs)
+        source = assemble(interpolate(self.sinusoid(), Vs))
         target = Cofunction(Vt.dual())
         project(function2cofunction(source), target)
         expected = function2cofunction(Function(Vt).project(source))
diff --git a/test/test_utility.py b/test/test_utility.py
index 6a9fa27b..eb4018c8 100644
--- a/test/test_utility.py
+++ b/test/test_utility.py
@@ -1,6 +1,7 @@
 """
 Test utility functions.
 """
+
 from firedrake import *
 from firedrake.norms import norm as fnorm
 from firedrake.norms import errornorm as ferrnorm
@@ -95,7 +96,7 @@ def setUp(self):
         self.mesh = uniform_mesh(2, 4)
         self.x, self.y = SpatialCoordinate(self.mesh)
         V = FunctionSpace(self.mesh, "CG", 1)
-        self.f = interpolate(self.x**2 + self.y, V)
+        self.f = assemble(interpolate(self.x**2 + self.y, V))
 
     def test_boundary_error(self):
         with self.assertRaises(NotImplementedError) as cm:
@@ -150,8 +151,7 @@ def test_consistency_firedrake(self, norm_type):
     def test_consistency_hdiv(self):
         V = VectorFunctionSpace(self.mesh, "CG", 1)
         x, y = SpatialCoordinate(self.mesh)
-        f = Function(V)
-        f.interpolate(as_vector([y * y, -x * x]))
+        f = assemble(interpolate(as_vector([y * y, -x * x]), V))
         expected = fnorm(f, norm_type="HDiv")
         got = norm(f, norm_type="HDiv")
         self.assertAlmostEqual(expected, got)
@@ -179,8 +179,8 @@ def setUp(self):
         self.mesh = uniform_mesh(2, 4)
         self.x, self.y = SpatialCoordinate(self.mesh)
         V = FunctionSpace(self.mesh, "CG", 1)
-        self.f = interpolate(self.x**2 + self.y, V)
-        self.g = interpolate(self.x + self.y**2, V)
+        self.f = assemble(interpolate(self.x**2 + self.y, V))
+        self.g = assemble(interpolate(self.x + self.y**2, V))
 
     def test_shape_error(self):
         with self.assertRaises(RuntimeError) as cm:
@@ -215,7 +215,7 @@ def test_zero_scalar(self, norm_type):
     def test_zero_hdiv(self):
         V = VectorFunctionSpace(self.mesh, "CG", 1)
         x, y = SpatialCoordinate(self.mesh)
-        f = interpolate(as_vector([y * y, -x * x]), V)
+        f = assemble(interpolate(as_vector([y * y, -x * x]), V))
         err = errornorm(f, f, norm_type="HDiv")
         self.assertAlmostEqual(err, 0.0)
 
@@ -228,8 +228,8 @@ def test_consistency_firedrake(self, norm_type):
     def test_consistency_hdiv(self):
         V = VectorFunctionSpace(self.mesh, "CG", 1)
         x, y = SpatialCoordinate(self.mesh)
-        f = interpolate(as_vector([y * y, -x * x]), V)
-        g = interpolate(as_vector([x * y, 1.0]), V)
+        f = assemble(interpolate(as_vector([y * y, -x * x]), V))
+        g = assemble(interpolate(as_vector([x * y, 1.0]), V))
         expected = ferrnorm(f, g, norm_type="HDiv")
         got = errornorm(f, g, norm_type="HDiv")
         self.assertAlmostEqual(expected, got)
diff --git a/test_adjoint/examples/burgers.py b/test_adjoint/examples/burgers.py
index fea98b83..497beb28 100644
--- a/test_adjoint/examples/burgers.py
+++ b/test_adjoint/examples/burgers.py
@@ -9,7 +9,9 @@
 Code here is based on that found at
     https://firedrakeproject.org/demos/burgers.py.html
 """
+
 from firedrake import *
+from firedrake.__future__ import interpolate
 
 
 # Problem setup
@@ -95,7 +97,7 @@ def get_initial_condition(self):
     """
     init_fs = self.function_spaces["uv_2d"][0]
     x, y = SpatialCoordinate(self.meshes[0])
-    return {"uv_2d": interpolate(as_vector([sin(pi * x), 0]), init_fs)}
+    return {"uv_2d": assemble(interpolate(as_vector([sin(pi * x), 0]), init_fs))}
 
 
 def get_qoi(self, sol, i):
diff --git a/test_adjoint/examples/steady_flow_past_cyl.py b/test_adjoint/examples/steady_flow_past_cyl.py
index 4f77ec59..7864a961 100644
--- a/test_adjoint/examples/steady_flow_past_cyl.py
+++ b/test_adjoint/examples/steady_flow_past_cyl.py
@@ -10,7 +10,9 @@
 Code here is based on that found at
     https://nbviewer.jupyter.org/github/firedrakeproject/firedrake/blob/master/docs/notebooks/06-pde-constrained-optimisation.ipynb
 """
+
 from firedrake import *
+from firedrake.__future__ import interpolate
 import os
 
 
@@ -63,7 +65,7 @@ def bcs(i):
         u_inflow = as_vector([y * (10 - y) / 25.0, 0])
         W = self.function_spaces["up"][i]
         noslip = DirichletBC(W.sub(0), (0, 0), (3, 5))
-        inflow = DirichletBC(W.sub(0), interpolate(u_inflow, W.sub(0)), 1)
+        inflow = DirichletBC(W.sub(0), assemble(interpolate(u_inflow, W.sub(0))), 1)
         return [inflow, noslip, DirichletBC(W.sub(0), 0, 4)]
 
     return bcs