From 25d6648fa8b08704f22d0b606d21c0fe2a7da5f8 Mon Sep 17 00:00:00 2001
From: Marcel Koch <marcel.koch@web.de>
Date: Fri, 7 Jul 2017 09:47:45 +0200
Subject: [PATCH] adds computation of quadrature point in cell
---
python/dune/perftool/__init__.py | 1 +
.../dune/perftool/blockstructured/__init__.py | 2 +
.../perftool/blockstructured/accumulation.py | 59 +++++++++++++++++++
.../perftool/blockstructured/quadrature.py | 40 +++++++++++++
python/dune/perftool/options.py | 2 +
python/dune/perftool/pdelab/geometry.py | 5 +-
python/dune/perftool/pdelab/localoperator.py | 5 +-
python/dune/perftool/pdelab/quadrature.py | 6 +-
test/CMakeLists.txt | 14 +++--
test/blockstructured/CMakeLists.txt | 1 +
test/blockstructured/poisson/CMakeLists.txt | 4 ++
test/blockstructured/poisson/poisson.mini | 19 ++++++
test/blockstructured/poisson/poisson.ufl | 13 ++++
13 files changed, 162 insertions(+), 9 deletions(-)
create mode 100644 python/dune/perftool/blockstructured/__init__.py
create mode 100644 python/dune/perftool/blockstructured/accumulation.py
create mode 100644 python/dune/perftool/blockstructured/quadrature.py
create mode 100644 test/blockstructured/CMakeLists.txt
create mode 100644 test/blockstructured/poisson/CMakeLists.txt
create mode 100644 test/blockstructured/poisson/poisson.mini
create mode 100644 test/blockstructured/poisson/poisson.ufl
diff --git a/python/dune/perftool/__init__.py b/python/dune/perftool/__init__.py
index f656fe92..3b55e444 100644
--- a/python/dune/perftool/__init__.py
+++ b/python/dune/perftool/__init__.py
@@ -5,3 +5,4 @@ import dune.perftool.loopy.symbolic
# to the selection mechanisms
import dune.perftool.pdelab
import dune.perftool.sumfact
+import dune.perftool.blockstructured
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/__init__.py b/python/dune/perftool/blockstructured/__init__.py
new file mode 100644
index 00000000..7208ad95
--- /dev/null
+++ b/python/dune/perftool/blockstructured/__init__.py
@@ -0,0 +1,2 @@
+import dune.perftool.blockstructured.accumulation
+import dune.perftool.blockstructured.quadrature
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/accumulation.py b/python/dune/perftool/blockstructured/accumulation.py
new file mode 100644
index 00000000..255eaaa0
--- /dev/null
+++ b/python/dune/perftool/blockstructured/accumulation.py
@@ -0,0 +1,59 @@
+from dune.perftool.generation import (backend,iname,domain)
+from dune.perftool.generation.loopy import transform
+from dune.perftool.pdelab.localoperator import (determine_accumulation_space,
+ boundary_predicates)
+
+@iname
+def sub_element_inames():
+ name = "subel"
+ domain(name, 4)
+ return (name,)
+
+
+@backend(interface="accum_insn", name='blockstructured')
+def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
+ # When we do not do sumfactorization we do not split the test function
+ assert(accterm.argument.expr is None)
+
+ # Do the tree traversal to get a pymbolic expression representing this expression
+ pymbolic_expr = visitor(accterm.term)
+
+ # It may happen that an entire accumulation term vanishes. We do nothing in that case
+ if pymbolic_expr == 0:
+ return
+
+ # Collect the lfs and lfs indices for the accumulate call
+ test_lfs = determine_accumulation_space(accterm.term, 0, measure)
+
+ # In the jacobian case, also determine the space for the ansatz space
+ ansatz_lfs = determine_accumulation_space(accterm.term, 1, measure)
+
+ # Collect the lfs and lfs indices for the accumulate call
+ from dune.perftool.pdelab.argument import name_accumulation_variable
+ accumvar = name_accumulation_variable((test_lfs.get_restriction() + ansatz_lfs.get_restriction()))
+
+ predicates = boundary_predicates(accterm.term, measure, subdomain_id, visitor)
+
+ rank = 1 if ansatz_lfs.lfs is None else 2
+
+ from dune.perftool.pdelab.argument import PDELabAccumulationFunction
+ from pymbolic.primitives import Call
+ expr = Call(PDELabAccumulationFunction(accumvar, rank),
+ (test_lfs.get_args() + ansatz_lfs.get_args() + (pymbolic_expr,))
+ )
+
+ from dune.perftool.generation import instruction
+ from dune.perftool.options import option_switch
+ quad_inames = visitor.interface.quadrature_inames()
+
+ instruction(assignees=(),
+ expression=expr,
+ forced_iname_deps=frozenset(visitor.inames).union(frozenset(quad_inames)),
+ forced_iname_deps_is_final=True,
+ predicates=predicates
+ )
+
+ subelem_inames = sub_element_inames()
+
+ from dune.perftool.sumfact.quadrature import nest_quadrature_loops
+ transform(nest_quadrature_loops, subelem_inames)
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/quadrature.py b/python/dune/perftool/blockstructured/quadrature.py
new file mode 100644
index 00000000..cb0d16ba
--- /dev/null
+++ b/python/dune/perftool/blockstructured/quadrature.py
@@ -0,0 +1,40 @@
+from dune.perftool.generation import (backend,
+ temporary_variable,
+ domain)
+from dune.perftool.pdelab.quadrature import (name_quadrature_points,
+ quadrature_iname
+ )
+from dune.perftool.pdelab.geometry import world_dimension
+from dune.perftool.blockstructured.accumulation import sub_element_inames
+from dune.perftool.generation.loopy import instruction
+from dune.perftool.options import get_option
+import pymbolic.primitives as prim
+
+
+@backend(interface="quad_pos", name='blockstructured')
+def pymbolic_quadrature_position():
+ quad_points = name_quadrature_points()
+ quad_iname = quadrature_iname()
+ from pymbolic.primitives import Subscript, Variable, Sum, Quotient
+
+ qp_in_micro = Subscript(Variable(quad_points), (Variable(quad_iname,)))
+
+ dim = world_dimension()
+ shape_impl = ('fv',)
+ qp_in_macro = 'macro_qp'
+ temporary_variable(qp_in_macro, shape_impl=shape_impl, shape=(dim,))
+ assignment_iname = 'assignment_iname'
+ domain(assignment_iname, dim)
+
+ expr = Subscript(qp_in_micro.aggregate, (qp_in_micro.index,) + (Variable(assignment_iname),))
+ expr = Sum((expr, Variable(assignment_iname)))
+ expr = Quotient(expr, get_option('number_of_blocks'))
+
+ instruction(assignee=Subscript(Variable(qp_in_macro), (Variable(assignment_iname),)),
+ expression=expr,
+ forced_iname_deps=frozenset((assignment_iname,)),
+ forced_iname_deps_is_final=True,
+ tags=frozenset({'quad'})
+ )
+
+ return Variable(qp_in_macro)
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 06670cd8..d3d361c8 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -66,6 +66,8 @@ class PerftoolOptionsArray(ImmutableRecord):
max_vector_width = PerftoolOption(default=256, helpstr=None)
unroll_dimension_loops = PerftoolOption(default=False, helpstr="whether loops over the gemetric dimension should be unrolled.")
precompute_quadrature_info = PerftoolOption(default=True, helpstr="whether loops over the gemetric dimension should be unrolled.")
+ blockstructured = PerftoolOption(default=False, helpstr="Use block structure")
+ number_of_blocks = PerftoolOption(default=1, helpstr="Number of sub blocks in one direction")
# Until more sophisticated logic is needed, we keep the actual option data in this module
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 65f40a10..c365e89e 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -368,7 +368,10 @@ def define_constant_jacobian_determinant(name):
def define_jacobian_determinant(name):
temporary_variable(name, shape=())
geo = name_geometry()
- pos = get_backend("quad_pos")()
+ if get_option('blockstructured'):
+ pos = get_backend("quad_pos", selector=lambda: 'blockstructured')()
+ else:
+ pos = get_backend("quad_pos")()
code = "{} = {}.integrationElement({});".format(name,
geo,
str(pos),
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 1e213f00..53580599 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -409,7 +409,10 @@ def visit_integrals(integrals):
from dune.perftool.ufl.visitor import UFL2LoopyVisitor
visitor = UFL2LoopyVisitor(interface, measure, indexmap)
- get_backend(interface="accum_insn")(visitor, accterm, measure, subdomain_id)
+ if get_option('blockstructured'):
+ get_backend(interface="accum_insn", selector=lambda: 'blockstructured')(visitor, accterm, measure, subdomain_id)
+ else:
+ get_backend(interface="accum_insn")(visitor, accterm, measure, subdomain_id)
def generate_kernel(integrals):
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 80fbf27a..1cd0d366 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -125,7 +125,11 @@ def pymbolic_quadrature_position():
@backend(interface="qp_in_cell")
def pymbolic_quadrature_position_in_cell(restriction):
from dune.perftool.pdelab.geometry import to_cell_coordinates
- return to_cell_coordinates(pymbolic_quadrature_position(), restriction)
+ if get_option("blockstructured"):
+ quad_pos = get_backend("quad_pos", selector=lambda: "blockstructured")()
+ else:
+ quad_pos = get_backend("quad_pos")()
+ return to_cell_coordinates(quad_pos, restriction)
@preamble
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f0812d3e..8a4cba8e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,7 +1,9 @@
-add_subdirectory(heatequation)
-add_subdirectory(laplace)
-add_subdirectory(nonlinear)
-add_subdirectory(poisson)
-add_subdirectory(stokes)
+#add_subdirectory(heatequation)
+#add_subdirectory(laplace)
+#add_subdirectory(nonlinear)
+#add_subdirectory(poisson)
+#add_subdirectory(stokes)
-add_subdirectory(sumfact)
+#add_subdirectory(sumfact)
+
+add_subdirectory(blockstructured)
\ No newline at end of file
diff --git a/test/blockstructured/CMakeLists.txt b/test/blockstructured/CMakeLists.txt
new file mode 100644
index 00000000..981f7857
--- /dev/null
+++ b/test/blockstructured/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(poisson)
\ No newline at end of file
diff --git a/test/blockstructured/poisson/CMakeLists.txt b/test/blockstructured/poisson/CMakeLists.txt
new file mode 100644
index 00000000..e8d17b90
--- /dev/null
+++ b/test/blockstructured/poisson/CMakeLists.txt
@@ -0,0 +1,4 @@
+dune_add_formcompiler_system_test(UFLFILE poisson.ufl
+ BASENAME blockstructured_poisson
+ INIFILE poisson.mini
+ )
\ No newline at end of file
diff --git a/test/blockstructured/poisson/poisson.mini b/test/blockstructured/poisson/poisson.mini
new file mode 100644
index 00000000..af1551d4
--- /dev/null
+++ b/test/blockstructured/poisson/poisson.mini
@@ -0,0 +1,19 @@
+__name = blockstructured_poisson_{__exec_suffix}
+__exec_suffix = numdiff, symdiff | expand num
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 10 10
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_ref
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-7
+blockstructured = 1
+number_of_blocks = 3
diff --git a/test/blockstructured/poisson/poisson.ufl b/test/blockstructured/poisson/poisson.ufl
new file mode 100644
index 00000000..a25c3a8d
--- /dev/null
+++ b/test/blockstructured/poisson/poisson.ufl
@@ -0,0 +1,13 @@
+cell = triangle
+
+x = SpatialCoordinate(cell)
+
+g = x[0]**2 + x[1]**2
+f = -4
+
+V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx]
--
GitLab