From 04e5a5309983658152e0d0b193a742524b50d290 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, 2 Feb 2018 17:40:55 +0100
Subject: [PATCH] Generate adjoint operators (only scalar controls!)
---
python/dune/perftool/options.py | 7 ++-
python/dune/perftool/pdelab/geometry.py | 7 ++-
python/dune/perftool/pdelab/localoperator.py | 55 ++++++++++++++++++--
python/dune/perftool/ufl/visitor.py | 7 +++
4 files changed, 69 insertions(+), 7 deletions(-)
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 067d6da6..50417725 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -91,7 +91,10 @@ class PerftoolFormOptionsArray(ImmutableRecord):
precompute_quadrature_info = PerftoolOption(default=True, helpstr="compute quadrature points and weights in the constructor of the local operator")
blockstructured = PerftoolOption(default=False, helpstr="Use block structure")
number_of_blocks = PerftoolOption(default=1, helpstr="Number of sub blocks in one direction")
-
+ adjoint = PerftoolOption(default=False, helpstr="Generate adjoint operator")
+ control = PerftoolOption(default=False, helpstr="Generate operator of derivative w.r.t. the control variable")
+ objective_function = PerftoolOption(default=None, helpstr="Name of form representing the objective function in UFL file")
+ control_variable = PerftoolOption(default=None, helpstr="Name of control variable in UFL file")
# Until more sophisticated logic is needed, we keep the actual option data in this module
_global_options = PerftoolGlobalOptionsArray()
@@ -164,7 +167,7 @@ def process_form_options(opt, form):
opt = opt.copy(form=form)
if opt.classname is None:
- opt = opt.copy(classname="{}Operator".format(opt.form))
+ opt = opt.copy(classname="{}Operator".format(form))
if opt.filename is None:
opt = opt.copy(filename="{}_{}_file.hh".format(get_option("target_name"), opt.classname))
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 6c400cc4..d5360273 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -231,8 +231,11 @@ def to_cell_coordinates(local, restriction):
def world_dimension():
- from dune.perftool.pdelab.driver import get_dimension
- return get_dimension()
+ data = get_global_context_value("data")
+ form = data.object_by_name[get_form_option("form")]
+ from dune.perftool.ufl.preprocess import preprocess_form
+ form = preprocess_form(form).preprocessed_form
+ return form.ufl_cell().geometric_dimension()
def intersection_dimension():
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 4873920b..cc8aa4e1 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -712,9 +712,58 @@ def generate_localoperator_kernels(operator):
data = get_global_context_value("data")
original_form = data.object_by_name[get_form_option("form")]
-
from dune.perftool.ufl.preprocess import preprocess_form
- form = preprocess_form(original_form).preprocessed_form
+
+ if get_form_option("adjoint"):
+ # Generate adjoint operator
+ #
+ # The jacobian of the adjoint form is just the jacobian of the
+ # original form with test and ansazt function swapped. A a
+ # linear form you have to subtract the derivative of the
+ # objective function w.r.t the ansatz function to get the
+ # final residual formulation of the adjoint.
+ #
+ # TODO: This might only be true for linear problems. Adjust
+ # documentation as knowledge improves ;)
+ assert get_form_option("objective_function") is not None
+ assert get_form_option("control") is False
+
+ from ufl import derivative, adjoint, action, replace
+ from ufl.classes import Coefficient
+
+ jacform = derivative(original_form, original_form.coefficients()[0])
+ adjoint_jacform = adjoint(jacform)
+ adjoint_form = action(adjoint_jacform, original_form.coefficients()[0])
+
+ objective = data.object_by_name[get_form_option("objective_function")]
+ jacobian_objective = derivative(objective, objective.coefficients()[0])
+
+ element = objective.coefficients()[0].ufl_element()
+ coeff = Coefficient(element)
+ jacobian_objective = replace(jacobian_objective, {objective.coefficients()[0]: coeff})
+
+ adjoint_form = adjoint_form + jacobian_objective
+
+ form = preprocess_form(adjoint_form).preprocessed_form
+
+ elif get_form_option("control"):
+ # Generate control operator
+ #
+ # This is the normal form derived w.r.t. the control
+ # variable. Right now this is olny implemented for scalar
+ # controls.
+ assert get_form_option("control_variable") is not None
+
+ control = data.object_by_name[get_form_option("control_variable")]
+ assert control.ufl_shape is ()
+
+ from ufl import diff
+ control_form = diff(original_form, control)
+ form = preprocess_form(control_form).preprocessed_form
+
+ else:
+ form = preprocess_form(original_form).preprocessed_form
+
# Reset the generation cache
from dune.perftool.generation import delete_cache_items
@@ -749,7 +798,7 @@ def generate_localoperator_kernels(operator):
# Set some options!
from dune.perftool.pdelab.driver import isQuadrilateral
- if isQuadrilateral(form.coefficients()[0].ufl_element().cell()):
+ if isQuadrilateral(form.arguments()[0].ufl_element().cell()):
from dune.perftool.options import set_form_option
# For Yasp Grids the jacobian of the transformation is diagonal and constant on each cell
set_form_option('diagonal_transformation_matrix', True)
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index aa9d7714..e208de92 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -169,6 +169,13 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return self.interface.pymbolic_gridfunction(o, restriction, self.grad)
+
+ def variable(self, o):
+ # Right now only scalar varibables are supported
+ assert o.ufl_shape is ()
+ return o.expression().value()
+
+
#
# Handlers for all indexing related stuff
#
--
GitLab