Skip to content

Commit

Permalink
correct test
Browse files Browse the repository at this point in the history
  • Loading branch information
@mscroggs
mscroggs committed Dec 20, 2024
1 parent 033475b commit 54ef0c4
Show file tree
Hide file tree
Showing 5 changed files with 76 additions and 59 deletions.
3 changes: 0 additions & 3 deletions ffcx/codegeneration/jit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,9 +370,6 @@ def _compile_objects(
logger.info("Calling JIT C compiler")
logger.info(79 * "#")

## TODO: remove
cffi_final_compile_args.remove("-Werror")

t0 = time.time()
f = io.StringIO()
# Temporarily set root logger handlers to string buffer only
Expand Down
4 changes: 3 additions & 1 deletion ffcx/ir/integral.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,9 @@ def compute_integral_ir(cell, integral_type, entity_type, integrands, argument_s
block_restrictions = tuple(block_restrictions)

# Check if each *each* factor corresponding to this argument is piecewise
all_factors_piecewise = all(F.nodes[ifi[0]]["status"] == "piecewise" for ifi in fi_ci)
all_factors_piecewise = all(
F.nodes[ifi[0]]["status"] == "piecewise" for ifi in fi_ci
)
block_is_permuted = False
for name in unames:
if tables[name].shape[0] > 1:
Expand Down
27 changes: 18 additions & 9 deletions ffcx/ir/representation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,8 +148,7 @@ def compute_ir(
ir_integrals = list(itertools.chain(*irs))

integral_domains = {
i.expression.name: set(j[0] for j in i.expression.integrand.keys())
for a in irs for i in a
i.expression.name: set(j[0] for j in i.expression.integrand.keys()) for a in irs for i in a
}

ir_forms = [
Expand Down Expand Up @@ -287,8 +286,10 @@ def _compute_integral_ir(
f"but requested degree is {degree}."
)
points = basix.cell.geometry(getattr(basix.CellType, cellname))
weights = np.array([1.0 / points.shape[0] / basix.cell.volume(
getattr(basix.CellType, cellname))] * points.shape[0])
weights = np.array(
[1.0 / points.shape[0] / basix.cell.volume(getattr(basix.CellType, cellname))]
* points.shape[0]
)
rules[cellname] = (points, weights, None)
else:
degree = md["quadrature_degree"]
Expand All @@ -300,7 +301,14 @@ def _compute_integral_ir(
form_data.argument_elements,
use_sum_factorization,
)
rules = {i: (points[i], weights[i], None if tensor_factors is None else tensor_factors[i]) for i in points}
rules = {
i: (
points[i],
weights[i],
None if tensor_factors is None else tensor_factors[i],
)
for i in points
}

for cell, (points, weights, tensor_factors) in rules.items():
points = np.asarray(points)
Expand All @@ -312,7 +320,9 @@ def _compute_integral_ir(
if rule not in grouped_integrands:
grouped_integrands[cell][rule] = []
grouped_integrands[cell][rule].append(integral.integrand())
sorted_integrals: dict[str, dict[QuadratureRule, Integral]] = {cell: {} for cell in grouped_integrands}
sorted_integrals: dict[str, dict[QuadratureRule, Integral]] = {
cell: {} for cell in grouped_integrands
}
for cell, integrands_by_cell in grouped_integrands.items():
for rule, integrands in integrands_by_cell.items():
integrands_summed = sorted_expr_sum(integrands)
Expand Down Expand Up @@ -359,9 +369,8 @@ def _compute_integral_ir(

# Create map from number of quadrature points -> integrand
integrand_map: dict[str, dict[QuadratureRule, ufl.core.expr.Expr]] = {
cell: {
rule: integral.integrand() for rule, integral in cell_integrals.items()
} for cell, cell_integrals in sorted_integrals.items()
cell: {rule: integral.integrand() for rule, integral in cell_integrals.items()}
for cell, cell_integrals in sorted_integrals.items()
}

# Build more specific intermediate representation
Expand Down
13 changes: 10 additions & 3 deletions ffcx/ir/representationutils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,13 +72,20 @@ def create_quadrature_points_and_weights(
pts["interval"] = np.array(
[tuple(i[0] for i in p) for p in itertools.product(*[f[0] for f in tensor_factors])]
)
wts["interval"] = np.array([np.prod(p) for p in itertools.product(*[f[1] for f in tensor_factors])])
wts["interval"] = np.array(
[np.prod(p) for p in itertools.product(*[f[1] for f in tensor_factors])]
)
else:
pts[cell.cellname()], wts[cell.cellname()] = create_quadrature(cell.cellname(), degree, rule, elements)
pts[cell.cellname()], wts[cell.cellname()] = create_quadrature(
cell.cellname(), degree, rule, elements
)
elif integral_type in ufl.measure.facet_integral_types:
for ft in cell.facet_types():
pts[ft.cellname()], wts[ft.cellname()] = create_quadrature(
ft.cellname(), degree, rule, elements,
ft.cellname(),
degree,
rule,
elements,
)
elif integral_type in ufl.measure.point_integral_types:
pts["vertex"], wts["vertex"] = create_quadrature("vertex", degree, rule, elements)
Expand Down
88 changes: 45 additions & 43 deletions test/test_jit_forms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1223,6 +1223,7 @@ def tabulate_tensor(ele_type, V_cell_type, W_cell_type, coeffs):

assert np.allclose(A, A_ref)


@pytest.mark.parametrize("dtype", ["float64"])
def test_ds_prism(compile_args, dtype):
element = basix.ufl.element("Lagrange", "prism", 1)
Expand Down Expand Up @@ -1276,10 +1277,17 @@ def test_ds_prism(compile_args, dtype):
entity_index = np.array([0], dtype=int)

xdtype = dtype_to_scalar_dtype(dtype)
coords = np.array([
[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0],
[0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0],
], dtype=xdtype)
coords = np.array(
[
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
[1.0, 0.0, 1.0],
[0.0, 1.0, 1.0],
],
dtype=xdtype,
)

c_type, c_xtype = dtype_to_c_type(dtype), dtype_to_c_type(xdtype)

Expand All @@ -1296,29 +1304,38 @@ def test_ds_prism(compile_args, dtype):

assert np.allclose(
A,
np.array([
[1 / 12, 1 / 24, 1 / 24, 0, 0, 0],
[1 / 24, 1 / 12, 1 / 24, 0, 0, 0],
[1 / 24, 1 / 24, 1 / 12, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
])
np.array(
[
[1 / 12, 1 / 24, 1 / 24, 0, 0, 0],
[1 / 24, 1 / 12, 1 / 24, 0, 0, 0],
[1 / 24, 1 / 24, 1 / 12, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
]
),
)

# Test integral over quadrilateral (facet 1)
A = np.zeros((6, 6), dtype=dtype)
entity_index = np.array([1], dtype=int)
entity_index = np.array([1], dtype=np.int64)

xdtype = dtype_to_scalar_dtype(dtype)
coords = np.array([
[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0],
[0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0],
], dtype=xdtype)
coords = np.array(
[
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
[1.0, 0.0, 1.0],
[0.0, 1.0, 1.0],
],
dtype=xdtype,
)

c_type, c_xtype = dtype_to_c_type(dtype), dtype_to_c_type(xdtype)

kernel = getattr(integral_tri, f"tabulate_tensor_{dtype}")
kernel = getattr(integral_quad, f"tabulate_tensor_{dtype}")

kernel(
ffi.cast(f"{c_type} *", A.ctypes.data),
Expand All @@ -1329,31 +1346,16 @@ def test_ds_prism(compile_args, dtype):
ffi.cast(f"uint8_t *", entity_perm.ctypes.data),
)

print(A)
for i, j in zip(
A,
np.array([
[1/9, 1/18, 0, 1/18, 1/36, 0],
[1/18, 1/9, 0, 1/36, 1/18, 0],
[0, 0, 0, 0, 0, 0],
[1/18, 1/36, 0, 1/9, 1/18, 0],
[1/36, 1/18, 0, 1/18, 1/9, 0],
[0, 0, 0, 0, 0, 0],
])):
for a, b in zip(i, j):
if a > 1/100:
from math import gcd
numerator = int(a * 9 * 32 / b + 0.01)
g = gcd(numerator, 9 * 32)
print(a, b, a / b, a * 9 * 32 / b, "-->", numerator//g, "/", 9*32//g)
assert np.allclose(
A,
np.array([
[1/9, 1/18, 0, 1/18, 1/36, 0],
[1/18, 1/9, 0, 1/36, 1/18, 0],
[0, 0, 0, 0, 0, 0],
[1/18, 1/36, 0, 1/9, 1/18, 0],
[1/36, 1/18, 0, 1/18, 1/9, 0],
[0, 0, 0, 0, 0, 0],
])
np.array(
[
[1 / 9, 1 / 18, 0, 1 / 18, 1 / 36, 0],
[1 / 18, 1 / 9, 0, 1 / 36, 1 / 18, 0],
[0, 0, 0, 0, 0, 0],
[1 / 18, 1 / 36, 0, 1 / 9, 1 / 18, 0],
[1 / 36, 1 / 18, 0, 1 / 18, 1 / 9, 0],
[0, 0, 0, 0, 0, 0],
]
),
)

0 comments on commit 54ef0c4

Please sign in to comment.