From 0e5b2a7ee372e40df070718fd40289310ceb6156 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Mon, 25 Apr 2016 17:12:46 +0200
Subject: [PATCH] Towards neumann boundaries: subdomain_data implemented
---
python/dune/perftool/compile.py | 11 ++--
python/dune/perftool/loopy/transformer.py | 48 ++++++++++++++---
python/dune/perftool/pdelab/localoperator.py | 54 +++++++++++---------
python/dune/perftool/pdelab/parameter.py | 8 +--
python/dune/perftool/ufl/execution.py | 4 +-
test/poisson/CMakeLists.txt | 12 +++++
test/poisson/poisson_neumann.ufl | 13 +++++
7 files changed, 111 insertions(+), 39 deletions(-)
create mode 100644 test/poisson/poisson_neumann.ufl
diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index 450e0b1e..83265ae7 100644
--- a/python/dune/perftool/compile.py
+++ b/python/dune/perftool/compile.py
@@ -33,7 +33,8 @@ def read_ufl(uflfile):
# Extract the form from the given data
data = interpret_ufl_namespace(namespace)
form = data.forms[0]
- form = compute_form_data(form).preprocessed_form
+ formdata = compute_form_data(form)
+ form = formdata.preprocessed_form
# We do not expect more than one form
assert len(data.forms) == 1
@@ -48,7 +49,7 @@ def read_ufl(uflfile):
form = transform_form(form, reindexing)
# form = transform_form(form, split_arguments)
- return form, data.object_names
+ return formdata, data.object_names
def generate_driver(form, filename):
@@ -78,16 +79,16 @@ def generate_driver(form, filename):
# This function is the entrypoint of the ufl2pdelab executable
def compile_form():
from dune.perftool.options import get_option
- form, namedata = read_ufl(get_option("uflfile"))
+ formdata, namedata = read_ufl(get_option("uflfile"))
from dune.perftool.generation import cache_context
if get_option("driver_file"):
with cache_context('driver', delete=True):
- generate_driver(form, get_option("driver_file"))
+ generate_driver(formdata.preprocessed_form, get_option("driver_file"))
if get_option("operator_file"):
from dune.perftool.pdelab.localoperator import generate_localoperator_kernels
- kernels = generate_localoperator_kernels(form, namedata)
+ kernels = generate_localoperator_kernels(formdata, namedata)
# TODO insert sophisticated analysis/feedback loops here
if get_option("interactive"):
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 1c0643dc..909ed752 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -85,7 +85,7 @@ def get_count():
return c
-def transform_accumulation_term(term):
+def transform_accumulation_term(term, measure, subdomain_id):
from dune.perftool.ufl.transformations.replace import ReplaceExpression
from pymbolic.primitives import Variable
@@ -180,11 +180,47 @@ def transform_accumulation_term(term):
from dune.perftool.pdelab.quadrature import name_factor
factor = name_factor()
- instruction(code="{}.accumulate({}, {}*{});".format(accumvar,
- ", ".join(accumargs),
- expr_tv_name,
- factor,
- ),
+
+ # Generate the code snippet for this accumulation instruction
+ code = "{}.accumulate({}, {}*{});".format(accumvar,
+ ", ".join(accumargs),
+ expr_tv_name,
+ factor,
+ )
+
+ # Maybe wrap this instruction into a condiditional. This mostly happens with mixed boundary conditions
+ if subdomain_id not in ['everywhere', 'otherwise']:
+ # We need to reconstruct the subdomain_data parameter of the measure
+ # I am *totally* confused as to why this information is not at hand anyway,
+ # but conversation with Martin pointed me to dolfin.fem.assembly where this
+ # is done in preprocessing with the limitation of only one possible type of
+ # modified measure per integral type.
+
+ # Get the original form and inspect the present measures
+ from dune.perftool.generation import get_generic_context_value
+ original_form = get_generic_context_value("formdata").original_form
+ sd = original_form.subdomain_data()
+ assert len(sd) == 1
+ subdomains, = list(sd.values())
+ domain, = list(sd.keys())
+ for k in list(subdomains.keys()):
+ if subdomains[k] is None:
+ del subdomains[k]
+
+ # Finally extract the original subdomain_data (which needs to be present!)
+ subdomain_data = subdomains[measure]
+
+ # Determine the name of the parameter function
+ name = get_generic_context_value("namedata")[id(subdomain_data)]
+
+ # Trigger the generation of code for this thing in the parameter class
+ from dune.perftool.pdelab.parameter import parameter_function
+ parameter_function(name, subdomain_data, t='int')
+
+ code = "if ({} == {})\n {}".format(name, subdomain_id, code)
+
+ # Finally, issue the instruction
+ instruction(code=code,
assignees=frozenset({accumvar}),
read_variables=frozenset({accumvar, factor, expr_tv_name}),
forced_iname_deps=acc_inames,
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 5289f829..8dc39a4a 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -135,8 +135,11 @@ def measure_specific_details(measure):
return ret
-def generate_kernel(integrand=None, measure=None):
- assert integrand and measure
+def generate_kernel(integral):
+ integrand = integral.integrand()
+ measure = integral.integral_type()
+ subdomain_id = integral.subdomain_id()
+ subdomain_data = integral.subdomain_data()
# Get the measure specifics
specifics = measure_specific_details(measure)
@@ -148,7 +151,7 @@ def generate_kernel(integrand=None, measure=None):
# Iterate over the terms and generate a kernel
for term in accterms:
from dune.perftool.loopy.transformer import transform_accumulation_term
- transform_accumulation_term(term)
+ transform_accumulation_term(term, measure, subdomain_id)
# Extract the information, which is needed to create a loopy kernel.
# First extracting it, might be useful to alter it before kernel generation.
@@ -204,7 +207,10 @@ def cgen_class_from_cache(tag, members=[]):
return Class(basename, base_classes=base_classes, members=members + pm, constructors=[constructor], tparam_decls=tparams)
-def generate_localoperator_kernels(form, namedata):
+def generate_localoperator_kernels(formdata, namedata):
+ # Extract the relevant attributes of the form data
+ form = formdata.preprocessed_form
+
# For the moment, I do assume that there is but one integral of each type. This might differ
# if you use different quadrature orders for different terms.
assert len(form.integrals()) == len(set(i.integral_type() for i in form.integrals()))
@@ -240,28 +246,30 @@ def generate_localoperator_kernels(form, namedata):
import functools
from dune.perftool.generation import generic_context
namedata_context = functools.partial(generic_context, "namedata")
+ formdata_context = functools.partial(generic_context, "formdata")
- # Generate the necessary residual methods
- for integral in form.integrals():
- with namedata_context(namedata):
- kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
- operator_kernels[(integral.integral_type(), 'residual')] = kernel
-
- # Generate the necessary jacobian methods
- from dune.perftool.options import get_option
- if get_option("numerical_jacobian"):
- include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
- else:
- from ufl import derivative
- from ufl.algorithms import expand_derivatives
- jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
-
- for integral in jacform.integrals():
+ with formdata_context(formdata):
+ # Generate the necessary residual methods
+ for integral in form.integrals():
with namedata_context(namedata):
- kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
- operator_kernels[(integral.integral_type(), 'jacobian')] = kernel
+ kernel = generate_kernel(integral)
+ operator_kernels[(integral.integral_type(), 'residual')] = kernel
- # TODO: JacobianApply for matrix-free computations.
+ # Generate the necessary jacobian methods
+ from dune.perftool.options import get_option
+ if get_option("numerical_jacobian"):
+ include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
+ else:
+ from ufl import derivative
+ from ufl.algorithms import expand_derivatives
+ jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
+
+ for integral in jacform.integrals():
+ with namedata_context(namedata):
+ kernel = generate_kernel(integrand=integral.integrand(), measure=integral.integral_type())
+ operator_kernels[(integral.integral_type(), 'jacobian')] = kernel
+
+ # TODO: JacobianApply for matrix-free computations.
# Return the set of generated kernels
return operator_kernels
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index 8abb856a..9b4885d3 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -39,9 +39,9 @@ def name_paramclass():
@class_member("parameterclass", access=AccessModifier.PUBLIC)
-def define_parameter_function_class_member(name, expr):
+def define_parameter_function_class_member(name, expr, t):
result = ["template<typename E, typename X>",
- "double {}(const E& e, const X& local) const".format(name),
+ "{} {}(const E& e, const X& local) const".format(t, name),
"{",
]
@@ -70,7 +70,7 @@ def evaluate_parameter_function(name):
)
-def parameter_function(name, expr):
+def parameter_function(name, expr, t='double'):
temporary_variable(name, shape=())
- define_parameter_function_class_member(name, expr)
+ define_parameter_function_class_member(name, expr, t)
evaluate_parameter_function(name)
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 09b66bf3..99114f50 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -44,6 +44,8 @@ def TrialFunctions(element):
return split(TrialFunction(element))
+_dg0 = FiniteElement("DG", "triangle", 0)
+
class Expression(Coefficient):
def __init__(self, expr, is_global=True):
assert isinstance(expr, str)
@@ -51,7 +53,7 @@ class Expression(Coefficient):
self.is_global = is_global
# Initialize a coefficient with a dummy finite element map.
- Coefficient.__init__(self, FiniteElement("DG", "triangle", 0))
+ Coefficient.__init__(self, _dg0)
class FiniteElement(ufl.FiniteElement):
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index 71fe4fcf..2ccd6cce 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -1,3 +1,4 @@
+# 1. Poisson Test Case: source f, bc g + numerical differentiation
add_generated_executable(UFLFILE poisson.ufl
TARGET poisson_numdiff
FORM_COMPILER_ARGS --numerical-jacobian
@@ -12,6 +13,7 @@ dune_add_system_test(TARGET poisson_numdiff
dune_symlink_to_source_files(FILES poisson_ref.vtu)
+# 2. Poisson Test Case: source f, bc g + symbolic differentiation
add_generated_executable(UFLFILE poisson.ufl
TARGET poisson_symdiff
)
@@ -20,3 +22,13 @@ dune_add_system_test(TARGET poisson_symdiff
INIFILE poisson.mini
SCRIPT dune_vtkcompare.py
)
+
+# 3. Poisson Test Case: source f, mixed neumann/dirichlet boundary + numerical differentiation
+add_generated_executable(UFLFILE poisson_neumann.ufl
+ TARGET poisson_neumann_numdiff
+ FORM_COMPILER_ARGS --numerical-jacobian
+ )
+
+dune_add_system_test(TARGET poisson_neumann_numdiff
+ INIFILE poisson.mini
+ )
diff --git a/test/poisson/poisson_neumann.ufl b/test/poisson/poisson_neumann.ufl
new file mode 100644
index 00000000..385a1bb8
--- /dev/null
+++ b/test/poisson/poisson_neumann.ufl
@@ -0,0 +1,13 @@
+f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
+g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
+j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return s*4*x[1]*c.two_norm2()*std::exp(-1.*c.two_norm2());")
+bctype = Expression("if ((x[1]<1e-8) || (x[1]>1-1e-8)) return 0; else return 1;")
+
+V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g, dirichlet_constraints=bctype)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# Define the boundary measure that knows where we are...
+ds = ds(subdomain_data=bctype)
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx - j*v*ds(0)]
--
GitLab