Skip to content

Commit

Permalink
Work mostly on Expression
Browse files Browse the repository at this point in the history
  • Loading branch information
chrisrichardson committed Aug 25, 2023
1 parent 3a81fe8 commit b85eb4f
Show file tree
Hide file tree
Showing 5 changed files with 52 additions and 70 deletions.
8 changes: 0 additions & 8 deletions ffcx/codegeneration/access.py
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,6 @@ def cell_coordinate(self, e, mt, tabledata, num_points):
raise RuntimeError("Expecting reference cell coordinate to be symbolically rewritten.")

def facet_coordinate(self, e, mt, tabledata, num_points):
L = self.language
if mt.global_derivatives:
raise RuntimeError("Not expecting derivatives of FacetCoordinate.")
if mt.local_derivatives:
Expand Down Expand Up @@ -185,15 +184,13 @@ def reference_cell_volume(self, e, mt, tabledata, access):
raise RuntimeError(f"Unhandled cell types {cellname}.")

def reference_facet_volume(self, e, mt, tabledata, access):
L = self.language
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("interval", "triangle", "tetrahedron", "quadrilateral", "hexahedron"):
return L.Symbol(f"{cellname}_reference_facet_volume", dtype=L.DataType.REAL)
else:
raise RuntimeError(f"Unhandled cell types {cellname}.")

def reference_normal(self, e, mt, tabledata, access):
L = self.language
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("interval", "triangle", "tetrahedron", "quadrilateral", "hexahedron"):
table = L.Symbol(f"{cellname}_reference_facet_normals", dtype=L.DataType.REAL)
Expand All @@ -203,7 +200,6 @@ def reference_normal(self, e, mt, tabledata, access):
raise RuntimeError(f"Unhandled cell types {cellname}.")

def cell_facet_jacobian(self, e, mt, tabledata, num_points):
L = self.language
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("triangle", "tetrahedron", "quadrilateral", "hexahedron"):
table = L.Symbol(f"{cellname}_reference_facet_jacobian", dtype=L.DataType.REAL)
Expand All @@ -215,7 +211,6 @@ def cell_facet_jacobian(self, e, mt, tabledata, num_points):
raise RuntimeError(f"Unhandled cell types {cellname}.")

def reference_cell_edge_vectors(self, e, mt, tabledata, num_points):
L = self.language
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("triangle", "tetrahedron", "quadrilateral", "hexahedron"):
table = L.Symbol(f"{cellname}_reference_edge_vectors", dtype=L.DataType.REAL)
Expand All @@ -226,7 +221,6 @@ def reference_cell_edge_vectors(self, e, mt, tabledata, num_points):
raise RuntimeError(f"Unhandled cell types {cellname}.")

def reference_facet_edge_vectors(self, e, mt, tabledata, num_points):
L = self.language
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname in ("tetrahedron", "hexahedron"):
table = L.Symbol(f"{cellname}_reference_edge_vectors", dtype=L.DataType.REAL)
Expand All @@ -240,7 +234,6 @@ def reference_facet_edge_vectors(self, e, mt, tabledata, num_points):
raise RuntimeError(f"Unhandled cell types {cellname}.")

def facet_orientation(self, e, mt, tabledata, num_points):
L = self.language
cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
if cellname not in ("interval", "triangle", "tetrahedron"):
raise RuntimeError(f"Unhandled cell types {cellname}.")
Expand Down Expand Up @@ -312,7 +305,6 @@ def cell_edge_vectors(self, e, mt, tabledata, num_points):
)

def facet_edge_vectors(self, e, mt, tabledata, num_points):
L = self.language

# Get properties of domain
domain = ufl.domain.extract_unique_domain(mt.terminal)
Expand Down
4 changes: 2 additions & 2 deletions ffcx/codegeneration/definitions.py
Original file line number Diff line number Diff line change
Expand Up @@ -96,8 +96,8 @@ def coefficient(self, t, mt, tabledata, quadrature_rule, access):

# If a map is necessary from stride 1 to bs, the code must be added before the quadrature loop.
if dof_access_map:
pre_code += [L.ArrayDecl(self.options["scalar_type"], dof_access.array, num_dofs)]
pre_body = L.Assign(dof_access, dof_access_map)
pre_code += [L.ArrayDecl(dof_access.array, sizes=num_dofs)]
pre_body = [L.Assign(dof_access, dof_access_map)]
pre_code += [L.ForRange(ic, 0, num_dofs, pre_body)]
else:
dof_access = self.symbols.coefficient_dof_access(mt.terminal, ic * bs + begin)
Expand Down
48 changes: 15 additions & 33 deletions ffcx/codegeneration/expression_generator.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@
import ffcx.codegeneration.lnodes as L
from ffcx.codegeneration.lnodes import LNode
from ffcx.ir.representation import ExpressionIR
from ffcx.naming import scalar_to_value_type

logger = logging.getLogger("ffcx")

Expand All @@ -28,22 +27,18 @@ def __init__(self, ir: ExpressionIR, backend: FFCXBackend):

self.ir = ir
self.backend = backend
self.ufl_to_language = L.UFL2LNodes()
self.scope: Dict[Any, LNode] = {}
self._ufl_names: Set[Any] = set()
self.symbol_counters: DefaultDict[Any, int] = collections.defaultdict(int)
self.shared_symbols: Dict[Any, Any] = {}
self.quadrature_rule = list(self.ir.integrand.keys())[0]

def generate(self):
L = self.backend.language

parts = []
scalar_type = self.backend.access.options["scalar_type"]
value_type = scalar_to_value_type(scalar_type)

parts += self.generate_element_tables(value_type)
parts += self.generate_element_tables()
# Generate the tables of geometry data that are needed
parts += self.generate_geometry_tables(value_type)
parts += self.generate_geometry_tables()
parts += self.generate_piecewise_partition()

all_preparts = []
Expand Down Expand Up @@ -78,11 +73,11 @@ def generate_geometry_tables(self):
parts = []
for i, cell_list in cells.items():
for c in cell_list:
parts.append(geometry.write_table(L, ufl_geometry[i], c))
parts.append(geometry.write_table(ufl_geometry[i], c))

return parts

def generate_element_tables(self, float_type: str):
def generate_element_tables(self):
"""Generate tables of FE basis evaluated at specified points."""
parts = []

Expand All @@ -91,7 +86,8 @@ def generate_element_tables(self, float_type: str):

for name in table_names:
table = tables[name]
decl = L.ArrayDecl(name, table)
symbol = L.Symbol(name, dtype=L.DataType.REAL)
decl = L.ArrayDecl(symbol, sizes=table.shape, values=table, const=True)
parts += [decl]

# Add leading comment if there are any tables
Expand All @@ -107,8 +103,6 @@ def generate_quadrature_loop(self):
In the context of expressions quadrature loop is not accumulated.
"""
L = self.backend.language

# Generate varying partition
body = self.generate_varying_partition()
body = L.commented_code_list(
Expand All @@ -133,25 +127,21 @@ def generate_quadrature_loop(self):

def generate_varying_partition(self):
"""Generate factors of blocks which are not cellwise constant."""
L = self.backend.language

# Get annotated graph of factorisation
F = self.ir.integrand[self.quadrature_rule]["factorization"]

arraysymbol = L.Symbol(f"sv_{self.quadrature_rule.id()}")
arraysymbol = L.Symbol(f"sv_{self.quadrature_rule.id()}", dtype=L.DataType.SCALAR)
parts = self.generate_partition(arraysymbol, F, "varying")
parts = L.commented_code_list(
parts, f"Unstructured varying computations for quadrature rule {self.quadrature_rule.id()}")
return parts

def generate_piecewise_partition(self):
"""Generate factors of blocks which are constant (i.e. do not depend on quadrature points)."""
L = self.backend.language

# Get annotated graph of factorisation
F = self.ir.integrand[self.quadrature_rule]["factorization"]

arraysymbol = L.Symbol("sp")
arraysymbol = L.Symbol("sp", dtype=L.DataType.SCALAR)
parts = self.generate_partition(arraysymbol, F, "piecewise")
parts = L.commented_code_list(parts, "Unstructured piecewise computations")
return parts
Expand Down Expand Up @@ -183,8 +173,6 @@ def generate_dofblock_partition(self):

def generate_block_parts(self, blockmap, blockdata):
"""Generate and return code parts for a given block."""
L = self.backend.language

# The parts to return
preparts = []
quadparts = []
Expand Down Expand Up @@ -283,8 +271,6 @@ def get_arg_factors(self, blockdata, block_rank, indices):
Indices used to index element tables
"""
L = self.backend.language

arg_factors = []
for i in range(block_rank):
mad = blockdata.ma_data[i]
Expand All @@ -304,21 +290,18 @@ def get_arg_factors(self, blockdata, block_rank, indices):

def new_temp_symbol(self, basename):
"""Create a new code symbol named basename + running counter."""
L = self.backend.language
name = "%s%d" % (basename, self.symbol_counters[basename])
self.symbol_counters[basename] += 1
return L.Symbol(name)
return L.Symbol(name, dtype=L.DataType.SCALAR)

def get_var(self, v):
if v._ufl_is_literal_:
return self.backend.ufl_to_language.get(v)
return self.ufl_to_language.get(v)
f = self.scope.get(v)
return f

def generate_partition(self, symbol, F, mode):
"""Generate computations of factors of blocks."""
L = self.backend.language

definitions = []
pre_definitions = dict()
intermediates = []
Expand All @@ -332,7 +315,7 @@ def generate_partition(self, symbol, F, mode):
mt = attr.get('mt')

if v._ufl_is_literal_:
vaccess = self.backend.ufl_to_language.get(v)
vaccess = self.ufl_to_language.get(v)
elif mt is not None:
# All finite element based terminals have table data, as well
# as some, but not all, of the symbolic geometric terminals
Expand Down Expand Up @@ -361,7 +344,7 @@ def generate_partition(self, symbol, F, mode):

# Mapping UFL operator to target language
self._ufl_names.add(v._ufl_handler_name_)
vexpr = self.backend.ufl_to_language.get(v, *vops)
vexpr = self.ufl_to_language.get(v, *vops)

# Create a new intermediate for each subexpression
# except boolean conditions and its childs
Expand All @@ -387,7 +370,7 @@ def generate_partition(self, symbol, F, mode):
intermediates.append(L.Assign(vaccess, vexpr))
else:
scalar_type = self.backend.access.options["scalar_type"]
vaccess = L.Symbol("%s_%d" % (symbol.name, j))
vaccess = L.Symbol("%s_%d" % (symbol.name, j), dtype=L.DataType.SCALAR)
intermediates.append(L.VariableDecl(f"const {scalar_type}", vaccess, vexpr))

# Store access node for future reference
Expand All @@ -405,7 +388,6 @@ def generate_partition(self, symbol, F, mode):

if intermediates:
if use_symbol_array:
scalar_type = self.backend.access.options["scalar_type"]
parts += [L.ArrayDecl(scalar_type, symbol, len(intermediates))]
parts += [L.ArrayDecl(symbol, sizes=len(intermediates))]
parts += intermediates
return parts
60 changes: 34 additions & 26 deletions ffcx/codegeneration/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,31 +5,31 @@
# SPDX-License-Identifier: LGPL-3.0-or-later

import numpy as np

import ffcx.codegeneration.lnodes as L
import basix


def write_table(L, tablename, cellname):
def write_table(tablename, cellname):
if tablename == "facet_edge_vertices":
return facet_edge_vertices(L, tablename, cellname)
return facet_edge_vertices(tablename, cellname)
if tablename == "reference_facet_jacobian":
return reference_facet_jacobian(L, tablename, cellname)
return reference_facet_jacobian(tablename, cellname)
if tablename == "reference_cell_volume":
return reference_cell_volume(L, tablename, cellname)
return reference_cell_volume(tablename, cellname)
if tablename == "reference_facet_volume":
return reference_facet_volume(L, tablename, cellname)
return reference_facet_volume(tablename, cellname)
if tablename == "reference_edge_vectors":
return reference_edge_vectors(L, tablename, cellname)
return reference_edge_vectors(tablename, cellname)
if tablename == "facet_reference_edge_vectors":
return facet_reference_edge_vectors(L, tablename, cellname)
return facet_reference_edge_vectors(tablename, cellname)
if tablename == "reference_facet_normals":
return reference_facet_normals(L, tablename, cellname)
return reference_facet_normals(tablename, cellname)
if tablename == "facet_orientation":
return facet_orientation(L, tablename, cellname)
return facet_orientation(tablename, cellname)
raise ValueError(f"Unknown geometry table name: {tablename}")


def facet_edge_vertices(L, tablename, cellname):
def facet_edge_vertices(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
topology = basix.topology(celltype)
triangle_edges = basix.topology(basix.CellType.triangle)[1]
Expand All @@ -48,40 +48,45 @@ def facet_edge_vertices(L, tablename, cellname):
raise ValueError("Only triangular and quadrilateral faces supported.")

out = np.array(edge_vertices, dtype=int)
return L.ArrayDecl("static const unsigned int", f"{cellname}_{tablename}", out.shape, out)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.INT)
return L.ArrayDecl(symbol, values=out, const=True)


def reference_facet_jacobian(L, tablename, cellname, type: str):
def reference_facet_jacobian(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
out = basix.cell.facet_jacobians(celltype)
return L.ArrayDecl(f"static const {type}", f"{cellname}_{tablename}", out.shape, out)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.ArrayDecl(symbol, values=out, const=True)


def reference_cell_volume(L, tablename, cellname, type: str):
def reference_cell_volume(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
out = basix.cell.volume(celltype)
return L.VariableDecl(f"static const {type}", f"{cellname}_{tablename}", out)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.VariableDecl(symbol, out)


def reference_facet_volume(L, tablename, cellname, type: str):
def reference_facet_volume(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
volumes = basix.cell.facet_reference_volumes(celltype)
for i in volumes[1:]:
if not np.isclose(i, volumes[0]):
raise ValueError("Reference facet volume not supported for this cell type.")
return L.VariableDecl(f"static const {type}", f"{cellname}_{tablename}", volumes[0])
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.VariableDecl(symbol, volumes[0])


def reference_edge_vectors(L, tablename, cellname, type: str):
def reference_edge_vectors(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
topology = basix.topology(celltype)
geometry = basix.geometry(celltype)
edge_vectors = [geometry[j] - geometry[i] for i, j in topology[1]]
out = np.array(edge_vectors[cellname])
return L.ArrayDecl(f"static const {type}", f"{cellname}_{tablename}", out.shape, out)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.ArrayDecl(symbol, values=out, const=True)


def facet_reference_edge_vectors(L, tablename, cellname, type: str):
def facet_reference_edge_vectors(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
topology = basix.topology(celltype)
geometry = basix.geometry(celltype)
Expand All @@ -101,16 +106,19 @@ def facet_reference_edge_vectors(L, tablename, cellname, type: str):
raise ValueError("Only triangular and quadrilateral faces supported.")

out = np.array(edge_vectors)
return L.ArrayDecl(f"static const {type}", f"{cellname}_{tablename}", out.shape, out)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.ArrayDecl(symbol, values=out, const=True)


def reference_facet_normals(L, tablename, cellname, type: str):
def reference_facet_normals(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
out = basix.cell.facet_outward_normals(celltype)
return L.ArrayDecl(f"static const {type}", f"{cellname}_{tablename}", out.shape, out)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.ArrayDecl(symbol, values=out, const=True)


def facet_orientation(L, tablename, cellname, type: str):
def facet_orientation(tablename, cellname):
celltype = getattr(basix.CellType, cellname)
out = basix.cell.facet_orientations(celltype)
return L.ArrayDecl(f"static const {type}", f"{cellname}_{tablename}", len(out), out)
symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
return L.ArrayDecl(symbol, values=out, const=True)
2 changes: 1 addition & 1 deletion ffcx/codegeneration/integral_generator.py
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,7 @@ def generate_geometry_tables(self):
parts = []
for i, cell_list in cells.items():
for c in cell_list:
parts.append(geometry.write_table(L, ufl_geometry[i], c))
parts.append(geometry.write_table(ufl_geometry[i], c))

return parts

Expand Down

0 comments on commit b85eb4f

Please sign in to comment.