Simplify interpolate
#4582
Simplify interpolate
#4582
Conversation
| def test_assemble_interp_operator(V2, f1): | ||
| # Check type | ||
| If1 = Interpolate(f1, V2) | ||
| If1 = Interpolate(f1, Argument(V2.dual(), 0)) |
There was a problem hiding this comment.
| If1 = Interpolate(f1, Argument(V2.dual(), 0)) | |
| If1 = Interpolate(f1, TestFunction(V2.dual())) |
leo-collins
force-pushed
the
leo/simplify-interpolate
branch
2 times, most recently
from
92645d1 to
aa30bf8
Compare
firedrake/interpolation.py
Outdated
| expr_args = extract_arguments(expr) | ||
| if isinstance(V, functionspaceimpl.WithGeometry): | ||
| # Need to create a Firedrake Argument so that it has a .function_space() method | ||
| expr_args = extract_arguments(ufl.as_ufl(expr)) |
There was a problem hiding this comment.
| expr_args = extract_arguments(ufl.as_ufl(expr)) | |
| if isinstance(expr, ufl.BaseForm): | |
| expr_args = [a for a in expr.arguments() if not isinstance(a, Coargument)] | |
| else: | |
| expr_args = extract_arguments(ufl.as_ufl(expr)) |
There was a problem hiding this comment.
Not sure if filtering all Coarguments is a sensible thing, but imagine if a Function were to advertise its Coargument (technically it should have one), we would have to discard it. An Interpolate is a BaseFormOperator that does advertise the Coargument, but here we should only view the operand as the output after assembling it.
There was a problem hiding this comment.
Maybe the more sensible thing is to drop the last Coargument?
| expr_args = extract_arguments(ufl.as_ufl(expr)) | ||
| is_adjoint = len(expr_args) and expr_args[0].number() == 0 | ||
| v = Argument(v.dual(), 1 if is_adjoint else 0) | ||
| V = Argument(V.dual(), 1 if is_adjoint else 0) |
There was a problem hiding this comment.
Put 0 if it is not there, otherwise keep increasing the count
| V = Argument(V.dual(), 1 if is_adjoint else 0) | |
| arg_numbers = [a.number() for a in expr_args] | |
| number = 0 if len(expr_args) == 0 or min(arg_numbers) > 0 else max(arg_numbers) + 1 | |
| V = Argument(V.dual(), number) |
There was a problem hiding this comment.
Similar fixes will be needed in UFL
There was a problem hiding this comment.
In UFL we're checking for contiguous argument numbers, I think this approach is more fool-proof: just get the smallest number that's missing
| V = Argument(V.dual(), 1 if is_adjoint else 0) | |
| arg_numbers = set(a.number() for a in expr_args) | |
| number = min(set(range(len(expr_args) + 1)) - arg_numbers) | |
| V = Argument(V.dual(), number) |
There was a problem hiding this comment.
expr_args can only contain zero or one arguments so I think this logic is more confusing than is_adjoint
| class Interpolate(ufl.Interpolate): | ||
|
|
||
| def __init__(self, expr, v, | ||
| def __init__(self, expr, V, |
There was a problem hiding this comment.
This is very unrelated, but I think that a much more friendly interface is to allow either or both left and right arguments to be a primal FunctionSpace.
Right now we do this under the hood
Interpolate(Function(V1), V2) -> Interpolate(Function(V1), Argument(V2.dual(), 0))
It'd be reasonable to have a similar shortcut for the adjoint. When the left argument is a FunctionSpace, we would then automatically create the Argument for it.
Interpolate(V1, Cofunction(V2.dual())) -> Interpolate(Argument(V1, 0), Cofunction(V2.dual()))
And supplying two FunctionSpaces is a perfectly natural interface:
Interpolate(V1, V2) -> Interpolate(Argument(V1, 1), Argument(V2.dual(), 0))
Of course we need to arbitrarily decide who gets the lowest number, the more intuitive numbering that produces the forward Interpolation is to go from right to left.
leo-collins
force-pushed
the
leo/simplify-interpolate
branch
from
4b21ac6 to
c8f2b37
Compare
test -> trial fix hypre-ads
cofunction docstring
fix
leo-collins
force-pushed
the
leo/simplify-interpolate
branch
from
10fb533 to
f8c2318
Compare
Simplify
interpolateandInterpolate. This removes the argument renumbering 0 -> 1 in the primal expression if V is a FunctionSpace. Hence this is an API change and should be merged after the next release. PR #4572 warns users about this.