diff --git a/.github/actionlint.yaml b/.github/actionlint.yml similarity index 100% rename from .github/actionlint.yaml rename to .github/actionlint.yml diff --git a/.github/actions/install/action.yml b/.github/actions/install/action.yml new file mode 100644 index 0000000000..8ba3dbbe57 --- /dev/null +++ b/.github/actions/install/action.yml @@ -0,0 +1,162 @@ +name: Install Firedrake + +inputs: + os: + description: Operating system ('linux' or 'macos') + type: string + required: true + source_ref: + description: Source branch + type: string + required: true + base_ref: + description: Target branch ('main' or 'release') + type: string + required: true + scalar_type: + description: Scalar type ('default' or 'complex') + type: string + default: default + gpu: + description: Target GPU platform + type: string + default: none + deps: + description: Optional dependencies to install along with Firedrake (e.g. 'check') + type: string + default: check + +runs: + using: composite + steps: + - name: Install system dependencies (1) + shell: bash + if: inputs.os == 'linux' + run: | + apt-get update + apt-get -y install curl git python3 python3-venv + + - name: Add NVIDIA CUDA deb repositories + if: inputs.gpu == 'cuda' + shell: bash + run: | + deburl=$(python3 ./firedrake-repo/scripts/firedrake-configure --show-extra-repo-pkg-url --gpu-arch cuda) + debfile=$( basename "${deburl}" ) + curl -fsSLO "${deburl}" + dpkg -i "${debfile}" + apt-get update + + - name: Install system dependencies (2) + shell: bash + run: | + if [ ${{ inputs.os }} = 'linux' ]; then + apt-get -y install \ + $(python3 ./firedrake-repo/scripts/firedrake-configure \ + --scalar-type ${{ inputs.scalar_type }} \ + --gpu-arch ${{ inputs.gpu }} \ + --show-system-packages) + : # Dependencies needed to run the test suite + apt-get -y install \ + fonts-dejavu \ + graphviz \ + graphviz-dev \ + parallel \ + poppler-utils + elif [ ${{ inputs.os }} = 'macos' ]; then + brew install \ + $(python3 ./firedrake-repo/scripts/firedrake-configure \ + --scalar-type ${{ inputs.scalar_type }} \ + --show-system-packages) + else + echo "Unrecognised 'os' input: '${{ inputs.os }}" + exit 1 + fi + + - name: Install PETSc + shell: bash + env: + EXTRA_OPTIONS: -use_gpu_aware_mpi 0 + run: | + : # Clone PETSc + if [ ${{ inputs.base_ref }} = 'main' ]; then + git clone --depth 1 https://gitlab.com/petsc/petsc.git + elif [ ${{ inputs.base_ref }} = 'release' ]; then + git clone --depth 1 \ + --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --show-petsc-version) \ + https://gitlab.com/petsc/petsc.git + else + echo "Unrecognised 'base_ref' input: '${{ inputs.base_ref }}" + exit 1 + fi + + if [ ${{ inputs.os }} = 'linux' ]; then + NPROCS=8 + elif [ ${{ inputs.os }} = 'macos' ]; then + NPROCS=4 + else + echo "Unrecognised 'os' input: '${{ inputs.os }}" + exit 1 + fi + + cd petsc + python3 ../firedrake-repo/scripts/firedrake-configure \ + --scalar-type ${{ inputs.scalar_type }} \ + --gpu-arch ${{ inputs.gpu }} \ + --show-petsc-configure-options | \ + xargs -L1 ./configure \ + --with-make-np=$NPROCS \ + --download-slepc + make + make check + + - name: Install Firedrake + shell: bash + run: | + export $(python3 ./firedrake-repo/scripts/firedrake-configure \ + --scalar-type ${{ inputs.scalar_type }} \ + --gpu-arch ${{ inputs.gpu }} \ + --show-env) + export SLEPC_DIR=$PETSC_DIR/$PETSC_ARCH + : # Disable compiler optimisation to speed up build (especially petsc4py) + export CFLAGS="-O0" + env + + python3 -m venv venv + . venv/bin/activate + + : # Empty the pip cache to ensure that everything is compiled from scratch + pip cache purge + + if [ ${{ inputs.base_ref }} = 'main' ]; then + : # Install build dependencies + pip install "$PETSC_DIR"/src/binding/petsc4py + pip install -r ./firedrake-repo/requirements-build.txt + + : # Install runtime dependencies that have been removed from the pyproject.toml + : # because they rely on non-PyPI versions of petsc4py. + pip install --no-build-isolation --no-deps \ + "$PETSC_DIR"/"$PETSC_ARCH"/externalpackages/git.slepc/src/binding/slepc4py + pip install --no-deps git+https://github.com/NGSolve/ngsPETSc.git netgen-mesher netgen-occt + + : # We have to pass '--no-build-isolation' to use a custom petsc4py + EXTRA_PIP_FLAGS='--no-build-isolation' + elif [ ${{ inputs.base_ref }} = 'release' ]; then + EXTRA_PIP_FLAGS='' + else + echo "Unrecognised 'base_ref' input: '${{ inputs.base_ref }}" + exit 1 + fi + + pip install --verbose $EXTRA_PIP_FLAGS \ + --no-binary h5py \ + --extra-index-url https://download.pytorch.org/whl/cpu \ + "./firedrake-repo[${{ inputs.deps }}]" + + firedrake-clean + pip list + + - name: Run firedrake-check + shell: bash + run: | + . venv/bin/activate + firedrake-check diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index 1391195e2c..67bde34979 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -11,10 +11,6 @@ on: description: The target branch (usually 'main' or 'release') type: string required: true - run_tests: - description: Whether to run the test suite - type: boolean - default: true build_docs: description: Whether to build the documentation type: boolean @@ -31,10 +27,6 @@ on: description: Whether to deploy the website type: boolean default: false - upload_pypi: - description: Whether to upload an sdist to PyPI - type: boolean - default: false workflow_dispatch: inputs: @@ -46,10 +38,6 @@ on: description: The target branch (usually 'main' or 'release') type: string required: true - run_tests: - description: Whether to run the test suite - type: boolean - default: true build_docs: description: Whether to build the documentation type: boolean @@ -66,10 +54,6 @@ on: description: Whether to deploy the website type: boolean default: false - upload_pypi: - description: Whether to upload an sdist to PyPI - type: boolean - default: false concurrency: # Cancel running jobs if new commits are pushed @@ -91,8 +75,6 @@ jobs: container: # TODO: set to 'ubuntu:latest' image: ubuntu:noble - outputs: - sdist_conclusion: ${{ steps.report_sdist.outputs.conclusion }} env: OMPI_ALLOW_RUN_AS_ROOT: 1 OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1 @@ -107,6 +89,7 @@ jobs: EXTRA_PYTEST_ARGS: --splitting-algorithm least_duration --timeout=600 --timeout-method=thread -o faulthandler_timeout=660 --durations-path=./firedrake-repo/tests/test_durations.json --durations=50 PYTEST_MPI_MAX_NPROCS: 8 steps: + # NOTE: is this needed any more? - name: Fix HOME # For unknown reasons GitHub actions overwrite HOME to /github/home # which will break everything unless fixed @@ -117,131 +100,24 @@ jobs: # Make sure the current directory is empty run: find . -delete - # Use a different mirror to fetch apt packages from to get around - # temporary outage. - # (https://askubuntu.com/questions/1549622/problem-with-archive-ubuntu-com-most-of-the-servers-are-not-responding) - # The mirror was chosen from https://launchpad.net/ubuntu/+archivemirrors. - - name: Configure apt - run: | - sed -i 's|http://archive.ubuntu.com/ubuntu|http://www.mirrorservice.org/sites/archive.ubuntu.com/ubuntu/|g' /etc/apt/sources.list.d/ubuntu.sources - apt-get update - - # Git is needed for actions/checkout and Python for firedrake-configure - - name: Install system dependencies (1) - run: apt-get -y install git python3 - - uses: actions/checkout@v5 with: path: firedrake-repo ref: ${{ inputs.source_ref }} - - name: Validate single source of truth - run: ./firedrake-repo/scripts/check-config - - # Check that the Dockerfile is using the latest Ubuntu version. - # The version is hardcoded into the Dockerfile so that the OS - # for each release is fixed. - - name: Check Dockerfile Ubuntu version - run: | - latest_version=$(grep "VERSION_ID=" /etc/os-release | cut -d '"' -f 2) - docker_version=$(grep FROM firedrake-repo/docker/Dockerfile.vanilla | cut -d ':' -f 2) - echo "Latest version: $latest_version" - echo "Docker version: $docker_version" - if [[ "$docker_version" != "$latest_version" ]]; then - echo "Ubuntu version ${docker_version} in Dockerfile is out of date with latest version ${latest_version}" - exit 1 - fi - - - name: Install system dependencies (2) - run: | - apt-get -y install \ - $(python3 ./firedrake-repo/scripts/firedrake-configure --arch ${{ matrix.arch }} --show-system-packages) - apt-get -y install python3-venv - : # Dependencies needed to run the test suite - apt-get -y install fonts-dejavu graphviz graphviz-dev parallel poppler-utils - - - name: Install PETSc - run: | - if [ ${{ inputs.target_branch }} = 'release' ]; then - git clone --depth 1 \ - --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --show-petsc-version) \ - https://gitlab.com/petsc/petsc.git - else - git clone --depth 1 https://gitlab.com/petsc/petsc.git - fi - cd petsc - python3 ../firedrake-repo/scripts/firedrake-configure \ - --arch ${{ matrix.arch }} --show-petsc-configure-options | \ - xargs -L1 ./configure --with-make-np=8 --download-slepc - make PETSC_DIR=/__w/firedrake/firedrake/petsc PETSC_ARCH=arch-firedrake-${{ matrix.arch }} - make check - { - echo "PETSC_DIR=/__w/firedrake/firedrake/petsc" - echo "PETSC_ARCH=arch-firedrake-${{ matrix.arch }}" - echo "SLEPC_DIR=/__w/firedrake/firedrake/petsc/arch-firedrake-${{ matrix.arch }}" - } >> "$GITHUB_ENV" - - name: Install Firedrake id: install - run: | - export $(python3 ./firedrake-repo/scripts/firedrake-configure --arch "${{ matrix.arch }}" --show-env) - python3 -m venv venv - . venv/bin/activate - - : # Empty the pip cache to ensure that everything is compiled from scratch - pip cache purge - - if [ ${{ inputs.target_branch }} = 'release' ]; then - EXTRA_BUILD_ARGS='' - EXTRA_PIP_FLAGS='' - else - : # Install build dependencies - pip install "$PETSC_DIR"/src/binding/petsc4py - pip install -r ./firedrake-repo/requirements-build.txt - - : # Install runtime dependencies that have been removed from the pyproject.toml - : # because they rely on non-PyPI versions of petsc4py. - pip install --no-build-isolation --no-deps \ - "$PETSC_DIR"/"$PETSC_ARCH"/externalpackages/git.slepc/src/binding/slepc4py - pip install --no-deps git+https://github.com/NGSolve/ngsPETSc.git netgen-mesher netgen-occt - - : # We have to pass '--no-build-isolation' to use a custom petsc4py - EXTRA_BUILD_ARGS='--no-isolation' - EXTRA_PIP_FLAGS='--no-build-isolation' - fi - - : # Install from an sdist so we can make sure that it is not ill-formed - pip install build - python -m build ./firedrake-repo --sdist $EXTRA_BUILD_ARGS - - pip install --verbose $EXTRA_PIP_FLAGS \ - --no-binary h5py \ - --extra-index-url https://download.pytorch.org/whl/cpu \ - "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]" - - firedrake-clean - pip list - - - name: Upload sdist (default ARCH only) - if: matrix.arch == 'default' - uses: actions/upload-artifact@v4 + uses: ./firedrake-repo/.github/actions/install with: - name: dist - path: firedrake-repo/dist/* - - - name: Report sdist build status - id: report_sdist - run: echo "conclusion=success" >> "$GITHUB_OUTPUT" - - - name: Run firedrake-check - run: | - . venv/bin/activate - firedrake-check - timeout-minutes: 5 + os: linux + source_ref: ${{ inputs.source_ref }} + base_ref: ${{ inputs.target_branch }} + scalar_type: ${{ matrix.arch }} + deps: ci - name: Run TSFC tests # Run even if earlier tests failed - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate : # Use pytest-xdist here so we can have a single collated output (not possible @@ -250,7 +126,7 @@ jobs: timeout-minutes: 10 - name: Run PyOP2 tests - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate : # Use pytest-xdist here so we can have a single collated output (not possible @@ -263,7 +139,7 @@ jobs: - name: Run Firedrake tests (nprocs = 1) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate : # Use pytest-xdist here so we can have a single collated output (not possible @@ -272,49 +148,49 @@ jobs: timeout-minutes: 90 - name: Run tests (nprocs = 2) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate firedrake-run-split-tests 2 4 "$EXTRA_PYTEST_ARGS" firedrake-repo/tests/firedrake timeout-minutes: 60 - name: Run tests (nprocs = 3) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate firedrake-run-split-tests 3 2 "$EXTRA_PYTEST_ARGS" firedrake-repo/tests/firedrake timeout-minutes: 60 - name: Run tests (nprocs = 4) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate firedrake-run-split-tests 4 2 "$EXTRA_PYTEST_ARGS" firedrake-repo/tests/firedrake timeout-minutes: 15 - name: Run tests (nprocs = 5) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate firedrake-run-split-tests 5 1 "$EXTRA_PYTEST_ARGS" firedrake-repo/tests/firedrake timeout-minutes: 15 - name: Run tests (nprocs = 6) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate firedrake-run-split-tests 6 1 "$EXTRA_PYTEST_ARGS" firedrake-repo/tests/firedrake timeout-minutes: 15 - name: Run tests (nprocs = 7) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate firedrake-run-split-tests 7 1 "$EXTRA_PYTEST_ARGS" firedrake-repo/tests/firedrake timeout-minutes: 15 - name: Run tests (nprocs = 8) - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' run: | . venv/bin/activate firedrake-run-split-tests 8 1 "$EXTRA_PYTEST_ARGS" firedrake-repo/tests/firedrake @@ -323,7 +199,6 @@ jobs: - name: Run Gusto smoke tests # Only test Gusto in real mode if: | - inputs.run_tests && (success() || steps.install.conclusion == 'success') && matrix.arch == 'default' run: | @@ -343,7 +218,6 @@ jobs: - name: Run Thetis smoke tests if: | - inputs.run_tests && (success() || steps.install.conclusion == 'success') && matrix.arch == 'default' run: | @@ -355,7 +229,6 @@ jobs: - name: Run spyro smoke tests if: | - inputs.run_tests && (success() || steps.install.conclusion == 'success') && matrix.arch == 'default' run: | @@ -367,7 +240,6 @@ jobs: - name: Run G-ADOPT smoke tests if: | - inputs.run_tests && (success() || steps.install.conclusion == 'success') && matrix.arch == 'default' run: | @@ -379,7 +251,7 @@ jobs: - name: Upload log files uses: actions/upload-artifact@v4 - if: inputs.run_tests && (success() || steps.install.conclusion == 'success') + if: success() || steps.install.conclusion == 'success' with: name: firedrake-logs-${{ matrix.arch }} path: pytest_*.log @@ -410,64 +282,14 @@ jobs: path: firedrake-repo ref: ${{ inputs.source_ref }} - - name: Install system dependencies - run: | - brew install $(python3 ./firedrake-repo/scripts/firedrake-configure --arch default --show-system-packages) - - - name: Install PETSc - run: | - if [ ${{ inputs.target_branch }} = 'release' ]; then - git clone --depth 1 \ - --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --show-petsc-version) \ - https://gitlab.com/petsc/petsc.git - else - git clone --depth 1 https://gitlab.com/petsc/petsc.git - fi - cd petsc - python3 ../firedrake-repo/scripts/firedrake-configure \ - --arch default --show-petsc-configure-options | \ - xargs -L1 ./configure --with-make-np=4 - make - make check - { - echo "PETSC_DIR=/Users/github/actions-runner/_work/firedrake/firedrake/petsc" - echo "PETSC_ARCH=arch-firedrake-default" - echo "SLEPC_DIR=/Users/github/actions-runner/_work/firedrake/firedrake/petsc/arch-firedrake-default" - } >> "$GITHUB_ENV" - - name: Install Firedrake id: install - run: | - export $(python3 ./firedrake-repo/scripts/firedrake-configure --arch default --show-env) - python3 -m venv venv - . venv/bin/activate - - : # Empty the pip cache to ensure that everything is compiled from scratch - pip cache purge - - if [ ${{ inputs.target_branch }} = 'release' ]; then - EXTRA_PIP_FLAGS='' - else - : # Install build dependencies - pip install "$PETSC_DIR"/src/binding/petsc4py - pip install -r ./firedrake-repo/requirements-build.txt - - : # We have to pass '--no-build-isolation' to use a custom petsc4py - EXTRA_PIP_FLAGS='--no-build-isolation' - fi - - pip install --verbose $EXTRA_PIP_FLAGS \ - --no-binary h5py \ - './firedrake-repo[check]' - - firedrake-clean - pip list - - - name: Run smoke tests - run: | - . venv/bin/activate - firedrake-check - timeout-minutes: 10 + uses: ./firedrake-repo/.github/actions/install + with: + os: macos + source_ref: ${{ inputs.source_ref }} + base_ref: ${{ inputs.target_branch }} + deps: check - name: Post-run cleanup if: always() @@ -510,15 +332,6 @@ jobs: # (https://github.com/actions/runner/issues/863) run: echo "HOME=/root" >> "$GITHUB_ENV" - - # Git is needed for actions/checkout and Python for firedrake-configure - # curl needed for adding new deb repositories to ubuntu - - name: Install system dependencies (1) - run: | - apt-get update - apt-get -y install git python3 curl - - - name: Pre-run cleanup # Make sure the current directory is empty run: find . -delete @@ -528,46 +341,8 @@ jobs: path: firedrake-repo ref: ${{ inputs.source_ref }} - - name: Add Nvidia CUDA deb repositories - run: | - deburl=$( python3 ./firedrake-repo/scripts/firedrake-configure --show-extra-repo-pkg-url --gpu-arch cuda ) - debfile=$( basename "${deburl}" ) - curl -fsSLO "${deburl}" - dpkg -i "${debfile}" - apt-get update - - - name: Install system dependencies (2) - run: | - apt-get -y install \ - $(python3 ./firedrake-repo/scripts/firedrake-configure --arch default --gpu-arch cuda --show-system-packages) - apt-get -y install python3-venv - : # Dependencies needed to run the test suite - apt-get -y install fonts-dejavu graphviz graphviz-dev parallel poppler-utils - - - name: Install PETSc - env: - EXTRA_OPTIONS: -use_gpu_aware_mpi 0 - run: | - if [ ${{ inputs.target_branch }} = 'release' ]; then - git clone --depth 1 \ - --branch $(python3 ./firedrake-repo/scripts/firedrake-configure --gpu-arch cuda --show-petsc-version) \ - https://gitlab.com/petsc/petsc.git - else - git clone --depth 1 https://gitlab.com/petsc/petsc.git - fi - cd petsc - python3 ../firedrake-repo/scripts/firedrake-configure \ - --arch default --gpu-arch cuda --show-petsc-configure-options | \ - xargs -L1 ./configure --with-make-np=4 - make - make check - { - echo "PETSC_DIR=/__w/firedrake/firedrake/petsc" - echo "PETSC_ARCH=arch-firedrake-default-cuda" - echo "SLEPC_DIR=/__w/firedrake/firedrake/petsc/arch-firedrake-default-cuda" - } >> "$GITHUB_ENV" - - name: Install Firedrake +<<<<<<< HEAD id: install run: | export $(python3 ./firedrake-repo/scripts/firedrake-configure --arch default --gpu-arch cuda --show-env) @@ -601,6 +376,15 @@ jobs: . venv/bin/activate firedrake-check timeout-minutes: 10 +======= + uses: ./firedrake-repo/.github/actions/install + with: + os: linux + source_ref: ${{ inputs.source_ref }} + base_ref: ${{ inputs.target_branch }} + gpu: cuda + deps: check +>>>>>>> release - name: Confirm GPU usage in OffloadPC test shell: bash @@ -630,6 +414,39 @@ jobs: run: | find . -delete + check: + name: Basic repository checks + runs-on: ubuntu-latest + container: + image: ubuntu:latest + steps: + - name: Install system dependencies + run: | + apt-get update + apt-get -y install python3 + + - uses: actions/checkout@v5 + with: + path: firedrake-repo + ref: ${{ inputs.source_ref }} + + - name: Validate single source of truth + run: ./firedrake-repo/scripts/check-config + + # Check that the Dockerfile is using the latest Ubuntu version. + # The version is hardcoded into the Dockerfile so that the OS + # for each release is fixed. + - name: Check Dockerfile Ubuntu version + run: | + latest_version=$(grep "VERSION_ID=" /etc/os-release | cut -d '"' -f 2) + docker_version=$(grep FROM firedrake-repo/docker/Dockerfile.vanilla | cut -d ':' -f 2) + echo "Latest version: $latest_version" + echo "Docker version: $docker_version" + if [[ "$docker_version" != "$latest_version" ]]; then + echo "Ubuntu version ${docker_version} in Dockerfile is out of date with latest version ${latest_version}" + exit 1 + fi + lint: name: Lint codebase runs-on: ubuntu-latest @@ -746,24 +563,3 @@ jobs: steps: - name: Deploy to GitHub Pages uses: actions/deploy-pages@v4 - - upload_pypi: - name: Upload to PyPI (optional) - needs: test_linux - if: | - always() && - inputs.upload_pypi && - inputs.target_branch == 'release' && - needs.test_linux.outputs.sdist_conclusion == 'success' - runs-on: ubuntu-latest - environment: - name: pypi - steps: - - uses: actions/download-artifact@v4 - with: - name: dist - path: dist - - name: Push to PyPI - uses: pypa/gh-action-pypi-publish@release/v1 - with: - password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/.github/workflows/pr.yml b/.github/workflows/pr.yml index 6b63f97cae..4ea32ccce6 100644 --- a/.github/workflows/pr.yml +++ b/.github/workflows/pr.yml @@ -2,7 +2,7 @@ name: Test pull request on: pull_request: - types: [ opened, synchronize, reopened, edited ] + types: [ opened, synchronize, reopened, labeled ] jobs: test: diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 3f2f11aa4d..9938ca54d9 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -45,15 +45,41 @@ jobs: fi pypi: - uses: ./.github/workflows/core.yml + name: Make PyPI release needs: check - with: - source_ref: ${{ inputs.branch }} - target_branch: release - run_tests: false - build_docs: false - upload_pypi: true - secrets: inherit + runs-on: [self-hosted, Linux] + container: + # TODO: set to 'ubuntu:latest' + image: ubuntu:noble + environment: + name: pypi + steps: + - uses: actions/checkout@v5 + with: + path: firedrake-repo + ref: ${{ inputs.branch }} + + - name: Install Firedrake + uses: ./firedrake-repo/.github/actions/install + with: + os: linux + source_ref: ${{ inputs.branch }} + base_ref: release + + - name: Create sdist + run: | + . venv/bin/activate + export $(python3 ./firedrake-repo/scripts/firedrake-configure --show-env) + env + pip install build + python -m build ./firedrake-repo --sdist + + - name: Upload to PyPI + uses: pypa/gh-action-pypi-publish@release/v1 + with: + packages-dir: firedrake-repo/dist/ + # TODO: Use trusted publishing + password: ${{ secrets.PYPI_API_TOKEN }} github_release: name: Create GitHub release diff --git a/AUTHORS.rst b/AUTHORS.rst index 67d18a36c9..8141939d59 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -108,7 +108,7 @@ Stephan C. Kramer............. Tuomas Kärnä -Michael Lange................. +Michael Lange Nicolas Loriant diff --git a/demos/demo_references.bib b/demos/demo_references.bib index 5792042258..9c84a461a5 100644 --- a/demos/demo_references.bib +++ b/demos/demo_references.bib @@ -1,3 +1,13 @@ +@Article{Egger:2018, + author = {Egger, Herbert and Fellner, Klemens and Pietschmann, Jan-Frederik and Tang, Bao Quoc}, + title = {{Analysis and numerical solution of coupled volume-surface reaction-diffusion systems with application to cell biology}}, + journal = {Applied Mathematics and Computation}, + year = {2018}, + volume = {336}, + pages = {351--367}, + doi = {10.1016/j.amc.2018.04.031}, +} + @Article{Shu:1988, author = {Shu, Chi-Wang and Osher, Stanley}, title = {{Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes}}, diff --git a/demos/submesh_reaction_diffusion/submesh_reaction_diffusion.py.rst b/demos/submesh_reaction_diffusion/submesh_reaction_diffusion.py.rst new file mode 100644 index 0000000000..87813bc032 --- /dev/null +++ b/demos/submesh_reaction_diffusion/submesh_reaction_diffusion.py.rst @@ -0,0 +1,263 @@ +Coupled volume-surface reaction-diffusion on a torus with submesh +================================================================= + +In many biological and physical processes, chemical species diffuse through +a volume and exchange with a surface species, coupled by an interfacial +transfer mechanism. A prototypical model of this class is the coupled +volume-surface reaction-diffusion system analysed by :cite:`Egger:2018`, which +we consider here on a solid torus :math:`\Omega \subset \mathbb{R}^3` with +surface :math:`\Gamma = \partial\Omega`. This demo illustrates how such +problems may be solved with the :func:`~.Submesh` functionality in Firedrake. + +The model +--------- + +Let :math:`L : \Omega \times (0, T] \to \mathbb{R}` denote the volume +concentration and :math:`\ell : \Gamma \times (0, T] \to \mathbb{R}` the +surface concentration. The evolution is governed by + +.. math:: + + \partial_t L - d_L \Delta L = 0 + \qquad \text{in } \Omega, + + \partial_t \ell - d_\ell \Delta_\Gamma \ell = \lambda L - \gamma \ell + \qquad \text{on } \Gamma, + +with the interface condition + +.. math:: + + d_L \,\partial_n L = \gamma \ell - \lambda L + \qquad \text{on } \Gamma. + +Here :math:`d_L` and :math:`d_\ell` are diffusion coefficients, while +:math:`\lambda` and :math:`\gamma` are positive transfer constants. The +interface condition states that the normal flux of :math:`L` out of the volume +equals the net transfer from surface to volume, so that mass is conserved globally: + +.. math:: + + M(t) = \int_\Omega L \,\mathrm{d}x + \int_\Gamma \ell \,\mathrm{d}s + = M(0) \qquad \text{for all } t > 0. + +Weak formulation +---------------- + +Multiplying the volume equation by :math:`v \in H^1(\Omega)`, integrating over +:math:`\Omega`, and using the interface condition in the boundary term gives + +.. math:: + + (\partial_t L, v)_\Omega + + d_L (\nabla L, \nabla v)_\Omega + + (\lambda L - \gamma \ell,\, v)_\Gamma = 0. + +Multiplying the surface equation by :math:`w \in H^1(\Gamma)` and integrating +over :math:`\Gamma` gives + +.. math:: + + (\partial_t \ell, w)_\Gamma + + d_\ell (\nabla_\Gamma \ell, \nabla_\Gamma w)_\Gamma + - (\lambda L - \gamma \ell,\, w)_\Gamma = 0. + +The coupling terms :math:`(\lambda L - \gamma \ell, v)_\Gamma` and +:math:`(\lambda L - \gamma \ell, w)_\Gamma` involve both the volume function +:math:`L` (restricted to :math:`\Gamma`) and the surface function :math:`\ell`, +integrated over the same surface. In Firedrake this is handled by a +*cross-mesh measure* on the submesh. + +Implementation +-------------- + +We begin by importing Firedrake and ngsPETSc to create the torus mesh. +:: + + from firedrake import * + try: + import netgen + except ImportError: + import sys + warning("Unable to import Netgen.") + sys.exit(0) + +Building the torus mesh with ngsPETSc +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +We construct a solid torus using Open CASCADE Technology via Netgen. The +generating circle has major radius :math:`R = 3` and minor radius :math:`r = 1` +and is swept around the :math:`z`-axis. :: + + from netgen.occ import * + + n_pts = 80 + def _curve(t): + return Pnt(0, 3 + cos(t), sin(t)) + + pnts = [_curve(2 * pi * k / n_pts) for k in range(n_pts + 1)] + spline = SplineApproximation(pnts) + face = Face(Wire(spline)) + torus = face.Revolve(Axis((0, 0, 0), Z), 360) + + ngmesh = OCCGeometry(torus).GenerateMesh(maxh=0.5) + base_v = Mesh(ngmesh) + +Mesh hierarchy and submesh hierarchy +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The submesh for the surface can be built from +either the :func:`~.Submesh` or :func:`~.SubmeshHierarchy` constructors, +but only the latter enables us to use a multigrid solver. +Geometric multigrid requires a hierarchy of uniformly refined meshes. We +build the volume and surface hierarchies together. + +Like :class:`~.DirichletBC`, :func:`~.Submesh` +takes in a subdomain id to indicate which part of the mesh should be +extracted. In this case we want the entire exterior facet mesh, which we +can specify directly. :: + + nref = 1 + mh_v = MeshHierarchy(base_v, nref) + mh_s = SubmeshHierarchy(mh_v, subdomain_id="on_boundary") + mesh_v = mh_v[-1] + mesh_s = mh_s[-1] + +Function spaces +~~~~~~~~~~~~~~~ + +We use continuous piecewise-linear elements on both the volume and the surface, +collected into a :func:`~.MixedFunctionSpace`. :: + + V_v = FunctionSpace(mesh_v, "CG", 1) + V_s = FunctionSpace(mesh_s, "CG", 1) + Z = V_v * V_s + +Initial data and time-stepping setup +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +We initialise with smooth functions on the torus. :: + + dt = Constant(0.1) + T = 1 + + z = Function(Z) # solution at t^{n+1} + z_old = Function(Z) # solution at t^n + + L, l = split(z) + L_old, l_old = split(z_old) + v, w = split(TestFunction(Z)) + + + X_v = SpatialCoordinate(mesh_v) + X_s = SpatialCoordinate(mesh_s) + + z_old.subfunctions[0].interpolate(2 + sin(X_v[0]) * cos(X_v[1])) + z_old.subfunctions[1].interpolate(2 + cos(X_s[2])) + z.assign(z_old) + +The model parameters are kept simple. :: + + d_L = Constant(1) # volume diffusion + d_ell = Constant(1) # surface diffusion + lam = Constant(1) # transfer rate: volume → surface + gam = Constant(1) # transfer rate: surface → volume + +Integration measures +~~~~~~~~~~~~~~~~~~~~~ + +Three measures are needed. ``dV`` integrates over :math:`\Omega`, ``dA`` over +:math:`\Gamma`, and ``dC`` is the *cross-mesh measure* that integrates over the +submesh but also queries degrees of freedom from the parent volume mesh. +Firedrake computes the intersection automatically with ``intersect_measures``. :: + + dV = dx(mesh_v) + dA = dx(mesh_s) + dC = Measure("dx", domain=mesh_s, intersect_measures=[ds(mesh_v)]) + +The variational form +~~~~~~~~~~~~~~~~~~~~~ + +We discretise in time with the L-stable backward Euler scheme. The cross-mesh +coupling terms use ``dC``: the first argument to ``inner`` may come from the +volume space ``V_v`` or the surface space ``V_s``, and the test functions +``v`` and ``w`` live on their respective spaces. :: + + transfer = lam * L - gam * l + + F = ( + inner((L - L_old) / dt, v) * dV + + d_L * inner(grad(L), grad(v)) * dV + + inner(transfer, v) * dC + + inner((l - l_old) / dt, w) * dA + + d_ell * inner(grad(l), grad(w)) * dA + - inner(transfer, w) * dC + ) + +Fieldsplit geometric-multigrid solver +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The system is not symmetric because :math:`\lambda \neq \gamma` in general, so +we use GMRES as the outer solver. (It could be made so by rescaling, but we do not +pursue this here.) The two diagonal blocks—a volume elliptic +problem and a surface elliptic problem—are each preconditioned independently +by a full-cycle geometric multigrid solver. (Although monolithic multigrid approaches on +the coupled system are also available.) Firedrake automatically +rediscretises the operators on each level using the mesh hierarchies we built +above. :: + + gmg_block = { + "ksp_type": "preonly", + "pc_type": "mg", + "pc_mg_type": "full", + "mg_levels_ksp_type": "chebyshev", + "mg_levels_pc_type": "jacobi", + "mg_levels_ksp_max_it": 2, + } + + sp = { + "snes_type": "ksponly", + "ksp_type": "gmres", + "ksp_monitor": None, + "ksp_rtol": 1e-8, + "pc_type": "fieldsplit", + "pc_fieldsplit_type": "additive", + "fieldsplit_0": gmg_block, + "fieldsplit_1": gmg_block, + } + +We create the problem and solver once, then step in time. :: + + problem = NonlinearVariationalProblem(F, z) + solver = NonlinearVariationalSolver(problem, solver_parameters=sp) + +Time loop +~~~~~~~~~ + +We advance the solution and write VTK output at each step. :: + + L_out, l_out = z_old.subfunctions + L_out.rename("Volume") + l_out.rename("Surface") + + pvd_v = VTKFile("output/volume.pvd") + pvd_s = VTKFile("output/surface.pvd") + + t = 0.0 + pvd_v.write(L_out, time=t) + pvd_s.write(l_out, time=t) + + while t < T - 1e-8: + solver.solve() + z_old.assign(z) + t += float(dt) + pvd_v.write(L_out, time=t) + pvd_s.write(l_out, time=t) + +A python script version of this demo can be found :demo:`here +`. + +.. rubric:: References + +.. bibliography:: demo_references.bib + :filter: docname in docnames diff --git a/docs/source/advanced_tut.rst b/docs/source/advanced_tut.rst index 5007d12f4b..38d9a3d849 100644 --- a/docs/source/advanced_tut.rst +++ b/docs/source/advanced_tut.rst @@ -35,3 +35,4 @@ element systems. Steady Boussinesq problem with integral constraints. Steady multicomponent flow -- microfluidic mixing of hydrocarbons. Deflation techniques for computing multiple solutions of nonlinear problems. + Coupled volume-surface reaction-diffusion on a torus using Submesh and geometric multigrid. diff --git a/docs/source/firedrake_19.rst b/docs/source/firedrake_19.rst index 2adf060d3a..b67ac14cf5 100644 --- a/docs/source/firedrake_19.rst +++ b/docs/source/firedrake_19.rst @@ -68,8 +68,8 @@ presenters. Workshop dinner ~~~~~~~~~~~~~~~ -There will be a workshop dinner in the evening of 26 September at `ASK -Italian `_ on +There will be a workshop dinner in the evening of 26 September at ASK +Italian on Durham's Millenium Square (easily accessible by foot from the conference location and town centre). diff --git a/docs/source/firedrake_26.rst b/docs/source/firedrake_26.rst index 1030179fa2..67cdb1e6f0 100644 --- a/docs/source/firedrake_26.rst +++ b/docs/source/firedrake_26.rst @@ -25,7 +25,7 @@ The PETSc meeting will run from 1-3 June, and the Firedrake meeting from 3-5 June, with overlapping sessions of common interest on Wednesday 3 June. Participants may register for either part of the week, but we hope that many will take the opportunity to learn more about these two related projects by -staying for the whole week. +staying for the whole week. Venue & registration -------------------- @@ -41,13 +41,13 @@ night for PETSc, Wednesday and Thursday night for Firedrake), 3 meals per day coffee breaks. There are also options without the room and breakfast component for those who live locally. -================= =================================== ====== ======= ======================= ========================== +================= =================================== ====== ======= ======================= ========================== Event Dates Full Student Full (no accommodation) Student (no accommodation) -================= =================================== ====== ======= ======================= ========================== -PETSc & Firedrake 1 June (10:00) - 5 June (16:00) £1000 £400 £420 £170 -PETSc only 1 June (10:00) - 3 June (lunchtime) £550 £220 £250 £100 -Firedrake only 3 June (lunchtime) - 5 June (16:00) £550 £220 £250 £100 -================= =================================== ====== ======= ======================= ========================== +================= =================================== ====== ======= ======================= ========================== +PETSc & Firedrake 1 June (10:00) - 5 June (16:00) £1000 £400 £420 £170 +PETSc only 1 June (10:00) - 3 June (lunchtime) £550 £220 £250 £100 +Firedrake only 3 June (lunchtime) - 5 June (16:00) £550 £220 £250 £100 +================= =================================== ====== ======= ======================= ========================== International attendees will likely need to book a hotel for the night before the conference. We recommend either booking an extra night at the `conference @@ -105,7 +105,7 @@ exact): 13:30-14:45 Firedrake tutorial 1 14:45-15:30 coffee break 15:30-16:45 Firedrake tutorial 2 - 16:45-17:30 poster session + 16:45-17:30 poster session and drinks reception, sponsored by `Devito Codes `__ =========== ==================== **Thursday** @@ -159,7 +159,7 @@ By air ~~~~~~ Oxford lies between London and Birmingham. The closest London airport is London -Heathrow. +Heathrow. From Heathrow, a `regular coach service `__ runs to Oxford. There @@ -189,7 +189,7 @@ Support These workshops are brought to you by `CCP-DCM `__, `CoSeC `__ and `STFC -`__ +`__ .. |STFC| image:: /images/stfc.* :height: 60px @@ -207,6 +207,20 @@ These workshops are brought to you by `CCP-DCM `__, +--------+---------+ +The poster session and drinks reception is sponsored by +`Devito Codes `__. + +.. |DevitoLogo| image:: /images/devito.svg + :height: 60px + :target: https://www.devitocodes.com + +.. only:: html + + +--------------+ + | |DevitoLogo| | + +--------------+ + + Further details --------------- diff --git a/docs/source/images/devito.svg b/docs/source/images/devito.svg new file mode 100644 index 0000000000..d19958e959 --- /dev/null +++ b/docs/source/images/devito.svg @@ -0,0 +1,28 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/team.ini b/docs/source/team.ini index 994cb27bcc..ce3e17bc1e 100644 --- a/docs/source/team.ini +++ b/docs/source/team.ini @@ -62,7 +62,7 @@ Tianjiao Sun: Andrew T. T. McRae: Fabio Luporini: https://www.imperial.ac.uk/people/f.luporini12 Alastair Gregory: -Michael Lange: https://www.linkedin.com/in/michael-lange-56675994/ +Michael Lange: Simon W. Funke: Florian Rathgeber: https://kynan.github.io Gheorghe-Teodor Bercea: diff --git a/firedrake/__init__.py b/firedrake/__init__.py index 1c18b04cd5..5652943c6f 100644 --- a/firedrake/__init__.py +++ b/firedrake/__init__.py @@ -99,7 +99,7 @@ def init_petsc(): ) from firedrake.mg import ( # noqa: F401 HierarchyBase, MeshHierarchy, ExtrudedMeshHierarchy, - NonNestedHierarchy, SemiCoarsenedExtrudedHierarchy, + NonNestedHierarchy, SemiCoarsenedExtrudedHierarchy, SubmeshHierarchy, prolong, restrict, inject, TransferManager, OpenCascadeMeshHierarchy, AdaptiveMeshHierarchy, AdaptiveTransferManager diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index a4d4d92460..785f91c51b 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -13,6 +13,8 @@ from firedrake.utils import IntType, ScalarType from libc.string cimport memset from libc.stdlib cimport qsort from finat.element_factory import as_fiat_cell +from numbers import Integral +from collections.abc import Sequence cimport numpy as np cimport mpi4py.MPI as MPI @@ -2203,11 +2205,11 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): def _get_periodicity(dm: PETSc.DM) -> tuple[tuple[bool, bool], ...]: """Return mesh periodicity information. - + This function returns a 2-tuple of bools per dimension where the first entry indicates whether the mesh is periodic in that dimension, and the second indicates whether the mesh is single-cell periodic in that dimension. - + """ cdef: const PetscReal *maxCell, *L @@ -3971,7 +3973,7 @@ def create_halo_exchange_sf(PETSc.DM dm): def submesh_create(PETSc.DM dm, PetscInt subdim, label_name, - PetscInt label_value, + subdomain_id, PetscBool ignore_label_halo, comm=None): """Create submesh. @@ -3984,8 +3986,8 @@ def submesh_create(PETSc.DM dm, Topological dimension of the submesh label_name : str Name of the label - label_value : int - Value in the label + subdomain_id : int | Sequence + Values in the label ignore_label_halo : bool If labeled points in the halo are ignored. comm : PETSc.Comm | None @@ -3993,20 +3995,36 @@ def submesh_create(PETSc.DM dm, """ cdef: + PETSc.IS points, subpoints PETSc.DMLabel label, temp_label char *temp_label_name = "firedrake_submesh_temp_label" - PetscInt pStart, pEnd, p, i, stratum_size - PETSc.PetscIS stratum_is = NULL + PetscInt pStart, pEnd, p, i, stratum_size = 0, label_value = 1 const PetscInt *stratum_indices = NULL + # Cast subdomain_id into an iterable + if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): + subdomain_id = (subdomain_id,) + # Take the union of the all the label values label = dm.getLabel(label_name) + points = PETSc.IS() + for sub in subdomain_id: + if isinstance(sub, Integral): + subpoints = label.getStratumIS(sub) + elif sub == "on_boundary": + subpoints = dm.getStratumIS("exterior_facets", 1) + else: + raise ValueError(f"Submesh construction got invalid subdomain_id {sub}.") + if points: + points = points.union(subpoints) + else: + points = subpoints # Create temp_label that contains no lower-dimensional points. dm.createLabel(temp_label_name) temp_label = dm.getLabel(temp_label_name) - CHKERR(DMLabelGetStratumSize(label.dmlabel, label_value, &stratum_size)) + if points: + CHKERR(ISGetSize(points.iset, &stratum_size)) if stratum_size > 0: - CHKERR(DMLabelGetStratumIS(label.dmlabel, label_value, &stratum_is)) - CHKERR(ISGetIndices(stratum_is, &stratum_indices)) + CHKERR(ISGetIndices(points.iset, &stratum_indices)) CHKERR(DMPlexGetDepthStratum(dm.dm, subdim, &pStart, &pEnd)) for i in range(stratum_size): p = stratum_indices[i] @@ -4014,8 +4032,7 @@ def submesh_create(PETSc.DM dm, # culling all lower-dimensional points. if pStart <= p < pEnd: CHKERR(DMLabelSetValue(temp_label.dmlabel, p, label_value)) - CHKERR(ISRestoreIndices(stratum_is, &stratum_indices)) - CHKERR(ISDestroy(&stratum_is)) + CHKERR(ISRestoreIndices(points.iset, &stratum_indices)) # Make submesh using temp_label. subdm, ownership_transfer_sf = dm.filter(label=temp_label, value=label_value, diff --git a/firedrake/dmhooks.py b/firedrake/dmhooks.py index 5b1562e84c..a2b60a1f50 100644 --- a/firedrake/dmhooks.py +++ b/firedrake/dmhooks.py @@ -417,9 +417,8 @@ def coarsen(dm, comm): """ from firedrake.mg.utils import get_level V = get_function_space(dm) - # TODO: Think harder. - m, = set(m_ for m_ in V.mesh()) - hierarchy, level = get_level(m) + mesh = V.mesh() + hierarchy, level = get_level(mesh) if level < 1: raise RuntimeError("Cannot coarsen coarsest DM") coarsen = get_ctx_coarsener(dm) diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py index b5b2e90a0a..36c2652615 100644 --- a/firedrake/functionspaceimpl.py +++ b/firedrake/functionspaceimpl.py @@ -926,7 +926,7 @@ def local_to_global_map(self, bcs, lgmap=None, mat_type=None): return PETSc.LGMap().create(indices, bsize=bsize, comm=lgmap.comm) def collapse(self): - return type(self)(self.mesh(), self.ufl_element()) + return type(self)(self.mesh(), self.ufl_element(), name=self.name) class RestrictedFunctionSpace(FunctionSpace): diff --git a/firedrake/logging.py b/firedrake/logging.py index 1841f7a4e4..7a773fde33 100644 --- a/firedrake/logging.py +++ b/firedrake/logging.py @@ -72,6 +72,9 @@ def set_log_handlers(handlers=None, comm=COMM_WORLD): for package in packages: logger = logging.getLogger(package) + # Avoids duplicating messages if someone directly invokes logging.info: + logger.propagate = False + for handler in logger.handlers: logger.removeHandler(handler) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 64c929c210..b6aa2e3f15 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3917,16 +3917,19 @@ def _pic_swarm_in_mesh( base_parent_cell_nums, extrusion_heights = _parent_extrusion_numbering( parent_cell_nums_local, parent_mesh.layers ) - # mesh.topology.cell_closure[:, -1] maps Firedrake cell numbers to plex - # numbers. - plex_parent_cell_nums = parent_mesh.topology.cell_closure[ - base_parent_cell_nums, -1 + # cell_closure[:, -1] maps Firedrake cell numbers to plex numbers. + # Index only visible rows: -1 sentinels crash on empty-rank arrays. + plex_parent_cell_nums = np.full_like(base_parent_cell_nums, -1) + plex_parent_cell_nums[visible_idxs] = parent_mesh.topology.cell_closure[ + base_parent_cell_nums[visible_idxs], -1 ] base_parent_cell_nums_visible = base_parent_cell_nums[visible_idxs] extrusion_heights_visible = extrusion_heights[visible_idxs] else: - plex_parent_cell_nums = parent_mesh.topology.cell_closure[ - parent_cell_nums_local, -1 + # Index only visible rows: -1 sentinels crash on empty-rank arrays. + plex_parent_cell_nums = np.full_like(parent_cell_nums_local, -1) + plex_parent_cell_nums[visible_idxs] = parent_mesh.topology.cell_closure[ + parent_cell_nums_local[visible_idxs], -1 ] base_parent_cell_nums_visible = None extrusion_heights_visible = None @@ -4869,10 +4872,11 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig subdim : int | None Topological dimension of the submesh. Defaults to ``mesh.topological_dimension``. - subdomain_id : int | None + subdomain_id : int | Sequence | None Subdomain ID representing the submesh. - If `None` the submesh will cover the entire domain. - This is useful to obtain a codim-1 submesh over all facets or + If multiple subdomain IDs are provided, their union is taken. + If `None` the submesh will cover the entire domain, + this is useful to obtain a codim-1 submesh over all facets or a submesh over a different communicator. label_name : str | None Name of the label to search ``subdomain_id`` in. @@ -4915,6 +4919,47 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig ridges to be contained in the quad mesh are shared by at most two facets to make the quad mesh orientation algorithm work. + Examples + -------- + >>> mesh = UnitSquareMesh(4, 4) + >>> x, y = SpatialCoordinate(mesh) + >>> DG = FunctionSpace(mesh, "DG", 0) + >>> DGT = FunctionSpace(mesh, "DGT", 0) + + Mark a cell subdomain and construct a codim-0 submesh from all cells in the subdomain + + >>> cell_marker = assemble(interpolate(conditional(lt(x, 0.5), 1, 0), DG)) + >>> mesh.mark_entities(cell_marker, 111) + >>> submesh = Submesh(mesh, subdomain_id=111) + + Mark a facet subdomain and construct a codim-1 submesh from all facets in the subdomain + + >>> facet_marker = assemble(interpolate(conditional(lt(abs(x-0.5), 1E-12), 1, 0), DGT)) + >>> mesh.mark_entities(facet_marker, 222) + >>> submesh = Submesh(mesh, subdim=mesh.topological_dimension-1, subdomain_id=222) + + Construct a codim-0 submesh of the union of multiple subdomains by passing a list + + >>> mesh.mark_entities(assemble(interpolate(conditional(lt(x, 0.5), 1, 0), DG)), 1) + >>> mesh.mark_entities(assemble(interpolate(conditional(lt(y, 0.5), 1, 0), DG)), 2) + >>> submesh = Submesh(mesh, subdomain_id=[1, 2]) + + Construct a codim-1 submesh of all the facets (the skeleton mesh) + + >>> submesh = Submesh(mesh, subdim=1) + + Construct a codim-1 submesh of the entire boundary + + >>> submesh = Submesh(mesh, subdomain_id="on_boundary") + + Construct a codim-1 submesh of the union of multiple boundaries + + >>> submesh = Submesh(mesh, subdim=mesh.topological_dimension-1, subdomain_id=[1, 2, 3]) + + Construct a codim-0 submesh of the part of the mesh owned by each MPI rank + + >>> submesh = Submesh(mesh, ignore_halo=True, comm=COMM_SELF) + """ if not isinstance(mesh, MeshGeometry): raise TypeError("Parent mesh must be a `MeshGeometry`") @@ -4922,6 +4967,18 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig raise NotImplementedError("Can not create a submesh of an ``ExtrudedMesh``") elif isinstance(mesh.topology, VertexOnlyMeshTopology): raise NotImplementedError("Can not create a submesh of a ``VertexOnlyMesh``") + + if subdomain_id == "on_boundary": + if subdim is None: + subdim = mesh.topological_dimension - 1 + elif subdim != mesh.topological_dimension - 1: + raise ValueError('subdomain_id="on_boundary" requires subdim=dim-1') + if label_name is None: + label_name = "exterior_facets" + elif label_name != "exterior_facets": + raise ValueError('subdomain_id="on_boundary" requires label_name="exterior_facets"') + subdomain_id = 1 + if subdim is None: subdim = mesh.topological_dimension plex = mesh.topology_dm @@ -4947,7 +5004,7 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig if subplex.getDimension() != subdim: raise RuntimeError(f"Found subplex dim ({subplex.getDimension()}) != expected ({subdim})") if reorder is None: - # Ideally we should set perm_is = mesh.dm_reordering[label_indices] + # Ideally we should set perm_is = mesh._dm_renumbering[label_indices] reorder = mesh._did_reordering submesh = Mesh( diff --git a/firedrake/mg/__init__.py b/firedrake/mg/__init__.py index 12f8ec7c36..c73e5c7849 100644 --- a/firedrake/mg/__init__.py +++ b/firedrake/mg/__init__.py @@ -1,6 +1,7 @@ from firedrake.mg.mesh import ( # noqa F401 HierarchyBase, MeshHierarchy, ExtrudedMeshHierarchy, - NonNestedHierarchy, SemiCoarsenedExtrudedHierarchy + NonNestedHierarchy, SemiCoarsenedExtrudedHierarchy, + SubmeshHierarchy, ) from firedrake.mg.interface import ( # noqa F401 prolong, restrict, inject diff --git a/firedrake/mg/embedded.py b/firedrake/mg/embedded.py index 85f3b719d1..0ff46b67c0 100644 --- a/firedrake/mg/embedded.py +++ b/firedrake/mg/embedded.py @@ -246,16 +246,22 @@ def op(self, source, target, transfer_op): Vt = target.function_space() source_element = Vs.ufl_element() target_element = Vt.ufl_element() + + # Recurse on sub-elements before any cache lookup: the mixed element may + # span multiple cell types (e.g. volume + surface submesh), which would + # cause cache() to call get_embedding_dg_element on a multi-cell element. + if type(source_element) is finat.ufl.MixedElement: + assert type(target_element) is finat.ufl.MixedElement + for source_, target_ in zip(source.subfunctions, target.subfunctions): + self.op(source_, target_, transfer_op=transfer_op) + return + if not self.requires_transfer(Vs, transfer_op, source, target): return gdim = Vt.mesh().geometric_dimension if self.is_native(target_element, gdim, transfer_op): self._native_transfer(target_element, gdim, transfer_op)(source, target) - elif type(source_element) is finat.ufl.MixedElement: - assert type(target_element) is finat.ufl.MixedElement - for source_, target_ in zip(source.subfunctions, target.subfunctions): - self.op(source_, target_, transfer_op=transfer_op) else: # Get some work vectors dgsource = self.DG_work(Vs) @@ -312,16 +318,22 @@ def restrict(self, source, target): Vt_star = target.function_space() source_element = Vs_star.ufl_element() target_element = Vt_star.ufl_element() + + # Recurse on sub-elements before any cache lookup: the mixed element may + # span multiple cell types (e.g. volume + surface submesh), which would + # cause cache() to call get_embedding_dg_element on a multi-cell element. + if type(source_element) is finat.ufl.MixedElement: + assert type(target_element) is finat.ufl.MixedElement + for source_, target_ in zip(source.subfunctions, target.subfunctions): + self.restrict(source_, target_) + return + if not self.requires_transfer(Vs_star, Op.RESTRICT, source, target): return gdim = Vs_star.mesh().geometric_dimension if self.is_native(source_element, gdim, Op.RESTRICT): self._native_transfer(source_element, gdim, Op.RESTRICT)(source, target) - elif type(source_element) is finat.ufl.MixedElement: - assert type(target_element) is finat.ufl.MixedElement - for source_, target_ in zip(source.subfunctions, target.subfunctions): - self.restrict(source_, target_) else: Vs = Vs_star.dual() Vt = Vt_star.dual() diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py index 2279cbc67d..fecbae149e 100644 --- a/firedrake/mg/mesh.py +++ b/firedrake/mg/mesh.py @@ -1,6 +1,7 @@ import numpy as np from fractions import Fraction from collections import defaultdict +from collections.abc import Sequence from pyop2.datatypes import IntType @@ -14,7 +15,7 @@ from .utils import set_level __all__ = ("HierarchyBase", "MeshHierarchy", "ExtrudedMeshHierarchy", "NonNestedHierarchy", - "SemiCoarsenedExtrudedHierarchy") + "SemiCoarsenedExtrudedHierarchy", "SubmeshHierarchy") class HierarchyBase(object): @@ -310,3 +311,79 @@ def SemiCoarsenedExtrudedHierarchy(base_mesh, height, nref=1, base_layer=-1, ref def NonNestedHierarchy(*meshes): return HierarchyBase(meshes, [None for _ in meshes], [None for _ in meshes], nested=False) + + +def SubmeshHierarchy(parent_hierarchy: HierarchyBase, + subdim: int | None = None, + subdomain_id: int | Sequence | None = None, + label_name: str | None = None, + name: str | None = None, + ignore_halo: bool = False, + reorder: bool | None = None, + comm=None): + """Build a hierarchy of submeshes from an existing mesh hierarchy. + + Parameters + ---------- + parent_hierarchy : + Parent mesh hierarchy. + subdim : + Topological dimension of the submesh. + Defaults to ``mesh.topological_dimension``. + subdomain_id : + Subdomain ID representing the submesh. + If `None` the submesh will cover the entire domain. + This is useful to obtain a codim-1 submesh over all facets or + a submesh over a different communicator. + label_name : + Name of the label to search ``subdomain_id`` in. + Defaults to ``'Cell Sets'`` or ``'Face Sets'`` depending on ``subdim``. + name : str | None + Name of the submesh. + Defaults to ``mesh.name + "_submesh"``· + ignore_halo : + Whether to exclude the halo from the submesh. + reorder : + Whether to reorder the mesh entities. By default, + the submesh will be reordered if the parent mesh was reordered. + comm : PETSc.Comm | None + An optional sub-communicator to define the submesh. + By default, the submesh is defined on `mesh.comm`. + + Returns + ------- + HierarchyBase + The submesh hierarchy. + + """ + meshes = [firedrake.Submesh(mesh, + subdim=subdim, subdomain_id=subdomain_id, + label_name=label_name, name=name, + ignore_halo=ignore_halo, reorder=reorder, + comm=comm) + for mesh in parent_hierarchy._meshes] + + lgmaps_with_overlap = [] + for i, m in enumerate(meshes): + lgmaps_with_overlap.append(impl.create_lgmap(m.topology_dm)) + m.topology_dm.setRefineLevel(i) + lgmaps_without_overlap = lgmaps_with_overlap + lgmaps = [ + (no, o) for no, o in zip(lgmaps_without_overlap, lgmaps_with_overlap) + ] + coarse_to_fine_cells = [] + fine_to_coarse_cells = [None] + for (coarse, fine), (clgmaps, flgmaps) in zip(zip(meshes[:-1], meshes[1:]), + zip(lgmaps[:-1], lgmaps[1:])): + c2f, f2c = impl.coarse_to_fine_cells(coarse, fine, clgmaps, flgmaps) + coarse_to_fine_cells.append(c2f) + fine_to_coarse_cells.append(f2c) + + refinements_per_level = parent_hierarchy.refinements_per_level + coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), c2f) + for i, c2f in enumerate(coarse_to_fine_cells)) + fine_to_coarse_cells = dict((Fraction(i, refinements_per_level), f2c) + for i, f2c in enumerate(fine_to_coarse_cells)) + return HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, + refinements_per_level=refinements_per_level, + nested=parent_hierarchy.nested) diff --git a/firedrake/mg/ufl_utils.py b/firedrake/mg/ufl_utils.py index d5a26a36e6..b40bfbd7ac 100644 --- a/firedrake/mg/ufl_utils.py +++ b/firedrake/mg/ufl_utils.py @@ -1,7 +1,7 @@ import ufl from ufl.corealg.map_dag import map_expr_dag from ufl.corealg.multifunction import MultiFunction -from ufl.domain import as_domain, extract_unique_domain +from ufl.domain import extract_unique_domain from ufl.duals import is_dual from functools import singledispatch, partial @@ -99,14 +99,18 @@ def coarsen_form(form, self, coefficient_mapping=None): integrals = [] for it in form.integrals(): integrand = map_expr_dag(mapper, it.integrand()) - mesh = as_domain(it) + mesh = it.ufl_domain() new_mesh = self(mesh, self) if isinstance(integrand, ufl.classes.Zero): continue if it.subdomain_data() is not None: raise CoarseningError("Don't know how to coarsen subdomain data") + # Coarsen secondary meshes in cross-mesh integrals (e.g. intersect_measures). + integral_type_map = {self(d, self): itype + for d, itype in it.extra_domain_integral_type_map().items()} new_itg = it.reconstruct(integrand=integrand, - domain=new_mesh) + domain=new_mesh, + extra_domain_integral_type_map=integral_type_map) integrals.append(new_itg) form = ufl.Form(integrals) return form @@ -278,47 +282,45 @@ def coarsen_snescontext(context, self, coefficient_mapping=None): # Assume not something that needs coarsening (e.g. float) new_appctx[k] = v + # Get options prefix for current level + parent_context = context + while parent_context._fine: + parent_context = parent_context._fine + parent_prefix = parent_context.options_prefix + opts = PETSc.Options(parent_prefix) + if opts.getString("snes_type", "") == "fas": + solver_prefix = "fas_" + else: + solver_prefix = "mg_" _, level = utils.get_level(problem.u_restrict.function_space().mesh()) if level == 0: - # Use different mat_type on coarsest level - opts = PETSc.Options(context.options_prefix) - if opts.getString("snes_type", "") == "fas": - solver_prefix = "fas_" - else: - solver_prefix = "mg_" - - coarse_mat_type = opts.getString(f"{solver_prefix}coarse_mat_type", "") - if coarse_mat_type == "": - coarse_mat_type = context.mat_type - default_pmat_type = context.pmat_type - else: - default_pmat_type = coarse_mat_type - coarse_pmat_type = opts.getString(f"{solver_prefix}coarse_pmat_type", - default_pmat_type) - - coarse_sub_mat_type = opts.getString(f"{solver_prefix}coarse_sub_mat_type", "") - if coarse_sub_mat_type == "": - coarse_sub_mat_type = context.sub_mat_type - default_sub_pmat_type = context.sub_pmat_type - else: - default_sub_pmat_type = coarse_sub_mat_type - coarse_sub_pmat_type = opts.getString(f"{solver_prefix}coarse_sub_pmat_type", - default_sub_pmat_type or "") or None + levels_prefix = f"{solver_prefix}coarse_" else: - coarse_mat_type = context.mat_type - coarse_pmat_type = context.pmat_type - coarse_sub_mat_type = context.sub_mat_type - coarse_sub_pmat_type = context.sub_pmat_type - - coarse = _SNESContext(problem, - mat_type=coarse_mat_type, - pmat_type=coarse_pmat_type, - sub_mat_type=coarse_sub_mat_type, - sub_pmat_type=coarse_sub_pmat_type, - appctx=new_appctx, - options_prefix=context.options_prefix, - transfer_manager=context.transfer_manager, - pre_apply_bcs=context.pre_apply_bcs) + levels_prefix = f"{solver_prefix}levels_" + current_level_prefix = f"{solver_prefix}levels_{level}_" + options_prefix = f"{parent_prefix}{current_level_prefix}" + + # Use different mat_type on each level + mat_type = None + pmat_type = None + sub_mat_type = None + sub_pmat_type = None + for prefix in (levels_prefix, current_level_prefix): + mat_type = opts.getString(f"{prefix}mat_type", "") or mat_type + pmat_type = opts.getString(f"{prefix}pmat_type", "") or pmat_type + sub_mat_type = opts.getString(f"{prefix}sub_mat_type", "") or sub_mat_type + sub_pmat_type = opts.getString(f"{prefix}sub_pmat_type", "") or sub_pmat_type + + pmat_type = pmat_type or mat_type + sub_pmat_type = sub_pmat_type or sub_mat_type + coarse = context.reconstruct(problem=problem, + mat_type=mat_type, + pmat_type=pmat_type, + sub_mat_type=sub_mat_type, + sub_pmat_type=sub_pmat_type, + appctx=new_appctx, + options_prefix=options_prefix, + ) coarse._coefficient_mapping = coefficient_mapping coarse._fine = context context._coarse = coarse diff --git a/firedrake/solving_utils.py b/firedrake/solving_utils.py index b1090d1111..db18865712 100644 --- a/firedrake/solving_utils.py +++ b/firedrake/solving_utils.py @@ -373,7 +373,7 @@ def split(self, fields): splits = [] problem = self._problem splitter = ExtractSubBlock() - for field in fields: + for field_num, field in enumerate(fields): F = splitter.split(problem.F, argument_indices=(field, )) J = splitter.split(problem.J, argument_indices=(field, field)) us = problem.u_restrict.subfunctions @@ -443,10 +443,10 @@ def split(self, fields): new_problem = NLVP(F, subu, bcs=bcs, J=J, Jp=Jp, is_linear=problem.is_linear, form_compiler_parameters=problem.form_compiler_parameters) new_problem._constant_jacobian = problem._constant_jacobian - splits.append(type(self)(new_problem, mat_type=self.mat_type, pmat_type=self.pmat_type, - appctx=self.appctx, - transfer_manager=self.transfer_manager, - pre_apply_bcs=self.pre_apply_bcs)) + name = V.name if len(V) == 1 else None + field_prefix = f"fieldsplit_{name or field_num}_" + options_prefix = f"{self.options_prefix}{field_prefix}" + splits.append(self.reconstruct(new_problem, options_prefix=options_prefix)) return self._splits.setdefault(tuple(fields), splits) @staticmethod diff --git a/pyop2/types/set.py b/pyop2/types/set.py index a913aab2b7..d2c28d8b49 100644 --- a/pyop2/types/set.py +++ b/pyop2/types/set.py @@ -1,3 +1,4 @@ +import copy import ctypes import numbers @@ -438,6 +439,9 @@ def _kernel_args_(self): def _argtypes_(self): return self._superset._argtypes_ + (ctypes.c_voidp, ) + def __deepcopy__(self, memo): + return type(self)(copy.deepcopy(self._superset, memo), self._indices.copy()) + # Look up any unspecified attributes on the _set. def __getattr__(self, name): """Returns a :class:`Set` specific attribute.""" diff --git a/pyproject.toml b/pyproject.toml index e3a63900e3..8811a56148 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,10 +18,17 @@ dependencies = [ "decorator<=4.4.2", "mpi4py>3; python_version >= '3.13'", "mpi4py; python_version < '3.13'", +<<<<<<< HEAD # TODO RELEASE "fenics-ufl @ git+https://github.com/FEniCS/ufl.git@main", # TODO RELEASE "firedrake-fiat @ git+https://github.com/firedrakeproject/fiat.git@main", +======= + # FEniCS may make new minor releases of UFL without warning so we pin each Firedrake + # release to a specific UFL minor version (e.g. 2025.3.x) + "fenics-ufl>=2025.3,<2025.4", + "firedrake-fiat>=2026.4", +>>>>>>> release "h5py>3.12.1", "firedrake-rtree", "immutabledict", diff --git a/scripts/firedrake-configure b/scripts/firedrake-configure index 865dd8eed6..7ab392d617 100755 --- a/scripts/firedrake-configure +++ b/scripts/firedrake-configure @@ -73,24 +73,25 @@ def main(): warn( "Passing '--no-package-manager' is deprecated, please set '--os unknown' instead" ) - args.os = OS.UNKNOWN - os_ = args.os if args.os else detect_os() + args.os = "unknown" - arch_key = (os_, args.scalar_type, args.gpu_arch) - arch_key_str = f"({', '.join(f'{k.__class__.__name__}: {k.name}' for k in arch_key)})" + os_ = OS(args.os) if args.os else detect_os() + scalar_type = ScalarType(args.scalar_type) + gpu_arch = GPUPlatform(args.gpu_arch) + + arch_key = (os_, scalar_type, gpu_arch) + arch_key_str = f"{{{', '.join(f'{k.__class__.__name__}: {k.value}' for k in arch_key)}}}" if arch_key in OFFICIAL_ARCHS: arch = OFFICIAL_ARCHS[arch_key] elif arch_key in COMMUNITY_ARCHS: arch = COMMUNITY_ARCHS[arch_key] - warn( - "You have selected a community-maintained Firedrake configuration, " - f"if you encounter issues please contact {arch.maintainer}" - ) + warn(f"""\ +You have selected a community-maintained Firedrake configuration, +if you encounter issues please contact {arch.maintainer}""") else: - error( - f"Firedrake configuration not found for option set '{arch_key_str}', " - "please consider passing '--os unknown' or building manually" - ) + error(f"""\ +Firedrake configuration not found for option set '{arch_key_str}', +please consider passing '--os unknown' or building manually""") if args.show_system_packages: try: @@ -132,8 +133,7 @@ Please see https://firedrakeproject.org/install for more information.""" package_manager_group.add_argument( "--os", "--package-manager", # old alias - choices=list(OS), - type=OS, + choices=[os_.value for os_ in OS], required=False, help="The operating system, if not provided 'firedrake-configure' " "will attempt to guess it.", @@ -149,16 +149,14 @@ Please see https://firedrakeproject.org/install for more information.""" parser.add_argument( "--scalar-type", "--arch", # old alias - choices=list(ScalarType), - type=ScalarType, - default=ScalarType.DEFAULT, + choices=[st.value for st in ScalarType], + default=ScalarType.DEFAULT.value, help="PETSc scalar type.", ) parser.add_argument( "--gpu-arch", - choices=list(GPUPlatform), - type=GPUPlatform, - default="none", + choices=[gpu_opt.value for gpu_opt in GPUPlatform], + default=GPUPlatform.NO_GPU.value, help=( "Target GPU architecture. WARNING: This is an experimental feature. " "GPU support in Firedrake is currently very limited." @@ -218,10 +216,10 @@ def detect_os() -> "OS": if platform.machine() == "arm64": return OS.MACOS_ARM64 - error( - "Cannot identify operating system, if you want to configure Firedrake " - "without any system packages please pass '--os unknown'" - ) + error("""\ +Cannot identify operating system, to configure Firedrake please specify the +operating system manually. Run 'python3 firedrake-configure --help' to see +the available options.""") class OS(enum.Enum): diff --git a/tests/firedrake/demos/test_demos_run.py b/tests/firedrake/demos/test_demos_run.py index f28371849f..bffd70ff51 100644 --- a/tests/firedrake/demos/test_demos_run.py +++ b/tests/firedrake/demos/test_demos_run.py @@ -52,7 +52,8 @@ Demo(("saddle_point_pc", "saddle_point_systems"), ["hypre", "mumps"]), Demo(("fast_diagonalisation", "fast_diagonalisation_poisson"), ["mumps"]), Demo(('vlasov_poisson_1d', 'vp1d'), []), - Demo(('shape_optimization', 'shape_optimization'), ["adjoint", "vtk"]) + Demo(('shape_optimization', 'shape_optimization'), ["adjoint", "vtk"]), + Demo(('submesh_reaction_diffusion', 'submesh_reaction_diffusion'), ["netgen", "vtk"]), ] PARALLEL_DEMOS = [ Demo(("full_waveform_inversion", "full_waveform_inversion"), ["adjoint"]), diff --git a/tests/firedrake/multigrid/test_options_prefix.py b/tests/firedrake/multigrid/test_options_prefix.py new file mode 100644 index 0000000000..204b584c7f --- /dev/null +++ b/tests/firedrake/multigrid/test_options_prefix.py @@ -0,0 +1,69 @@ +from firedrake import * +import pytest + + +@pytest.mark.parametrize("named", [False, True], ids=["unnamed", "named"]) +def test_fieldsplit_mg_options_prefix(named): + if named: + names = ("V1", "V2") + else: + names = (None, None) + refine = 2 + base = UnitIntervalMesh(10) + mh = MeshHierarchy(base, refine) + mesh = mh[-1] + V1 = FunctionSpace(mesh, 'CG', 1, name=names[0]) + V2 = FunctionSpace(mesh, 'CG', 2, name=names[1]) + W = MixedFunctionSpace([V1, V2]) + params0 = { + "pc_type": "mg", + "mg_levels_ksp_type": "chebyshev", + "mg_levels_pc_type": "jacobi", + "mg_coarse_mat_type": "aij", + "mg_coarse_ksp_type": "preonly", + "mg_coarse_pc_type": "lu", + } + params1 = { + "pc_type": "mg", + "mg_levels_ksp_type": "chebyshev", + "mg_levels_pc_type": "jacobi", + "mg_levels_0_mat_type": "aij", + "mg_levels_0_ksp_type": "preonly", + "mg_levels_0_pc_type": "lu", + } + sp = { + "mat_type": "matfree", + "ksp_rtol": 1E-12, + "ksp_max_it": 12, + "ksp_type": "cg", + "pc_type": "fieldsplit", + "pc_fieldsplit_type": "additive", + "fieldsplit_ksp_type": "preonly", + f"fieldsplit_{names[0] or 0}": params0, + f"fieldsplit_{names[1] or 1}": params1, + } + + x = SpatialCoordinate(mesh) + z_exact = as_vector([x[0], x[0]**2]) + z = Function(W) + w = TrialFunction(W) + v = TestFunction(W) + a = inner(grad(w), grad(v)) * dx + L = a(v, z_exact) + bcs = [DirichletBC(W.sub(i), z_exact[i], "on_boundary") for i in range(len(W))] + problem = LinearVariationalProblem(a, L, z, bcs=bcs) + solver = LinearVariationalSolver(problem, solver_parameters=sp) + solver.solve() + assert errornorm(z_exact, z) / norm(z_exact) < 1E-12 + + assert solver.snes.ksp.pc.getOperators()[0].getType() == "python" + fsplit = solver.snes.ksp.pc.getFieldSplitSubKSP() + assert len(fsplit) == len(W) + for ksp in fsplit: + coarse = ksp.pc.getMGSmoother(0) + assert coarse.pc.getType() == "lu" + assert coarse.getOperators()[0].getType().endswith("aij") + for l in range(1, len(mh)): + level = ksp.pc.getMGSmoother(l) + assert level.pc.getType() == "jacobi" + assert level.getOperators()[0].getType() == "python" diff --git a/tests/firedrake/multigrid/test_submesh_mg.py b/tests/firedrake/multigrid/test_submesh_mg.py new file mode 100644 index 0000000000..5b24f4da6f --- /dev/null +++ b/tests/firedrake/multigrid/test_submesh_mg.py @@ -0,0 +1,170 @@ +"""Tests for multigrid on coupled volume+surface PDE problems using Submesh. + +Three things are tested: + 1. Explicit submesh hierarchy construction: SubmeshHierarchy(mh, ...) should + tag all levels with __level_info__ so that MixedFunctionSpace([V, sV]) + no longer crashes. + 2. GMG on individual fieldsplit blocks (Stage 2). + 3. Monolithic GMG on the full mixed system (Stage 3). + +Problem: coupled Poisson–Helmholtz on Ω = [0,1]² with a 1D surface Γ = {y=0}. + + ∫_Ω ∇u·∇v dV + ∫_Γ (u−λ) v dC = ∫_Ω f_u v dV + ∫_Γ ∇λ·∇μ dA − ∫_Γ (u−λ) μ dC = ∫_Γ f_λ μ dA + +with dC = dx(smesh, intersect_measures=(ds(mesh),)) — a cross-mesh measure +that exercises the extra_domain_integral_type_map coarsening path. + +Manufactured solution: + u = cos(π x) cos(π y), + λ = cos(π x) on Γ, +""" + +import pytest +from firedrake import * +from firedrake.mg.utils import has_level, get_level + + +def build_problem(base_n=4, nref=1): + """Return the problem objects for the coupled Poisson–Helmholtz system.""" + base = UnitSquareMesh(base_n, base_n) + mh = MeshHierarchy(base, nref) + # marker 3 on UnitSquareMesh is the y=0 edge + smh = SubmeshHierarchy(mh, subdim=1, subdomain_id=3) + + mesh = mh[-1] + smesh = smh[-1] + + V = FunctionSpace(mesh, "CG", 1) + sV = FunctionSpace(smesh, "CG", 1) + Z = MixedFunctionSpace([V, sV]) + + dV = dx(mesh) + dA = dx(smesh) + # Cross-mesh measure: integrates over smesh intersected with all exterior + # facets of mesh. Since smesh IS a subset of ∂mesh, this picks out Γ. + dC = Measure("dx", smesh, intersect_measures=[ds(mesh)]) + + x, y = SpatialCoordinate(mesh) + xs = SpatialCoordinate(smesh)[0] # x-coordinate along Γ = {y=0} + + u, lam = TrialFunctions(Z) + v, mu = TestFunctions(Z) + a = (inner(grad(u), grad(v)) * dV + + inner(u - lam, v) * dC + + inner(grad(lam), grad(mu)) * dA + - inner(u - lam, mu) * dC) + + u_exact = cos(pi * x) * cos(pi * y) + lam_exact = cos(pi * xs) + + test, trial = a.arguments() + L = a(test, as_vector([u_exact, lam_exact])) + + bcs = [DirichletBC(Z.sub(0), u_exact, (1, 2, 4)), + DirichletBC(Z.sub(1), lam_exact, "on_boundary")] + + return mesh, smesh, Z, a, L, bcs, u_exact, lam_exact + + +# --------------------------------------------------------------------------- +# Stage 1: hierarchy construction +# --------------------------------------------------------------------------- + +def test_submesh_hierarchy_construction(): + """MeshHierarchy(submesh, nref) tags all levels with __level_info__.""" + base = UnitSquareMesh(4, 4) + mh = MeshHierarchy(base, 1) + + smh = SubmeshHierarchy(mh, subdim=1, subdomain_id=3) + assert len(mh) == len(smh) + + mesh = mh[-1] + smesh = smh[-1] + + assert has_level(smesh), "Fine submesh should have level info" + hierarchy, level = get_level(smesh) + assert level == 1, f"Fine submesh should be at level 1, got {level}" + assert len(hierarchy) == 2, f"Hierarchy should have 2 levels, got {len(hierarchy)}" + + smesh_coarse = hierarchy[0] + _, coarse_level = get_level(smesh_coarse) + assert coarse_level == 0, f"Coarse submesh should be at level 0, got {coarse_level}" + + assert smesh.submesh_parent is not None + assert smesh.submesh_parent is mesh + + # coarse_to_fine and fine_to_coarse maps must be non-None + assert hierarchy.coarse_to_fine_cells[0] is not None + assert hierarchy.fine_to_coarse_cells[1] is not None + + # MixedFunctionSpace construction must not crash (this was the original bug) + V = FunctionSpace(mesh, "CG", 1) + sV = FunctionSpace(smesh, "CG", 1) + Z = MixedFunctionSpace([V, sV]) + assert tuple(Z.mesh()) == (mesh, smesh) + + +# --------------------------------------------------------------------------- +# Stages 2 & 3: GMG solver tests +# --------------------------------------------------------------------------- + +def fieldsplit_gmg_params(): + """GMG on individual fieldsplit blocks.""" + return { + "ksp_type": "cg", + "ksp_monitor": None, + "ksp_rtol": 1e-10, + "ksp_atol": 0, + "pc_type": "fieldsplit", + "pc_fieldsplit_type": "additive", + "fieldsplit_ksp_type": "preonly", + "fieldsplit_pc_type": "mg", + "fieldsplit_mg_levels_ksp_type": "chebyshev", + "fieldsplit_mg_levels_pc_type": "jacobi", + "fieldsplit_mg_levels_ksp_max_it": 2, + "fieldsplit_mg_coarse_pc_type": "lu", + "fieldsplit_mg_coarse_pc_factor_mat_solver_type": "mumps", + } + + +def monolithic_gmg_params(): + """Monolithic GMG on the full coupled system.""" + return { + "ksp_type": "cg", + "ksp_monitor": None, + "ksp_rtol": 1e-10, + "ksp_atol": 0, + "pc_type": "mg", + "mg_levels_ksp_type": "chebyshev", + "mg_levels_pc_type": "jacobi", + "mg_levels_ksp_max_it": 2, + "mg_coarse_pc_type": "lu", + "mg_coarse_pc_factor_mat_solver_type": "mumps", + } + + +@pytest.mark.parallel([1, 4]) +@pytest.mark.parametrize("solver_type", ["fieldsplit_gmg", "monolithic_gmg"]) +def test_submesh_gmg(solver_type): + """GMG converges in O(1) iterations and recovers the correct solution.""" + mesh, smesh, Z, a, L, bcs, u_exact, lam_exact = build_problem(base_n=4, nref=2) + + params = fieldsplit_gmg_params() if solver_type == "fieldsplit_gmg" else monolithic_gmg_params() + + z = Function(Z) + problem = LinearVariationalProblem(a, L, z, bcs=bcs) + solver = LinearVariationalSolver(problem, solver_parameters=params) + solver.solve() + + ksp_its = solver.snes.ksp.getIterationNumber() + assert ksp_its < 15, ( + f"Expected < 15 KSP iterations with {solver_type}, got {ksp_its}. " + "This suggests the multigrid hierarchy or preconditioner is broken." + ) + + u_h, lam_h = z.subfunctions + err_u = errornorm(u_exact, u_h) / norm(u_exact) + err_lam = errornorm(lam_exact, lam_h) / norm(lam_exact) + assert err_u < 2e-2, f"Volume error too large: {err_u}" + assert err_lam < 4e-3, f"Surface error too large: {err_lam}" diff --git a/tests/firedrake/regression/test_assemble.py b/tests/firedrake/regression/test_assemble.py index 9972715746..6b5018b157 100644 --- a/tests/firedrake/regression/test_assemble.py +++ b/tests/firedrake/regression/test_assemble.py @@ -415,3 +415,41 @@ def test_assemble_tensor_empty_shape(mesh): v = Function(V).assign(1) expected = assemble(inner(v, v)*dx) assert np.allclose(result, expected) + + +@pytest.mark.parametrize("coefficient", [False, True], ids=["Expr", "Function"]) +def test_cell_avg(coefficient): + mesh = UnitSquareMesh(3, 3) + x = SpatialCoordinate(mesh) + expr = dot(x, x) ** 2 + if coefficient: + V = FunctionSpace(mesh, "CG", 4) + expr = Function(V).interpolate(expr) + + result = assemble(inner(cell_avg(expr), expr) * dx) + + Q = FunctionSpace(mesh, "DG", 0) + p = Function(Q) + p.project(expr) + expect = assemble(inner(p, expr) * dx) + assert np.isclose(result, expect) + + +def test_cell_avg_mfs(): + mesh = UnitSquareMesh(3, 3) + V = VectorFunctionSpace(mesh, "CG", 2) + Q = FunctionSpace(mesh, "DG", 0) + Z = MixedFunctionSpace([V, Q]) + z = Function(Z) + usub, psub = z.subfunctions + usub.interpolate(SpatialCoordinate(mesh)) + psub.interpolate(Constant(1)) + + expect = 6 + result1 = assemble(inner(cell_avg(div(usub)), psub + div(usub))*dx) + assert np.isclose(result1, expect) + + # This fails if do_replace_functions=True in entity_avg in tscf/ufl_utils.py + u, p = split(z) + result2 = assemble(inner(cell_avg(div(u)), p + div(u))*dx) + assert np.isclose(result2, expect) diff --git a/tests/firedrake/regression/test_interpolate.py b/tests/firedrake/regression/test_interpolate.py index 66d733c57c..f8c3c528e8 100644 --- a/tests/firedrake/regression/test_interpolate.py +++ b/tests/firedrake/regression/test_interpolate.py @@ -741,3 +741,17 @@ def test_interpolate_form_mixed(): res3 = assemble(inner(u, q) * dx) # V x W -> R assert mat_equals(res1, res3) + + +def test_interpolation_cell_subset(): + mesh = UnitIntervalMesh(2) + x = SpatialCoordinate(mesh) + + ind = Function(FunctionSpace(mesh, "DG", 0)).interpolate(conditional(x[0] > 0.5, 1, 0)) + + mesh = RelabeledMesh(mesh, [ind], [100]) + + a = Constant(1) + u = Function(FunctionSpace(mesh, "DG", 0)).interpolate(-a*-a, subset=mesh.cell_subset(100)) + + assert assemble((u-a/2)*dx) < 1e-14 diff --git a/tests/firedrake/submesh/test_submesh_interface.py b/tests/firedrake/submesh/test_submesh_interface.py new file mode 100644 index 0000000000..394738fab8 --- /dev/null +++ b/tests/firedrake/submesh/test_submesh_interface.py @@ -0,0 +1,45 @@ +import pytest +import numpy as np +from firedrake import * + + +def test_submesh_subdomain_id_union(): + mesh = UnitSquareMesh(4, 4) + x, y = SpatialCoordinate(mesh) + M = FunctionSpace(mesh, "DG", 0) + m1 = Function(M).interpolate(conditional(lt(x, 0.5), 1, 0)) + m2 = Function(M).interpolate(conditional(lt(y, 0.5), 1, 0)) + mesh.mark_entities(m1, 111) + mesh.mark_entities(m2, 222) + + subdomain_id = [111, 222] + submesh1 = Submesh(mesh, mesh.topological_dimension, subdomain_id=subdomain_id) + + m3 = Function(M).interpolate(m1 + m2 - m1 * m2) + expected = assemble(m3*dx) + assert abs(assemble(1*dx(domain=submesh1)) - expected) < 1E-12 + + mesh.mark_entities(m3, 333) + submesh2 = Submesh(mesh, mesh.topological_dimension, 333) + assert submesh2.cell_set.size == submesh1.cell_set.size + assert np.allclose(submesh2.coordinates.dat.data, submesh1.coordinates.dat.data) + + +@pytest.mark.parametrize("subdomain_id", ["on_boundary", (1, 3, 6)]) +def test_submesh_facet_subdomain_id_union(subdomain_id): + mesh = UnitCubeMesh(2, 2, 2) + submesh1 = Submesh(mesh, mesh.topological_dimension - 1, subdomain_id=subdomain_id) + if subdomain_id == "on_boundary": + area = assemble(1*ds(domain=mesh)) + else: + area = assemble(1*ds(subdomain_id, domain=mesh)) + assert abs(assemble(1*dx(domain=submesh1)) - area) < 1E-12 + + DGT = FunctionSpace(mesh, "DGT", 0) + facet_function = Function(DGT) + DirichletBC(DGT, 1, subdomain_id).apply(facet_function) + facet_value = 999 + rmesh = RelabeledMesh(mesh, [facet_function], [facet_value]) + submesh2 = Submesh(rmesh, mesh.topological_dimension - 1, facet_value) + assert submesh2.cell_set.size == submesh1.cell_set.size + assert np.allclose(submesh2.coordinates.dat.data, submesh1.coordinates.dat.data) diff --git a/tests/firedrake/submesh/test_submesh_interpolate.py b/tests/firedrake/submesh/test_submesh_interpolate.py index 39d7f8d70d..c3b084a878 100644 --- a/tests/firedrake/submesh/test_submesh_interpolate.py +++ b/tests/firedrake/submesh/test_submesh_interpolate.py @@ -346,3 +346,51 @@ def test_submesh_interpolate_adjoint(fe_fesub): # Test 0-form result_0 = assemble(interpolate(u1, ustar2, allow_missing_dofs=True)) assert np.isclose(result_0, expected) + + +@pytest.mark.parallel(nprocs=8) +def test_submesh_interpolate_3Dcell_2Dfacet_empty_rank_8_processes(): + # Regression test: ranks with zero Submesh cells crashed with IndexError + # in _pic_swarm_in_mesh (-1 sentinels invalid on empty-rank cell_closure). + mesh = UnitCubeMesh(2, 2, 2) + x, y, z = SpatialCoordinate(mesh) + V_marker = FunctionSpace(mesh, "HDiv Trace", 0) + facet_indicator = Function(V_marker).interpolate( + conditional(x > 0.999, 1.0, 0.0) + ) + facet_value = 999 + mesh = RelabeledMesh(mesh, [facet_indicator], [facet_value]) + subm = Submesh(mesh, mesh.topological_dimension - 1, facet_value) + subm.tolerance = 0.1 + x, y, z = SpatialCoordinate(mesh) + xs, ys, zs = SpatialCoordinate(subm) + V_sub = FunctionSpace(subm, "CG", 1) + V_parent = FunctionSpace(mesh, "CG", 1) + f_sub = Function(V_sub).interpolate(ys + zs) + f_parent = Function(V_parent) + f_parent.interpolate(f_sub, allow_missing_dofs=True) + expected = Function(V_parent).interpolate( + conditional(x > 0.999, y + z, 0.0) + ) + assert np.allclose( + f_parent.dat.data_with_halos, expected.dat.data_with_halos + ) + + +@pytest.mark.parallel(nprocs=8) +def test_submesh_interpolate_3Dcell_extruded_empty_rank_8_processes(): + # Regression test for the same IndexError in the extruded branch of + # _pic_swarm_in_mesh: ExtrudedMesh(UnitSquareMesh(1,1), layers=3) has 6 + # cells, so with nprocs=8 at least two ranks own zero extruded cells. + base = UnitSquareMesh(1, 1) + ext = ExtrudedMesh(base, layers=3) + xe, ye, ze = SpatialCoordinate(ext) + V_ext = FunctionSpace(ext, "CG", 1) + f_ext = Function(V_ext).interpolate(xe + ye + ze) + mesh = UnitCubeMesh(2, 2, 2) + xt, yt, zt = SpatialCoordinate(mesh) + V = FunctionSpace(mesh, "CG", 1) + f = Function(V) + f.interpolate(f_ext, allow_missing_dofs=True) + expected = Function(V).interpolate(xt + yt + zt) + assert np.allclose(f.dat.data_with_halos, expected.dat.data_with_halos) diff --git a/tsfc/ufl_utils.py b/tsfc/ufl_utils.py index c26febd68e..d16b5b93f5 100644 --- a/tsfc/ufl_utils.py +++ b/tsfc/ufl_utils.py @@ -43,6 +43,7 @@ def compute_form_data(form, do_apply_default_restrictions=True, do_apply_restrictions=True, do_estimate_degrees=True, + do_replace_functions=True, coefficients_to_split=None, complex_mode=False): """Preprocess UFL form in a format suitable for TSFC. Return @@ -52,6 +53,9 @@ def compute_form_data(form, kwargs overriden in the way TSFC needs it and is provided for other form compilers based on TSFC. """ + # Multidomain problems require further index simplifications to ensure + # that unwanted quantities do not appear inside single-domain integrals. + do_remove_component_tensors = len(form.ufl_domains()) > 1 fd = ufl_compute_form_data( form, do_apply_function_pullbacks=do_apply_function_pullbacks, @@ -61,9 +65,10 @@ def compute_form_data(form, do_apply_default_restrictions=do_apply_default_restrictions, do_apply_restrictions=do_apply_restrictions, do_estimate_degrees=do_estimate_degrees, - do_replace_functions=True, + do_replace_functions=do_replace_functions, coefficients_to_split=coefficients_to_split, - complex_mode=complex_mode + complex_mode=complex_mode, + do_remove_component_tensors=do_remove_component_tensors, ) constants = extract_firedrake_constants(form) fd.constants = constants @@ -108,7 +113,9 @@ def entity_avg(integrand, measure, argument_multiindices): degree = estimate_total_polynomial_degree(integrand) form = integrand * measure fd = compute_form_data(form, do_estimate_degrees=False, - do_apply_function_pullbacks=False) + do_apply_function_pullbacks=False, + do_replace_functions=False, + ) itg_data, = fd.integral_data integral, = itg_data.integrals integrand = integral.integrand()