Dear all, I'm having a variable Sw, which I calculate pointwise in a routine calc_Sw (see attachment). This variable depends on my primary unknown h: Sw = (h-hn)/(hw-hn) (hn,hw known values). But now it seems like automatic differentiation for this routine does not work. I'm getting zero entries for the Jacobian, whereas d(Sw)/dh = 1.0 should hold. Do I have to write this routine with conditionals? I tried this, but failed. How would my routine calc_Sw have to look like with conditionals? Thank you! Henrik -- Dipl.-Math. Henrik Büsing Institute for Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University ------------------------------------------------------ Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889 ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing@eonerc.rwth-aachen.de ------------------------------------------------------
On 5 Nov 2015, at 09:49, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de> wrote:
Dear all,
I’m having a variable Sw, which I calculate pointwise in a routine calc_Sw (see attachment). This variable depends on my primary unknown h: Sw = (h-hn)/(hw-hn) (hn,hw known values).
But now it seems like automatic differentiation for this routine does not work. I’m getting zero entries for the Jacobian, whereas d(Sw)/dh = 1.0 should hold.
Yes, this is because the AD doesn't know about the relationship between Sw and dh. UFL has a facility for this, but I notice we don't expose it in firedrake (however, it is straightforward to do): We want derivative(F, u) But F contains a coefficient, S, whose derivative wrt u is 1.0 (however, they are not symbolically related in a way UFL understands). So we build a mapping from this coefficient to its derivative wrt u: coefficient_derivatives = {S: 1.0} and then pass this additional information to the derivative call. derivative(F, u, coefficient_derivatives=coefficient_derivatives) Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through. If this works for you, do you want to propose a patch that adds this functionality? Cheers, Lawrence
On 05/11/15 10:08, Lawrence Mitchell wrote:
On 5 Nov 2015, at 09:49, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de> wrote:
Dear all,
I’m having a variable Sw, which I calculate pointwise in a routine calc_Sw (see attachment). This variable depends on my primary unknown h: Sw = (h-hn)/(hw-hn) (hn,hw known values).
But now it seems like automatic differentiation for this routine does not work. I’m getting zero entries for the Jacobian, whereas d(Sw)/dh = 1.0 should hold.
Yes, this is because the AD doesn't know about the relationship between Sw and dh. UFL has a facility for this, but I notice we don't expose it in firedrake (however, it is straightforward to do):
We want
derivative(F, u)
But F contains a coefficient, S, whose derivative wrt u is 1.0 (however, they are not symbolically related in a way UFL understands). So we build a mapping from this coefficient to its derivative wrt u:
coefficient_derivatives = {S: 1.0}
and then pass this additional information to the derivative call.
derivative(F, u, coefficient_derivatives=coefficient_derivatives)
Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through. If this works for you, do you want to propose a patch that adds this functionality?
Cheers,
Lawrence
Are you sure the problem isn't that he needs bounds for his solve? AFAICS dS/dh=0 for h<hw or h>hn. Otherwise shouldn't rewriting this with conditionals have worked? Henrik, I presume you are looking for a solution hw<h<hn ? Cheers Stephan
On 5 Nov 2015, at 10:26, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
On 05/11/15 10:08, Lawrence Mitchell wrote:
On 5 Nov 2015, at 09:49, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de> wrote:
Dear all,
I’m having a variable Sw, which I calculate pointwise in a routine calc_Sw (see attachment). This variable depends on my primary unknown h: Sw = (h-hn)/(hw-hn) (hn,hw known values).
But now it seems like automatic differentiation for this routine does not work. I’m getting zero entries for the Jacobian, whereas d(Sw)/dh = 1.0 should hold.
Yes, this is because the AD doesn't know about the relationship between Sw and dh. UFL has a facility for this, but I notice we don't expose it in firedrake (however, it is straightforward to do):
We want
derivative(F, u)
But F contains a coefficient, S, whose derivative wrt u is 1.0 (however, they are not symbolically related in a way UFL understands). So we build a mapping from this coefficient to its derivative wrt u:
coefficient_derivatives = {S: 1.0}
and then pass this additional information to the derivative call.
derivative(F, u, coefficient_derivatives=coefficient_derivatives)
Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through. If this works for you, do you want to propose a patch that adds this functionality?
Cheers,
Lawrence
Are you sure the problem isn't that he needs bounds for his solve? AFAICS dS/dh=0 for h<hw or h>hn. Otherwise shouldn't rewriting this with conditionals have worked? Henrik, I presume you are looking for a solution hw<h<hn ?
I didn't think hard about this, I simply considered that there's a functional relationship which UFL doesn't know about if you just pass in unrelated coefficients. Lawrence
-----Ursprüngliche Nachricht----- Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake- bounces@imperial.ac.uk] Im Auftrag von Stephan Kramer Gesendet: 05 November 2015 11:27 An: firedrake@imperial.ac.uk Betreff: Re: [firedrake] Problem with Jacobian
On 05/11/15 10:08, Lawrence Mitchell wrote:
On 5 Nov 2015, at 09:49, Buesing, Henrik <HBuesing@eonerc.rwth-
aachen.de> wrote:
Dear all,
I'm having a variable Sw, which I calculate pointwise in a routine
calc_Sw (see attachment). This variable depends on my primary unknown h: Sw = (h-hn)/(hw-hn) (hn,hw known values).
But now it seems like automatic differentiation for this routine does
not work. I'm getting zero entries for the Jacobian, whereas d(Sw)/dh = 1.0 should hold.
Yes, this is because the AD doesn't know about the relationship between Sw and dh. UFL has a facility for this, but I notice we don't expose it in firedrake (however, it is straightforward to do):
We want
derivative(F, u)
But F contains a coefficient, S, whose derivative wrt u is 1.0 (however, they are not symbolically related in a way UFL understands). So we build a mapping from this coefficient to its derivative wrt u:
coefficient_derivatives = {S: 1.0}
and then pass this additional information to the derivative call.
derivative(F, u, coefficient_derivatives=coefficient_derivatives)
Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through. If this works for you, do you want to propose a patch that adds this functionality?
Cheers,
Lawrence
Are you sure the problem isn't that he needs bounds for his solve? AFAICS dS/dh=0 for h<hw or h>hn. Otherwise shouldn't rewriting this with conditionals have worked? [Buesing, Henrik] I failed to rewrite my function with conditionals. Just because I didn't know how to properly write down the conditional() stuff and how to set Sw. If I could rewrite this function, I assume this would work!
Henrik, I presume you are looking for a solution hw<h<hn ? [Buesing, Henrik] Yes.
Cheers Stephan
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 06/11/15 12:40, Buesing, Henrik wrote:
Are you sure the problem isn't that he needs bounds for his solve? AFAICS dS/dh=0 for h<hw or h>hn. Otherwise shouldn't rewriting this with conditionals have worked?
@Stephan: So how can I rewrite the attached saturation.py with conditionals?
Sw=conditional(hin>hn, Constant(0.0), conditional(hin<hw, Constant(1.0), (hin-hn)/(hw-hn))) This gives you a UFL expression that you can use in a larger system, or project to some function space. Cheers Stephan
On 06/11/15 12:40, Buesing, Henrik wrote:
Are you sure the problem isn't that he needs bounds for his solve? AFAICS dS/dh=0 for h<hw or h>hn. Otherwise shouldn't rewriting this with conditionals have worked?
@Stephan: So how can I rewrite the attached saturation.py with conditionals?
Sw=conditional(hin>hn, Constant(0.0), conditional(hin<hw, Constant(1.0), (hin-hn)/(hw-hn)))
This gives you a UFL expression that you can use in a larger system, or project to some function space.
@Stephan: Thank you! Is there also a way to write the attached enthn.py in UFL notation? Basically I failed to do the int(dm), if dm is a UFL expression. And probably the array access might be a problem? Thank you! Henrik
On 06/11/15 14:23, Buesing, Henrik wrote:
On 06/11/15 12:40, Buesing, Henrik wrote:
Are you sure the problem isn't that he needs bounds for his solve? AFAICS dS/dh=0 for h<hw or h>hn. Otherwise shouldn't rewriting this with conditionals have worked?
@Stephan: So how can I rewrite the attached saturation.py with conditionals?
Sw=conditional(hin>hn, Constant(0.0), conditional(hin<hw, Constant(1.0), (hin-hn)/(hw-hn)))
This gives you a UFL expression that you can use in a larger system, or project to some function space.
@Stephan: Thank you! Is there also a way to write the attached enthn.py in UFL notation? Basically I failed to do the int(dm), if dm is a UFL expression. And probably the array access might be a problem?
Yes, looks like a piecewise linear interpolation from an arbitrary lookup table? I don't see how to do this in UFL indeed. You could approximate it with a polynomial but that might not be accurate enough and have nasty under/over shoots if the data is not smooth. So the best way to do this might indeed be to keep your point-wise evaluation and add a hook for the evaluation of its derivative Cheers Stephan
On 06/11/15 14:23, Buesing, Henrik wrote:
On 06/11/15 12:40, Buesing, Henrik wrote:
Are you sure the problem isn't that he needs bounds for his solve? AFAICS dS/dh=0 for h<hw or h>hn. Otherwise shouldn't rewriting this with conditionals have worked?
@Stephan: So how can I rewrite the attached saturation.py with conditionals?
Sw=conditional(hin>hn, Constant(0.0), conditional(hin<hw, Constant(1.0), (hin-hn)/(hw-hn)))
This gives you a UFL expression that you can use in a larger system, or project to some function space.
@Stephan: Thank you! Is there also a way to write the attached enthn.py in UFL notation? Basically I failed to do the int(dm), if dm is a UFL expression. And probably the array access might be a problem?
Yes, looks like a piecewise linear interpolation from an arbitrary lookup table?
Yes, exactly!
I don't see how to do this in UFL indeed. You could approximate it with a polynomial but that might not be accurate enough and have nasty under/over shoots if the data is not smooth.
So the best way to do this might indeed be to keep your point-wise evaluation and add a hook for the evaluation of its derivative
Ok. I will see if I get the "coefficient_derivatives" argument to work. Else I just might pass the whole Jacobian. Thank you! Henrik
-----Ursprüngliche Nachricht----- Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake- bounces@imperial.ac.uk] Im Auftrag von Lawrence Mitchell Gesendet: 05 November 2015 11:08 An: firedrake@imperial.ac.uk Betreff: Re: [firedrake] Problem with Jacobian
On 5 Nov 2015, at 09:49, Buesing, Henrik <HBuesing@eonerc.rwth- aachen.de> wrote:
Dear all,
I’m having a variable Sw, which I calculate pointwise in a routine calc_Sw (see attachment). This variable depends on my primary unknown h: Sw = (h-hn)/(hw-hn) (hn,hw known values).
But now it seems like automatic differentiation for this routine does not work. I’m getting zero entries for the Jacobian, whereas d(Sw)/dh = 1.0 should hold.
Yes, this is because the AD doesn't know about the relationship between Sw and dh. UFL has a facility for this, but I notice we don't expose it in firedrake (however, it is straightforward to do):
We want
derivative(F, u)
But F contains a coefficient, S, whose derivative wrt u is 1.0 (however, they are not symbolically related in a way UFL understands). So we build a mapping from this coefficient to its derivative wrt u:
coefficient_derivatives = {S: 1.0}
and then pass this additional information to the derivative call.
derivative(F, u, coefficient_derivatives=coefficient_derivatives)
Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through. [Buesing, Henrik] So I just add
def derivative(form, u, du=None,coefficient_derivatives=None): "" return ufl.derivative(form, u, du, coefficient_derivatives) or how would this work? And then I have to reinstall firedrake?
If this works for you, do you want to propose a patch that adds this functionality?
[Buesing, Henrik] How would I do this?
Cheers,
Lawrence
Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through.
[Buesing, Henrik] @Lawrence: So how would I pass the "coefficient_derivatives" field through?
On 6 Nov 2015, at 12:37, Buesing, Henrik <HBuesing@eonerc.rwth-aachen.de> wrote:
Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through.
[Buesing, Henrik]
@Lawrence: So how would I pass the "coefficient_derivatives" field through?
If you change the definition of derivative in firedrake/ufl_expr.py to take an optional coefficient_derivatives argument: def derivative(F, u, du=None, coefficient_derivatives=None): ... And just pass that value through to the ufl.derivative call: ufl.derivative(..., coefficient_derivatives=coefficient_derivatives) If you do this in the firedrake source directly, you will be able to commit your change in the local git repository. You can send us these changes by forking firedrake on github and proposing a pull-request that way. Cheers, Lawrence
Firedrake uses the UFL derivative function, but does not expose this extra argument in the interface. It is straightforward to alter the definition in firedrake/ufl_expr.py to take this extra argument and pass it through.
[Buesing, Henrik]
@Lawrence: So how would I pass the "coefficient_derivatives" field through?
If you change the definition of derivative in firedrake/ufl_expr.py to take an optional coefficient_derivatives argument:
def derivative(F, u, du=None, coefficient_derivatives=None): ...
And just pass that value through to the ufl.derivative call:
ufl.derivative(..., coefficient_derivatives=coefficient_derivatives)
If you do this in the firedrake source directly, you will be able to commit your change in the local git repository. You can send us these changes by forking firedrake on github and proposing a pull-request that way.
Hmm... Okay I did this. Then I did a "make" in the firedrake directory. But then I'm still getting a TypeError: derivative() takes at most 3 arguments (4 given) So how can I test this? Thank you! Henrik
participants (3)
-
Buesing, Henrik
-
Lawrence Mitchell
-
Stephan Kramer