Compatibility fix for Firedrake 2024-12 Docker image (#24)

* Remove import of firedrake_configuration package
- Monkey-patch u_restrict property required by newer PETSc TS
- Package is no longer available due to pip install migration

* Fix unit test failures
- Migrate adjoint.py to newer interpolate API
- Set max SNES failures to make TS retry steps

* Migrate examples to newer VTKFile output API
- Silences deprecation warning emitted from Firedrake

* Adapt DAEProblem class to use to restricted solution
- Duck typing compatible with VariationalProblem class
- Use the new member everywhere

* Remove unused & forgotten LOC

---

This PR removes the (actually unused) import of `firedrake_configuration` package that is no longer present in recent (2 weeks or so) versions of [Firedrake Docker Hub image](https://hub.docker.com/r/firedrakeproject/firedrake/tags).

Root cause of the issue is the migration of Firedrake to be `pip`-installable. There's details about the topic [here](firedrakeproject/firedrake#3877). Specifically, the problem started with [the PR that migrates _PyOP2_ into the main _Firedrake_ repo](firedrakeproject/firedrake#3817).

---

It seems like the obvious fix, i.e. just remove the import, worked. The `firedrake_configuration` module wasn't used for anything directly within the `firedrake_ts` package. Besides this, the PR contains a few small fixes to make the unit tests and examples work.

It seems like PETSc 3.22 upgrage that's bundled in the Firedrake 2024-12 image doesn't set the [`ts_max_snes_failures`](https://petsc.org/release/manualpages/TS/TSSetMaxSNESFailures/) option anymore. Somehow the default value causes the inner SNES solver never to retry any linear solves. So the TS solver has no other option but to fail whenever it happens instead of retrying.

I also migrated the example applications to use the newer interpolate and VTK output Firedrake API.

Author: vojtamolda

examples/adjoint.py

@@ -14,7 +14,8 @@
 
 x = ufl.SpatialCoordinate(mesh)
 expr = ufl.as_vector([ufl.sin(2 * ufl.pi * x[0]), ufl.cos(2 * ufl.pi * x[1])])
-u = fd.interpolate(expr, V)
+u = fd.Function(V)
+u.interpolate(expr)
 
 u_dot = fd.Function(V)
 v = fd.TestFunction(V)
@@ -35,7 +36,7 @@
 end = 0.1
 tspan = (t, end)
 
-state_out = fd.File("result/state.pvd")
+state_out = fd.output.VTKFile("result/state.pvd")
 
 
 def ts_monitor(ts, steps, time, X):
@@ -81,5 +82,5 @@ def ts_monitor(ts, steps, time, X):
 print(f"Norm of dJdu after the adjoint solve: {dJdu_vec.norm()=} {fd.norm(fdJdu)=}")
 
 
-adj_out = fd.File("result/adj.pvd")
+adj_out = fd.output.VTKFile("result/adj.pvd")
 adj_out.write(fdJdu, time=0)

examples/burgers.py

@@ -24,7 +24,7 @@
     inner(u_t, v) + inner(dot(u, nabla_grad(u)), v) + nu * inner(grad(u), grad(v))
 ) * dx
 
-outfile = File("result/burgers.pvd")
+outfile = output.VTKFile("result/burgers.pvd")
 outfile.write(project(u, V_out, name="Velocity"), time=0.0)
 
 

firedrake_ts/solving_utils.py

@@ -3,7 +3,6 @@
 import numpy
 
 from pyop2 import op2
-from firedrake_configuration import get_config
 from firedrake import function, cofunction, dmhooks
 from firedrake.exceptions import ConvergenceError
 from firedrake.petsc import PETSc

firedrake_ts/ts_solver.py

@@ -81,15 +81,16 @@ def __init__(
         is_form_consistent(self.is_linear, self.bcs)
         self.Jp_eq_J = Jp is None
 
-        self.u = u
+        self.u_restrict = u
         self.udot = udot
         self.tspan = tspan
         self.F = F
         self.G = G
         self.Jp = Jp
-        if not isinstance(self.u, function.Function):
+
+        if not isinstance(self.u_restrict, function.Function):
             raise TypeError(
-                "Provided solution is a '%s', not a Function" % type(self.u).__name__
+                "Provided solution is a '%s', not a Function" % type(self.u_restrict).__name__
             )
         if not isinstance(self.udot, function.Function):
             raise TypeError(
@@ -123,7 +124,7 @@ def dirichlet_bcs(self):
 
     @utils.cached_property
     def dm(self):
-        return self.u.function_space().dm
+        return self.u_restrict.function_space().dm
 
 
 class DAESolver(OptionsManager):
@@ -227,14 +228,15 @@ def update_diffusivity(current_solution, current_time_derivative):
             # Mixed problem, use jacobi pc if user has not supplied
             # one.
             self.set_default_parameter("pc_type", "jacobi")
+        self.set_default_parameter("ts_max_snes_failures", -1)
 
         self.ts = PETSc.TS().create(comm=problem.dm.comm)
         self.snes = self.ts.getSNES()
 
         self._problem = problem
 
         self._ctx = ctx
-        self._work = problem.u.dof_dset.layout_vec.duplicate()
+        self._work = problem.u_restrict.dof_dset.layout_vec.duplicate()
 
         self.ts.setDM(problem.dm)
         self.ts.setMonitor(monitor_callback)
@@ -338,15 +340,15 @@ def solve(self, bounds=None):
         # Make sure appcontext is attached to the DM before we solve.
         dm = self.ts.getDM()
         for dbc in self._problem.dirichlet_bcs():
-            dbc.apply(self._problem.u)
+            dbc.apply(self._problem.u_restrict)
 
         if bounds is not None:
             lower, upper = bounds
             with lower.dat.vec_ro as lb, upper.dat.vec_ro as ub:
                 self.snes.setVariableBounds(lb, ub)
 
         work = self._work
-        with self._problem.u.dat.vec as u:
+        with self._problem.u_restrict.dat.vec as u:
             u.copy(work)
             with ExitStack() as stack:
                 # Ensure options database has full set of options (so monitors