Merge branch 'main' into dokken/group_integrals_v2

Author: mscroggs

ffcx/analysis.py

@@ -181,8 +181,8 @@ def _analyze_form(form: ufl.form.Form, options: typing.Dict) -> ufl.algorithms.f
                 assert np.allclose(e._points, custom_q[0])
                 assert np.allclose(e._weights, custom_q[1])
 
-    # Determine unique quadrature degree, quadrature scheme and
-    # precision per each integral data
+    # Determine unique quadrature degree and quadrature scheme
+    # per each integral data
     for id, integral_data in enumerate(form_data.integral_data):
         # Iterate through groups of integral data. There is one integral
         # data for all integrals with same domain, itype, subdomain_id
@@ -192,22 +192,6 @@ def _analyze_form(form: ufl.form.Form, options: typing.Dict) -> ufl.algorithms.f
         # all integrals in this integral data group, i.e. must be the
         # same for for the same (domain, itype, subdomain_id)
 
-        # Extract precision
-        p_default = -1
-        precisions = set([integral.metadata().get("precision", p_default)
-                          for integral in integral_data.integrals])
-        precisions.discard(p_default)
-
-        if len(precisions) == 1:
-            p = precisions.pop()
-        elif len(precisions) == 0:
-            # Default precision
-            p = None
-        else:
-            raise RuntimeError("Only one precision allowed within integrals grouped by subdomain.")
-
-        integral_data.metadata["precision"] = p
-
         qd_default = -1
         qr_default = "default"
 
@@ -228,12 +212,11 @@ def _analyze_form(form: ufl.form.Form, options: typing.Dict) -> ufl.algorithms.f
                 logger.info(f"Integral {i}, integral group {id}:")
                 logger.info(f"--- quadrature rule: {qr}")
                 logger.info(f"--- quadrature degree: {qd}")
-                logger.info(f"--- precision: {p}")
 
-                metadata.update({"quadrature_degree": qd, "quadrature_rule": qr, "precision": p})
+                metadata.update({"quadrature_degree": qd, "quadrature_rule": qr})
             else:
                 metadata.update({"quadrature_points": custom_q[0], "quadrature_weights": custom_q[1],
-                                 "quadrature_rule": "custom", "precision": p})
+                                 "quadrature_rule": "custom"})
 
             integral_data.integrals[i] = integral.reconstruct(metadata=metadata)
 

ffcx/codegeneration/C/c_implementation.py

@@ -6,8 +6,7 @@
 
 import warnings
 import ffcx.codegeneration.lnodes as L
-from ffcx.codegeneration.utils import scalar_to_value_type, cdtype_to_numpy
-import numpy as np
+from ffcx.codegeneration.utils import scalar_to_value_type
 
 math_table = {
     "double": {
@@ -139,22 +138,16 @@
 
 
 class CFormatter(object):
-    def __init__(self, scalar, precision=None) -> None:
+    def __init__(self, scalar) -> None:
         self.scalar_type = scalar
         self.real_type = scalar_to_value_type(scalar)
-        if precision is None:
-            np_type = cdtype_to_numpy(self.real_type)
-            self.precision = np.finfo(np_type).precision + 1
-        else:
-            assert isinstance(precision, int)
-            self.precision = precision
 
     def _format_number(self, x):
-        p = self.precision
+        # Use 16sf for precision (good for float64 or less)
         if isinstance(x, complex):
-            return f"({x.real:.{p}}+I*{x.imag:.{p}})"
+            return f"({x.real:.16}+I*{x.imag:.16})"
         elif isinstance(x, float):
-            return f"{x:.{p}}"
+            return f"{x:.16}"
         return str(x)
 
     def _build_initializer_lists(self, values):

ffcx/codegeneration/C/integrals.py

@@ -36,7 +36,7 @@ def generator(ir, options):
     parts = ig.generate()
 
     # Format code as string
-    CF = CFormatter(options["scalar_type"], ir.precision)
+    CF = CFormatter(options["scalar_type"])
     body = CF.c_format(parts)
 
     # Generate generic FFCx code snippets and add specific parts

ffcx/ir/representation.py

@@ -134,7 +134,6 @@ class IntegralIR(typing.NamedTuple):
     unique_table_types: typing.Dict[str, str]
     integrand: typing.Dict[QuadratureRule, dict]
     name: str
-    precision: int
     needs_facet_permutations: bool
     coordinate_element: str
 
@@ -485,7 +484,6 @@ def _compute_integral_ir(form_data, form_index, element_numbers, integral_names,
             _offset += np.prod(constant.ufl_shape, dtype=int)
 
         ir["original_constant_offsets"] = original_constant_offsets
-        ir["precision"] = itg_data.metadata["precision"]
 
         # Create map from number of quadrature points -> integrand
         integrands = {rule: integral.integrand() for rule, integral in sorted_integrals.items()}