diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py index 067d6da638b6c42101ace4ea2910bd66c0bd659b..5041772548f6646a9c1de40f0d6e929078900faa 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 6c400cc4db61472abef9e0d39e953eb9e1e8d1e0..d536027362ee2ef347ad1847934c50fb94cf3fc4 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 4873920b2901f812bbe06e01f233c0f9d3423630..cc8aa4e117c263fbceb2a8f310f558b3f8f2fbca 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 aa9d7714a2ae60da6fd37911eced5f15af259d46..e208de92b907519763f4131bcfacfe6bee796192 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 #