diff --git a/dune/perftool/matrixfree.hh b/dune/perftool/matrixfree.hh
index 3ffa62f1a9fbba76168e55ce0c8dbfdb43301dd7..4ab94377e2aa7d8790360042c89a9423bf4a6352 100644
--- a/dune/perftool/matrixfree.hh
+++ b/dune/perftool/matrixfree.hh
@@ -3,12 +3,14 @@
#include <iostream>
-#include <dune/common/timer.hh>
-#include <dune/common/parametertree.hh>
+// #include <dune/common/timer.hh>
+// #include <dune/common/parametertree.hh>
-#include <dune/pdelab/backend/interface.hh>
-#include <dune/pdelab/constraints/common/constraints.hh>
-#include <dune/pdelab/backend/solver.hh>
+// #include <dune/pdelab/backend/interface.hh>
+// #include <dune/pdelab/constraints/common/constraints.hh>
+// #include <dune/pdelab/backend/solver.hh>
+
+#include <dune/pdelab/backend/istl/seqistlsolverbackend.hh>
namespace Dune{
namespace PDELab{
@@ -18,7 +20,7 @@ namespace Dune{
typedef Dune::PDELab::OnTheFlyOperator<V,V,GO> ISTLOnTheFlyOperator;
ISTLOnTheFlyOperator opb(go);
Dune::Richardson<V,V> richardson(1.0);
- Dune::CGSolver<V> solverb(opb,richardson,1E-10,5000,2);
+ 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;
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 54938c547e06ccf743fdcac29594cd1eed7396da..3cad73b4eb16076948c536db2951ea486916c4b2 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, testfunction=False, trialfunction=True)
+ apply_ma = extract_modified_arguments(o, testfunction=False, applyfunction=True)
# 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/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 6ff8444d27df79f4822629de12d15310422448e6..818075addc1ac1e04f95191007361e6c2517991f 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -47,6 +47,22 @@ def name_trialfunction(element, 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)
+ evaluate_trialfunction_gradient(element, name, 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)
+ evaluate_trialfunction(element, name, restriction, component)
+ return name
+
+
@symbol
def name_testfunctionspace(restriction):
return restricted_name("lfsv", restriction)
@@ -73,6 +89,12 @@ 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):
# TODO introduce a proper type for local function spaces!
@@ -101,7 +123,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 +144,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 +165,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 +175,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/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 3b89f76e513eb8b1a26c2e82eafcd7fccb1ebef9..9b2ad5fa68a1f6f2b5b3733d342c95012e0737cd 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
@@ -365,7 +377,27 @@ 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"):
+ from dune.perftool.ufl.execution import ApplyCoefficient
+ apply_coefficient = ApplyCoefficient(form.coefficients()[0].ufl_element(),1)
+ from ufl import action
+ jac_apply_form = action(expand_derivatives(derivative(form, form.coefficients()[0])),apply_coefficient)
+
+ 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..45745524d9063c8783dc704ff62e357981d5eca5 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,114 @@ 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()
+ cc = name_coefficientcontainer(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,
+ cc,
+ 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()
+ cc = name_coefficientcontainer(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,
+ cc,
+ 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()
+ cc_s = name_coefficientcontainer(Restriction.NEGATIVE)
+ cc_n = name_coefficientcontainer(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,
+ cc_s,
+ lfsvt,
+ lfsv_s,
+ lfsut,
+ lfsu_n,
+ cct,
+ cc_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..d420edf2e66904c49f51282366a52e0682b99a45 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -26,13 +26,22 @@ class TrialFunction(ufl.Coefficient):
ufl.Coefficient.__init__(self, element, count=0)
+class ApplyCoefficient(ufl.Coefficient):
+ """ A coefficient that always takes the reserved index 1 """
+ def __init__(self, element, count=None):
+ if count and count is not 1:
+ raise ValueError("The jacobian apply coefficient must be of index 1 in uflpdelab")
+ ufl.Coefficient.__init__(self, element, count=1)
+
class Coefficient(ufl.Coefficient):
""" A coefficient that honors the reserved index 0. """
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..88eed9589636b1fcbd40e163b7889bcc98b1dc05 100644
--- a/python/dune/perftool/ufl/modified_terminals.py
+++ b/python/dune/perftool/ufl/modified_terminals.py
@@ -103,10 +103,11 @@ 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, testfunction=True, trialfunction=False, applyfunction=False):
self.argnumber = argnumber
self.trialfunction = trialfunction
self.testfunction = testfunction
+ self.applyfunction = applyfunction
self.modified_arguments = set()
ret = self.call(o)
if ret:
@@ -140,6 +141,9 @@ class _ModifiedArgumentExtractor(MultiFunction):
if self.trialfunction:
if o.count() == 0:
return o
+ if self.applyfunction:
+ if o.count() == 1:
+ return o
call = MultiFunction.__call__
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index 3f95382ff6ec91b8f3fd8f3c1f61d231a7ea02b4..c096c63a9da4dfa8a3bc70e50a8d599b55a5ee38 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -88,7 +88,7 @@ dune_add_system_test(TARGET poisson_dg_neumann_symdiff
INIFILE poisson_dg_neumann_symdiff.mini
)
-# 9. Poisson Test Case: source f, bc g + numerical differentiation
+# 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
@@ -99,6 +99,38 @@ dune_add_system_test(TARGET poisson_numdiff_matrix_free
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
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_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