firedrake-adjoint and quadrature degree

[danshapero]

firedrake-adjoint seems to use the (very conservative) default quadrature degree estimation for adjoint solves, even if you specify a much lower quadrature degree for both the forward solver and the assembly of the quantity of interest. This causes some distracting warnings and reduced performance.

It would be nice if firedrake-adjoint could apply the same form compiler from the adjoint problem to the adjoint problem. If there's no way to recover the form compiler, there should at least be a way to specify form compiler parameters for the adjoint problem.

The code below demonstrates the problem; it's reduced from a demo I was making with icepack. I can make a more reduced example that doesn't use icepack if that helps with testing.

import firedrake
from firedrake import inner, ds
from firedrake_adjoint import Control, ReducedFunctional
import icepack
from icepack.constants import (
    ice_density as rho_I,
    water_density as rho_W,
    gravity as g,
    glen_flow_law as n,
)

nx, ny = 32, 32
Lx, Ly = 20e3, 20e3
mesh = firedrake.RectangleMesh(nx, ny, Lx, Ly)
Q = firedrake.FunctionSpace(mesh, "CG", 1)
V = firedrake.VectorFunctionSpace(mesh, "CG", 1)

x = firedrake.SpatialCoordinate(mesh)

# An exact solution for ice shelf flow
u_in = 100.0
h_in = 500.0
dh = 100.0
T = firedrake.Constant(260.0)
A = icepack.rate_factor(T)

rho = rho_I * (1 - rho_I / rho_W)
expr = u_in + (
    A * (rho * g * h_in / 4) ** n *
    (1 - (1 - dh / h_in) * (x[0] / Lx)) ** (n + 1) *
    Lx * (h_in / dh) / (n + 1)
)
u_exact = firedrake.as_vector((expr, 0))
u = firedrake.interpolate(u_exact, V)
h = firedrake.interpolate(h_in - dh * x[0] / Lx, Q)

model = icepack.models.IceShelf()
J = model.action(
    velocity=u,
    thickness=h,
    fluidity=A,
    ice_front_ids=[2],
    side_wall_ids=[3, 4],
)
F = firedrake.derivative(J, u)
inflow_ids = [1]
bc = firedrake.DirichletBC(V, u_exact, inflow_ids)
fc_params = {"form_compiler_parameters": {"quadrature_degree": 4}}
problem = firedrake.NonlinearVariationalProblem(F, u, bc, **fc_params)
solver = firedrake.NonlinearVariationalSolver(
    problem, solver_parameters={"snes_type": "newtontr"}
)
solver.solve()

nu = firedrake.FacetNormal(mesh)
flux = firedrake.assemble(h * inner(u, nu) * ds((2,)), **fc_params)
print(f"flux out of terminus: {float(flux) / 1e9} km^3 / year")
F = ReducedFunctional(flux, Control(h))
dF_dh = F.derivative()

[danshapero]

I'm closing this issue because I managed to fix it on the icepack side. The change was to attach the quadrature degree to dx, ds, etc. This seems to propagate through to the adjoint.