diff --git a/.github/workflows/build-mac.yml b/.github/workflows/build-mac.yml index 51aa9f5104..f8530311df 100644 --- a/.github/workflows/build-mac.yml +++ b/.github/workflows/build-mac.yml @@ -50,6 +50,9 @@ jobs: run: | . ../firedrake_venv/bin/activate python -m pytest -v tests/firedrake/regression/ -k "poisson_strong or stokes_mini or dg_advection" + # also test for 'problem libraries' (spatialindex and libsupermesh) + python -m pytest -v tests/firedrake/regression/test_locate_cell.py + python -m pytest -v tests/firedrake/supermesh/test_assemble_mixed_mass_matrix.py timeout-minutes: 30 - name: Post-run cleanup if: ${{ always() }} diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 909db6990a..0eb616c24d 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -66,6 +66,7 @@ jobs: --mpicxx="$MPICH_DIR"/mpicxx \ --mpif90="$MPICH_DIR"/mpif90 \ --mpiexec="$MPICH_DIR"/mpiexec \ + --mpihome="$MPICH_DIR"/.. \ --venv-name firedrake_venv \ --no-package-manager \ --disable-ssh \ diff --git a/.github/workflows/pip-mac.yml b/.github/workflows/pip-mac.yml index c8d1fbfe2e..54e2ada7cc 100644 --- a/.github/workflows/pip-mac.yml +++ b/.github/workflows/pip-mac.yml @@ -68,25 +68,6 @@ jobs: --download-superlu_dist make - - name: Install libsupermesh - run: | - source pip_venv/bin/activate - python -m pip install 'rtree>=1.2' - cd pip_venv/src - git clone https://github.com/firedrakeproject/libsupermesh.git - mkdir -p libsupermesh/build - cd libsupermesh/build - cmake .. \ - -DBUILD_SHARED_LIBS=ON \ - -DCMAKE_INSTALL_PREFIX="$VIRTUAL_ENV" \ - -DMPI_C_COMPILER=/opt/homebrew/bin/mpicc \ - -DMPI_CXX_COMPILER=/opt/homebrew/bin/mpicxx \ - -DMPI_Fortran_COMPILER=/opt/homebrew/bin/mpif90 \ - -DCMAKE_Fortran_COMPILER=/opt/homebrew/bin/mpif90 \ - -DMPIEXEC_EXECUTABLE=/opt/homebrew/bin/mpiexec - make - make install - - uses: actions/checkout@v4 with: path: pip_venv/src/firedrake @@ -118,6 +99,9 @@ jobs: cd pip_venv/src/firedrake python -m pytest --timeout=1800 -v tests/firedrake/regression \ -k "poisson_strong or stokes_mini or dg_advection" + # also test for 'problem libraries' (spatialindex and libsupermesh) + python -m pytest -v tests/firedrake/regression/test_locate_cell.py + python -m pytest -v tests/firedrake/supermesh/test_assemble_mixed_mass_matrix.py timeout-minutes: 30 - name: Cleanup (post) diff --git a/.github/workflows/pip.yml b/.github/workflows/pip.yml index c358b70141..868c71cf9c 100644 --- a/.github/workflows/pip.yml +++ b/.github/workflows/pip.yml @@ -56,25 +56,6 @@ jobs: with: path: src/firedrake - - name: Install libsupermesh - run: | - source pip_venv/bin/activate - python -m pip install 'rtree>=1.2' - cd pip_venv/src - git clone https://github.com/firedrakeproject/libsupermesh.git - mkdir -p libsupermesh/build - cd libsupermesh/build - cmake .. \ - -DBUILD_SHARED_LIBS=ON \ - -DCMAKE_INSTALL_PREFIX="$VIRTUAL_ENV" \ - -DMPI_C_COMPILER="$MPICH_DIR/mpicc" \ - -DMPI_CXX_COMPILER="$MPICH_DIR/mpicxx" \ - -DMPI_Fortran_COMPILER="$MPICH_DIR/mpif90" \ - -DCMAKE_Fortran_COMPILER="$MPICH_DIR/mpif90" \ - -DMPIEXEC_EXECUTABLE="$MPICH_DIR/mpiexec" - make - make install - - name: Pip install run: | source pip_venv/bin/activate diff --git a/.github/workflows/pyop2.yml b/.github/workflows/pyop2.yml index b5e750242a..5104284d05 100644 --- a/.github/workflows/pyop2.yml +++ b/.github/workflows/pyop2.yml @@ -35,17 +35,6 @@ jobs: id: setup-python with: python-version: ${{ matrix.python-version }} - # By default setup-python pollutes the environment in such a way that virtual - # environments cannot be used. This prevents us from building libsupermesh because - # it relies on having rtree installed into a venv. - # https://github.com/actions/setup-python/issues/851 - # https://github.com/actions/setup-python/blob/main/docs/advanced-usage.md#using-update-environment-flag - update-environment: false - - - name: Create virtual environment - shell: bash - run: | - ${{ steps.setup-python.outputs.python-path }} -m venv venv - name: Clone PETSc uses: actions/checkout@v4 @@ -65,25 +54,6 @@ jobs: --with-fortran-bindings=0 make - - name: Install libsupermesh - shell: bash - run: | - source venv/bin/activate - python -m pip install 'rtree>=1.2' - git clone https://github.com/firedrakeproject/libsupermesh.git - mkdir -p libsupermesh/build - cd libsupermesh/build - cmake .. \ - -DBUILD_SHARED_LIBS=ON \ - -DCMAKE_INSTALL_PREFIX="$VIRTUAL_ENV" \ - -DMPI_C_COMPILER=mpicc \ - -DMPI_CXX_COMPILER=mpicxx \ - -DMPI_Fortran_COMPILER=mpif90 \ - -DCMAKE_Fortran_COMPILER=mpif90 \ - -DMPIEXEC_EXECUTABLE=mpiexec - make - make install - - name: Checkout Firedrake uses: actions/checkout@v4 with: @@ -93,7 +63,6 @@ jobs: shell: bash working-directory: firedrake run: | - source ../venv/bin/activate python -m pip install -U pip python -m pip install -U pytest-timeout @@ -101,7 +70,6 @@ jobs: shell: bash working-directory: firedrake run: | - source ../venv/bin/activate export CC=mpicc export HDF5_DIR="$PETSC_DIR/$PETSC_ARCH" export HDF5_MPI=ON @@ -111,7 +79,6 @@ jobs: shell: bash working-directory: firedrake run: | - source ../venv/bin/activate pytest --tb=native --timeout=480 --timeout-method=thread -o faulthandler_timeout=540 -v tests/tsfc timeout-minutes: 10 @@ -119,7 +86,6 @@ jobs: shell: bash working-directory: firedrake run: | - source ../venv/bin/activate pytest -m "not parallel" --tb=native --timeout=480 --timeout-method=thread -o faulthandler_timeout=540 -v tests/pyop2 mpiexec -n 2 --oversubscribe pytest -m "parallel[2]" --tb=native --timeout=480 --timeout-method=thread -o faulthandler_timeout=540 -v tests/pyop2 mpiexec -n 3 --oversubscribe pytest -m "parallel[3]" --tb=native --timeout=480 --timeout-method=thread -o faulthandler_timeout=540 -v tests/pyop2 diff --git a/docker/Dockerfile.complex b/docker/Dockerfile.complex index 1528f5c4ac..6aa73f280c 100644 --- a/docker/Dockerfile.complex +++ b/docker/Dockerfile.complex @@ -18,5 +18,6 @@ RUN bash -c "python3 firedrake-install \ --mpicxx=$MPICH_DIR/mpicxx \ --mpif90=$MPICH_DIR/mpif90 \ --mpiexec=$MPICH_DIR/mpiexec \ + --mpihome=$MPICH_DIR/.. \ --slepc \ --documentation-dependencies" diff --git a/docker/Dockerfile.vanilla b/docker/Dockerfile.vanilla index 2de26b1916..eee3ef53e3 100644 --- a/docker/Dockerfile.vanilla +++ b/docker/Dockerfile.vanilla @@ -16,4 +16,5 @@ RUN bash -c "python3 firedrake-install \ --mpicc=$MPICH_DIR/mpicc \ --mpicxx=$MPICH_DIR/mpicxx \ --mpif90=$MPICH_DIR/mpif90 \ - --mpiexec=$MPICH_DIR/mpiexec" + --mpiexec=$MPICH_DIR/mpiexec \ + --mpihome=$MPICH_DIR/.." diff --git a/docs/source/download.rst b/docs/source/download.rst index 74645fc682..5f8f478864 100644 --- a/docs/source/download.rst +++ b/docs/source/download.rst @@ -241,7 +241,6 @@ Requirements * An activated virtual environment. * All the system requirements listed in :ref:`system-requirements`. * A Firedrake-compatible PETSc installation (using our `fork of PETSc `_). The set of flags passed to PETSc can be retrieved by passing the command ``--show-petsc-configure-options`` to ``firedrake-install``. -* `libsupermesh `_ to be installed inside the virtual environment (see `here `_ for an example of how to do this). * The following environment variables to be set: * ``PETSC_DIR`` and ``PETSC_ARCH`` to point to the correct location for the PETSc installation. diff --git a/docs/source/element_list.py b/docs/source/element_list.py index 2c2d048f01..6e8c6ab981 100644 --- a/docs/source/element_list.py +++ b/docs/source/element_list.py @@ -1,6 +1,6 @@ from finat.ufl.elementlist import ufl_elements # ~ from ufl.finiteelement.elementlist import ufl_elements -from tsfc.finatinterface import supported_elements +from finat.element_factory import supported_elements import csv shape_names = { diff --git a/docs/source/firedrake_usa_25.rst b/docs/source/firedrake_usa_25.rst index e1d7a2b293..8cac5f4a9c 100644 --- a/docs/source/firedrake_usa_25.rst +++ b/docs/source/firedrake_usa_25.rst @@ -27,13 +27,13 @@ The conference will begin with a tutorial session on the morning of 28 February Conference venue ---------------- -The conference will take place in the Bill Daniel Student Center in Room 202 in the heart of the Baylor campus. +The conference will take place in the Bill Daniel Student Center in Room 202 in the heart of the Baylor campus. A `campus map `__ is available online, and both Apple Maps and Google Maps accurately locate the Bill Daniel Student Center. Accommodation ------------- -We will be reserving a room block in a hotel near campus. More information to follow. +We have reserved a block of hotel rooms at SpringHill Suites Waco, about a 15 minute walk from the Bill Daniel Student Center. Follow `this link `__ to hold your room. The group rate is only available until 20 January 2025. Conference dinner @@ -53,9 +53,9 @@ The registration fees are as follows: :widths: 25 50 :header-rows: 0 - * - Student + * - Student: - $50 - * - Non-student + * - Non-student: - $200 The `SIAM Texas-Louisiana Section `__ is providing some support for students currently attending universities in Texas or Louisiana to attend. @@ -66,7 +66,7 @@ Conference registration is coming soon. Abstract submission ------------------- -Abstract submission will open soon via Easy Chair. +Abstracts can be submitted `via EasyChair `__. @@ -79,17 +79,11 @@ The conference has been kindly supported by the SIAM TX-LA Section and EPSRC. Travel to Waco -------------- -* By air +* By air: Waco has a small airport. There is daily service between Waco and Dallas/Fort Worth International Airport via American Airlines. - Waco has a small airport. There is daily service between Waco and Dallas/Fort Worth International Airport via American Airlines. +* By ground: We are less than two hours by car from the Dallas and Austin airports, and just under three hours from Bush Intercontinental Airport in Houston. Additionally, Waco is reachable by bus services such as Greyhound and FlixBus. -* By ground - - We are less than two hours by car from the Dallas and Austin airports, and just under three hours from Bush Intercontinental Airport in Houston. Additionally, Waco is reachable by bus services such as Greyhound and FlixBus. - -* Parking on campus - - Baylor has plenty of visitor parking for your personal or rental vehicle, but conslut `these instructions __` and make sure to `register your vehicle __`. +* Parking on campus: Baylor has plenty of visitor parking for your personal or rental vehicle, but consult `these instructions `__ and make sure to `register your vehicle `__. diff --git a/docs/source/parallelism.rst b/docs/source/parallelism.rst index 70fc119e98..d71b3d9305 100644 --- a/docs/source/parallelism.rst +++ b/docs/source/parallelism.rst @@ -45,11 +45,12 @@ firedrake installer to use it, by running: .. code-block:: shell - python3 firedrake-install --mpiexec=mpiexec --mpicc=mpicc --mpicxx=mpicxx --mpif90=mpif90 + python3 firedrake-install --mpiexec=mpiexec --mpicc=mpicc --mpicxx=mpicxx --mpif90=mpif90 --mpihome mpihome where ``mpiexec``, ``mpicc``, ``mpicxx``, and ``mpif90`` are the commands to run an MPI job and to compile C, C++, and Fortran 90 code, -respectively. +respectively. ``mpihome`` is an extra variable that must point to the +root directory of the MPI installation (e.g. ``/usr`` or ``/opt/mpich``). Printing in parallel ==================== diff --git a/firedrake/assemble.py b/firedrake/assemble.py index f451b3f596..f3049ae01c 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -12,7 +12,7 @@ import numpy from pyadjoint.tape import annotate_tape from tsfc import kernel_args -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from tsfc.ufl_utils import extract_firedrake_constants import ufl import finat.ufl diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 8b8e1da170..b5e0777dfe 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -11,7 +11,7 @@ from mpi4py import MPI from firedrake.utils import IntType, ScalarType from libc.string cimport memset from libc.stdlib cimport qsort -from tsfc.finatinterface import as_fiat_cell +from finat.element_factory import as_fiat_cell cimport numpy as np cimport mpi4py.MPI as MPI diff --git a/firedrake/cython/extrusion_numbering.pyx b/firedrake/cython/extrusion_numbering.pyx index 278af86c75..86167ed20b 100644 --- a/firedrake/cython/extrusion_numbering.pyx +++ b/firedrake/cython/extrusion_numbering.pyx @@ -193,7 +193,7 @@ from mpi4py.libmpi cimport (MPI_Op_create, MPI_OP_NULL, MPI_Op_free, MPI_User_function) from pyop2 import op2 from firedrake.utils import IntType -from tsfc.finatinterface import as_fiat_cell +from finat.element_factory import as_fiat_cell cimport numpy cimport mpi4py.MPI as MPI diff --git a/firedrake/extrusion_utils.py b/firedrake/extrusion_utils.py index b038d904af..4298a21d7b 100644 --- a/firedrake/extrusion_utils.py +++ b/firedrake/extrusion_utils.py @@ -8,7 +8,7 @@ from pyop2.caching import serial_cache from firedrake.petsc import PETSc from firedrake.utils import IntType, RealType, ScalarType -from tsfc.finatinterface import create_element +from finat.element_factory import create_element import loopy as lp from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_2 # noqa: F401 from firedrake.parameters import target diff --git a/firedrake/formmanipulation.py b/firedrake/formmanipulation.py index 35a6789107..c3a3166ae8 100644 --- a/firedrake/formmanipulation.py +++ b/firedrake/formmanipulation.py @@ -2,13 +2,17 @@ import numpy import collections -from ufl import as_vector -from ufl.classes import Zero, FixedIndex, ListTensor +from ufl import as_vector, FormSum, Form, split +from ufl.classes import Zero, FixedIndex, ListTensor, ZeroBaseForm from ufl.algorithms.map_integrands import map_integrand_dags from ufl.corealg.map_dag import MultiFunction, map_expr_dags +from pyop2 import MixedDat + from firedrake.petsc import PETSc from firedrake.ufl_expr import Argument +from firedrake.cofunction import Cofunction +from firedrake.functionspace import FunctionSpace, MixedFunctionSpace, DualSpace class ExtractSubBlock(MultiFunction): @@ -85,9 +89,8 @@ def coefficient_derivative(self, o, expr, coefficients, arguments, cds): @PETSc.Log.EventDecorator() def argument(self, o): - from ufl import split - from firedrake import MixedFunctionSpace, FunctionSpace V = o.function_space() + if len(V) == 1: # Not on a mixed space, just return ourselves. return o @@ -98,15 +101,11 @@ def argument(self, o): V_is = V.subfunctions indices = self.blocks[o.number()] - try: - indices = tuple(indices) - nidx = len(indices) - except TypeError: - # Only one index provided. + # Only one index provided. + if isinstance(indices, int): indices = (indices, ) - nidx = 1 - if nidx == 1: + if len(indices) == 1: W = V_is[indices[0]] W = FunctionSpace(W.mesh(), W.ufl_element()) a = (Argument(W, o.number(), part=o.part()), ) @@ -127,6 +126,50 @@ def argument(self, o): for j in numpy.ndindex(V_is[i].value_shape)] return self._arg_cache.setdefault(o, as_vector(args)) + def cofunction(self, o): + V = o.function_space() + + # Not on a mixed space, just return ourselves. + if len(V) == 1: + return o + + # We only need the test space for Cofunction + indices = self.blocks[0] + V_is = V.subfunctions + + # Only one index provided. + if isinstance(indices, int): + indices = (indices, ) + + # for two-forms, the cofunction should only + # be returned for the diagonal blocks, so + # if we are asked for an off-diagonal block + # then we return a zero form, analogously to + # the off components of arguments. + if len(self.blocks) == 2: + itest, itrial = self.blocks + on_diag = (itest == itrial) + else: + on_diag = True + + # if we are on the diagonal, then return a Cofunction + # in the relevant subspace that points to the data in + # the full space. This means that the right hand side + # of the fieldsplit problem will be correct. + if on_diag: + if len(indices) == 1: + i = indices[0] + W = V_is[i] + W = DualSpace(W.mesh(), W.ufl_element()) + c = Cofunction(W, val=o.subfunctions[i].dat) + else: + W = MixedFunctionSpace([V_is[i] for i in indices]) + c = Cofunction(W, val=MixedDat(o.dat[i] for i in indices)) + else: + c = ZeroBaseForm(o.arguments()) + + return c + SplitForm = collections.namedtuple("SplitForm", ["indices", "form"]) @@ -168,7 +211,20 @@ def split_form(form, diagonal=False): assert len(shape) == 2 for idx in numpy.ndindex(shape): f = splitter.split(form, idx) - if len(f.integrals()) > 0: + + # does f actually contain anything? + if isinstance(f, Cofunction): + flen = 1 + elif isinstance(f, FormSum): + flen = len(f.components()) + elif isinstance(f, Form): + flen = len(f.integrals()) + else: + raise ValueError( + "ExtractSubBlock.split should have returned an instance of " + "either Form, FormSum, or Cofunction") + + if flen > 0: if diagonal: i, j = idx if i != j: diff --git a/firedrake/functionspacedata.py b/firedrake/functionspacedata.py index be00e3ffff..a1b6190cfb 100644 --- a/firedrake/functionspacedata.py +++ b/firedrake/functionspacedata.py @@ -20,7 +20,7 @@ from decorator import decorator from functools import partial -from tsfc.finatinterface import create_element as _create_element +from finat.element_factory import create_element as _create_element from pyop2 import op2 from firedrake.utils import IntType diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index e8555ba7dd..f27024ff28 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -14,7 +14,7 @@ from pyop2 import op2 from pyop2.caching import memory_and_disk_cache -from tsfc.finatinterface import create_element, as_fiat_cell +from finat.element_factory import create_element, as_fiat_cell from tsfc import compile_expression_dual_evaluation from tsfc.ufl_utils import extract_firedrake_constants @@ -1092,8 +1092,7 @@ def _interpolator(V, tensor, expr, subset, arguments, access, bcs=None): # interpolation) we have to pass the finat element we construct # here. Ideally we would only pass the UFL element through. kernel = compile_expression(cell_set.comm, expr, to_element, V.ufl_element(), - domain=source_mesh, parameters=parameters, - log=PETSc.Log.isActive()) + domain=source_mesh, parameters=parameters) ast = kernel.ast oriented = kernel.oriented needs_cell_sizes = kernel.needs_cell_sizes @@ -1221,10 +1220,9 @@ def _interpolator(V, tensor, expr, subset, arguments, access, bcs=None): f"firedrake-tsfc-expression-kernel-cache-uid{os.getuid()}") -def _compile_expression_key(comm, expr, to_element, ufl_element, domain, parameters, log): +def _compile_expression_key(comm, expr, to_element, ufl_element, domain, parameters): """Generate a cache key suitable for :func:`tsfc.compile_expression_dual_evaluation`.""" - key = hash_expr(expr), hash(ufl_element), utils.tuplify(parameters), log - return key + return (hash_expr(expr), hash(ufl_element), utils.tuplify(parameters)) @memory_and_disk_cache( diff --git a/firedrake/mesh.py b/firedrake/mesh.py index bef4d38bf3..4e209340c7 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -45,7 +45,7 @@ ngsPETSc = None # Only for docstring import mpi4py # noqa: F401 -from tsfc.finatinterface import as_fiat_cell +from finat.element_factory import as_fiat_cell __all__ = [ diff --git a/firedrake/mg/kernels.py b/firedrake/mg/kernels.py index 5405b6d726..f892f6260c 100644 --- a/firedrake/mg/kernels.py +++ b/firedrake/mg/kernels.py @@ -30,7 +30,7 @@ from tsfc.driver import TSFCIntegralDataInfo from tsfc.kernel_interface.common import lower_integral_type from tsfc.parameters import default_parameters -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from finat.quadrature import make_quadrature from firedrake.pointquery_utils import dX_norm_square, X_isub_dX, init_X, inside_check, is_affine, celldist_l1_c_expr from firedrake.pointquery_utils import to_reference_coords_newton_step as to_reference_coords_newton_step_body @@ -45,7 +45,7 @@ def to_reference_coordinates(ufl_coordinate_element, parameters=None): parameters = _ # Create FInAT element - element = tsfc.finatinterface.create_element(ufl_coordinate_element) + element = finat.element_factory.create_element(ufl_coordinate_element) gdim, = ufl_coordinate_element.reference_value_shape cell = ufl_coordinate_element.cell diff --git a/firedrake/output/paraview_reordering.py b/firedrake/output/paraview_reordering.py index 20fbe50099..8b6edb147d 100644 --- a/firedrake/output/paraview_reordering.py +++ b/firedrake/output/paraview_reordering.py @@ -1,4 +1,4 @@ -from tsfc.finatinterface import create_base_element +from finat.element_factory import create_base_element import numpy as np from pyop2.utils import as_tuple diff --git a/firedrake/parloops.py b/firedrake/parloops.py index 08cb57e3d8..0a33cd4ae5 100644 --- a/firedrake/parloops.py +++ b/firedrake/parloops.py @@ -171,7 +171,7 @@ def par_loop(kernel, measure, args, kernel_kwargs=None, **kwargs): domain = '{[i]: 0 <= i < A.dofs}' instructions = ''' for i - A[i] = max(A[i], B[0]) + A[i] = fmax(A[i], B[0]) end ''' par_loop((domain, instructions), dx, {'A' : (A, RW), 'B': (B, READ)}) diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 4793af0bfe..2308616622 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -221,7 +221,7 @@ def compile_coordinate_element(mesh: MeshGeometry, contains_eps: float, paramete ufl_coordinate_element = mesh.ufl_coordinate_element() # Create FInAT element - element = tsfc.finatinterface.create_element(ufl_coordinate_element) + element = finat.element_factory.create_element(ufl_coordinate_element) code = { "geometric_dimension": mesh.geometric_dimension(), diff --git a/firedrake/preconditioners/fdm.py b/firedrake/preconditioners/fdm.py index 1696801cb4..e75172dc8e 100644 --- a/firedrake/preconditioners/fdm.py +++ b/firedrake/preconditioners/fdm.py @@ -19,7 +19,7 @@ from firedrake_citations import Citations from ufl.algorithms.ad import expand_derivatives from ufl.algorithms.expand_indices import expand_indices -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from pyop2.compilation import load from pyop2.mpi import COMM_SELF from pyop2.sparsity import get_preallocation diff --git a/firedrake/preconditioners/pmg.py b/firedrake/preconditioners/pmg.py index 251f585cbe..48df37d615 100644 --- a/firedrake/preconditioners/pmg.py +++ b/firedrake/preconditioners/pmg.py @@ -9,7 +9,7 @@ from firedrake.solving_utils import _SNESContext from firedrake.tsfc_interface import extract_numbered_coefficients from firedrake.utils import ScalarType_c, IntType_c, cached_property -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from tsfc import compile_expression_dual_evaluation from pyop2 import op2 from pyop2.caching import serial_cache @@ -530,7 +530,7 @@ def coarsen_bc_value(self, bc, cV): def prolongation_transfer_kernel_action(Vf, expr): to_element = create_element(Vf.ufl_element()) - kernel = compile_expression_dual_evaluation(expr, to_element, Vf.ufl_element(), log=PETSc.Log.isActive()) + kernel = compile_expression_dual_evaluation(expr, to_element, Vf.ufl_element()) coefficients = extract_numbered_coefficients(expr, kernel.coefficient_numbers) if kernel.needs_external_coords: coefficients = [Vf.mesh().coordinates] + coefficients diff --git a/firedrake/slate/slac/compiler.py b/firedrake/slate/slac/compiler.py index 7e1b14281c..d9e7bced02 100644 --- a/firedrake/slate/slac/compiler.py +++ b/firedrake/slate/slac/compiler.py @@ -237,4 +237,4 @@ def gem_to_loopy(gem_expr, var2terminal, scalar_type): # Part B: impero_c to loopy output_arg = OutputKernelArg(output_loopy_arg) - return generate_loopy(impero_c, args, scalar_type, "slate_loopy", [], log=PETSc.Log.isActive()), output_arg + return generate_loopy(impero_c, args, scalar_type, "slate_loopy", []), output_arg diff --git a/firedrake/slate/slac/kernel_builder.py b/firedrake/slate/slac/kernel_builder.py index cbc9b6fed2..419931232f 100644 --- a/firedrake/slate/slac/kernel_builder.py +++ b/firedrake/slate/slac/kernel_builder.py @@ -14,7 +14,7 @@ from firedrake.slate.slac.tsfc_driver import compile_terminal_form from tsfc import kernel_args -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from tsfc.loopy import create_domains, assign_dtypes from pytools import UniqueNameGenerator diff --git a/firedrake/supermeshing.py b/firedrake/supermeshing.py index a1ce2cde17..ee576ea6e5 100644 --- a/firedrake/supermeshing.py +++ b/firedrake/supermeshing.py @@ -1,7 +1,8 @@ # Code for projections and other fun stuff involving supermeshes. import firedrake import ctypes -import sys +import pathlib +import libsupermesh from firedrake.cython.supermeshimpl import assemble_mixed_mass_matrix as ammm, intersection_finder from firedrake.mg.utils import get_level from firedrake.petsc import PETSc @@ -428,7 +429,8 @@ def likely(cell_A): "complex_mode": 1 if complex_mode else 0 } - dirs = get_petsc_dir() + (sys.prefix, ) + libsupermesh_dir = pathlib.Path(libsupermesh.get_include()).parent.absolute() + dirs = get_petsc_dir() + (libsupermesh_dir,) includes = ["-I%s/include" % d for d in dirs] libs = ["-L%s/lib" % d for d in dirs] libs = libs + ["-Wl,-rpath,%s/lib" % d for d in dirs] + ["-lpetsc", "-lsupermesh"] diff --git a/firedrake/tsfc_interface.py b/firedrake/tsfc_interface.py index a4a57ae0cb..ba10d79507 100644 --- a/firedrake/tsfc_interface.py +++ b/firedrake/tsfc_interface.py @@ -53,8 +53,8 @@ ) -def tsfc_compile_form_hashkey(form, prefix, parameters, interface, diagonal, log): - # Drop prefix as it's only used for naming and log +def tsfc_compile_form_hashkey(form, prefix, parameters, interface, diagonal): + # Drop prefix as it's only used for naming return default_parallel_hashkey(form.signature(), prefix, parameters, interface, diagonal) @@ -94,7 +94,7 @@ def __init__( """ tree = tsfc_compile_form(form, prefix=name, parameters=parameters, interface=interface, - diagonal=diagonal, log=PETSc.Log.isActive()) + diagonal=diagonal) kernels = [] for kernel in tree: # Individual kernels do not have to use all of the coefficients diff --git a/firedrake/utility_meshes.py b/firedrake/utility_meshes.py index 4f2a0f34f1..99d687449e 100644 --- a/firedrake/utility_meshes.py +++ b/firedrake/utility_meshes.py @@ -24,6 +24,7 @@ as_tensor, dot, And, + Or, sin, cos, real @@ -1800,6 +1801,8 @@ def PeriodicBoxMesh( Lx, Ly, Lz, + directions=(True, True, True), + hexahedral=False, reorder=None, distribution_parameters=None, comm=COMM_WORLD, @@ -1809,105 +1812,185 @@ def PeriodicBoxMesh( ): """Generate a periodic mesh of a 3D box. - :arg nx: The number of cells in the x direction - :arg ny: The number of cells in the y direction - :arg nz: The number of cells in the z direction - :arg Lx: The extent in the x direction - :arg Ly: The extent in the y direction - :arg Lz: The extent in the z direction - :kwarg reorder: (optional), should the mesh be reordered? - :kwarg distribution_parameters: options controlling mesh - distribution, see :func:`.Mesh` for details. - :kwarg comm: Optional communicator to build the mesh on. - :kwarg name: Optional name of the mesh. - :kwarg distribution_name: the name of parallel distribution used - when checkpointing; if `None`, the name is automatically - generated. - :kwarg permutation_name: the name of entity permutation (reordering) used - when checkpointing; if `None`, the name is automatically - generated. + Parameters + ---------- + nx : int + Number of cells in the x direction. + ny : int + Number of cells in the y direction. + nz : int + Number of cells in the z direction. + Lx : float + Extent in the x direction. + Ly : float + Extent in the y direction. + Lz : float + Extent in the z direction. + directions : list or tuple + Directions of periodicity. + hexahedral : bool + Whether to make hexahedral mesh or not. + reorder : bool or None + Whether to reorder the mesh. + distribution_parameters : dict or None + Options controlling mesh distribution, see :func:`.Mesh` for details. + comm : + Communicator to build the mesh on. + name : str + Name of the mesh. + distribution_name : str or None + Name of parallel distribution used when checkpointing; + if `None`, the name is automatically generated. + permutation_name : str or None + Name of entity permutation (reordering) used when checkpointing; + if `None`, the name is automatically generated. + + Returns + ------- + MeshGeometry + The mesh. + + Notes + ----- + + The boundary surfaces are numbered as follows: + + * 1: plane x == 0 + * 2: plane x == 1 + * 3: plane y == 0 + * 4: plane y == 1 + * 5: plane z == 0 + * 6: plane z == 1 + + where periodic surfaces are regarded as interior, for which dS integral is to be used. + """ for n in (nx, ny, nz): if n < 3: raise ValueError( "3D periodic meshes with fewer than 3 cells are not currently supported" ) + if hexahedral: + if len(directions) != 3: + raise ValueError(f"directions must have exactly dim (=3) elements : Got {directions}") + plex = PETSc.DMPlex().createBoxMesh( + (nx, ny, nz), + lower=(0., 0., 0.), + upper=(Lx, Ly, Lz), + simplex=False, + periodic=directions, + interpolate=True, + sparseLocalize=False, + comm=comm, + ) + m = mesh.Mesh( + plex, + reorder=reorder, + distribution_parameters=distribution_parameters, + name=name, + distribution_name=distribution_name, + permutation_name=permutation_name, + comm=comm) + x, y, z = SpatialCoordinate(m) + V = FunctionSpace(m, "Q", 2) + eps = min([Lx / nx, Ly / ny, Lz / nz]) / 1000. + if directions[0]: # x + fx0 = Function(V).interpolate(conditional(Or(x < eps, x > Lx - eps), 1., 0.)) + fx1 = fx0 + else: + fx0 = Function(V).interpolate(conditional(x < eps, 1., 0.)) + fx1 = Function(V).interpolate(conditional(x > Lx - eps, 1., 0.)) + if directions[1]: # y + fy0 = Function(V).interpolate(conditional(Or(y < eps, y > Ly - eps), 1., 0.)) + fy1 = fy0 + else: + fy0 = Function(V).interpolate(conditional(y < eps, 1., 0.)) + fy1 = Function(V).interpolate(conditional(y > Ly - eps, 1., 0.)) + if directions[2]: # z + fz0 = Function(V).interpolate(conditional(Or(z < eps, z > Lz - eps), 1., 0.)) + fz1 = fz0 + else: + fz0 = Function(V).interpolate(conditional(z < eps, 1., 0.)) + fz1 = Function(V).interpolate(conditional(z > Lz - eps, 1., 0.)) + return mesh.RelabeledMesh(m, [fx0, fx1, fy0, fy1, fz0, fz1], [1, 2, 3, 4, 5, 6], name=name) + else: + if tuple(directions) != (True, True, True): + raise NotImplementedError("Can only specify directions with hexahedral = True") + xcoords = np.arange(0.0, Lx, Lx / nx, dtype=np.double) + ycoords = np.arange(0.0, Ly, Ly / ny, dtype=np.double) + zcoords = np.arange(0.0, Lz, Lz / nz, dtype=np.double) + coords = ( + np.asarray(np.meshgrid(xcoords, ycoords, zcoords)).swapaxes(0, 3).reshape(-1, 3) + ) + i, j, k = np.meshgrid( + np.arange(nx, dtype=np.int32), + np.arange(ny, dtype=np.int32), + np.arange(nz, dtype=np.int32), + ) + v0 = k * nx * ny + j * nx + i + v1 = k * nx * ny + j * nx + (i + 1) % nx + v2 = k * nx * ny + ((j + 1) % ny) * nx + i + v3 = k * nx * ny + ((j + 1) % ny) * nx + (i + 1) % nx + v4 = ((k + 1) % nz) * nx * ny + j * nx + i + v5 = ((k + 1) % nz) * nx * ny + j * nx + (i + 1) % nx + v6 = ((k + 1) % nz) * nx * ny + ((j + 1) % ny) * nx + i + v7 = ((k + 1) % nz) * nx * ny + ((j + 1) % ny) * nx + (i + 1) % nx - xcoords = np.arange(0.0, Lx, Lx / nx, dtype=np.double) - ycoords = np.arange(0.0, Ly, Ly / ny, dtype=np.double) - zcoords = np.arange(0.0, Lz, Lz / nz, dtype=np.double) - coords = ( - np.asarray(np.meshgrid(xcoords, ycoords, zcoords)).swapaxes(0, 3).reshape(-1, 3) - ) - i, j, k = np.meshgrid( - np.arange(nx, dtype=np.int32), - np.arange(ny, dtype=np.int32), - np.arange(nz, dtype=np.int32), - ) - v0 = k * nx * ny + j * nx + i - v1 = k * nx * ny + j * nx + (i + 1) % nx - v2 = k * nx * ny + ((j + 1) % ny) * nx + i - v3 = k * nx * ny + ((j + 1) % ny) * nx + (i + 1) % nx - v4 = ((k + 1) % nz) * nx * ny + j * nx + i - v5 = ((k + 1) % nz) * nx * ny + j * nx + (i + 1) % nx - v6 = ((k + 1) % nz) * nx * ny + ((j + 1) % ny) * nx + i - v7 = ((k + 1) % nz) * nx * ny + ((j + 1) % ny) * nx + (i + 1) % nx - - cells = [ - [v0, v1, v3, v7], - [v0, v1, v7, v5], - [v0, v5, v7, v4], - [v0, v3, v2, v7], - [v0, v6, v4, v7], - [v0, v2, v6, v7], - ] - cells = np.asarray(cells).reshape(-1, ny, nx, nz).swapaxes(0, 3).reshape(-1, 4) - plex = mesh.plex_from_cell_list( - 3, cells, coords, comm, mesh._generate_default_mesh_topology_name(name) - ) - m = mesh.Mesh( - plex, - reorder=reorder_noop, - distribution_parameters=distribution_parameters_no_overlap, - name=name, - distribution_name=distribution_name, - permutation_name=permutation_name, - comm=comm, - ) + cells = [ + [v0, v1, v3, v7], + [v0, v1, v7, v5], + [v0, v5, v7, v4], + [v0, v3, v2, v7], + [v0, v6, v4, v7], + [v0, v2, v6, v7], + ] + cells = np.asarray(cells).reshape(-1, ny, nx, nz).swapaxes(0, 3).reshape(-1, 4) + plex = mesh.plex_from_cell_list( + 3, cells, coords, comm, mesh._generate_default_mesh_topology_name(name) + ) + m = mesh.Mesh( + plex, + reorder=reorder_noop, + distribution_parameters=distribution_parameters_no_overlap, + name=name, + distribution_name=distribution_name, + permutation_name=permutation_name, + comm=comm, + ) - new_coordinates = Function( - VectorFunctionSpace( - m, FiniteElement("DG", tetrahedron, 1, variant="equispaced") - ), - name=mesh._generate_default_mesh_coordinates_name(name), - ) - new_coordinates.interpolate(m.coordinates) + new_coordinates = Function( + VectorFunctionSpace( + m, FiniteElement("DG", tetrahedron, 1, variant="equispaced") + ), + name=mesh._generate_default_mesh_coordinates_name(name), + ) + new_coordinates.interpolate(m.coordinates) - coords_by_cell = new_coordinates.dat.data.reshape((-1, 4, 3)).transpose(1, 0, 2) + coords_by_cell = new_coordinates.dat.data.reshape((-1, 4, 3)).transpose(1, 0, 2) - # ensure we really got a view: - assert coords_by_cell.base is new_coordinates.dat.data.base + # ensure we really got a view: + assert coords_by_cell.base is new_coordinates.dat.data.base - # Find the cells that are too big in each direction because they are - # wrapped. - cell_is_wrapped = ( - (coords_by_cell.max(axis=0) - coords_by_cell.min(axis=0)) - / (Lx/nx, Ly/ny, Lz/nz) > 1.1 - ) + # Find the cells that are too big in each direction because they are + # wrapped. + cell_is_wrapped = ( + (coords_by_cell.max(axis=0) - coords_by_cell.min(axis=0)) + / (Lx/nx, Ly/ny, Lz/nz) > 1.1 + ) - # Move wrapped coordinates to the other end of the domain. - for i, extent in enumerate((Lx, Ly, Lz)): - coords = coords_by_cell[:, :, i] - coords[np.logical_and(cell_is_wrapped[:, i], np.isclose(coords, 0))] \ - = extent + # Move wrapped coordinates to the other end of the domain. + for i, extent in enumerate((Lx, Ly, Lz)): + coords = coords_by_cell[:, :, i] + coords[np.logical_and(cell_is_wrapped[:, i], np.isclose(coords, 0))] \ + = extent - return _postprocess_periodic_mesh(new_coordinates, - comm, - distribution_parameters, - reorder, - name, - distribution_name, - permutation_name) + return _postprocess_periodic_mesh(new_coordinates, + comm, + distribution_parameters, + reorder, + name, + distribution_name, + permutation_name) @PETSc.Log.EventDecorator() @@ -1915,6 +1998,8 @@ def PeriodicUnitCubeMesh( nx, ny, nz, + directions=(True, True, True), + hexahedral=False, reorder=None, distribution_parameters=None, comm=COMM_WORLD, @@ -1924,20 +2009,52 @@ def PeriodicUnitCubeMesh( ): """Generate a periodic mesh of a unit cube - :arg nx: The number of cells in the x direction - :arg ny: The number of cells in the y direction - :arg nz: The number of cells in the z direction - :kwarg reorder: (optional), should the mesh be reordered? - :kwarg distribution_parameters: options controlling mesh - distribution, see :func:`.Mesh` for details. - :kwarg comm: Optional communicator to build the mesh on. - :kwarg name: Optional name of the mesh. - :kwarg distribution_name: the name of parallel distribution used - when checkpointing; if `None`, the name is automatically - generated. - :kwarg permutation_name: the name of entity permutation (reordering) used - when checkpointing; if `None`, the name is automatically - generated. + Parameters + ---------- + nx : int + Number of cells in the x direction. + ny : int + Number of cells in the y direction. + nz : int + Number of cells in the z direction. + directions : list or tuple + Directions of periodicity. + hexahedral : bool + Whether to make hexahedral mesh or not. + reorder : bool or None + Should the mesh be reordered? + distribution_parameters : dict or None + Options controlling mesh distribution, see :func:`.Mesh` for details. + comm : + Communicator to build the mesh on. + name : str + Name of the mesh. + distribution_name : str or None + Name of parallel distribution used when checkpointing; + if `None`, the name is automatically generated. + permutation_name : str or None + Name of entity permutation (reordering) used when checkpointing; + if `None`, the name is automatically generated. + + Returns + ------- + MeshGeometry + The mesh. + + Notes + ----- + + The boundary surfaces are numbered as follows: + + * 1: plane x == 0 + * 2: plane x == 1 + * 3: plane y == 0 + * 4: plane y == 1 + * 5: plane z == 0 + * 6: plane z == 1 + + where periodic surfaces are regarded as interior, for which dS integral is to be used. + """ return PeriodicBoxMesh( nx, @@ -1946,6 +2063,8 @@ def PeriodicUnitCubeMesh( 1.0, 1.0, 1.0, + directions=directions, + hexahedral=hexahedral, reorder=reorder, distribution_parameters=distribution_parameters, comm=comm, diff --git a/pyproject.toml b/pyproject.toml index fe15d8e7c1..bb934e4182 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,7 @@ dependencies = [ "fenics-fiat @ git+https://github.com/firedrakeproject/fiat.git", "pyadjoint-ad @ git+https://github.com/dolfin-adjoint/pyadjoint.git", "loopy @ git+https://github.com/firedrakeproject/loopy.git@main", + "libsupermesh @ git+https://github.com/firedrakeproject/libsupermesh.git", ] classifiers = [ "Development Status :: 5 - Production/Stable", @@ -98,6 +99,7 @@ requires = [ "mpi4py; python_version < '3.13'", "petsc4py", "rtree>=1.2", + "libsupermesh @ git+https://github.com/firedrakeproject/libsupermesh.git", ] build-backend = "setuptools.build_meta" diff --git a/requirements-git.txt b/requirements-git.txt index 93e6df0f47..b616f15bbc 100644 --- a/requirements-git.txt +++ b/requirements-git.txt @@ -4,3 +4,4 @@ git+https://github.com/dolfin-adjoint/pyadjoint.git#egg=pyadjoint-ad git+https://github.com/firedrakeproject/loopy.git@main#egg=loopy git+https://github.com/firedrakeproject/pytest-mpi.git@main#egg=pytest-mpi git+https://github.com/firedrakeproject/petsc.git@firedrake#egg=petsc +git+https://github.com/firedrakeproject/libsupermesh.git#egg=libsupermesh diff --git a/scripts/firedrake-install b/scripts/firedrake-install index f4406a6ffb..17a8d7fdbc 100755 --- a/scripts/firedrake-install +++ b/scripts/firedrake-install @@ -16,10 +16,21 @@ from glob import iglob from itertools import chain import re import importlib + + +class InstallError(Exception): + # Exception for generic install problems. + pass + + try: from pkg_resources.extern.packaging.version import Version, InvalidVersion except ModuleNotFoundError: - from packaging.version import Version, InvalidVersion + try: + from packaging.version import Version, InvalidVersion + except ModuleNotFoundError: + raise InstallError("Neither setuptools or packaging found. Please " + "install one of these packages before trying again.") osname = platform.uname().system arch = platform.uname().machine @@ -52,11 +63,6 @@ firedrake_apps = { } -class InstallError(Exception): - # Exception for generic install problems. - pass - - class FiredrakeConfiguration(dict): """A dictionary extended to facilitate the storage of Firedrake configuration information.""" @@ -76,7 +82,7 @@ class FiredrakeConfiguration(dict): self["options"][o] = args.__dict__[o] _persistent_options = ["package_manager", - "minimal_petsc", "mpicc", "mpicxx", "mpif90", "mpiexec", "disable_ssh", + "minimal_petsc", "mpicc", "mpicxx", "mpif90", "mpiexec", "mpihome", "disable_ssh", "honour_petsc_dir", "with_parmetis", "slepc", "packages", "honour_pythonpath", "opencascade", "torch", "jax", @@ -306,6 +312,9 @@ honoured.""", parser.add_argument("--mpiexec", type=str, action="store", default=None, help="MPI launcher. If not set, MPICH will be downloaded and used.") + parser.add_argument("--mpihome", type=str, + action="store", default=None, + help="Location of MPI files. If not set, MPICH will be downloaded and used.") parser.add_argument("--mpi4py-version", help="Specify an exact version of mpi4py to install") parser.add_argument("--show-petsc-configure-options", action="store_true", help="Print out the configure options passed to PETSc and exit") @@ -339,9 +348,9 @@ honoured.""", args = parser.parse_args() # If the user has set any MPI info, they must set them all - if args.mpicc or args.mpicxx or args.mpif90 or args.mpiexec: - if not (args.mpicc and args.mpicxx and args.mpif90 and args.mpiexec): - log.error("If you set any MPI information, you must set all of {mpicc, mpicxx, mpif90, mpiexec}.") + if args.mpicc or args.mpicxx or args.mpif90 or args.mpiexec or args.mpihome: + if not (args.mpicc and args.mpicxx and args.mpif90 and args.mpiexec and args.mpihome): + log.error("If you set any MPI information, you must set all of {mpicc, mpicxx, mpif90, mpiexec, mpihome}.") sys.exit(1) if args.package_branch: @@ -1252,36 +1261,6 @@ def build_and_install_h5py(): log.info("No need to rebuild h5py") -def build_and_install_libsupermesh(cc, cxx, f90, mpiexec): - log.info("Installing libsupermesh") - url = "git+https://github.com/firedrakeproject/libsupermesh.git" - if os.path.exists("libsupermesh"): - changed = git_update("libsupermesh", url) - else: - git_clone(url) - changed = True - if changed: - with directory("libsupermesh"): - check_call(["git", "reset", "--hard"]) - check_call(["git", "clean", "-f", "-x", "-d"]) - check_call(["mkdir", "-p", "build"]) - with directory("build"): - cmd = [ - "cmake", "..", "-DBUILD_SHARED_LIBS=ON", - "-DCMAKE_INSTALL_PREFIX=" + firedrake_env, - "-DMPI_C_COMPILER=" + cc, - "-DMPI_CXX_COMPILER=" + cxx, - "-DMPI_Fortran_COMPILER=" + f90, - "-DCMAKE_Fortran_COMPILER=" + f90, - "-DMPIEXEC_EXECUTABLE=" + mpiexec, - ] - check_call(cmd) - check_call(["make"]) - check_call(["make", "install"]) - else: - log.info("No need to rebuild libsupermesh") - - def build_and_install_pythonocc(): log.info("Installing pythonocc-core") url = "git+https://github.com/tpaviot/pythonocc-core.git@595b0a4e8e60e8d6011bea0cdb54ac878efcfcd2" @@ -1420,7 +1399,7 @@ if args.rebuild_script: sys.exit(0) -def create_compiler_env(cc, cxx, f90): +def create_compiler_env(cc, cxx, f90, mpihome): env = dict() if cc: env["MPICC"] = cc @@ -1439,8 +1418,10 @@ def create_compiler_env(cc, cxx, f90): env["CXX"] = cxx if f90: env["MPIF90"] = f90 - env["MPI_C_COMPILER"] = f90 + env["MPI_Fortran_COMPILER"] = f90 env["F90"] = f90 + if mpihome: + env["MPI_HOME"] = mpihome return env @@ -1688,6 +1669,8 @@ run_pip(["install", "-U", "hatch"]) run_pip(["install", "-U", "editables"]) run_pip(["install", "-U", "pip"]) run_pip(["install", "-U", "wheel"]) +run_pip(["install", "-U", "scikit-build-core"]) +run_pip(["install", "-U", "rtree>=1.2"]) # Extra numpy dependendencies, see # https://github.com/numpy/numpy/blob/main/pyproject.toml @@ -1733,7 +1716,18 @@ cc = options["mpicc"] cxx = options["mpicxx"] f90 = options["mpif90"] mpiexec = options["mpiexec"] -compiler_env = create_compiler_env(cc, cxx, f90) +try: + mpihome = options["mpihome"] +except KeyError: + if cc: + try: + mpihome = os.environ["MPI_HOME"] + except KeyError: + raise InstallError("MPI compilers specified but mpihome is missing from " + "configuration. Please set MPI_HOME to continue.") + else: + mpihome = None +compiler_env = create_compiler_env(cc, cxx, f90, mpihome) os.chdir(firedrake_env) # Dict to store BLAS and OPENBLAS environment variables @@ -1830,7 +1824,8 @@ if mode == "install": cxx = os.path.join(compilerbin, "mpicxx") f90 = os.path.join(compilerbin, "mpif90") mpiexec = os.path.join(compilerbin, "mpiexec") - compiler_env = create_compiler_env(cc, cxx, f90) + mpihome = os.path.join(compilerbin, "..") + compiler_env = create_compiler_env(cc, cxx, f90, mpihome) # Make sure that we link against the right MPI and PETSc shared libraries link_env = { @@ -1877,9 +1872,6 @@ if mode == "install": install(p+"/") sys.path.append(os.getcwd() + "/" + p) - with environment(**compiler_env): - build_and_install_libsupermesh(cc, cxx, f90, mpiexec) - with pipargs("--no-deps"), environment(**compiler_env, **link_env): install("firedrake/") @@ -1986,7 +1978,8 @@ else: cxx = os.path.join(compilerbin, "mpicxx") f90 = os.path.join(compilerbin, "mpif90") mpiexec = os.path.join(compilerbin, "mpiexec") - compiler_env = create_compiler_env(cc, cxx, f90) + mpihome = os.path.join(compilerbin, "..") + compiler_env = create_compiler_env(cc, cxx, f90, mpihome) # Make sure that we link against the right MPI and PETSc shared libraries link_env = { @@ -2059,9 +2052,6 @@ Please consider updating your PETSc manually. else: install(p+"/") - with environment(**compiler_env): - build_and_install_libsupermesh(cc, cxx, f90, mpiexec) - with pipargs("--no-deps"), environment(**compiler_env, **link_env): install("firedrake/") diff --git a/setup.py b/setup.py index cf970c1712..52112e1500 100644 --- a/setup.py +++ b/setup.py @@ -6,6 +6,7 @@ import pybind11 import petsc4py import rtree +import libsupermesh import pkgconfig from dataclasses import dataclass, field from setuptools import setup, find_packages, Extension @@ -145,6 +146,9 @@ def __getitem__(self, key): # In the next 2 linkages we are using `site.getsitepackages()[0]`, which isn't # guaranteed to be the correct place we could also use "$ORIGIN/../../lib_dir", # but that definitely doesn't work with editable installs. +# This is necessary because Python build isolation means that the compile-time +# library dirs (in the isolated build env) are different to the run-time +# library dirs (in the venv). # libspatialindex # example: @@ -159,15 +163,15 @@ def __getitem__(self, key): # libsupermesh # example: -# gcc -I/supermesh/include -# gcc /supermesh/supermesh.cpython-311-x86_64-linux-gnu.so \ +# gcc -Ipath/to/libsupermesh/include +# gcc path/to/libsupermesh/libsupermesh.cpython-311-x86_64-linux-gnu.so \ # -lsupermesh \ -# -Wl,-rpath,$ORIGIN/../../supermesh -supermesh_ = ExternalDependency( - include_dirs=[f"{sys.prefix}/include"], - library_dirs=[f"{sys.prefix}/lib"], +# -Wl,-rpath,$ORIGIN/../../libsupermesh +libsupermesh_ = ExternalDependency( + include_dirs=[libsupermesh.get_include()], + library_dirs=[str(Path(libsupermesh.get_library()).parent)], + runtime_library_dirs=[os.path.join(site.getsitepackages()[0], "libsupermesh", "lib")], libraries=["supermesh"], - runtime_library_dirs=[f"{sys.prefix}/lib"], ) # The following extensions need to be linked accordingly: @@ -221,7 +225,7 @@ def extensions(): name="firedrake.cython.supermeshimpl", language="c", sources=[os.path.join("firedrake", "cython", "supermeshimpl.pyx")], - **(petsc_ + numpy_ + supermesh_) + **(petsc_ + numpy_ + libsupermesh_) )) # pyop2/sparsity.pyx: petsc, numpy, cython_list.append(Extension( diff --git a/tests/firedrake/regression/test_interpolate_p3intmoments.py b/tests/firedrake/regression/test_interpolate_p3intmoments.py index a56cf13ad9..6ef0ef34d2 100644 --- a/tests/firedrake/regression/test_interpolate_p3intmoments.py +++ b/tests/firedrake/regression/test_interpolate_p3intmoments.py @@ -9,7 +9,7 @@ from FIAT.quadrature import make_quadrature from FIAT.polynomial_set import ONPolynomialSet from finat.fiat_elements import ScalarFiatElement -from tsfc.finatinterface import convert, as_fiat_cell +from finat.element_factory import convert, as_fiat_cell import finat.ufl ufcint = UFCInterval() @@ -89,7 +89,7 @@ def __init__(self, cell, degree): super().__init__(P3IntMoments(cell, degree)) -# Replace the old tsfc.finatinterface.convert dispatch with a new one that +# Replace the old finat.element_factory.convert dispatch with a new one that # gives the the new FInAT element for P3 on an interval with variant # "interior-moment" old_convert = convert.dispatch(finat.ufl.FiniteElement) diff --git a/tests/firedrake/regression/test_linesmoother.py b/tests/firedrake/regression/test_linesmoother.py index ba1def7c33..049421ddcc 100644 --- a/tests/firedrake/regression/test_linesmoother.py +++ b/tests/firedrake/regression/test_linesmoother.py @@ -43,7 +43,8 @@ def backend(request): return request.param -def test_linesmoother(mesh, S1family, expected, backend): +@pytest.mark.parametrize("rhs", ["form_rhs", "cofunc_rhs"]) +def test_linesmoother(mesh, S1family, expected, backend, rhs): base_cell = mesh._base_mesh.ufl_cell() S2family = "DG" if base_cell.is_simplex() else "DQ" DGfamily = "DG" if mesh.ufl_cell().is_simplex() else "DQ" @@ -86,6 +87,10 @@ def test_linesmoother(mesh, S1family, expected, backend): f = exp(-rsq) L = inner(f, q)*dx(degree=2*(degree+1)) + if rhs == 'cofunc_rhs': + L = assemble(L) + elif rhs != 'form_rhs': + raise ValueError("Unknown right hand side type") w0 = Function(W) problem = LinearVariationalProblem(a, L, w0, bcs=bcs, aP=aP, form_compiler_parameters={"mode": "vanilla"}) diff --git a/tests/firedrake/regression/test_matrix_free.py b/tests/firedrake/regression/test_matrix_free.py index 6cc6b207f7..f8528a1ed5 100644 --- a/tests/firedrake/regression/test_matrix_free.py +++ b/tests/firedrake/regression/test_matrix_free.py @@ -130,6 +130,7 @@ def test_matrixfree_action(a, V, bcs): @pytest.mark.parametrize("preassembled", [False, True], ids=["variational", "preassembled"]) +@pytest.mark.parametrize("rhs", ["form_rhs", "cofunc_rhs"]) @pytest.mark.parametrize("parameters", [{"ksp_type": "preonly", "pc_type": "python", @@ -168,7 +169,7 @@ def test_matrixfree_action(a, V, bcs): "fieldsplit_1_fieldsplit_1_pc_type": "python", "fieldsplit_1_fieldsplit_1_pc_python_type": "firedrake.AssembledPC", "fieldsplit_1_fieldsplit_1_assembled_pc_type": "lu"}]) -def test_fieldsplitting(mesh, preassembled, parameters): +def test_fieldsplitting(mesh, preassembled, parameters, rhs): V = FunctionSpace(mesh, "CG", 1) P = FunctionSpace(mesh, "DG", 0) Q = VectorFunctionSpace(mesh, "DG", 1) @@ -185,6 +186,10 @@ def test_fieldsplitting(mesh, preassembled, parameters): a = inner(u, v)*dx L = inner(expect, v)*dx + if rhs == 'cofunc_rhs': + L = assemble(L) + elif rhs != 'form_rhs': + raise ValueError("Unknown right hand side type") f = Function(W) diff --git a/tests/firedrake/regression/test_mesh_generation.py b/tests/firedrake/regression/test_mesh_generation.py index 637b0a1502..d2efe2a8b4 100644 --- a/tests/firedrake/regression/test_mesh_generation.py +++ b/tests/firedrake/regression/test_mesh_generation.py @@ -476,6 +476,28 @@ def test_boxmesh_kind(kind, num_cells): assert m.num_cells() == num_cells +@pytest.mark.parallel(nprocs=2) +def test_periodic_unit_cube_hex_cell(): + mesh = PeriodicUnitCubeMesh(3, 3, 3, directions=[True, True, False], hexahedral=True) + x, y, z = SpatialCoordinate(mesh) + V = FunctionSpace(mesh, "CG", 3) + expr = (1 - x) * x + (1 - y) * y + z + f = Function(V).interpolate(expr) + error = assemble((f - expr) ** 2 * dx) + assert error < 1.e-30 + + +@pytest.mark.parallel(nprocs=4) +def test_periodic_unit_cube_hex_facet(): + mesh = PeriodicUnitCubeMesh(3, 3, 3, directions=[True, False, False], hexahedral=True) + for subdomain_id in [1, 2]: + area = assemble(Constant(1.) * dS(domain=mesh, subdomain_id=subdomain_id)) + assert abs(area - 1.0) < 1.e-15 + for subdomain_id in [3, 4, 5, 6]: + area = assemble(Constant(1.) * ds(domain=mesh, subdomain_id=subdomain_id)) + assert abs(area - 1.0) < 1.e-15 + + @pytest.mark.parallel(nprocs=4) def test_split_comm_dm_mesh(): nspace = 2 diff --git a/tests/firedrake/regression/test_nullspace.py b/tests/firedrake/regression/test_nullspace.py index bdef9fd6a3..e4f436ff8d 100644 --- a/tests/firedrake/regression/test_nullspace.py +++ b/tests/firedrake/regression/test_nullspace.py @@ -291,7 +291,8 @@ def test_nullspace_mixed_multiple_components(): @pytest.mark.parallel(nprocs=2) @pytest.mark.parametrize("aux_pc", [False, True], ids=["PC(mu)", "PC(DG0-mu)"]) -def test_near_nullspace_mixed(aux_pc): +@pytest.mark.parametrize("rhs", ["form_rhs", "cofunc_rhs"]) +def test_near_nullspace_mixed(aux_pc, rhs): # test nullspace and nearnullspace for a mixed Stokes system # this is tested on the SINKER case of May and Moresi https://doi.org/10.1016/j.pepi.2008.07.036 # fails in parallel if nullspace is copied to fieldsplit_1_Mp_ksp solve (see PR #3488) @@ -323,6 +324,10 @@ def test_near_nullspace_mixed(aux_pc): f = as_vector((0, -9.8*conditional(inside_box, 2, 1))) L = inner(f, v)*dx + if rhs == 'cofunc_rhs': + L = assemble(L) + elif rhs != 'form_rhs': + raise ValueError("Unknown right hand side type") bcs = [DirichletBC(W[0].sub(0), 0, (1, 2)), DirichletBC(W[0].sub(1), 0, (3, 4))] diff --git a/tests/tsfc/test_create_fiat_element.py b/tests/tsfc/test_create_fiat_element.py deleted file mode 100644 index f8a7d6efc4..0000000000 --- a/tests/tsfc/test_create_fiat_element.py +++ /dev/null @@ -1,150 +0,0 @@ -import pytest - -import FIAT -from FIAT.discontinuous_lagrange import DiscontinuousLagrange as FIAT_DiscontinuousLagrange - -import ufl -import finat.ufl -from tsfc.finatinterface import create_element as _create_element - - -supported_elements = { - # These all map directly to FIAT elements - "Brezzi-Douglas-Marini": FIAT.BrezziDouglasMarini, - "Brezzi-Douglas-Fortin-Marini": FIAT.BrezziDouglasFortinMarini, - "Lagrange": FIAT.Lagrange, - "Nedelec 1st kind H(curl)": FIAT.Nedelec, - "Nedelec 2nd kind H(curl)": FIAT.NedelecSecondKind, - "Raviart-Thomas": FIAT.RaviartThomas, - "Regge": FIAT.Regge, -} -"""A :class:`.dict` mapping UFL element family names to their -FIAT-equivalent constructors.""" - - -def create_element(ufl_element): - """Create a FIAT element given a UFL element.""" - finat_element = _create_element(ufl_element) - return finat_element.fiat_equivalent - - -@pytest.fixture(params=["BDM", - "BDFM", - "Lagrange", - "N1curl", - "N2curl", - "RT", - "Regge"]) -def triangle_names(request): - return request.param - - -@pytest.fixture -def ufl_element(triangle_names): - return finat.ufl.FiniteElement(triangle_names, ufl.triangle, 2) - - -def test_triangle_basic(ufl_element): - element = create_element(ufl_element) - assert isinstance(element, supported_elements[ufl_element.family()]) - - -@pytest.fixture(params=["CG", "DG", "DG L2"], scope="module") -def tensor_name(request): - return request.param - - -@pytest.fixture(params=[ufl.interval, ufl.triangle, - ufl.quadrilateral], - ids=lambda x: x.cellname(), - scope="module") -def ufl_A(request, tensor_name): - return finat.ufl.FiniteElement(tensor_name, request.param, 1) - - -@pytest.fixture -def ufl_B(tensor_name): - return finat.ufl.FiniteElement(tensor_name, ufl.interval, 1) - - -def test_tensor_prod_simple(ufl_A, ufl_B): - tensor_ufl = finat.ufl.TensorProductElement(ufl_A, ufl_B) - - tensor = create_element(tensor_ufl) - A = create_element(ufl_A) - B = create_element(ufl_B) - - assert isinstance(tensor, FIAT.TensorProductElement) - - assert tensor.A is A - assert tensor.B is B - - -@pytest.mark.parametrize(('family', 'expected_cls'), - [('P', FIAT.GaussLobattoLegendre), - ('DP', FIAT.GaussLegendre), - ('DP L2', FIAT.GaussLegendre)]) -def test_interval_variant_default(family, expected_cls): - ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3) - assert isinstance(create_element(ufl_element), expected_cls) - - -@pytest.mark.parametrize(('family', 'variant', 'expected_cls'), - [('P', 'equispaced', FIAT.Lagrange), - ('P', 'spectral', FIAT.GaussLobattoLegendre), - ('DP', 'equispaced', FIAT_DiscontinuousLagrange), - ('DP', 'spectral', FIAT.GaussLegendre), - ('DP L2', 'equispaced', FIAT_DiscontinuousLagrange), - ('DP L2', 'spectral', FIAT.GaussLegendre)]) -def test_interval_variant(family, variant, expected_cls): - ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3, variant=variant) - assert isinstance(create_element(ufl_element), expected_cls) - - -def test_triangle_variant_spectral(): - ufl_element = finat.ufl.FiniteElement('DP', ufl.triangle, 2, variant='spectral') - create_element(ufl_element) - - -def test_triangle_variant_spectral_l2(): - ufl_element = finat.ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral') - create_element(ufl_element) - - -def test_quadrilateral_variant_spectral_q(): - element = create_element(finat.ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral')) - assert isinstance(element.element.A, FIAT.GaussLobattoLegendre) - assert isinstance(element.element.B, FIAT.GaussLobattoLegendre) - - -def test_quadrilateral_variant_spectral_dq(): - element = create_element(finat.ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='spectral')) - assert isinstance(element.element.A, FIAT.GaussLegendre) - assert isinstance(element.element.B, FIAT.GaussLegendre) - - -def test_quadrilateral_variant_spectral_dq_l2(): - element = create_element(finat.ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral')) - assert isinstance(element.element.A, FIAT.GaussLegendre) - assert isinstance(element.element.B, FIAT.GaussLegendre) - - -def test_quadrilateral_variant_spectral_rtcf(): - element = create_element(finat.ufl.FiniteElement('RTCF', ufl.quadrilateral, 2, variant='spectral')) - assert isinstance(element.element._elements[0].A, FIAT.GaussLobattoLegendre) - assert isinstance(element.element._elements[0].B, FIAT.GaussLegendre) - assert isinstance(element.element._elements[1].A, FIAT.GaussLegendre) - assert isinstance(element.element._elements[1].B, FIAT.GaussLobattoLegendre) - - -def test_cache_hit(ufl_element): - A = create_element(ufl_element) - B = create_element(ufl_element) - - assert A is B - - -if __name__ == "__main__": - import os - import sys - pytest.main(args=[os.path.abspath(__file__)] + sys.argv[1:]) diff --git a/tests/tsfc/test_create_finat_element.py b/tests/tsfc/test_create_finat_element.py deleted file mode 100644 index c0f34c292c..0000000000 --- a/tests/tsfc/test_create_finat_element.py +++ /dev/null @@ -1,138 +0,0 @@ -import pytest - -import ufl -import finat.ufl -import finat -from tsfc.finatinterface import create_element, supported_elements - - -@pytest.fixture(params=["BDM", - "BDFM", - "Lagrange", - "N1curl", - "N2curl", - "RT", - "Regge"]) -def triangle_names(request): - return request.param - - -@pytest.fixture -def ufl_element(triangle_names): - return finat.ufl.FiniteElement(triangle_names, ufl.triangle, 2) - - -def test_triangle_basic(ufl_element): - element = create_element(ufl_element) - assert isinstance(element, supported_elements[ufl_element.family()]) - - -@pytest.fixture -def ufl_vector_element(triangle_names): - return finat.ufl.VectorElement(triangle_names, ufl.triangle, 2) - - -def test_triangle_vector(ufl_element, ufl_vector_element): - scalar = create_element(ufl_element) - vector = create_element(ufl_vector_element) - - assert isinstance(vector, finat.TensorFiniteElement) - assert scalar == vector.base_element - - -@pytest.fixture(params=["CG", "DG", "DG L2"]) -def tensor_name(request): - return request.param - - -@pytest.fixture(params=[ufl.interval, ufl.triangle, - ufl.quadrilateral], - ids=lambda x: x.cellname()) -def ufl_A(request, tensor_name): - return finat.ufl.FiniteElement(tensor_name, request.param, 1) - - -@pytest.fixture -def ufl_B(tensor_name): - return finat.ufl.FiniteElement(tensor_name, ufl.interval, 1) - - -def test_tensor_prod_simple(ufl_A, ufl_B): - tensor_ufl = finat.ufl.TensorProductElement(ufl_A, ufl_B) - - tensor = create_element(tensor_ufl) - A = create_element(ufl_A) - B = create_element(ufl_B) - - assert isinstance(tensor, finat.TensorProductElement) - - assert tensor.factors == (A, B) - - -@pytest.mark.parametrize(('family', 'expected_cls'), - [('P', finat.GaussLobattoLegendre), - ('DP', finat.GaussLegendre), - ('DP L2', finat.GaussLegendre)]) -def test_interval_variant_default(family, expected_cls): - ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3) - assert isinstance(create_element(ufl_element), expected_cls) - - -@pytest.mark.parametrize(('family', 'variant', 'expected_cls'), - [('P', 'equispaced', finat.Lagrange), - ('P', 'spectral', finat.GaussLobattoLegendre), - ('DP', 'equispaced', finat.DiscontinuousLagrange), - ('DP', 'spectral', finat.GaussLegendre), - ('DP L2', 'equispaced', finat.DiscontinuousLagrange), - ('DP L2', 'spectral', finat.GaussLegendre)]) -def test_interval_variant(family, variant, expected_cls): - ufl_element = finat.ufl.FiniteElement(family, ufl.interval, 3, variant=variant) - assert isinstance(create_element(ufl_element), expected_cls) - - -def test_triangle_variant_spectral(): - ufl_element = finat.ufl.FiniteElement('DP', ufl.triangle, 2, variant='spectral') - create_element(ufl_element) - - -def test_triangle_variant_spectral_l2(): - ufl_element = finat.ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral') - create_element(ufl_element) - - -def test_quadrilateral_variant_spectral_q(): - element = create_element(finat.ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral')) - assert isinstance(element.product.factors[0], finat.GaussLobattoLegendre) - assert isinstance(element.product.factors[1], finat.GaussLobattoLegendre) - - -def test_quadrilateral_variant_spectral_dq(): - element = create_element(finat.ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='spectral')) - assert isinstance(element.product.factors[0], finat.GaussLegendre) - assert isinstance(element.product.factors[1], finat.GaussLegendre) - - -def test_quadrilateral_variant_spectral_dq_l2(): - element = create_element(finat.ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral')) - assert isinstance(element.product.factors[0], finat.GaussLegendre) - assert isinstance(element.product.factors[1], finat.GaussLegendre) - - -def test_cache_hit(ufl_element): - A = create_element(ufl_element) - B = create_element(ufl_element) - - assert A is B - - -def test_cache_hit_vector(ufl_vector_element): - A = create_element(ufl_vector_element) - B = create_element(ufl_vector_element) - - assert A is B - - -if __name__ == "__main__": - import os - import sys - pytest.main(args=[os.path.abspath(__file__)] + sys.argv[1:]) diff --git a/tests/tsfc/test_dual_evaluation.py b/tests/tsfc/test_dual_evaluation.py index b4f6e9770a..85f1617678 100644 --- a/tests/tsfc/test_dual_evaluation.py +++ b/tests/tsfc/test_dual_evaluation.py @@ -1,7 +1,7 @@ import pytest import ufl import finat.ufl -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from tsfc import compile_expression_dual_evaluation diff --git a/tests/tsfc/test_interpolation_factorisation.py b/tests/tsfc/test_interpolation_factorisation.py index b3d4e3288b..4355c24b1f 100644 --- a/tests/tsfc/test_interpolation_factorisation.py +++ b/tests/tsfc/test_interpolation_factorisation.py @@ -7,7 +7,7 @@ from finat.ufl import FiniteElement, VectorElement, TensorElement from tsfc import compile_expression_dual_evaluation -from tsfc.finatinterface import create_element +from finat.element_factory import create_element @pytest.fixture(params=[interval, quadrilateral, hexahedron], diff --git a/tests/tsfc/test_tsfc_274.py b/tests/tsfc/test_tsfc_274.py index 453d8746e8..39da190496 100644 --- a/tests/tsfc/test_tsfc_274.py +++ b/tests/tsfc/test_tsfc_274.py @@ -2,7 +2,7 @@ import numpy from finat.point_set import PointSet from gem.interpreter import evaluate -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from ufl import quadrilateral from finat.ufl import FiniteElement, RestrictedElement diff --git a/tsfc/driver.py b/tsfc/driver.py index 6e3c3baaf3..38a780a0ad 100644 --- a/tsfc/driver.py +++ b/tsfc/driver.py @@ -47,14 +47,13 @@ """ -def compile_form(form, prefix="form", parameters=None, interface=None, diagonal=False, log=False): +def compile_form(form, prefix="form", parameters=None, interface=None, diagonal=False): """Compiles a UFL form into a set of assembly kernels. :arg form: UFL form :arg prefix: kernel name will start with this string :arg parameters: parameters object :arg diagonal: Are we building a kernel for the diagonal of a rank-2 element tensor? - :arg log: bool if the Kernel should be profiled with Log events :returns: list of kernels """ cpu_time = time.time() @@ -71,7 +70,7 @@ def compile_form(form, prefix="form", parameters=None, interface=None, diagonal= kernels = [] for integral_data in fd.integral_data: start = time.time() - kernel = compile_integral(integral_data, fd, prefix, parameters, interface=interface, diagonal=diagonal, log=log) + kernel = compile_integral(integral_data, fd, prefix, parameters, interface=interface, diagonal=diagonal) if kernel is not None: kernels.append(kernel) logger.info(GREEN % "compile_integral finished in %g seconds.", time.time() - start) @@ -80,7 +79,7 @@ def compile_form(form, prefix="form", parameters=None, interface=None, diagonal= return kernels -def compile_integral(integral_data, form_data, prefix, parameters, interface, *, diagonal=False, log=False): +def compile_integral(integral_data, form_data, prefix, parameters, interface, *, diagonal=False): """Compiles a UFL integral into an assembly kernel. :arg integral_data: UFL integral data @@ -89,7 +88,6 @@ def compile_integral(integral_data, form_data, prefix, parameters, interface, *, :arg parameters: parameters object :arg interface: backend module for the kernel interface :arg diagonal: Are we building a kernel for the diagonal of a rank-2 element tensor? - :arg log: bool if the Kernel should be profiled with Log events :returns: a kernel constructed by the kernel interface """ parameters = preprocess_parameters(parameters) @@ -137,7 +135,7 @@ def compile_integral(integral_data, form_data, prefix, parameters, interface, *, integrand_exprs = builder.compile_integrand(integrand, params, ctx) integral_exprs = builder.construct_integrals(integrand_exprs, params) builder.stash_integrals(integral_exprs, params, ctx) - return builder.construct_kernel(kernel_name, ctx, log) + return builder.construct_kernel(kernel_name, ctx, parameters["add_petsc_events"]) def preprocess_parameters(parameters): @@ -157,7 +155,7 @@ def preprocess_parameters(parameters): def compile_expression_dual_evaluation(expression, to_element, ufl_element, *, domain=None, interface=None, - parameters=None, log=False): + parameters=None): """Compile a UFL expression to be evaluated against a compile-time known reference element's dual basis. Useful for interpolating UFL expressions into e.g. N1curl spaces. @@ -168,7 +166,6 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *, :arg domain: optional UFL domain the expression is defined on (required when expression contains no domain). :arg interface: backend module for the kernel interface :arg parameters: parameters object - :arg log: bool if the Kernel should be profiled with Log events :returns: Loopy-based ExpressionKernel object. """ if parameters is None: @@ -267,7 +264,7 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *, builder.register_requirements([evaluation]) builder.set_output(return_var) # Build kernel tuple - return builder.construct_kernel(impero_c, index_names, needs_external_coords, log=log) + return builder.construct_kernel(impero_c, index_names, needs_external_coords, parameters["add_petsc_events"]) class DualEvaluationCallable(object): diff --git a/tsfc/fem.py b/tsfc/fem.py index 99251ed0c6..a5fa2009b1 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -15,6 +15,7 @@ PhysicalGeometry) from finat.point_set import PointSet, PointSingleton from finat.quadrature import make_quadrature +from finat.element_factory import as_fiat_cell, create_element from gem.node import traversal from gem.optimise import constant_fold_zero, ffc_rounding from gem.unconcatenate import unconcatenate @@ -32,7 +33,6 @@ from ufl.domain import extract_unique_domain from tsfc import ufl2gem -from tsfc.finatinterface import as_fiat_cell, create_element from tsfc.kernel_interface import ProxyKernelInterface from tsfc.modified_terminals import (analyse_modified_terminal, construct_modified_terminal) @@ -270,9 +270,10 @@ def get_quadrature_rule(fiat_cell, integration_dim, quadrature_degree, scheme): def make_basis_evaluation_key(ctx, finat_element, mt, entity_id): + ufl_element = mt.terminal.ufl_element() domain = extract_unique_domain(mt.terminal) coordinate_element = domain.ufl_coordinate_element() - return (finat_element, mt.local_derivatives, ctx.point_set, ctx.integration_dim, entity_id, coordinate_element, mt.restriction) + return (ufl_element, mt.local_derivatives, ctx.point_set, ctx.integration_dim, entity_id, coordinate_element, mt.restriction) class PointSetContext(ContextBase): @@ -697,8 +698,7 @@ def take_singleton(xs): for alpha, tables in per_derivative.items()} # Coefficient evaluation - ctx.index_cache.setdefault(terminal.ufl_element(), element.get_indices()) - beta = ctx.index_cache[terminal.ufl_element()] + beta = ctx.index_cache.setdefault(terminal.ufl_element(), element.get_indices()) zeta = element.get_value_indices() vec_beta, = gem.optimise.remove_componenttensors([gem.Indexed(vec, beta)]) value_dict = {} diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py deleted file mode 100644 index b7e3d0ad72..0000000000 --- a/tsfc/finatinterface.py +++ /dev/null @@ -1,365 +0,0 @@ -# This file was modified from FFC -# (http://bitbucket.org/fenics-project/ffc), copyright notice -# reproduced below. -# -# Copyright (C) 2009-2013 Kristian B. Oelgaard and Anders Logg -# -# This file is part of FFC. -# -# FFC is free software: you can redistribute it and/or modify -# it under the terms of the GNU Lesser General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# FFC is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Lesser General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with FFC. If not, see . - -import weakref -from functools import singledispatch - -import FIAT -import finat -import finat.ufl -import ufl - -__all__ = ("as_fiat_cell", "create_base_element", - "create_element", "supported_elements") - - -supported_elements = { - # These all map directly to FInAT elements - "Bernstein": finat.Bernstein, - "Bernardi-Raugel": finat.BernardiRaugel, - "Bernardi-Raugel Bubble": finat.BernardiRaugelBubble, - "Brezzi-Douglas-Marini": finat.BrezziDouglasMarini, - "Brezzi-Douglas-Fortin-Marini": finat.BrezziDouglasFortinMarini, - "Bubble": finat.Bubble, - "FacetBubble": finat.FacetBubble, - "Crouzeix-Raviart": finat.CrouzeixRaviart, - "Discontinuous Lagrange": finat.DiscontinuousLagrange, - "Discontinuous Raviart-Thomas": lambda c, d: finat.DiscontinuousElement(finat.RaviartThomas(c, d)), - "Discontinuous Taylor": finat.DiscontinuousTaylor, - "Gauss-Legendre": finat.GaussLegendre, - "Gauss-Lobatto-Legendre": finat.GaussLobattoLegendre, - "HDiv Trace": finat.HDivTrace, - "Hellan-Herrmann-Johnson": finat.HellanHerrmannJohnson, - "Johnson-Mercier": finat.JohnsonMercier, - "Nonconforming Arnold-Winther": finat.ArnoldWintherNC, - "Conforming Arnold-Winther": finat.ArnoldWinther, - "Hu-Zhang": finat.HuZhang, - "Hermite": finat.Hermite, - "Kong-Mulder-Veldhuizen": finat.KongMulderVeldhuizen, - "Argyris": finat.Argyris, - "Hsieh-Clough-Tocher": finat.HsiehCloughTocher, - "QuadraticPowellSabin6": finat.QuadraticPowellSabin6, - "QuadraticPowellSabin12": finat.QuadraticPowellSabin12, - "Reduced-Hsieh-Clough-Tocher": finat.ReducedHsiehCloughTocher, - "Mardal-Tai-Winther": finat.MardalTaiWinther, - "Alfeld-Sorokina": finat.AlfeldSorokina, - "Arnold-Qin": finat.ArnoldQin, - "Reduced-Arnold-Qin": finat.ReducedArnoldQin, - "Christiansen-Hu": finat.ChristiansenHu, - "Guzman-Neilan 1st kind H1": finat.GuzmanNeilanFirstKindH1, - "Guzman-Neilan 2nd kind H1": finat.GuzmanNeilanSecondKindH1, - "Guzman-Neilan Bubble": finat.GuzmanNeilanBubble, - "Guzman-Neilan H1(div)": finat.GuzmanNeilanH1div, - "Morley": finat.Morley, - "Bell": finat.Bell, - "Lagrange": finat.Lagrange, - "Nedelec 1st kind H(curl)": finat.Nedelec, - "Nedelec 2nd kind H(curl)": finat.NedelecSecondKind, - "Raviart-Thomas": finat.RaviartThomas, - "Regge": finat.Regge, - "Gopalakrishnan-Lederer-Schoberl 1st kind": finat.GopalakrishnanLedererSchoberlFirstKind, - "Gopalakrishnan-Lederer-Schoberl 2nd kind": finat.GopalakrishnanLedererSchoberlSecondKind, - "BDMCE": finat.BrezziDouglasMariniCubeEdge, - "BDMCF": finat.BrezziDouglasMariniCubeFace, - # These require special treatment below - "DQ": None, - "Q": None, - "RTCE": None, - "RTCF": None, - "NCE": None, - "NCF": None, - "Real": finat.Real, - "DPC": finat.DPC, - "S": finat.Serendipity, - "SminusF": finat.TrimmedSerendipityFace, - "SminusDiv": finat.TrimmedSerendipityDiv, - "SminusE": finat.TrimmedSerendipityEdge, - "SminusCurl": finat.TrimmedSerendipityCurl, - "DPC L2": finat.DPC, - "Discontinuous Lagrange L2": finat.DiscontinuousLagrange, - "Gauss-Legendre L2": finat.GaussLegendre, - "DQ L2": None, - "Direct Serendipity": finat.DirectSerendipity, -} -"""A :class:`.dict` mapping UFL element family names to their -FInAT-equivalent constructors. If the value is ``None``, the UFL -element is supported, but must be handled specially because it doesn't -have a direct FInAT equivalent.""" - - -def as_fiat_cell(cell): - """Convert a ufl cell to a FIAT cell. - - :arg cell: the :class:`ufl.Cell` to convert.""" - if not isinstance(cell, ufl.AbstractCell): - raise ValueError("Expecting a UFL Cell") - return FIAT.ufc_cell(cell) - - -@singledispatch -def convert(element, **kwargs): - """Handler for converting UFL elements to FInAT elements. - - :arg element: The UFL element to convert. - - Do not use this function directly, instead call - :func:`create_element`.""" - if element.family() in supported_elements: - raise ValueError("Element %s supported, but no handler provided" % element) - raise ValueError("Unsupported element type %s" % type(element)) - - -cg_interval_variants = { - "fdm": finat.FDMLagrange, - "fdm_ipdg": finat.FDMLagrange, - "fdm_quadrature": finat.FDMQuadrature, - "fdm_broken": finat.FDMBrokenH1, - "fdm_hermite": finat.FDMHermite, -} - - -dg_interval_variants = { - "fdm": finat.FDMDiscontinuousLagrange, - "fdm_quadrature": finat.FDMDiscontinuousLagrange, - "fdm_ipdg": lambda *args: finat.DiscontinuousElement(finat.FDMLagrange(*args)), - "fdm_broken": finat.FDMBrokenL2, -} - - -# Base finite elements first -@convert.register(finat.ufl.FiniteElement) -def convert_finiteelement(element, **kwargs): - cell = as_fiat_cell(element.cell) - if element.family() == "Quadrature": - degree = element.degree() - scheme = element.quadrature_scheme() - if degree is None or scheme is None: - raise ValueError("Quadrature scheme and degree must be specified!") - - return finat.make_quadrature_element(cell, degree, scheme), set() - lmbda = supported_elements[element.family()] - if element.family() == "Real" and element.cell.cellname() in {"quadrilateral", "hexahedron"}: - lmbda = None - element = finat.ufl.FiniteElement("DQ", element.cell, 0) - if lmbda is None: - if element.cell.cellname() == "quadrilateral": - # Handle quadrilateral short names like RTCF and RTCE. - element = element.reconstruct(cell=quadrilateral_tpc) - elif element.cell.cellname() == "hexahedron": - # Handle hexahedron short names like NCF and NCE. - element = element.reconstruct(cell=hexahedron_tpc) - else: - raise ValueError("%s is supported, but handled incorrectly" % - element.family()) - finat_elem, deps = _create_element(element, **kwargs) - return finat.FlattenedDimensions(finat_elem), deps - - finat_kwargs = {} - kind = element.variant() - if kind is None: - kind = 'spectral' # default variant - - if element.family() == "Lagrange": - if kind == 'spectral': - lmbda = finat.GaussLobattoLegendre - elif element.cell.cellname() == "interval" and kind in cg_interval_variants: - lmbda = cg_interval_variants[kind] - elif any(map(kind.startswith, ['integral', 'demkowicz', 'fdm'])): - lmbda = finat.IntegratedLegendre - finat_kwargs["variant"] = kind - elif kind in ['mgd', 'feec', 'qb', 'mse']: - degree = element.degree() - shift_axes = kwargs["shift_axes"] - restriction = kwargs["restriction"] - deps = {"shift_axes", "restriction"} - return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction), deps - else: - # Let FIAT handle the general case - lmbda = finat.Lagrange - finat_kwargs["variant"] = kind - - elif element.family() in ["Discontinuous Lagrange", "Discontinuous Lagrange L2"]: - if kind == 'spectral': - lmbda = finat.GaussLegendre - elif element.cell.cellname() == "interval" and kind in dg_interval_variants: - lmbda = dg_interval_variants[kind] - elif any(map(kind.startswith, ['integral', 'demkowicz', 'fdm'])): - lmbda = finat.Legendre - finat_kwargs["variant"] = kind - elif kind in ['mgd', 'feec', 'qb', 'mse']: - degree = element.degree() - shift_axes = kwargs["shift_axes"] - restriction = kwargs["restriction"] - deps = {"shift_axes", "restriction"} - return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction, continuous=False), deps - else: - # Let FIAT handle the general case - lmbda = finat.DiscontinuousLagrange - finat_kwargs["variant"] = kind - - elif element.variant() is not None: - finat_kwargs["variant"] = element.variant() - - return lmbda(cell, element.degree(), **finat_kwargs), set() - - -# Element modifiers and compound element types -@convert.register(finat.ufl.BrokenElement) -def convert_brokenelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.DiscontinuousElement(finat_elem), deps - - -@convert.register(finat.ufl.EnrichedElement) -def convert_enrichedelement(element, **kwargs): - elements, deps = zip(*[_create_element(elem, **kwargs) - for elem in element._elements]) - return finat.EnrichedElement(elements), set.union(*deps) - - -@convert.register(finat.ufl.NodalEnrichedElement) -def convert_nodalenrichedelement(element, **kwargs): - elements, deps = zip(*[_create_element(elem, **kwargs) - for elem in element._elements]) - return finat.NodalEnrichedElement(elements), set.union(*deps) - - -@convert.register(finat.ufl.MixedElement) -def convert_mixedelement(element, **kwargs): - elements, deps = zip(*[_create_element(elem, **kwargs) - for elem in element.sub_elements]) - return finat.MixedElement(elements), set.union(*deps) - - -@convert.register(finat.ufl.VectorElement) -@convert.register(finat.ufl.TensorElement) -def convert_tensorelement(element, **kwargs): - inner_elem, deps = _create_element(element.sub_elements[0], **kwargs) - shape = element.reference_value_shape - shape = shape[:len(shape) - len(inner_elem.value_shape)] - shape_innermost = kwargs["shape_innermost"] - return (finat.TensorFiniteElement(inner_elem, shape, not shape_innermost), - deps | {"shape_innermost"}) - - -@convert.register(finat.ufl.TensorProductElement) -def convert_tensorproductelement(element, **kwargs): - cell = element.cell - if type(cell) is not ufl.TensorProductCell: - raise ValueError("TensorProductElement not on TensorProductCell?") - shift_axes = kwargs["shift_axes"] - dim_offset = 0 - elements = [] - deps = set() - for elem in element.sub_elements: - kwargs["shift_axes"] = shift_axes + dim_offset - dim_offset += elem.cell.topological_dimension() - finat_elem, ds = _create_element(elem, **kwargs) - elements.append(finat_elem) - deps.update(ds) - return finat.TensorProductElement(elements), deps - - -@convert.register(finat.ufl.HDivElement) -def convert_hdivelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.HDivElement(finat_elem), deps - - -@convert.register(finat.ufl.HCurlElement) -def convert_hcurlelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.HCurlElement(finat_elem), deps - - -@convert.register(finat.ufl.WithMapping) -def convert_withmapping(element, **kwargs): - return _create_element(element.wrapee, **kwargs) - - -@convert.register(finat.ufl.RestrictedElement) -def convert_restrictedelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps - - -hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval) -quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval) -_cache = weakref.WeakKeyDictionary() - - -def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None): - """Create a FInAT element (suitable for tabulating with) given a UFL element. - - :arg ufl_element: The UFL element to create a FInAT element from. - :arg shape_innermost: Vector/tensor indices come after basis function indices - :arg restriction: cell restriction in interior facet integrals - (only for runtime tabulated elements) - """ - finat_element, deps = _create_element(ufl_element, - shape_innermost=shape_innermost, - shift_axes=shift_axes, - restriction=restriction) - return finat_element - - -def _create_element(ufl_element, **kwargs): - """A caching wrapper around :py:func:`convert`. - - Takes a UFL element and an unspecified set of parameter options, - and returns the converted element with the set of keyword names - that were relevant for conversion. - """ - # Look up conversion in cache - try: - cache = _cache[ufl_element] - except KeyError: - _cache[ufl_element] = {} - cache = _cache[ufl_element] - - for key, finat_element in cache.items(): - # Cache hit if all relevant parameter values match. - if all(kwargs[param] == value for param, value in key): - return finat_element, set(param for param, value in key) - - # Convert if cache miss - if ufl_element.cell is None: - raise ValueError("Don't know how to build element when cell is not given") - - finat_element, deps = convert(ufl_element, **kwargs) - - # Store conversion in cache - key = frozenset((param, kwargs[param]) for param in deps) - cache[key] = finat_element - - # Forward result - return finat_element, deps - - -def create_base_element(ufl_element, **kwargs): - """Create a "scalar" base FInAT element given a UFL element. - Takes a UFL element and an unspecified set of parameter options, - and returns the converted element. - """ - finat_element = create_element(ufl_element, **kwargs) - if isinstance(finat_element, finat.TensorFiniteElement): - finat_element = finat_element.base_element - return finat_element diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index f757af951f..09bb881678 100644 --- a/tsfc/kernel_interface/common.py +++ b/tsfc/kernel_interface/common.py @@ -15,7 +15,7 @@ from gem.optimise import remove_componenttensors as prune from numpy import asarray from tsfc import fem, ufl_utils -from tsfc.finatinterface import as_fiat_cell, create_element +from finat.element_factory import as_fiat_cell, create_element from tsfc.kernel_interface import KernelInterface from tsfc.logging import logger from ufl.utils.sequences import max_degree @@ -394,7 +394,7 @@ def lower_integral_type(fiat_cell, integral_type): elif integral_type == 'exterior_facet_top': entity_ids = [1] else: - entity_ids = list(range(len(fiat_cell.get_topology()[integration_dim]))) + entity_ids = list(fiat_cell.get_topology()[integration_dim]) return integration_dim, entity_ids diff --git a/tsfc/kernel_interface/firedrake_loopy.py b/tsfc/kernel_interface/firedrake_loopy.py index 0969854cd3..6ace74a0e1 100644 --- a/tsfc/kernel_interface/firedrake_loopy.py +++ b/tsfc/kernel_interface/firedrake_loopy.py @@ -12,7 +12,7 @@ import loopy as lp from tsfc import kernel_args, fem -from tsfc.finatinterface import create_element +from finat.element_factory import create_element from tsfc.kernel_interface.common import KernelBuilderBase as _KernelBuilderBase, KernelBuilderMixin, get_index_names, check_requirements, prepare_coefficient, prepare_arguments, prepare_constant from tsfc.loopy import generate as generate_loopy diff --git a/tsfc/parameters.py b/tsfc/parameters.py index 1277713ad5..af44ce0cd4 100644 --- a/tsfc/parameters.py +++ b/tsfc/parameters.py @@ -20,6 +20,9 @@ # So that tests pass (needs to match scalar_type) "scalar_type_c": "double", + + # Whether to wrap the generated kernels in a PETSc event + "add_petsc_events": False, }