...
 
Commits (5)
......@@ -7,6 +7,7 @@ from pymbolic.mapper.substitutor import make_subst_func
import loopy as lp
import pymbolic.primitives as prim
from pytools import memoize_method
from six.moves import intern
......@@ -152,9 +153,9 @@ def get_dependencies(expr):
# This class assumes that no assignments writes to its indices
class Assignment(lp.Assignment):
def assignee_subscript_deps(self):
if isinstance(self.assignee, prim.Subscript):
from loopy.symbolic import get_dependencies
return (get_dependencies(self.assignee.aggregate),)
else:
return frozenset()
@memoize_method
def write_dependency_names(self):
"""Return a set of dependencies of the left hand side of the
assignments performed by this instruction, (NOT!!!!) including indices.
"""
return frozenset(self.assignee_var_names())
......@@ -116,6 +116,10 @@ use_sde:
type: boolean
default: False
helpstr: "Use sde instead of own performance measurements."
use_extruded_mesh:
type: boolean
default: False
helpstr: ""
with_mpi:
type: boolean
default: True
......
......@@ -251,7 +251,10 @@ def generate_driver():
name_driver_block()
# Entrypoint for driver generation
if get_option("opcounter") or get_option("performance_measuring"):
if get_option("use_extruded_mesh"):
from dune.codegen.pdelab.driver.extrude import extruded_test
extruded_test()
elif get_option("opcounter") or get_option("performance_measuring"):
if get_option("performance_measuring"):
assert(not get_option("opcounter"))
assert(isQuadrilateral(get_cell()))
......
from dune.codegen.pdelab.driver.gridfunctionspace import name_test_gfs
from dune.codegen.pdelab.driver.solve import type_vector, name_vector
from dune.codegen.pdelab.driver import get_form_ident
from dune.codegen.pdelab.driver.gridoperator import name_gridoperator
from dune.codegen.generation import preamble
@preamble(section='solver')
def extruded_test():
vector_type = type_vector(get_form_ident())
vector = name_vector(get_form_ident())
gfs = name_test_gfs()
go = name_gridoperator(get_form_ident())
return ["{0} y_{1}({2}); y_{1} = 0.;".format(vector_type, get_form_ident(), gfs),
"{}.residual({}, y_{});".format(go, vector, get_form_ident())]
......@@ -60,6 +60,8 @@ def main_type_range():
@preamble(section="grid", kernel="main")
def typedef_grid(name):
dim = get_dimension()
if get_option("use_extruded_mesh"):
dim = dim - 1
if isQuadrilateral(get_trial_element().cell()):
range_type = main_type_range()
if get_option("grid_unstructured"):
......
......@@ -78,7 +78,7 @@ def blockstructured_boundary_predicated(visitor, measure, subdomain_id):
if subdomain_id not in ['everywhere', 'otherwise']:
if isinstance(visitor.refinement, HyperCubeRefinement):
subelem_inames = visitor.cell_inames # TODO: maybe these inames should be accessible elsewhere?
subelem_inames = visitor.cell_inames
def iname_equals(iname, i):
return prim.Comparison(prim.Variable(iname), "==", i)
......
......@@ -2,18 +2,39 @@ import pymbolic.primitives as prim
from dune.codegen.generation import (geometry_mixin,
temporary_variable,
instruction,
get_global_context_value,
domain,
kernel_cached,
transform)
from dune.codegen.loopy.symbolic import substitute
from dune.codegen.pdelab.geometry import (world_dimension,
name_cell_geometry,
GeometryMixinBase)
from ufl import VectorElement
from dune.codegen.tools import get_pymbolic_basename, maybe_wrap_subscript
from dune.codegen.ufl.modified_terminals import Restriction
from loopy.match import Writes
from meshstructure import Point
def pymbolic_basis(dim, gradient=False):
_basis = (1 - prim.Variable('x0'), prim.Variable('x0'))
_grad = ((-1,), (1,))
def substitute_dimension(dim):
return tuple(substitute(p, {'x0': prim.Variable('x{}'.format(dim))}) for p in _basis)
coordinate_grad = _grad
coordinate_basis = _basis
for i in range(1, dim):
coordinate_grad = tuple(tuple(gi * p_1d for gi in g) + (b * g_1d[0],)
for p_1d, g_1d in zip(substitute_dimension(i), _grad)
for g, b in zip(coordinate_grad, coordinate_basis))
coordinate_basis = tuple(p * p_1d for p_1d in substitute_dimension(i) for p in coordinate_basis)
if gradient:
return coordinate_grad
else:
return coordinate_basis
@geometry_mixin("structured_refinement")
class StructuredRefinementGeometryMixin(GeometryMixinBase):
......@@ -26,24 +47,13 @@ class StructuredRefinementGeometryMixin(GeometryMixinBase):
super().__init__(*args, **kwargs)
self._add_transformation(self.quadrature_inames())
# TODO: currently only VectorElement Q1 allowed
_basis = (1 - prim.Variable('x0'), prim.Variable('x0'))
_grad = ((-1,), (1,))
def substitute_dimension(dim):
return tuple(substitute(p, {'x0': prim.Variable('x{}'.format(dim))}) for p in _basis)
self.coordinate_grad = _grad
self.coordinate_basis = _basis
for i in range(1, world_dimension()):
self.coordinate_grad = tuple(tuple(gi * p_1d for gi in g) + (b * g_1d[0],)
for p_1d, g_1d in zip(substitute_dimension(i), _grad)
for g, b in zip(self.coordinate_grad, self.coordinate_basis))
self.coordinate_basis = tuple(p * p_1d for p_1d in substitute_dimension(i) for p in self.coordinate_basis)
from dune.codegen.structured_refinement.tools import construct_geometry
self.implementation = construct_geometry(self.refinement,
VectorElement("Lagrange", self.refinement.cell, 1,
dim=self.refinement.cell.geometric_dimension()))
def geometry_dofs(self):
raise NotImplementedError
return self.implementation.geometry_dofs(Point(self.source_set.indices, self.source_set))
@kernel_cached
def _init_transformation(self, name, basis, qp):
......@@ -69,7 +79,7 @@ class StructuredRefinementGeometryMixin(GeometryMixinBase):
**insn_kwargs)
def spatial_coordinate(self, o):
basis = self.coordinate_grad if self.reference_grad else self.coordinate_basis
basis = pymbolic_basis(self.implementation.geometric_dimension, self.reference_grad)
measure = self.measure
......@@ -81,29 +91,3 @@ class StructuredRefinementGeometryMixin(GeometryMixinBase):
self._init_transformation(name, basis, qp)
return prim.Variable(name)
@kernel_cached
def define_element_corners(name, cell):
n_corners = cell.num_vertices()
temporary_variable(name, shape_impl=('fv', 'fv'), shape=(n_corners, cell.geometric_dimension()))
iname = "i_corner"
domain(iname, n_corners)
kernel = get_global_context_value("pdelab_kernel")
restriction = {'volume': Restriction.NONE,
'skeleton': Restriction.POSITIVE,
'boundary': Restriction.POSITIVE}
instruction(code="{0}[{1}] = {2}.corner({3});".format(name, iname, name_cell_geometry(restriction[kernel]), iname),
assignees=frozenset({name}),
within_inames=frozenset({iname}), within_inames_is_final=True,
tags=frozenset({'element_corners'}))
def name_element_corners(cell):
name = "corners"
define_element_corners(name, cell)
return name
......@@ -5,38 +5,43 @@ from dune.codegen.options import get_form_option
from dune.codegen.pdelab.geometry import world_dimension, local_dimension, enforce_boundary_restriction, \
name_in_cell_geometry
from dune.codegen.pdelab.quadrature import quadrature_preamble
from dune.codegen.structured_refinement.geometry import name_element_corners, StructuredRefinementGeometryMixin
from dune.codegen.structured_refinement.geometry import StructuredRefinementGeometryMixin, pymbolic_basis
import numpy as np
from dune.codegen.tools import get_pymbolic_basename
from dune.codegen.ufl.modified_terminals import Restriction
from loopy.match import Writes
from .structured import StructuredMeshGeometry
@geometry_mixin("blockstructured_refinement")
class BlockstructuredRefinementSymbolicGeometryMixin(StructuredRefinementGeometryMixin):
def geometry_dofs(self):
dimension = world_dimension()
vertices_name = name_element_corners(self.refinement.base.cell)
macro_vertices = prim.Variable(vertices_name)
class HyperCubeGeometry(StructuredMeshGeometry):
def geometry_dofs(self, point):
macro_vertices = self.base_geometry.geometry_dofs(point)
inames = self.cell_inames
inames = point.multiindex
reference_corners = [[0], [1]]
for _ in range(1, dimension):
for _ in range(1, self.geometric_dimension):
reference_corners = [r + r_1d for r_1d in [[0], [1]] for r in reference_corners]
micro_vertices = [[(ci + prim.Variable(iname)) / get_form_option('number_of_blocks')
micro_vertices = [[(ci + iname) / get_form_option('number_of_blocks')
for ci, iname in zip(c, inames)] for c in reference_corners]
basis_at_micro_vertices = [[substitute(b, dict(list(zip(['x{}'.format(i) for i in range(dimension)], v))))
for b in self.coordinate_basis] for v in micro_vertices]
coordinate_basis = pymbolic_basis(self.geometric_dimension, gradient=False)
basis_at_micro_vertices = [[substitute(b, dict(list(zip(['x{}'.format(i)
for i in range(self.geometric_dimension)], v))))
for b in coordinate_basis] for v in micro_vertices]
vertices = [[prim.CommonSubexpression(sum(prim.Subscript(macro_vertices, (i, j)) * b for i, b in enumerate(basis)),)
for j in range(dimension)] for basis in basis_at_micro_vertices]
vertices = [[prim.CommonSubexpression(sum(macro_vertices[i][j] * b
for i, b in enumerate(basis)),)
for j in range(self.geometric_dimension)] for basis in basis_at_micro_vertices]
return vertices
@geometry_mixin("blockstructured_refinement")
class BlockstructuredRefinementSymbolicGeometryMixin(StructuredRefinementGeometryMixin):
def to_cell(self, local):
if local_dimension() < world_dimension():
return self._to_cell(local)
......
......@@ -5,13 +5,16 @@ from dune.codegen.generation import geometry_mixin, temporary_variable, instruct
get_counted_variable, kernel_cached, get_global_context_value
from dune.codegen.pdelab.geometry import enforce_boundary_restriction, world_dimension, \
local_dimension
from dune.codegen.structured_refinement.geometry import StructuredRefinementGeometryMixin, name_element_corners
from dune.codegen.structured_refinement.geometry import StructuredRefinementGeometryMixin
import pymbolic.primitives as prim
import numpy as np
from dune.codegen.tools import get_pymbolic_basename, maybe_wrap_subscript
from meshstructure import Point, TensorProductEntitySet
from .structured import StructuredMeshGeometry
@kernel_cached
def copy_base(src, dst, bound, additional_inames):
......@@ -28,29 +31,24 @@ def copy_base(src, dst, bound, additional_inames):
depends_on=frozenset({Writes(get_pymbolic_basename(src))}))
@geometry_mixin("extruded_refinement")
class ExtrudedRefinementGeometryMixin(StructuredRefinementGeometryMixin):
# TODO: cannot cache since self.cell_inames might change within volume kernel
def geometry_dofs(self):
vertices_name = name_element_corners(self.refinement.base.cell)
macro_vertices = prim.Variable(vertices_name)
class MeshExtrusionGeometry(StructuredMeshGeometry):
def geometry_dofs(self, point=None):
base_point = Point(point.multiindex[:-1], TensorProductEntitySet(*point.entity_set.factors[:-1],
variant_tag=point.entity_set.variant_tag))
macro_vertices = self.base_geometry.geometry_dofs(base_point)
nlevel = self.refinement.nlevel
nlevel = self.topology.nlevel
if self.measure == 'bottom_facet':
level = 0
elif self.measure == 'top_facet':
level = nlevel - 1
else:
level = prim.Variable(self.cell_inames[-1]) # TODO: always assume extrusion iname is the last one
level = point.multiindex[-1] # TODO: always assume extrusion iname is the last one
lower_dofs = [[prim.Subscript(macro_vertices, (i, j)) for j in range(world_dimension() - 1)] +
[get_height() / nlevel * level] for i in range(2)]
upper_dofs = [[prim.Subscript(macro_vertices, (i, j)) for j in range(world_dimension() - 1)] +
[get_height() / nlevel * (level + 1)] for i in range(2)]
lower_dofs = [vertex + [get_height() / nlevel * level] for vertex in macro_vertices]
upper_dofs = [vertex + [get_height() / nlevel * (level + 1)] for vertex in macro_vertices]
dofs = lower_dofs + upper_dofs
return dofs
@geometry_mixin("extruded_refinement")
class ExtrudedRefinementGeometryMixin(StructuredRefinementGeometryMixin):
def to_cell(self, local):
return local
......
from ufl import VectorElement
from meshstructure import MeshGeometry
from meshstructure.topology import StructuredMeshTopology
class StructuredMeshGeometry(MeshGeometry):
def __init__(self, topology, element):
assert isinstance(topology, StructuredMeshTopology)
from dune.codegen.structured_refinement.tools import construct_geometry
base = topology.base
self.base_geometry = construct_geometry(base,
VectorElement("Lagrange", base.cell, 1,
dim=base.cell.geometric_dimension())
)
super().__init__(topology, element)
\ No newline at end of file
from dune.codegen.generation import kernel_cached, temporary_variable, domain, get_global_context_value, instruction
from dune.codegen.pdelab.geometry import name_cell_geometry
from dune.codegen.ufl.modified_terminals import Restriction
from meshstructure.geometry import MeshGeometry
import pymbolic.primitives as prim
class UnstructuredGeometry(MeshGeometry):
def geometry_dofs(self, point):
vertices_name = name_element_corners(self.topology.cell)
vertex_set, = self.topology.entity_variants(codimension=self.topological_dimension)
return [[prim.Subscript(prim.Variable(vertices_name), (i, j)) for j in range(self.geometric_dimension)]
for i in range(vertex_set.size)]
@kernel_cached
def define_element_corners(name, cell):
n_corners = cell.num_vertices()
temporary_variable(name, shape_impl=('fv', 'fv'), shape=(n_corners, cell.geometric_dimension()))
iname = "i_corner"
domain(iname, n_corners)
kernel = get_global_context_value("pdelab_kernel")
restriction = {'volume': Restriction.NONE,
'skeleton': Restriction.POSITIVE,
'boundary': Restriction.POSITIVE}
instruction(code="{0}[{1}] = {2}.corner({3});".format(name, iname, name_cell_geometry(restriction[kernel]), iname),
assignees=frozenset({name}),
within_inames=frozenset({iname}), within_inames_is_final=True,
tags=frozenset({'element_corners'}))
def name_element_corners(cell):
name = "corners"
define_element_corners(name, cell)
return name
\ No newline at end of file
......@@ -29,7 +29,6 @@ class StructuredRefinementQuadratureMixin(GenericQuadratureMixin):
source_set, = self.refinement.entity_variants(codimension=0)
else:
if local_dimension() < kernel_dimension():
# TODO: here I need the interior of the face set
source_set_unconstrained = self.refinement.entity_variants(codimension=1)[0]
factors = source_set_unconstrained.factors
extrusion_factor = factors[-1]
......
from dune.codegen.generation import (globalarg, preamble, get_global_context_value)
from dune.codegen.pdelab.geometry import world_dimension, local_dimension
from dune.codegen.options import get_form_option
from dune.codegen.structured_refinement.implementations.blockstructured_geometry import HyperCubeGeometry
from dune.codegen.structured_refinement.implementations.extruded_geometry import MeshExtrusionGeometry
from dune.codegen.structured_refinement.implementations.unstructured import UnstructuredGeometry
from meshstructure import MeshExtrusion, HyperCubeRefinement
from meshstructure.topology import UnstructuredMeshTopology
def kernel_dimension():
......@@ -37,6 +43,17 @@ def name_container_alias(container, lfs, element):
return name
def construct_geometry(topology, element):
if isinstance(topology, UnstructuredMeshTopology):
return UnstructuredGeometry(topology, element)
elif isinstance(topology, MeshExtrusion):
return MeshExtrusionGeometry(topology, element)
elif isinstance(topology, HyperCubeRefinement):
return HyperCubeGeometry(topology, element)
else:
raise NotImplementedError()
def construct_visitor():
basis_mixins = ["structured_refinement"]
accumulation_mixins = ["structured_refinement"]
......
......@@ -11,7 +11,7 @@ from dune.codegen.pdelab.localoperator import generate_kernel, generate_kernels_
from dune.codegen.pdelab.signatures import (assembly_routine_args,
assembly_routine_signature,
kernel_name,
)
assembler_routine_name)
from dune.codegen.options import get_form_option, get_option, form_option_context
from dune.codegen.cgen.clazz import ClassMember
......@@ -36,7 +36,8 @@ def sumfact_generate_kernels_per_integral(integrals):
# Generate all necessary kernels
for facedir in range(dim):
for facemod in range(2):
with global_context(facedir_s=facedir, facemod_s=facemod):
with global_context(facedir_s=facedir, facemod_s=facemod,
kernel=get_kernel_name(facedir, facemod)):
yield generate_kernel(integrals)
# Generate switch statement
......@@ -49,7 +50,8 @@ def sumfact_generate_kernels_per_integral(integrals):
for facedir_n in range(dim):
for facemod_n in range(2):
if decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n):
with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n,
kernel=get_kernel_name(facedir_s, facemod_s, facedir_n, facemod_n)):
yield generate_kernel(integrals)
# Generate switch statement
......@@ -57,8 +59,15 @@ def sumfact_generate_kernels_per_integral(integrals):
def get_kernel_name(facedir_s=None, facemod_s=None, facedir_n=None, facemod_n=None):
with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n):
return kernel_name()
arn = assembler_routine_name()
suffix = "{}{}{}{}".format("_facedirs{}".format(facedir_s) if facedir_s is not None else "",
"_facedirn{}".format(facedir_n) if facedir_n is not None else "",
"_facemods{}".format(facemod_s) if facemod_s is not None else "",
"_facemodn{}".format(facemod_n) if facemod_n is not None else "",
)
return "{}{}".format(arn, suffix)
def decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
......
__name = extruded_lineartransport
cells = 20 20
extension = 1. 1.
cells = 20
extension = 1.
[instat]
T = 1.5
......
......@@ -7,8 +7,9 @@ x = SpatialCoordinate(cell)
dirichlet = conditional(x[0] + x[1] < 1.0, 1, 0)
line = UnstructuredHyperCube(1)
structured = HyperCubeRefinement(line, 4)
metadata = {"refinement": MeshExtrusion(line, 10), "levels": 10, "height": 1., "quadrature_order": 2, "apply_cse": True}
metadata = {"refinement": MeshExtrusion(structured, 10), "levels": 10, "height": 1., "quadrature_order": 2, "apply_cse": True}
dx = dx(metadata=metadata)
dS = dS(metadata=metadata)
......