From 318fb56e9bda97f8bf7cae0dc8a210b1bdede380 Mon Sep 17 00:00:00 2001
From: Marcel Koch <marcel.koch@web.de>
Date: Mon, 10 Jul 2017 14:42:16 +0200
Subject: [PATCH] adds computation of index in macro cell
---
.../dune/perftool/blockstructured/__init__.py | 3 +-
.../perftool/blockstructured/accumulation.py | 82 +++++++++++--------
.../dune/perftool/blockstructured/argument.py | 23 ++++++
.../perftool/blockstructured/quadrature.py | 30 ++++---
python/dune/perftool/blockstructured/tools.py | 22 +++++
python/dune/perftool/pdelab/__init__.py | 1 +
python/dune/perftool/pdelab/argument.py | 2 +
python/dune/perftool/pdelab/basis.py | 8 +-
python/dune/perftool/pdelab/driver.py | 12 ++-
python/dune/perftool/pdelab/geometry.py | 5 +-
python/dune/perftool/pdelab/quadrature.py | 5 +-
python/loopy | 2 +-
test/blockstructured/poisson/poisson.mini | 8 +-
test/blockstructured/poisson/poisson.ufl | 2 +-
14 files changed, 142 insertions(+), 63 deletions(-)
create mode 100644 python/dune/perftool/blockstructured/argument.py
create mode 100644 python/dune/perftool/blockstructured/tools.py
diff --git a/python/dune/perftool/blockstructured/__init__.py b/python/dune/perftool/blockstructured/__init__.py
index 7208ad95..85836626 100644
--- a/python/dune/perftool/blockstructured/__init__.py
+++ b/python/dune/perftool/blockstructured/__init__.py
@@ -1,2 +1,3 @@
import dune.perftool.blockstructured.accumulation
-import dune.perftool.blockstructured.quadrature
\ No newline at end of file
+import dune.perftool.blockstructured.quadrature
+import dune.perftool.blockstructured.argument
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/accumulation.py b/python/dune/perftool/blockstructured/accumulation.py
index 255eaaa0..a09ca22e 100644
--- a/python/dune/perftool/blockstructured/accumulation.py
+++ b/python/dune/perftool/blockstructured/accumulation.py
@@ -1,13 +1,29 @@
-from dune.perftool.generation import (backend,iname,domain)
+from dune.perftool.generation import (backend,
+ iname,
+ domain,
+ global_context,
+ get_global_context_value)
+from dune.perftool.options import get_option
from dune.perftool.generation.loopy import transform
+from dune.perftool.generation.counter import get_counted_variable
from dune.perftool.pdelab.localoperator import (determine_accumulation_space,
boundary_predicates)
+from dune.perftool.pdelab.geometry import world_dimension
+from dune.perftool.blockstructured.tools import micro_index_to_macro_index
+import pymbolic.primitives as prim
@iname
def sub_element_inames():
name = "subel"
- domain(name, 4)
- return (name,)
+ dim = world_dimension()
+ inames = tuple()
+ for i in range(dim):
+ name_counted = get_counted_variable(name)
+ inames = inames + (name_counted,)
+ domain(name_counted, get_option("number_of_blocks"))
+ return inames
+
+
@backend(interface="accum_insn", name='blockstructured')
@@ -16,44 +32,42 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
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
+ with global_context(subelem_inames=sub_element_inames()):
+ pymbolic_expr = visitor(accterm.term)
- # Collect the lfs and lfs indices for the accumulate call
- test_lfs = determine_accumulation_space(accterm.term, 0, measure)
+ pymbolic_expr = prim.Quotient(pymbolic_expr, 16)
- # In the jacobian case, also determine the space for the ansatz space
- ansatz_lfs = determine_accumulation_space(accterm.term, 1, measure)
+ # 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
- from dune.perftool.pdelab.argument import name_accumulation_variable
- accumvar = name_accumulation_variable((test_lfs.get_restriction() + ansatz_lfs.get_restriction()))
+ # Collect the lfs and lfs indices for the accumulate call
+ test_lfs = determine_accumulation_space(accterm.term, 0, measure)
- predicates = boundary_predicates(accterm.term, measure, subdomain_id, visitor)
+ # In the jacobian case, also determine the space for the ansatz space
+ ansatz_lfs = determine_accumulation_space(accterm.term, 1, measure)
- rank = 1 if ansatz_lfs.lfs is None else 2
+ # 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()))
- 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,))
- )
+ predicates = boundary_predicates(accterm.term, measure, subdomain_id, visitor)
- from dune.perftool.generation import instruction
- from dune.perftool.options import option_switch
- quad_inames = visitor.interface.quadrature_inames()
+ rank = 1 if ansatz_lfs.lfs is None else 2
- instruction(assignees=(),
- expression=expr,
- forced_iname_deps=frozenset(visitor.inames).union(frozenset(quad_inames)),
- forced_iname_deps_is_final=True,
- predicates=predicates
- )
+ from dune.perftool.pdelab.argument import PDELabAccumulationFunction
+ from pymbolic.primitives import Call
+ expr = Call(PDELabAccumulationFunction(accumvar, rank),
+ ((test_lfs.get_args()[0], micro_index_to_macro_index(test_lfs.get_args()[1])) + ansatz_lfs.get_args() + (pymbolic_expr,))
+ )
- subelem_inames = sub_element_inames()
+ from dune.perftool.generation import instruction
+ from dune.perftool.options import option_switch
+ quad_inames = visitor.interface.quadrature_inames()
- from dune.perftool.sumfact.quadrature import nest_quadrature_loops
- transform(nest_quadrature_loops, subelem_inames)
\ No newline at end of file
+ instruction(assignees=(),
+ expression=expr,
+ forced_iname_deps=frozenset(visitor.inames).union(frozenset(quad_inames)),
+ forced_iname_deps_is_final=True,
+ predicates=predicates
+ )
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/argument.py b/python/dune/perftool/blockstructured/argument.py
new file mode 100644
index 00000000..7d24a656
--- /dev/null
+++ b/python/dune/perftool/blockstructured/argument.py
@@ -0,0 +1,23 @@
+from dune.perftool.generation import (backend,
+ kernel_cached,
+ valuearg,get_global_context_value)
+from dune.perftool.pdelab.argument import CoefficientAccess
+from dune.perftool.blockstructured.tools import micro_index_to_macro_index
+from loopy.types import NumpyType
+from pymbolic.primitives import Variable, Call
+import pymbolic.primitives as prim
+
+# TODO rename prim
+@kernel_cached
+@backend(interface="pymbolic_coefficient", name="blockstructured")
+def pymbolic_coefficient(container, lfs, index):
+ # TODO introduce a proper type for local function spaces!
+ if isinstance(lfs, str):
+ valuearg(lfs, dtype=NumpyType("str"))
+
+ # If the LFS is not yet a pymbolic expression, make it one
+ from pymbolic.primitives import Expression
+ if not isinstance(lfs, Expression):
+ lfs = Variable(lfs)
+
+ return Call(CoefficientAccess(container), (lfs, micro_index_to_macro_index(index),))
diff --git a/python/dune/perftool/blockstructured/quadrature.py b/python/dune/perftool/blockstructured/quadrature.py
index cb0d16ba..36fd0acf 100644
--- a/python/dune/perftool/blockstructured/quadrature.py
+++ b/python/dune/perftool/blockstructured/quadrature.py
@@ -1,6 +1,8 @@
from dune.perftool.generation import (backend,
temporary_variable,
- domain)
+ domain,
+ get_backend,
+ get_global_context_value)
from dune.perftool.pdelab.quadrature import (name_quadrature_points,
quadrature_iname
)
@@ -26,15 +28,21 @@ def pymbolic_quadrature_position():
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'))
+ subelem_inames = get_global_context_value("subelem_inames")
+ for i in range(dim):
+ expr = prim.Subscript(qp_in_micro.aggregate, (qp_in_micro.index,) + (i,))
+ expr = prim.Sum((expr, prim.Variable(subelem_inames[i]),))
+ expr = prim.Quotient(expr, get_option('number_of_blocks'))
+ instruction(assignee=Subscript(Variable(qp_in_macro), (i,)),
+ expression=expr,
+ forced_iname_deps_is_final=True,
+ tags=frozenset({'quad'})
+ )
+ return Variable(qp_in_macro)
- 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)
+@backend(interface="qp_in_cell", name="blockstructured")
+def pymbolic_quadrature_position_in_cell(restriction):
+ from dune.perftool.pdelab.geometry import to_cell_coordinates
+ quad_pos = get_backend("quad_pos", selector=lambda: "blockstructured")()
+ return to_cell_coordinates(quad_pos, restriction)
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/tools.py b/python/dune/perftool/blockstructured/tools.py
new file mode 100644
index 00000000..f43e0df1
--- /dev/null
+++ b/python/dune/perftool/blockstructured/tools.py
@@ -0,0 +1,22 @@
+from dune.perftool.generation import get_global_context_value
+from dune.perftool.options import get_option
+import pymbolic.primitives as prim
+
+
+# TODO 3d
+def micro_index_to_macro_index(index):
+ if not isinstance(index, prim.Variable):
+ index = prim.Variable(index)
+
+ subelem_inames = get_global_context_value('subelem_inames')
+
+ k = get_option("number_of_blocks")
+
+ modified_index = prim.Product(((k+1), prim.Variable(subelem_inames[1])))
+ modified_index = prim.Sum((modified_index, prim.Variable(subelem_inames[0]), index))
+ index_div_2 = prim.FloorDiv(index, 2)
+ index_div_2 = prim.Product((index_div_2, k-1))
+ modified_index = prim.Sum((modified_index, index_div_2))
+
+ return modified_index
+
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 33571d68..4d7c050d 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -6,6 +6,7 @@ from dune.perftool.pdelab.argument import (pymbolic_apply_function,
pymbolic_apply_function_gradient,
pymbolic_trialfunction,
pymbolic_trialfunction_gradient,
+ pymbolic_coefficient
)
from dune.perftool.pdelab.basis import (pymbolic_basis,
pymbolic_reference_gradient,
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index f9906902..0d6d4ea7 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -13,6 +13,7 @@ from dune.perftool.generation import (domain,
valuearg,
get_global_context_value,
kernel_cached,
+ backend
)
from dune.perftool.pdelab.index import name_index
from dune.perftool.pdelab.basis import (evaluate_coefficient,
@@ -132,6 +133,7 @@ def name_applycontainer(restriction):
@kernel_cached
+@backend(interface="pymbolic_coefficient")
def pymbolic_coefficient(container, lfs, index):
# TODO introduce a proper type for local function spaces!
if isinstance(lfs, str):
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index f79f4df9..ccbb97f7 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -10,6 +10,7 @@ from dune.perftool.generation import (backend,
preamble,
temporary_variable,
)
+from dune.perftool.options import get_option
from dune.perftool.pdelab.spaces import (lfs_child,
lfs_iname,
lfs_inames,
@@ -32,6 +33,7 @@ from dune.perftool.tools import (get_pymbolic_basename,
from dune.perftool.pdelab.driver import FEM_name_mangling
from dune.perftool.pdelab.restriction import restricted_name
from pymbolic.primitives import Product, Subscript, Variable
+import pymbolic.primitives as prim
from loopy import Reduction
@@ -192,7 +194,10 @@ def evaluate_coefficient_gradient(element, name, container, restriction, tree_pa
lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]), symmetry=element.symmetry())
from dune.perftool.pdelab.argument import pymbolic_coefficient
- coeff = pymbolic_coefficient(container, lfs, index)
+ if get_option('blockstructured'):
+ coeff = get_backend("pymbolic_coefficient", selector=lambda: "blockstructured")(container, lfs, index)
+ else:
+ coeff = get_backend("pymbolic_coefficient")(container, lfs, index)
assignee = Subscript(Variable(name), tuple(Variable(i) for i in idims))
reduction_expr = Product((coeff, Subscript(Variable(get_pymbolic_basename(basis)), basis.index + (Variable(idims[-1]),))))
@@ -201,4 +206,5 @@ def evaluate_coefficient_gradient(element, name, container, restriction, tree_pa
assignee=assignee,
forced_iname_deps=frozenset(get_backend("quad_inames")()).union(frozenset(idims)),
forced_iname_deps_is_final=True,
+ tags=frozenset({"quad"})
)
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 16c50f88..7f6b8b8d 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -305,20 +305,24 @@ def typedef_fem(expr, name):
df = type_domainfield()
r = type_range()
dim = name_dimension()
+ if get_option("blockstructured"):
+ degree = expr._degree + get_option("number_of_blocks")
+ else:
+ degree = expr._degree
if isPk(expr):
include_file("dune/pdelab/finiteelementmap/pkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, expr._degree)
+ return "using {} = Dune::PDELab::PkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, degree)
if isQk(expr):
include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, expr._degree)
+ return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, degree)
if isDG(expr):
if isQuadrilateral(expr):
include_file("dune/pdelab/finiteelementmap/qkdg.hh", filetag="driver")
# TODO allow switching the basis here!
- return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, df, r, expr._degree, dim)
+ return "using {} = Dune::PDELab::QkDGLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, df, r, degree, dim)
if isSimplical(expr):
include_file("dune/pdelab/finiteelementmap/opbfem.hh", filetag="driver")
- return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;".format(name, df, r, expr._degree, dim)
+ return "using {} = Dune::PDELab::OPBLocalFiniteElementMap<{}, {}, {}, {}, Dune::GeometryType::simplex>;".format(name, df, r, degree, dim)
raise NotImplementedError("Geometry type not known in code generation")
raise NotImplementedError("FEM not implemented in dune-perftool")
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index c365e89e..9e4c7482 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -333,7 +333,10 @@ def define_jacobian_inverse_transposed(name, restriction):
dim = world_dimension()
temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
geo = name_cell_geometry(restriction)
- pos = get_backend("qp_in_cell")(restriction)
+ if get_option('blockstructured'):
+ pos = get_backend("qp_in_cell", selector=lambda : "blockstructured")(restriction)
+ else:
+ pos = get_backend("qp_in_cell")(restriction)
return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
geo,
str(pos),
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 1cd0d366..8356ecd2 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -125,10 +125,7 @@ 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
- if get_option("blockstructured"):
- quad_pos = get_backend("quad_pos", selector=lambda: "blockstructured")()
- else:
- quad_pos = get_backend("quad_pos")()
+ quad_pos = get_backend("quad_pos")()
return to_cell_coordinates(quad_pos, restriction)
diff --git a/python/loopy b/python/loopy
index 2f0cac0b..b6400046 160000
--- a/python/loopy
+++ b/python/loopy
@@ -1 +1 @@
-Subproject commit 2f0cac0bc927fe2b3394801419f88a781c2150e9
+Subproject commit b64000460e99bc7411abf3ca6c616ef5ddc5567d
diff --git a/test/blockstructured/poisson/poisson.mini b/test/blockstructured/poisson/poisson.mini
index af1551d4..adb4a949 100644
--- a/test/blockstructured/poisson/poisson.mini
+++ b/test/blockstructured/poisson/poisson.mini
@@ -1,10 +1,8 @@
__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
+cells = 1 1
+extension = 1. 1.
[wrapper.vtkcompare]
name = {__name}
@@ -16,4 +14,4 @@ numerical_jacobian = 1, 0 | expand num
exact_solution_expression = g
compare_l2errorsquared = 1e-7
blockstructured = 1
-number_of_blocks = 3
+number_of_blocks = 4
diff --git a/test/blockstructured/poisson/poisson.ufl b/test/blockstructured/poisson/poisson.ufl
index a25c3a8d..128cd990 100644
--- a/test/blockstructured/poisson/poisson.ufl
+++ b/test/blockstructured/poisson/poisson.ufl
@@ -1,4 +1,4 @@
-cell = triangle
+cell = quadrilateral
x = SpatialCoordinate(cell)
--
GitLab