diff --git a/python/dune/perftool/pdelab/adjoint.py b/python/dune/perftool/pdelab/adjoint.py
index 27786e3017a582a16d24ec1cdda08086ffd98ddb..cbba97984658c8e1169cbd2c546122ad1a4d78f0 100644
--- a/python/dune/perftool/pdelab/adjoint.py
+++ b/python/dune/perftool/pdelab/adjoint.py
@@ -48,14 +48,20 @@ def define_dJdm_member(name):
def generate_accumulation_instruction(expr, visitor, accumulation_index, number_of_controls):
+ # Create class member dJdm for accumulating
accumvar = "dJdm"
- assignee = prim.Subscript(prim.Variable(accumvar), accumulation_index)
shape = (number_of_controls,)
define_dJdm_member(accumvar)
+
+ # Tell loopy about
globalarg(accumvar, shape=shape)
- quad_inames = visitor.interface.quadrature_inames()
- from dune.perftool.generation import instruction
+ assignee = prim.Subscript(prim.Variable(accumvar), accumulation_index)
+
+ # We need to accumulate
expr = prim.Sum((assignee, expr))
+
+ from dune.perftool.generation import instruction
+ quad_inames = visitor.interface.quadrature_inames()
instruction(assignee=assignee,
expression=expr,
forced_iname_deps=frozenset(quad_inames),
@@ -65,8 +71,21 @@ def generate_accumulation_instruction(expr, visitor, accumulation_index, number_
def list_accumulation_infos(expr, visitor):
return ["control",]
-class AdjointInterface(PDELabInterface):
+class ControlInterface(PDELabInterface):
+ """Interface for generating the control localoperator
+
+ In this case we will not accumulate in the residual vector but use
+ a class member representing dJdm instead.
+
+ """
def __init__(self, accumulation_index, number_of_controls):
+ """Create ControlInterface
+
+ Arguments:
+ ----------
+ accumulation_index: In which component of the dJdm should be accumulated.
+ number_of_controls: Number of components of dJdm. Needed for creating the member variable.
+ """
self.accumulation_index = accumulation_index
self.number_of_controls = number_of_controls
@@ -81,7 +100,7 @@ class AdjointInterface(PDELabInterface):
def get_visitor(measure, subdomain_id, accumulation_index, number_of_controls):
- interface = AdjointInterface(accumulation_index, number_of_controls)
+ interface = ControlInterface(accumulation_index, number_of_controls)
from dune.perftool.ufl.visitor import UFL2LoopyVisitor
return UFL2LoopyVisitor(interface, measure, subdomain_id)
@@ -91,12 +110,18 @@ def visit_integral(integral, accumulation_index, number_of_controls):
measure = integral.integral_type()
subdomain_id = integral.subdomain_id()
- # Start the visiting process!
+ # The visitor needs to know about the current index and the number
+ # of controls in order to generate the accumulation instruction
visitor = get_visitor(measure, subdomain_id, accumulation_index, number_of_controls)
+
+ # Start the visiting process!
visitor.accumulate(integrand)
def generate_kernel(forms):
+ # Similar to the standard residual generation, except:
+ # - Have multiple forms
+ # - Pass index and number of forms along
logger = logging.getLogger(__name__)
# Visit all integrals once to collect information (dry-run)!
@@ -132,7 +157,7 @@ def generate_kernel(forms):
# @backend(interface="generate_kernels_per_integral")
-def adjoint_generate_kernels_per_integral(forms):
+def control_generate_kernels_per_integral(forms):
"""For the control problem forms will have one form for every
measure. Every form will only contain integrals of one type.
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index ce685263ec2b639d45dea108dc1604886a6aabec..31cc63ee395e79e57644d0abbd85f9b9bba60039 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -879,6 +879,9 @@ def generate_jacobian_kernels(form, original_form):
def generate_control_kernels(forms):
+ # All forms will we written in the residual method and
+ # accumulation will be done in a class member instead of the
+ # residual.
logger = logging.getLogger(__name__)
with global_context(form_type='residual'):
operator_kernels = {}
@@ -893,11 +896,12 @@ def generate_control_kernels(forms):
from dune.perftool.pdelab.signatures import assembler_routine_name
with global_context(kernel=assembler_routine_name()):
- # TODO: make backend switch for PDELab/sumfact
- from dune.perftool.pdelab.adjoint import adjoint_generate_kernels_per_integral
+ # TODO: Sumfactorization not yet implemented
+ assert not get_form_option('sumfact')
+ from dune.perftool.pdelab.adjoint import control_generate_kernels_per_integral
forms_measure = [form.integrals_by_type(measure) for form in forms]
- kernel = [k for k in adjoint_generate_kernels_per_integral(forms_measure)]
+ kernel = [k for k in control_generate_kernels_per_integral(forms_measure)]
operator_kernels[(measure, 'residual')] = kernel
@@ -954,23 +958,33 @@ def generate_localoperator_kernels(operator):
# Generate control operator
#
# This is the normal form derived w.r.t. the control
- # variable.
- assert get_form_option("control_variable") is not None
- controls = [data.object_by_name[ctrl.strip()] for ctrl in get_form_option("control_variable").split(",")]
- assert len(controls) == 1
-
+ # variable. We generate a form for every row of:
+ #
+ # \nabla \hat{J}(m) = (\nabla R(z,m))^T \lambda + \nabla_m J(z,m)
+ #
+ # These forms will not depend on the test function anymore and
+ # will need special treatment for the accumulation process.
from ufl import action, diff
from ufl.classes import Coefficient
- # We need to transform numpy ints to python native ints
+ # Get control variables
+ assert get_form_option("control_variable") is not None
+ controls = [data.object_by_name[ctrl.strip()] for ctrl in get_form_option("control_variable").split(",")]
+
+ # Transoform flat index to multiindex. Wrapper around numpy
+ # unravel since we need to transform numpy ints to native
+ # ints.
def _unravel(flat_index, shape):
multi_index = np.unravel_index(flat_index, shape)
multi_index = tuple(int(i) for i in multi_index)
return multi_index
- forms = []
+ # Will be used to replace ansatz function with adjoint function
element = original_form.coefficients()[0].ufl_element()
coeff = Coefficient(element, count=3)
+
+ # Store a form for every control
+ forms = []
for control in controls:
shape = control.ufl_shape
flat_length = int(np.prod(shape))