I am trying to solve an optimization problem on a compound space Z. I have some functions defined on a 2D domain \Omega and one more function (the control variable) defined on \partial \Omega.
I have defined my problem in FEniCS using MixedFunctionSpace from mixed-dimensional branch:
Z = MixedFunctionSpace(V, VT, U)
Then I define a variational problem:
z = Function(Z)
(v, vt, u) = z.split()
L = inner(nabla_grad(v), nabla_grad(vt))*dx + ...
Now I want to calcualate the derivatives of L w.r.t. all variables. By analogy to the stable version of FEniCS, I would write it like this:
You should be able to define the NonlinearVariationalProblem as:
z = Function(Z)
dz=TrialFunction(Z)
Jac = derivative(L, z, dz) #Define the Jacobian
problem = NonlinearVariationalProblem(L, z, bcs,Jac) #Define the Problem
solver = NonlinearVariationalSolver(problem)
Indeed, derivative doesn’t work on a tuple.
The system you are solving is block shaped, and derivative should be used per block, w.r.t to only one variable.
When using MixedNonlinearVariationalProblem, both the form and the Jacobian are defined as lists corresponding to the decomposition per block.
To obtain the corresponding list of subforms, you can do
L_blocks = extract_blocks(L)
and derivative can then be used on each subforms to build the list of blocks for the Jacobian :
Js = []
for Li in L_blocks :
for zi in z.split():
d = derivative(Li, zi)
Js.append(d)
You don’t need any extra modules. My last commit in the mixed-dimensional branch should fix this, so you just need to update the branch.
I am not quite sure I’ve fully understood which problem you are trying to solve but as far as I can tell, there is a MixedNonlinearVariationalSolver you can use with the previously defined MixedNonlinearVariationalProblem
It seems to me that you are not using your variable problem = MixedNonlinearVariationalProblem(...).
When you use solve(...) without specifying the solver, it will build the appropriate problem and solver objects depending on the type of form you give. For mixed-dimensional problems whose function space is a MixedFunctionSpace, the block extraction is done seamlessly from the given form composed of all the blocks. I guess you could sum the subforms of Js list and give it to solve(...), but I would rather recommend to use MixedNonlinearVariationalSolver from your problem
Indeed, the last build failed due to recent changes in UFL, but it does not impact the framework.
However, it did prevent the update of the Docker container (so my very last changes had not been included). I have just updated the container with the latest version of the mixed-dimensional branches.
(Sorry for the incovenience)
@ cdaversin thanks for updating the container, now MixedNonlinearVariationalProblem is visible. Further, I passed the Jacobian and the Hessian into it but ran into the following error:
L_blocks = extract_blocks(L)
Js = []
for Li in L_blocks:
for zi in z.split():
d = derivative(Li, zi)
Js.append(d)
Hes = []
for Js_i in Js:
for zi in z.split():
d = derivative(Js_i, zi)
Hes.append(d)
problem = MixedNonlinearVariationalProblem(Js, z.split(), bcs, Hes)
Hello s_ony and cdaversin. Can you explain why a Function is used for the derivative and not the Test and Trial functions? I.e. Why is it not with Test for the Js and Trial for the Hes ? eg like below
ztest=TestFunctions(Z)
ztrial=TrialFunctions(Z)
L_blocks = extract_blocks(L)
Js =sum([ [derivative(Li, zi) for zi in ztest.split() ] for Li in L_blocks], [])
Hes=sum([ [derivative(Li, zi) for zi in ztrial.split() ] for Js_i in Js], [])
I think it’s essential to define w.r.t. which coefficient (or tuple of coefficients) the derivative of form is taken. Additionally one can explicitly pass test and trial function as the third parameter in derivative to provide the direction in which the derivative should be taken. However, if it’s not provided, the new arguments are going to be created for you.
I would suggest taking a look at Automatic functional differentiation section in UFL Documentation.
I see now that the “derivative” function must always be told what coefficient to differentiate w.r.t - because a form (e.g. the one in you example) can use more than one Function object, but for the direction it just uses the ufl convention for “order” of Arguments (Test then Trial) to work out which is the lowest order of Argument missing and then uses that direction.
Thanks also for the pointer to the manual, I also found looking at the definition of “derivative” and “_handle_derivative_arguments” in uf/formoperators.py helped - the code shows that the function is able to deduce an appropriate direction from the form if no argument is given .
I can’t quite see how for a mixedmesh “_handle_derivative_arguments” manages to match the correct sub function of the test or trial function with the appropriate sub function of the “z” Coefficient, but I guess its not necessary; the code just differentiates w.r.t to the Coefficient in the direction of the the total Argument and that gives the same end result.
So I thinking that these two version below are equivalent
residuals =sum([ [df.derivative(Li,zi, zdir) for zi,zdir in zip(z.split(),ztest.split()) ] for Li in L_blocks], [])
jacobians =sum([ [df.derivative(Ri,zi, zdir) for zi,zdir in zip(z.split() ztrial.split()) ] for Ri in residuals], [])
residuals =sum([ [df.derivative(Li,zi) for zi in z.split() ] for Li in L_blocks], [])
jacobians =sum([ [df.derivative(Ri,zi) for zi in z.split() ] for Ri in residuals], [])
You don’t need to reply to that - you have given me enough information to work it out by doing a few tests.