domain-less Constants and the adjoint

[stephankramer]

If I understand correctly, since two weeks there appears to be a (total?) ban on using Constants without a domain in the adjoint. This breaks the adjoint of Thetis in a way that is quite hard to deal with. The Thetis code is littered with Constants, most of them just arbitrary parameters, rescaling factors, physical constants, etc., and mostly these are created in a context where no domain is available, e.g. options classes. The only way I can think of fixing this is to have some utility function that attaches a missing domain (i.e. re-creating Constants but now with a domain) before it's being used or modified, but figuring this out proves to be a little involved. I also have to make sure I don't recreate Constants that have been user-provided, as these may be modified by the user later on. This basically means the Thetis user also has to set a domain in any constant they use in their setup, even if they don't care about the adjoint.

I also have to carefully make sure that constants that are used on different meshes, e.g. on a normal FEM mesh and on a VertexOnlyMesh, get different copies of each. This appears to also not always be handled correctly by the adjoint, e.g. in the following code with a functional on the vertex-only mesh:

from firedrake import *
from firedrake_adjoint import *
mesh = UnitSquareMesh(10,10)
mesh0 = VertexOnlyMesh(mesh, [[.5,.5]])
c = Constant(1.0, domain=mesh)
d = Constant(2.0, domain=mesh0)
d.assign(c)
J = assemble(d*dx(domain=mesh0))
rf = ReducedFunctional(J, Control(c))
rf(Constant(3.0, domain=mesh))

The initial forward model works as expected, copying the value of c on mesh into a constant d on mesh0 which is then integrated over mesh0. In the replay however, the checkpointed value of d after assignment appears to be on mesh which then causes the assembly of the functional to fail.

My questions are:

• any suggestions how to best fix Thetis to deal with this
• is there any chance domain-less Constants could be supported again in the adjoint? I feel this change kind of breaks with the philosophy that the adjoint should be automatic regardless of the forward model. If it is too hard to do so, why support domain-less Constants at all? Not that I would be very happy if they got dropped altogether, but at least it would tell me that the forward model I develop including all its functionality that I've tested, should still run in principle with the adjoint/tlm/etc.
• currently I'm just adding a domain=... argument when creating Constants, but that now gives a user warning "FutureWarning: Giving Constants a domain is no longer supported. Instead please create a Function in the Real space" - so should I be doing that instead? This seems even more intrusive - and I don't really understand why it couldn't just wrap "Constant(value, domain=mesh)" to return a Function on the Real space. Of course I would have to change my isinstance(..., Constant) instances in that case.

[connorjward]

• is there any chance domain-less Constants could be supported again in the adjoint? I feel this change kind of breaks with the philosophy that the adjoint should be automatic regardless of the forward model. If it is too hard to do so, why support domain-less Constants at all? Not that I would be very happy if they got dropped altogether, but at least it would tell me that the forward model I develop including all its functionality that I've tested, should still run in principle with the adjoint/tlm/etc.

I think the main argument against Constants without a domain is that they are poorly mathematically defined. For things like computing derivatives involving Constants you really want to have a domain/function space. This leads to lot of special-purpose code that would go away if functions in the real space were used instead.

Unfortunately one instance where we cannot currently use real-space functions is with vector-valued constants. This would be quite involved to fix but I think we know what needs to be done.

The way I see it is that a domain-less Constant is effectively an incomplete real-space function. As such they are effectively syntactic sugar. Perhaps we could keep this syntax sugar by parsing forms to replace Constants with the right Functions when a user calls solve, assemble etc. Alternatively, a domain-less Constant is just a modifiable parameter (ufl.ConstantValue) in which case you shouldn't be able to differentiate w.r.t. it.

• currently I'm just adding a domain=... argument when creating Constants, but that now gives a user warning "FutureWarning: Giving Constants a domain is no longer supported. Instead please create a Function in the Real space" - so should I be doing that instead? This seems even more intrusive - and I don't really understand why it couldn't just wrap "Constant(value, domain=mesh)" to return a Function on the Real space. Of course I would have to change my isinstance(..., Constant) instances in that case.

We are already returning a function in the real space when you call Constant(x, domain=mesh).

[dham]

That assign issue is a bug which should be fixed.

The problem we have with Constants with no domain is that their parallel semantics don't make sense. What we were doing was just assuming that they were collective over COMM_WORLD, but that breaks under ensemble parallelism. So we need Constants to at least know what communicator they are on.

Now, if Constants are just user parameters and don't turn up as results of computations then I think the adjoint ought to be OK. Can you specify what the failure mode is?

[stephankramer]

Thanks guys. I've submitted the assign issue separately as #3009

The errors so far have mostly been copying Constants and in constant-to-constant assignments. The failure mode is that in the replay it hits a check in isconstant() (https://github.com/dolfin-adjoint/pyadjoint/blob/0378c81517626b8e256e9ac6dfd98c1925c14701/dolfin_adjoint_common/compat.py#L167) which get striggered anytime a Constant (which are now indeed just the domain-less Constants) is the source or target of an assignment.

The thing is we have always used Constants in Thetis as a way to specify parameters which may or may not be changed by the user through callbacks (e.g. to have a time-varying boundary condition). To deal with this in timesteppers we make copies of each of the specified constants and retain a "beginning of timestep" value, which gets reassigned from the potentially updated Constant at the beginning of the timestep. So probably what I should do is to tell users that if they use Constants that they want to update/reassign values to, they have to use Constants with a domain - aka Functions on the R space. Any other values, user set or defaults, that don't change can just remain Constants, and I should avoid making copies of these.

I'll try that approach for now and see where it gets me. I'm chasing some (possibly unrelated) other adjoint bugs as well at the same time....