diff --git a/dune/perftool/matrixfree.hh b/dune/perftool/matrixfree.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b5bf65ff1faacc89b4a1fd46b17bfec121572ed7
--- /dev/null
+++ b/dune/perftool/matrixfree.hh
@@ -0,0 +1,30 @@
+#ifndef DUNE_PERFTOOL_MATRIXFREE_HH
+#define DUNE_PERFTOOL_MATRIXFREE_HH
+
+#include <iostream>
+
+#include <dune/pdelab/backend/istl/seqistlsolverbackend.hh>
+
+namespace Dune{
+ namespace PDELab{
+
+ template <typename GO, typename V>
+ void solveMatrixFree(GO& go, V& x ){
+ typedef Dune::PDELab::OnTheFlyOperator<V,V,GO> ISTLOnTheFlyOperator;
+ ISTLOnTheFlyOperator opb(go);
+ Dune::Richardson<V,V> richardson(1.0);
+ Dune::BiCGSTABSolver<V> solverb(opb,richardson,1E-10,5000,2);
+ Dune::InverseOperatorResult stat;
+ // evaluate residual w.r.t initial guess
+ using TrialGridFunctionSpace = typename GO::Traits::TrialGridFunctionSpace;
+ using W = Dune::PDELab::Backend::Vector<TrialGridFunctionSpace,typename V::ElementType>;
+ W r(go.testGridFunctionSpace(),0.0);
+ go.residual(x,r);
+ // solve the jacobian system
+ V z(go.trialGridFunctionSpace(),0.0);
+ solverb.apply(z,r,stat);
+ x -= z;
+ }
+ }
+}
+#endif
diff --git a/python/dune/perftool/loopy/functions.py b/python/dune/perftool/loopy/functions.py
index 58bd1d1f71b537c753cb6c56c2d2e2a7580352b4..2338e33f19bbb022f5ef30654b36b94bfb0913e1 100644
--- a/python/dune/perftool/loopy/functions.py
+++ b/python/dune/perftool/loopy/functions.py
@@ -23,16 +23,15 @@ def lfs_child_mangler(target, func, dtypes):
class CoefficientAccess(FunctionIdentifier):
- def __init__(self, restriction):
- self.restriction = restriction
+ def __init__(self, container):
+ self.container = container
def __getinitargs__(self):
- return (self.restriction,)
+ return (self.container,)
@property
def name(self):
- from dune.perftool.pdelab.argument import name_coefficientcontainer
- return name_coefficientcontainer(self.restriction)
+ return self.container
def coefficient_mangler(target, func, dtypes):
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 54938c547e06ccf743fdcac29594cd1eed7396da..630ed6e871963f0f2483282865586550ed955946 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -50,12 +50,13 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
# Initialize the local function spaces that we might need for this term
# We therefore gather a list of modified trial functions too.
from dune.perftool.ufl.modified_terminals import extract_modified_arguments
- test_ma = extract_modified_arguments(o)
- trial_ma = extract_modified_arguments(o, trialfunction=True, testfunction=False)
+ test_ma = extract_modified_arguments(o,)
+ trial_ma = extract_modified_arguments(o, coeffcount=0)
+ apply_ma = extract_modified_arguments(o, coeffcount=1)
# All test and trial functions on boundary integrals are technically restricted
import itertools
- for ma in itertools.chain(test_ma, trial_ma):
+ for ma in itertools.chain(test_ma, trial_ma, apply_ma):
if self.measure == 'exterior_facet':
ma.restriction = Restriction.NEGATIVE
@@ -209,21 +210,29 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
return Subscript(Variable(name_testfunction(leaf_element, restriction)), (Variable(iname),))
def coefficient(self, o):
- # If this is a trialfunction, we do something entirely different
- if o.count() == 0:
+ # Do something different for trial function and coefficients from jacobian apply
+ if o.count() == 0 or o.count() == 1:
# Correct the restriction on boundary integrals
restriction = self.restriction
if self.measure == 'exterior_facet':
restriction = Restriction.NEGATIVE
if self.grad:
- from dune.perftool.pdelab.argument import name_trialfunction_gradient
- return Variable(name_trialfunction_gradient(o.ufl_element(), restriction, self.component))
+ if o.count() == 0:
+ from dune.perftool.pdelab.argument import name_trialfunction_gradient
+ return Variable(name_trialfunction_gradient(o.ufl_element(), restriction, self.component))
+ else:
+ from dune.perftool.pdelab.argument import name_apply_function_gradient
+ return Variable(name_apply_function_gradient(o.ufl_element(), restriction, self.component))
else:
- from dune.perftool.pdelab.argument import name_trialfunction
- return Variable(name_trialfunction(o.ufl_element(), restriction, self.component))
-
- # so this is a parameter function
+ if o.count() == 0:
+ from dune.perftool.pdelab.argument import name_trialfunction
+ return Variable(name_trialfunction(o.ufl_element(), restriction, self.component))
+ else:
+ from dune.perftool.pdelab.argument import name_apply_function
+ return Variable(name_apply_function(o.ufl_element(), restriction, self.component))
+
+ # Check if this is a parameter function
else:
# We expect all coefficients to be of type Expression or VectorExpression!
from dune.perftool.ufl.execution import Expression, VectorExpression
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index ff563cfb2bd1a50af135571fc43a248009b26018..bbea79cfe4e129bc19fc3e96563dd2c34cc2ff43 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -13,6 +13,7 @@ def get_form_compiler_arguments():
parser.add_argument("--operator-file", type=str, help="The filename for the generated local operator header")
parser.add_argument('--version', action='version', version='%(prog)s 0.1')
parser.add_argument("--numerical-jacobian", action="store_true", help="use numerical jacobians (only makes sense, if uflpdelab for some reason fails to generate analytic jacobians)")
+ parser.add_argument("--matrix-free", action="store_true", help="use iterative solver with matrix free jacobian application")
parser.add_argument(
"--exact-solution-expression",
type=str,
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 6ff8444d27df79f4822629de12d15310422448e6..dc3a508b0da97591af582b75119cfaca0134e161 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -5,8 +5,8 @@ from dune.perftool.ufl.modified_terminals import ModifiedArgumentDescriptor
from dune.perftool.pdelab import (name_index,
restricted_name,
)
-from dune.perftool.pdelab.basis import (evaluate_trialfunction,
- evaluate_trialfunction_gradient,
+from dune.perftool.pdelab.basis import (evaluate_coefficient,
+ evaluate_coefficient_gradient,
lfs_iname,
name_basis,
name_basis_gradient,
@@ -35,7 +35,8 @@ def name_testfunction(element, restriction):
def name_trialfunction_gradient(element, restriction, component):
rawname = "gradu" + "_".join(str(c) for c in component)
name = restricted_name(rawname, restriction)
- evaluate_trialfunction_gradient(element, name, restriction, component)
+ container = name_coefficientcontainer(restriction)
+ evaluate_coefficient_gradient(element, name, container, restriction, component)
return name
@@ -43,7 +44,26 @@ def name_trialfunction_gradient(element, restriction, component):
def name_trialfunction(element, restriction, component):
rawname = "u" + "_".join(str(c) for c in component)
name = restricted_name(rawname, restriction)
- evaluate_trialfunction(element, name, restriction, component)
+ container = name_coefficientcontainer(restriction)
+ evaluate_coefficient(element, name, container, restriction, component)
+ return name
+
+
+@symbol
+def name_apply_function_gradient(element, restriction, component):
+ rawname = "gradz_func" + "_".join(str(c) for c in component)
+ name = restricted_name(rawname, restriction)
+ container = name_applycontainer(restriction)
+ evaluate_coefficient_gradient(element, name, container, restriction, component)
+ return name
+
+
+@symbol
+def name_apply_function(element, restriction, component):
+ rawname = "z_func" + "_".join(str(c) for c in component)
+ name = restricted_name(rawname, restriction)
+ container = name_applycontainer(restriction)
+ evaluate_coefficient(element, name, container, restriction, component)
return name
@@ -73,8 +93,14 @@ def name_coefficientcontainer(restriction):
return name
+@symbol
+def name_applycontainer(restriction):
+ name = restricted_name("z", restriction)
+ return name
+
+
@pymbolic_expr
-def pymbolic_coefficient(lfs, index, restriction):
+def pymbolic_coefficient(container, lfs, index):
# TODO introduce a proper type for local function spaces!
if isinstance(lfs, str):
valuearg(lfs, dtype=loopy.types.NumpyType("str"))
@@ -85,7 +111,7 @@ def pymbolic_coefficient(lfs, index, restriction):
lfs = Variable(lfs)
from dune.perftool.loopy.functions import CoefficientAccess
- return Call(CoefficientAccess(restriction), (lfs, Variable(index),))
+ return Call(CoefficientAccess(container), (lfs, Variable(index),))
@symbol
@@ -101,7 +127,8 @@ def name_argumentspace(ma):
if ma.argexpr.number() == 1:
return name_trialfunctionspace(ma.restriction)
if isinstance(ma.argexpr, Coefficient):
- assert ma.argexpr.count() == 0
+ # Index 0 is reserved for trialfunction, index 1 is reserved for jacobian apply function
+ assert ma.argexpr.count() < 2
return name_trialfunctionspace(ma.restriction)
# We should never encounter an argument other than 0 or 1
assert False
@@ -121,7 +148,6 @@ def name_jacobian(restriction1, restriction2):
# Restrictions may only differ if NONE
if (restriction1 == Restriction.NONE) or (restriction2 == Restriction.NONE):
assert restriction1 == restriction2
-
return restricted_name(restricted_name("jac", restriction1), restriction2)
@@ -143,7 +169,7 @@ def type_residual():
def name_accumulation_variable(restrictions=(Restriction.NONE,)):
from dune.perftool.generation import get_global_context_value
ft = get_global_context_value("form_type")
- if ft == 'residual':
+ if ft == 'residual' or ft == 'jacobian_apply':
return name_residual(*restrictions)
if ft == 'jacobian':
return name_jacobian(*restrictions)
@@ -153,7 +179,7 @@ def name_accumulation_variable(restrictions=(Restriction.NONE,)):
def type_accumulation_variable():
from dune.perftool.generation import get_global_context_value
ft = get_global_context_value("form_type")
- if ft == 'residual':
+ if ft == 'residual' or ft == 'jacobian_apply':
return type_residual()
if ft == 'jacobian':
return type_jacobian()
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index cacdb4cd07b1c8dccb01880b2b406d33d4e74b67..d6adc144be19f513a6f27830a9fa5483d9ce84e7 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -323,7 +323,7 @@ def shape_as_pymbolic(shape):
@cached
-def evaluate_trialfunction(element, name, restriction, component):
+def evaluate_coefficient(element, name, container, restriction, component):
from ufl.functionview import select_subelement
sub_element = select_subelement(element, component)
@@ -348,7 +348,7 @@ def evaluate_trialfunction(element, name, restriction, component):
lfs = lfs_child(lfs, idims, shape=shape_as_pymbolic(shape))
from dune.perftool.pdelab.argument import pymbolic_coefficient
- coeff = pymbolic_coefficient(lfs, index, restriction)
+ coeff = pymbolic_coefficient(container, lfs, index)
assignee = Subscript(Variable(name), tuple(Variable(i) for i in idims))
reduction_expr = Product((coeff, Subscript(Variable(basis), Variable(index))))
@@ -360,7 +360,7 @@ def evaluate_trialfunction(element, name, restriction, component):
@cached
-def evaluate_trialfunction_gradient(element, name, restriction, component):
+def evaluate_coefficient_gradient(element, name, container, restriction, component):
# First we determine the rank of the tensor we are talking about
from ufl.functionview import select_subelement
sub_element = select_subelement(element, component)
@@ -387,10 +387,11 @@ def evaluate_trialfunction_gradient(element, name, restriction, component):
lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]))
from dune.perftool.pdelab.argument import pymbolic_coefficient
- coeff = pymbolic_coefficient(lfs, index, restriction)
+ coeff = pymbolic_coefficient(container, lfs, index)
assignee = Subscript(Variable(name), tuple(Variable(i) for i in idims))
reduction_expr = Product((coeff, Subscript(Variable(basis), (Variable(index), Variable(idims[-1])))))
+
instruction(expression=Reduction("sum", index, reduction_expr, allow_simultaneous=True),
assignee=assignee,
forced_iname_deps=frozenset({quadrature_iname()}).union(frozenset(idims)),
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 9940aba70a83b17e606b3912bda72aa55ef4dd9a..ff87bd114c7ee54107f25ece4e5d0d83531e3c87 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -887,6 +887,7 @@ def define_stationarynonlinearproblemsolver(name):
ls = name_linearsolver()
return "{} {}({}, {}, {});".format(snptype, name, go, x, ls)
+
@symbol
def name_stationarylinearproblemsolver():
define_stationarylinearproblemsolver("slp")
@@ -903,11 +904,20 @@ def name_stationarynonlinearproblemsolver():
def dune_solve():
# Test if form is linear in ansatzfunction
if is_linear(_form):
- slp = name_stationarylinearproblemsolver()
- return "{}.apply();".format(slp)
+ if get_option("matrix_free"):
+ go = name_gridoperator()
+ x = name_vector()
+ include_file("dune/perftool/matrixfree.hh", filetag="driver")
+ return "solveMatrixFree({},{});".format(go, x)
+ else:
+ slp = name_stationarylinearproblemsolver()
+ return "{}.apply();".format(slp)
else:
- snp = name_stationarynonlinearproblemsolver()
- return "{}.apply();".format(snp)
+ if get_option("matrix_free"):
+ raise NotImplementedError("Nonlinear and matrix free is not yet implemented")
+ else:
+ snp = name_stationarynonlinearproblemsolver()
+ return "{}.apply();".format(snp)
@preamble
@@ -1058,7 +1068,7 @@ def vtkoutput():
if get_option("exact_solution_expression"):
from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(_form.coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
+ if isinstance(_form.coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
if get_option("compare_dofs"):
@@ -1066,6 +1076,5 @@ def vtkoutput():
if get_option("compare_l2errorsquared"):
compare_L2_squared()
-
return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vec, predicate),
"{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index bf4f93e0ede645c726d630841c1b76f6f3b20696..b539e14fdce6942cf288968e47cc2effc7928a65 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -36,7 +36,8 @@ def lop_template_gfs(ma):
if ma.argexpr.number() == 1:
return lop_template_ansatz_gfs()
if isinstance(ma.argexpr, Coefficient):
- assert ma.argexpr.count() == 0
+ # Index 0 is reserved for trialfunction, index 1 is reserved for jacobian apply function
+ assert ma.argexpr.count() < 2
return lop_template_ansatz_gfs()
assert False
@@ -150,6 +151,17 @@ def assembly_routine_signature():
from dune.perftool.pdelab.signatures import jacobian_skeleton_signature
return jacobian_skeleton_signature()
+ if form_type == 'jacobian_apply':
+ if integral_type == 'cell':
+ from dune.perftool.pdelab.signatures import jacobian_apply_volume_signature
+ return jacobian_apply_volume_signature()
+ if integral_type == 'exterior_facet':
+ from dune.perftool.pdelab.signatures import jacobian_apply_boundary_signature
+ return jacobian_apply_boundary_signature()
+ if integral_type == 'interior_facet':
+ from dune.perftool.pdelab.signatures import jacobian_apply_skeleton_signature
+ return jacobian_apply_skeleton_signature()
+
assert False
@@ -319,16 +331,30 @@ def generate_localoperator_kernels(formdata, data):
_, loptype = class_type_from_cache("operator")
which = ufl_measure_to_pdelab_measure(measure)
+
+ # Numerical jacobian methods
base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")
# Add the initializer list for that base class
ini = name_initree_member()
ini_constructor = name_initree_constructor_param()
+
+ # Numerical jacobian methods
initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
classtag="operator",
)
+ if get_option("matrix_free"):
+ # Numeical jacobian apply base class
+ base_class("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype), classtag="operator")
+
+ # Numerical jacobian apply initializer list
+ initializer_list("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype),
+ ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
+ classtag="operator",
+ )
+
operator_kernels[(measure, 'residual')] = kernel
# Generate the necessary jacobian methods
@@ -351,7 +377,30 @@ def generate_localoperator_kernels(formdata, data):
for it in alpha_measures - jacobian_measures:
operator_kernels[(it, 'jacobian')] = None
- # TODO: JacobianApply for matrix-free computations.
+ # Jacobian apply methods for matrix-free computations
+ if get_option("matrix_free"):
+ # The apply vector has reserved index 1 so we directly use Coefficient class from ufl
+ from ufl import Coefficient
+ apply_coefficient = Coefficient(form.coefficients()[0].ufl_element(), 1)
+
+ # Create application of jacobian on vector
+ from ufl import action
+ jac_apply_form = action(jacform, apply_coefficient)
+
+ # Create kernel for jacobian application
+ with global_context(form_type="jacobian_apply"):
+ for measure in set(i.integral_type() for i in jac_apply_form.integrals()):
+ with global_context(integral_type=measure):
+ kernel = generate_kernel(jac_apply_form.integrals_by_type(measure))
+ operator_kernels[(measure, 'jacobian_apply')] = kernel
+
+ # Generate dummy functions for those kernels, that vanished in the differentiation process
+ # We *could* solve this problem by using lambda_* terms but we do not really want that, so
+ # we use empty jacobian assembly methods instead
+ alpha_measures = set(i.integral_type() for i in form.integrals())
+ jacobian_apply_measures = set(i.integral_type() for i in jac_apply_form.integrals())
+ for it in alpha_measures - jacobian_apply_measures:
+ operator_kernels[(it, 'jacobian_apply')] = None
# Return the set of generated kernels
return operator_kernels
diff --git a/python/dune/perftool/pdelab/signatures.py b/python/dune/perftool/pdelab/signatures.py
index b0e02593ba8975387c243ea999edec99fe183310..e778b5952ac299f852c4bb071c8e933a00a74207 100644
--- a/python/dune/perftool/pdelab/signatures.py
+++ b/python/dune/perftool/pdelab/signatures.py
@@ -9,6 +9,7 @@ from dune.perftool.pdelab.argument import (name_accumulation_variable,
type_accumulation_variable,
name_coefficientcontainer,
type_coefficientcontainer,
+ name_applycontainer,
name_testfunctionspace,
type_testfunctionspace,
name_trialfunctionspace,
@@ -236,3 +237,111 @@ def jacobian_skeleton_signature():
av_nn,
)
]
+
+
+@symbol
+def jacobian_apply_volume_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu = name_trialfunctionspace(Restriction.NONE)
+ lfsvt = type_testfunctionspace()
+ lfsv = name_testfunctionspace(Restriction.NONE)
+ cct = type_coefficientcontainer()
+ ac = name_applycontainer(Restriction.NONE)
+ avt = type_accumulation_variable()
+ av = name_accumulation_variable((Restriction.NONE,))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void jacobian_apply_volume(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(
+ geot,
+ geo,
+ lfsut,
+ lfsu,
+ cct,
+ ac,
+ lfsvt,
+ lfsv,
+ avt,
+ av,)
+ ]
+
+
+@symbol
+def jacobian_apply_boundary_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu = name_trialfunctionspace(Restriction.NEGATIVE)
+ lfsvt = type_testfunctionspace()
+ lfsv = name_testfunctionspace(Restriction.NEGATIVE)
+ cct = type_coefficientcontainer()
+ ac = name_applycontainer(Restriction.NEGATIVE)
+ avt = type_accumulation_variable()
+ av = name_accumulation_variable((Restriction.NEGATIVE,))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void jacobian_apply_boundary(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(
+ geot,
+ geo,
+ lfsut,
+ lfsu,
+ cct,
+ ac,
+ lfsvt,
+ lfsv,
+ avt,
+ av,)
+ ]
+
+
+@symbol
+def jacobian_apply_skeleton_signature():
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+ lfsut = type_trialfunctionspace()
+ lfsu_s = name_trialfunctionspace(Restriction.NEGATIVE)
+ lfsu_n = name_trialfunctionspace(Restriction.POSITIVE)
+ lfsvt = type_testfunctionspace()
+ lfsv_s = name_testfunctionspace(Restriction.NEGATIVE)
+ lfsv_n = name_testfunctionspace(Restriction.POSITIVE)
+ cct = type_coefficientcontainer()
+ ac_s = name_applycontainer(Restriction.NEGATIVE)
+ ac_n = name_applycontainer(Restriction.POSITIVE)
+ avt = type_accumulation_variable()
+ av_s = name_accumulation_variable((Restriction.NEGATIVE,))
+ av_n = name_accumulation_variable((Restriction.POSITIVE,))
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void jacobian_apply_skeleton(const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}, {}& {}) const'.format(
+ geot,
+ geo,
+ lfsut,
+ lfsu_s,
+ cct,
+ ac_s,
+ lfsvt,
+ lfsv_s,
+ lfsut,
+ lfsu_n,
+ cct,
+ ac_n,
+ lfsvt,
+ lfsv_n,
+ avt,
+ av_s,
+ avt,
+ av_n,)
+ ]
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 8018e8593cbc3cffb304a6aea1fa62f8423846ff..5fe6975b233a73a2629186489742f5e897b6d7f2 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -31,8 +31,10 @@ class Coefficient(ufl.Coefficient):
def __init__(self, element, count=None):
if count and count is 0:
raise ValueError("The coefficient of index 0 is reserved for the trial function in uflpdelab")
- if not count and ufl.Coefficient._globalcount is 0:
- count = 1
+ if count and count is 1:
+ raise ValueError("The coefficient of index 1 is reserved for the jacobian apply vector in uflpdelab")
+ if not count and ufl.Coefficient._globalcount < 2:
+ count = 2
ufl.Coefficient.__init__(self, element, count)
diff --git a/python/dune/perftool/ufl/modified_terminals.py b/python/dune/perftool/ufl/modified_terminals.py
index 27543551b6b218c299e4b82502a4e00c2cf1af47..82bb8eee94b2c0ec9e3a46835faae0e0dc610e70 100644
--- a/python/dune/perftool/ufl/modified_terminals.py
+++ b/python/dune/perftool/ufl/modified_terminals.py
@@ -103,10 +103,9 @@ class ModifiedArgumentDescriptor(MultiFunction):
class _ModifiedArgumentExtractor(MultiFunction):
""" A multifunction that extracts and returns the set of modified arguments """
- def __call__(self, o, argnumber=None, testfunction=True, trialfunction=False):
+ def __call__(self, o, argnumber=None, coeffcount=None):
self.argnumber = argnumber
- self.trialfunction = trialfunction
- self.testfunction = testfunction
+ self.coeffcount = coeffcount
self.modified_arguments = set()
ret = self.call(o)
if ret:
@@ -132,14 +131,12 @@ class _ModifiedArgumentExtractor(MultiFunction):
function_view = pass_on
def argument(self, o):
- if self.testfunction:
- if self.argnumber is None or o.number() == self.argnumber:
- return o
+ if self.argnumber is None or o.number() == self.argnumber:
+ return o
def coefficient(self, o):
- if self.trialfunction:
- if o.count() == 0:
- return o
+ if o.count() == self.coeffcount:
+ return o
call = MultiFunction.__call__
diff --git a/test/nonlinear/CMakeLists.txt b/test/nonlinear/CMakeLists.txt
index a78cccd05c560d4ed69edf860d2cbbc25f3928b4..cd481d317a73263b0df08383db3e1ab6f6f2dc1d 100644
--- a/test/nonlinear/CMakeLists.txt
+++ b/test/nonlinear/CMakeLists.txt
@@ -3,7 +3,6 @@ add_generated_executable(UFLFILE nonlinear.ufl
FORM_COMPILER_ARGS
)
-
dune_add_system_test(TARGET nonlinear_symdiff
INIFILE nonlinear_symdiff.mini
SCRIPT dune_vtkcompare.py
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index 88f33d3a7d5a8727eabec2abdf3d286bff882680..1a7aa445b5fc8cc706650cea3830dbce78b43559 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -88,6 +88,50 @@ dune_add_system_test(TARGET poisson_dg_neumann_symdiff
INIFILE poisson_dg_neumann_symdiff.mini
)
+# 9. Poisson Test Case: matrix free, numeric differentiation
+add_generated_executable(UFLFILE poisson.ufl
+ TARGET poisson_numdiff_matrix_free
+ FORM_COMPILER_ARGS --numerical-jacobian --matrix-free
+ )
+
+dune_add_system_test(TARGET poisson_numdiff_matrix_free
+ INIFILE poisson_numdiff_matrix_free.mini
+ SCRIPT dune_vtkcompare.py
+ )
+
+# 10. Poisson Test Case: matrix free, symbolic differentiation
+add_generated_executable(UFLFILE poisson.ufl
+ TARGET poisson_symdiff_matrix_free
+ FORM_COMPILER_ARGS --matrix-free
+ )
+
+dune_add_system_test(TARGET poisson_symdiff_matrix_free
+ INIFILE poisson_symdiff_matrix_free.mini
+ SCRIPT dune_vtkcompare.py
+ )
+
+# # 11. Poisson Test Case: dg, matrix free, numeric differentiation
+# add_generated_executable(UFLFILE poisson_dg.ufl
+# TARGET poisson_dg_numdiff_matrix_free
+# FORM_COMPILER_ARGS --numerical-jacobian --matrix-free
+# )
+
+# dune_add_system_test(TARGET poisson_dg_numdiff_matrix_free
+# INIFILE poisson_dg_numdiff_matrix_free.mini
+# SCRIPT dune_vtkcompare.py
+# )
+
+# # 12. Poisson Test Case: dg, matrix free, symbolic differentiation
+# add_generated_executable(UFLFILE poisson_dg.ufl
+# TARGET poisson_dg_symdiff_matrix_free
+# FORM_COMPILER_ARGS --matrix-free
+# )
+
+# dune_add_system_test(TARGET poisson_dg_symdiff_matrix_free
+# INIFILE poisson_dg_symdiff_matrix_free.mini
+# SCRIPT dune_vtkcompare.py
+# )
+
# Add the reference code with the PDELab localoperator that produced
# the reference vtk file
add_executable(poisson_dg_ref reference_main.cc)
diff --git a/test/poisson/dimension-grid-variations/CMakeLists.txt b/test/poisson/dimension-grid-variations/CMakeLists.txt
index 13b52b6dc48b7b6b17bdd5a5a4b723f91de2954b..9f2229d084feb7e05baa3111ce4b6d4879269a5f 100644
--- a/test/poisson/dimension-grid-variations/CMakeLists.txt
+++ b/test/poisson/dimension-grid-variations/CMakeLists.txt
@@ -13,7 +13,7 @@ dune_add_system_test(TARGET poisson_1d_cg_interval
add_generated_executable(UFLFILE poisson_1d_dg_interval.ufl
TARGET poisson_1d_dg_interval
- FORM_COMPILER_ARGS --exact-solution-expression g --compare-dofs 1e-2 --compare-l2errorsquared 1e-7
+ FORM_COMPILER_ARGS --exact-solution-expression g --compare-dofs 1e-2 --compare-l2errorsquared 1e-5
)
dune_add_system_test(TARGET poisson_1d_dg_interval
diff --git a/test/poisson/poisson_dg_numdiff_matrix_free.mini b/test/poisson/poisson_dg_numdiff_matrix_free.mini
new file mode 100644
index 0000000000000000000000000000000000000000..b45d4a8654e11107f8bbf111db0f996461085e2a
--- /dev/null
+++ b/test/poisson/poisson_dg_numdiff_matrix_free.mini
@@ -0,0 +1,11 @@
+__name = poisson_dg_numdiff_matrix_free
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_dg_ref
+extension = vtu
diff --git a/test/poisson/poisson_dg_symdiff_matrix_free.mini b/test/poisson/poisson_dg_symdiff_matrix_free.mini
new file mode 100644
index 0000000000000000000000000000000000000000..c3fcaf7fcf1576850a8d1a99f4a1d8b07d0ddcd8
--- /dev/null
+++ b/test/poisson/poisson_dg_symdiff_matrix_free.mini
@@ -0,0 +1,11 @@
+__name = poisson_dg_symdiff_matrix_free
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_dg_ref
+extension = vtu
diff --git a/test/poisson/poisson_numdiff_matrix_free.mini b/test/poisson/poisson_numdiff_matrix_free.mini
new file mode 100644
index 0000000000000000000000000000000000000000..4662d83881c20c13ffe49055338e1fc775020ab3
--- /dev/null
+++ b/test/poisson/poisson_numdiff_matrix_free.mini
@@ -0,0 +1,11 @@
+__name = poisson_numdiff_matrix_free
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_ref
+extension = vtu
diff --git a/test/poisson/poisson_symdiff_matrix_free.mini b/test/poisson/poisson_symdiff_matrix_free.mini
new file mode 100644
index 0000000000000000000000000000000000000000..8d288ef3358fdda09706b80874e444d0ff5c8b2c
--- /dev/null
+++ b/test/poisson/poisson_symdiff_matrix_free.mini
@@ -0,0 +1,11 @@
+__name = poisson_symdiff_matrix_free
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_ref
+extension = vtu