Adjoint of a form with TrialFunction & LinearSolver for a form with Function instead

variational formulation
1

Hello,

I aim to compute the sensitivity of a magnetostatic problem using the adjoint method.

Basically, my plan is to merge those 2 tutorials:

However, if I define my problem using a TrialFunction (magnetostatic) instead of a Function (adjoint)

J = Function(DG0)
mu = Function(DG0)
lmbda = Function(P1)
Az = TrialFunction(P1)
v = TestFunction(P1)
R = (1./mu) * dot( grad(Az), grad(v) ) * dx - J*v*dx

I get an error when I compute the adjoint

    dCostdmu = form(action(adjoint(derivative(R, mu)), lmbda))
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/path/2/ufl/formoperators.py", line 176, in adjoint
    return compute_form_adjoint(form, reordered_arguments)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/path/2/ufl/algorithms/formtransformations.py", line 458, in compute_form_adjoint
    arguments = form.arguments()
                ^^^^^^^^^^^^^^^^
  File "/path/2/ufl/form.py", line 102, in arguments
    self._analyze_form_arguments()
  File "/path/2/ufl/form.py", line 618, in _analyze_form_arguments
    arguments, coefficients = extract_arguments_and_coefficients(self)
                              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/path/2/ufl/algorithms/analysis.py", line 218, in extract_arguments_and_coefficients
    raise ValueError(
ValueError: Found different Arguments with same number and part.
Did you combine test or trial functions from different spaces?
The Arguments found are:
  v_1
  v_1
  v_0

If I use Function, it works but then I cannot solve the forward problem R using a LinearSolver(magnetostatic), but I have to define a NonlinearProblem (adjoint), despite the fact that R is linear.

Can you please explain me how to proceed such that I can solve my (linear) forward problem using a linear solver instead of a Newton scheme?

2

Please provide a minimal example that can reproduce the error, as it is quite some work to guess how you have exactly modified the code.

Please note that if you write the problem as a bilinear form,
where you have

a = (1./mu) * dot( grad(Az), grad(v) ) * dx
L = J*v*dx

you should create a residual form R after solving to compute the sensitivity, along the lines of

Azh = dolfinx.fem.Function(P1)
problem = dolfinx.fem.petsc.LinearProblem(a, L, u=Azh, ....)
R = ufl.replace(a-L, {Az: Azh})

and then use the same code as provided in

To get the correct residual derivatives.

3

Thanks, it does the job :slight_smile:

Btw, can you explain me when one should use TrialFunction and Function?
Or give me some reference explaining those (I haven’t found any yet).

Also, when should we use derivative(fun, var, trial) instead of derivative(fun, var)?

Thank you very much for your patience.

4

See: fem.Function vs ufl.TrialFunctions - #2 by dokken

It depends on what you want to achieve.
The third argument is what is inserted as the perturbation variable after differentiation.

5

Thanks for the info, I understand now the use of Function vs TrialFunction.

I guess the mathematical sense of derivative(fun, var, trial), but what puzzles me is that the tutorial about optimal control seems to behave the same w/ or w/o the third argument in derivative(...).

Here is the output w/o

$ python optimal_control.py 
J: 0.00022291602064715074
J: 9.216229129558542e-05
J: 9.008635501931439e-05
         Current function value: 0.000090
         Iterations: 3
         Function evaluations: 58
         Gradient evaluations: 46
opt_sol= message: Desired error not necessarily achieved due to precision loss.
 success: False
  status: 2
     fun: 9.008635501931439e-05
       x: [ 2.003e-02  2.023e-02 ...  2.023e-02  2.003e-02]
     nit: 3
     jac: [ 1.001e-09  9.999e-10 ...  9.999e-10  1.001e-09]
    nfev: 58
    njev: 46
Error: 0.00
(2) Error: [1.3508246588849946e-07, 3.377070371463124e-08, 8.442719550053301e-09, 2.1107017052109507e-09, 5.276863450313427e-10]
(2) Perturbance: [1.0, 0.5, 0.25, 0.125, 0.0625]
(2) Convergence rate: [1.9999962729679857, 1.9999925459405403, 1.9999850872121, 1.999970147875279]

As a reminder, here is what I got w/

$ python optimal_control.py 
J: 0.00025586196254962846
J: 0.00012317361886764953
J: 9.965292792495079e-05
J: 9.197593208180279e-05
J: 9.058358918936414e-05
         Current function value: 0.000091
         Iterations: 5
         Function evaluations: 76
         Gradient evaluations: 76
opt_sol= message: Desired error not necessarily achieved due to precision loss.
 success: False
  status: 2
     fun: 9.058358918936414e-05
       x: [ 5.439e-02  5.496e-02 ...  5.496e-02  5.439e-02]
     nit: 5
     jac: [ 8.990e-09  8.962e-09 ...  8.962e-09  8.990e-09]
    nfev: 76
    njev: 76
Error: 0.00
(2) Error: [1.3498198327760653e-07, 3.3745463081967856e-08, 8.436349415029447e-09, 2.1090791812234977e-09, 5.272657078802013e-10]
(2) Perturbance: [1.0, 0.5, 0.25, 0.125, 0.0625]
(2) Convergence rate: [2.000001399598824, 2.0000027969350995, 2.0000055903312397, 2.0000111838979926]
6

See the definition of the derivative:
https://docs.fenicsproject.org/ufl/2024.2.0/api-doc/ufl.html#ufl.derivative

If the argument is omitted, a new Argument is created in the same space as the coefficient, with argument number one higher than the highest one in the form.

7

I then assume the tutorial is verbose for completeness.

Thank you for your time.