propagate dtype through finat element construction - #263
Conversation
hardik-corintis
changed the title
53f7b0e to
80c771b
Compare
|
How are you constructing the reference element in firedrake? |
How can this be? Aren't you setting the precision when you construct the cell? |
|
Yes, just propagate the dtype at cell construction |
|
Arguably the cell should be a global constant, we should not be creating multiple instances in firedrake. Reducing the number of cell construction calls is a separate PR |
Ok, so if you agree, I can revert the changes and repurpose this PR to pass the Lines 113 to 120 in f16c2df |
Sure that sounds fine to me |
hardik-corintis
changed the title
Threads dtype=RealType (or the numpy-derived real dtype from scalar_type in TSFC) through every as_fiat_cell/create_element call site in firedrake and tsfc, using FIAT PR firedrakeproject#260 and FInAT's new dtype support (firedrakeproject/fiat#263). Removes the finat.element_factory.ufc_cell monkeypatch in firedrake/utils.py, which is no longer needed now that callers declare their working precision explicitly instead of relying on a globally-overridden default. Also fixes petsc_sparse's hardcoded rtol=1E-10 drop tolerance (firedrake/preconditioners/fdm.py), which is too tight for single precision: reference-cell geometry feeding into these matrices is only accurate to float32 round-off (~1e-7), so entries that should be exactly zero (e.g. in a Nedelec discrete gradient matrix) were leaking into the sparsity pattern above the fp64-tuned threshold, corrupting the exact combinatorial structure PETSc's PCBDDCNedelecSupport expects and causing test_bddc_aij_simplex[N1curl-3-False] to fail.
c70b2ef to
d9559f1
Compare
…tions Rank- and nullspace-determining linear algebra during macro-element construction (compute_normal, spanning_basis, AlfeldSorokinaSpace, hdiv_conforming_coefficients) needs float64 regardless of the working precision, since float32 round-off in vertex coordinates can flip the computed rank.
d9559f1 to
6ec3724
Compare
|
|
||
| if len(rows) > 0: | ||
| dual_mat = numpy.vstack(rows) | ||
| interior_facets = ref_complex.get_interior_facets(sd-1) |
There was a problem hiding this comment.
This patching is unacceptable. This means that storing the vertices in single precision was the wrong decision
| rows_fp64.append(numpy.tensordot(weights_fp64, njump_fp64, axes=(ax, ax))) | ||
|
|
||
| dual_mat = numpy.vstack(rows_fp64) | ||
| nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) |
There was a problem hiding this comment.
Why can't we just set the tolerance based on the working precision here?
| def _reconstruct_split_complex_fp64( | ||
| ref_el: SimplicialComplex) -> SimplicialComplex: | ||
| """Reconstruct a split complex from float64 root geometry. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| ref_el : SimplicialComplex | ||
| The reference complex to reconstruct. | ||
|
|
||
| Returns | ||
| ------- | ||
| SimplicialComplex | ||
| The reconstructed reference complex. | ||
| """ | ||
| if not isinstance(ref_el, SplitSimplicialComplex): | ||
| ref_el_fp64 = copy.copy(ref_el) | ||
| ref_el_fp64.vertices = reference_element.cast_vertices( | ||
| ref_el.vertices, numpy.float64) | ||
| ref_el_fp64._split_cache = {} | ||
| return ref_el_fp64 | ||
|
|
||
| parent_fp64 = _reconstruct_split_complex_fp64(ref_el._split_parent) | ||
| return ref_el.reconstruct(parent_fp64) |
There was a problem hiding this comment.
We need to find a better solution that does not involve going back and forth between precisions.
More fundamentally, we need to decide whether Firedrake float32 even needs a FIAT reimplementation. For most elements we can get away by casting the numerical results at the very end. The main issue was that symbolic FIAT tabulations need to be casted differently.
Tell Claude to find the path of least resistance (whichever is easier to maintain): to let FIAT compute everything for a single dtype or to compute everything in double and cast the final result to the target dtype.
2 participants
follow-up to #260
as_fiat_cell/convert_finiteelement/create_elementinfinat/element_factory.pynow take adtypeand forward it toufc_cell, soref_el.verticesis built at the caller's actual precision.firedrake side: firedrakeproject/firedrake#5033