FEM capabilities in the problem/data classes

.github/workflows/pytest.yml

@@ -0,0 +1,129 @@
+name: Pytest
+
+on:
+  # Trigger the workflow on push or pull request,
+  # but only for the master branch
+  push:
+    branches:
+      - master
+  pull_request:
+    branches:
+      - master
+
+jobs:
+  examples:
+    name: Pytest
+    runs-on: ubuntu-latest
+
+    env:
+      DEVITO_COMPILER: gcc
+      DEVITO_LANGUAGE: openmp
+      PYTHON_VERSION: 3.11.12
+
+    strategy:
+      # Prevent all build to stop if a single one fails
+      fail-fast: false
+
+    steps:
+    - name: Checkout stride
+      uses: actions/checkout@v3
+      with:
+        path: stride
+
+    - name: Checkout devito
+      uses: actions/checkout@v3
+      with:
+        repository: devitocodes/devito
+        path: devito
+
+    - name: Setup conda
+      uses: mamba-org/setup-micromamba@v2
+      with:
+        environment-file: stride/environment.yml
+        init-shell: bash
+        cache-environment: true
+        post-cleanup: 'all'
+
+    - name: Install dependencies
+      shell: bash -l {0}
+      run: |
+        cd stride
+        pip install -e .
+
+    - name: Install devito
+      shell: bash -l {0}
+      run: |
+        cd devito
+        pip install -e .
+
+    - name: Pytest
+      shell: bash -l {0}
+      run: |
+        cd stride
+        python -m pytest stride/tests/ -v
+
+  pytest-dolfinx:
+    name: Pytest DOLFINx
+    runs-on: ubuntu-latest
+
+    env:
+      DEVITO_COMPILER: gcc
+      DEVITO_LANGUAGE: openmp
+      PYTHON_VERSION: 3.11.12
+
+    strategy:
+      # Prevent all build to stop if a single one fails
+      fail-fast: false
+
+    steps:
+    - name: Checkout stride
+      uses: actions/checkout@v3
+      with:
+        path: stride
+
+    - name: Checkout devito
+      uses: actions/checkout@v3
+      with:
+        repository: devitocodes/devito
+        path: devito
+
+    - name: Setup conda
+      uses: mamba-org/setup-micromamba@v2
+      with:
+        environment-file: stride/environment.yml
+        init-shell: bash
+        cache-environment: true
+        post-cleanup: 'all'
+
+    # applied over the base environment rather than replacing it, the same way the file is meant
+    # to be used locally. The complex build is the one that covers both cases: it runs the
+    # resistive tests as well, whereas a real build has to skip the capacitive ones
+    - name: Install FEM dependencies
+      shell: bash -l {0}
+      run: |
+        micromamba install -y -n stride -f stride/environment-fem-complex.yml
+
+    # a missing or real-scalar build would turn the FEM tests into skips rather than failures,
+    # so the suite would go green having run none of them
+    - name: Check the build is complex
+      shell: bash -l {0}
+      run: |
+        python -c "import dolfinx, numpy; assert numpy.dtype(dolfinx.default_scalar_type).kind == 'c'"
+
+    - name: Install dependencies
+      shell: bash -l {0}
+      run: |
+        cd stride
+        pip install -e .
+
+    - name: Install devito
+      shell: bash -l {0}
+      run: |
+        cd devito
+        pip install -e .
+
+    - name: Pytest
+      shell: bash -l {0}
+      run: |
+        cd stride
+        python -m pytest stride/tests/ -v

environment-fem-complex.yml

@@ -0,0 +1,34 @@
+# Finite-element extras for stride, complex-scalar PETSc build.
+#
+# This file is ADDITIVE: it lists only what the FEM problem types need on top of the base
+# environment, and is applied over an existing environment rather than creating one:
+#
+#     conda env create -n stride-dolfinx-complex -f environment.yml
+#     conda env update -n stride-dolfinx-complex -f environment-fem-complex.yml
+#
+# Copying the base dependency list into this file instead is what lets the two drift apart, so
+# please keep it additive.
+#
+# DOLFINx cannot be installed with pip -- there are no wheels -- so these dependencies are
+# conda-only and deliberately not in requirements-optional.txt. Everything in stride works
+# without them; only the FEM problem types need them.
+#
+# PETSc is built for either real or complex scalars and the two cannot coexist. Use this file
+# for the complex build, environment-fem-real.yml for the real one.
+#
+# The MPI implementation is pinned here rather than in the base environment, which needs no MPI
+# at all. Mixing implementations between the base and this overlay causes hangs at exit inside
+# MPI finalisation, so if the base environment ever grows an MPI dependency, the two pins must
+# agree.
+name: stride-dolfinx-complex
+channels:
+  - conda-forge
+dependencies:
+  - fenics-dolfinx==0.9.0
+  - petsc=*=*complex*
+  - mpi=*=*openmpi*
+  - openmpi
+  - petsc4py
+  - mpi4py
+  - gmsh            # mesh generation, plus its Python bindings
+  - python-gmsh

environment-fem-real.yml

@@ -0,0 +1,34 @@
+# Finite-element extras for stride, real-scalar PETSc build.
+#
+# This file is ADDITIVE: it lists only what the FEM problem types need on top of the base
+# environment, and is applied over an existing environment rather than creating one:
+#
+#     conda env create -n stride-dolfinx-real -f environment.yml
+#     conda env update -n stride-dolfinx-real -f environment-fem-real.yml
+#
+# Copying the base dependency list into this file instead is what lets the two drift apart, so
+# please keep it additive.
+#
+# DOLFINx cannot be installed with pip -- there are no wheels -- so these dependencies are
+# conda-only and deliberately not in requirements-optional.txt. Everything in stride works
+# without them; only the FEM problem types need them.
+#
+# PETSc is built for either real or complex scalars and the two cannot coexist. Use this file
+# for the real build, environment-fem-complex.yml for the complex one.
+#
+# The MPI implementation is pinned here rather than in the base environment, which needs no MPI
+# at all. Mixing implementations between the base and this overlay causes hangs at exit inside
+# MPI finalisation, so if the base environment ever grows an MPI dependency, the two pins must
+# agree.
+name: stride-dolfinx-real
+channels:
+  - conda-forge
+dependencies:
+  - fenics-dolfinx==0.9.0
+  - petsc=*=*real*
+  - mpi=*=*openmpi*
+  - openmpi
+  - petsc4py
+  - mpi4py
+  - gmsh            # mesh generation, plus its Python bindings
+  - python-gmsh

stride/problem/base.py

@@ -1,7 +1,7 @@
 
 from mosaic import h5
 
-from .domain import Space, Time, SlowTime, Grid
+from .domain import Space, MeshedSpace, Time, SlowTime, Grid
 
 
 __all__ = ['Gridded', 'Saved', 'GriddedSaved', 'ProblemBase']
@@ -245,12 +245,41 @@ def load(self, *args, **kwargs):
         with h5.HDF5(*args, **kwargs, mode='r') as file:
             description = file.load(filter=kwargs.pop('filter', None), only=kwargs.pop('only', None))
 
-            # TODO If there's already a grid and they don't match, resample instead
             if 'space' in description and self._grid.space is None:
-                space = Space(shape=description.space.shape,
-                              spacing=description.space.spacing,
-                              extra=description.space.extra,
-                              absorbing=description.space.absorbing)
+                # NOTE the discriminator has to be the keys of description.space, not those of
+                # description: StructuredData.__get_desc__ also writes a top-level 'shape', so
+                # testing the outer description would send every field down the structured branch
+                space_description = description.space
+
+                if 'shape' in space_description:
+                    space = Space(shape=space_description.shape,
+                                  spacing=space_description.spacing,
+                                  extra=space_description.extra,
+                                  absorbing=space_description.absorbing)
+
+                elif 'nodes' in space_description:
+                    # the cell type is stored as a string, so HDF5 hands it back as bytes, and
+                    # the degree comes back as a numpy integer, which is not an int as far as
+                    # isinstance is concerned. Both have to be normalised before construction
+                    cell_type = self._materialise(space_description.get('cell_type', None))
+
+                    if isinstance(cell_type, bytes):
+                        cell_type = cell_type.decode()
+
+                    geometry_degree = self._materialise(
+                        space_description.get('geometry_degree', 1))
+
+                    space = MeshedSpace(
+                        nodes=self._materialise(space_description.nodes),
+                        cells=self._materialise(space_description.get('cells', None)),
+                        cell_tags=self._materialise(space_description.get('cell_tags', None)),
+                        cell_type=cell_type,
+                        geometry_degree=int(geometry_degree),
+                    )
+
+                else:
+                    raise ValueError('Unrecognised space description with keys %s'
+                                     % sorted(space_description.keys()))
 
                 self._grid.space = space
 
@@ -275,6 +304,31 @@ def load(self, *args, **kwargs):
             kwargs['filename'] = kwargs.pop('filename', file.filename)
             self.__set_desc__(description, **kwargs)
 
+    @staticmethod
+    def _materialise(value):
+        """
+        Read a lazily-loaded description entry into memory.
+
+        Loading a description defaults to being lazy, which yields the open dataset with a
+        ``load`` attribute attached rather than an ndarray. Anything that has to outlive the
+        open file must be materialised explicitly.
+
+        Parameters
+        ----------
+        value : object
+            Entry of a loaded description.
+
+        Returns
+        -------
+        object
+            The entry, read into memory if it was lazy.
+
+        """
+        if hasattr(value, 'load'):
+            return value.load()
+
+        return value
+
     def grid_description(self):
         """
         Get a description of the grid of the object.
@@ -289,12 +343,38 @@ def grid_description(self):
 
         if self.space is not None:
             space = self.space
-            grid_description['space'] = {
-                'shape': space.shape,
-                'spacing': space.spacing,
-                'extra': space.extra,
-                'absorbing': space.absorbing,
-            }
+            if isinstance(space, Space):
+                grid_description['space'] = {
+                    'shape': space.shape,
+                    'spacing': space.spacing,
+                    'extra': space.extra,
+                    'absorbing': space.absorbing,
+                }
+
+            elif isinstance(space, MeshedSpace):
+                # the connectivity travels with the nodes: a node table on its own is not a mesh
+                # and cannot be handed back to a solver. Facet tags are deliberately left out,
+                # being meaningless without the facet connectivity, which is not stored either
+                space_description = {'nodes': space.nodes}
+
+                if space.cells is not None:
+                    space_description['cells'] = space.cells
+
+                if space.cell_tags is not None:
+                    space_description['cell_tags'] = space.cell_tags
+
+                # the discretisation, without which the mesh cannot be rebuilt: a node and cell
+                # table alone does not say what shape a cell is or to what degree it is mapped
+                if space.cell_type is not None:
+                    space_description['cell_type'] = space.cell_type
+
+                space_description['geometry_degree'] = space.geometry_degree
+
+                grid_description['space'] = space_description
+
+            else:
+                raise TypeError('Cannot serialise a grid with space of type %s'
+                                % type(space).__name__)
 
         if self.time is not None:
             time = self.time