From 7868fe686ce109395fa2be0076571c500334bf44 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Fri, 9 Feb 2018 11:20:37 +0100
Subject: [PATCH] Fix bug in derivative of adjoint form

---
 python/dune/perftool/pdelab/localoperator.py | 26 ++++++++++++++------
 1 file changed, 19 insertions(+), 7 deletions(-)

diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index cc8aa4e1..8b3e012b 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -731,19 +731,26 @@ def generate_localoperator_kernels(operator):
         from ufl import derivative, adjoint, action, replace
         from ufl.classes import Coefficient
 
+        # Jacobian of the adjoint form
         jacform = derivative(original_form, original_form.coefficients()[0])
         adjoint_jacform = adjoint(jacform)
-        adjoint_form = action(adjoint_jacform, original_form.coefficients()[0])
 
+        # Derivative of objective function w.r.t. state
         objective = data.object_by_name[get_form_option("objective_function")]
-        jacobian_objective = derivative(objective, objective.coefficients()[0])
+        objective_jacobian = derivative(objective, objective.coefficients()[0])
 
-        element = objective.coefficients()[0].ufl_element()
-        coeff = Coefficient(element)
-        jacobian_objective = replace(jacobian_objective, {objective.coefficients()[0]: coeff})
+        # Replace coefficient belonging to ansatz function with new coefficient
+        element = adjoint_jacform.coefficients()[0].ufl_element()
+        coeff = Coefficient(element, count=3)
+        adjoint_jacform = replace(adjoint_jacform, {adjoint_jacform.coefficients()[0]: coeff})
+        objective_jacobian = replace(objective_jacobian, {objective.coefficients()[0]: coeff})
 
-        adjoint_form = adjoint_form + jacobian_objective
+        # Residual of the adjoint form
+        adjoint_form = action(adjoint_jacform, original_form.coefficients()[0])
+        adjoint_form = adjoint_form + objective_jacobian
 
+        # Update form and original_form
+        original_form = adjoint_form
         form = preprocess_form(adjoint_form).preprocessed_form
 
     elif get_form_option("control"):
@@ -757,8 +764,13 @@ def generate_localoperator_kernels(operator):
         control = data.object_by_name[get_form_option("control_variable")]
         assert control.ufl_shape is ()
 
-        from ufl import diff
+        from ufl import diff, replace
+        from ufl.classes import Coefficient
         control_form = diff(original_form, control)
+        # element = control_form.coefficients()[0].ufl_element()
+        # coeff = Coefficient(element, count=3)
+        # control_form = replace(control_form, {control_form.coefficients()[0]: coeff})
+        original_form = control_form
         form = preprocess_form(control_form).preprocessed_form
 
     else:
-- 
GitLab