diff --git a/.github/actions/install-dolfinx/action.yaml b/.github/actions/install-dolfinx/action.yaml index b39de196..694ffd24 100644 --- a/.github/actions/install-dolfinx/action.yaml +++ b/.github/actions/install-dolfinx/action.yaml @@ -25,7 +25,7 @@ inputs: petsc_arch: # id of input description: 'PETSc Arch' required: false - default: 'linux-gnu-real-32' + default: 'linux-gnu-real64-32' runs: using: composite @@ -33,14 +33,14 @@ runs: steps: - name: Get Basix - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: path: ./basix repository: FEniCS/basix ref: ${{ inputs.basix }} - name: Get DOLFINx - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: path: ./dolfinx repository: FEniCS/dolfinx diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml index 4ae02c5b..9dc3af4b 100644 --- a/.github/workflows/build_docs.yml +++ b/.github/workflows/build_docs.yml @@ -10,11 +10,11 @@ on: jobs: build: runs-on: ubuntu-22.04 - container: ghcr.io/fenics/dolfinx/dev-env:nightly + container: ghcr.io/fenics/dolfinx/dev-env:v0.7.0-mpich env: # Directory that will be published on github pages - PUBLISH_DIR: ./docs/_build/html - PETSC_ARCH: "linux-gnu-real-32" + PUBLISH_DIR: ./_build/html + PETSC_ARCH: "linux-gnu-real64-32" steps: - uses: actions/checkout@v4 @@ -32,7 +32,7 @@ jobs: run: python3 -m pip -v install python/[docs] - name: Build docs - run: jupyter book build docs + run: jupyter book build . - name: Upload documentation as artifact uses: actions/upload-artifact@v3 diff --git a/.github/workflows/deploy-pages.yml b/.github/workflows/deploy-pages.yml index 489787d3..0fd63fde 100644 --- a/.github/workflows/deploy-pages.yml +++ b/.github/workflows/deploy-pages.yml @@ -46,7 +46,7 @@ jobs: path: "./public/code-coverage-report" - name: Upload artifact - uses: actions/upload-pages-artifact@v1 + uses: actions/upload-pages-artifact@v2 with: path: "./public" @@ -54,8 +54,8 @@ jobs: uses: actions/checkout@v4 - name: Setup Pages - uses: actions/configure-pages@v2 + uses: actions/configure-pages@v3 - name: Deploy to GitHub Pages id: deployment - uses: actions/deploy-pages@v1 + uses: actions/deploy-pages@v2 diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml index fd2bdc84..9407a66f 100644 --- a/.github/workflows/docker.yml +++ b/.github/workflows/docker.yml @@ -23,13 +23,13 @@ jobs: uses: actions/checkout@v4 - name: Set up QEMU - uses: docker/setup-qemu-action@v2 + uses: docker/setup-qemu-action@v3 - name: Set up Docker Buildx - uses: docker/setup-buildx-action@v2 + uses: docker/setup-buildx-action@v3 - name: Log in to the Container registry - uses: docker/login-action@v2 + uses: docker/login-action@v3 with: registry: ${{ env.REGISTRY }} username: ${{ github.actor }} @@ -37,12 +37,12 @@ jobs: - name: Extract metadata (tags, labels) for Docker id: meta - uses: docker/metadata-action@v4 + uses: docker/metadata-action@v5 with: images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} - name: Build and push Docker image - uses: docker/build-push-action@v3 + uses: docker/build-push-action@v5 with: context: . push: true diff --git a/.github/workflows/sonarcloud.yml b/.github/workflows/sonarcloud.yml index 8e2b652d..83529a19 100644 --- a/.github/workflows/sonarcloud.yml +++ b/.github/workflows/sonarcloud.yml @@ -7,7 +7,6 @@ on: branches: - main - jobs: build: name: Build @@ -16,15 +15,15 @@ jobs: # avoid running on pull requests from forks # which don't have access to secrets if: (github.event_name != 'pull_request' || github.event.pull_request.head.repo.full_name == github.repository) - container: ghcr.io/fenics/test-env:nightly-mpich + container: ghcr.io/fenics/dolfinx/dev-env:v0.7.0 env: - SONAR_SCANNER_VERSION: - 4.8.0.2856 # Find the latest version in at: - # https://github.com/SonarSource/sonar-scanner-cli/tags + SONAR_SCANNER_VERSION: 5.0.1.3006 # Find the latest version in at: https://github.com/SonarSource/sonar-scanner-cli/tags SONAR_SERVER_URL: "https://sonarcloud.io" BUILD_WRAPPER_OUT_DIR: build_wrapper_output_directory # Directory where build-wrapper output will be placed - PETSC_ARCH: linux-gnu-real-32 + PETSC_ARCH: linux-gnu-real64-32 + + PETSC_DIR: /usr/local/petsc steps: - uses: actions/checkout@v4 with: @@ -35,12 +34,12 @@ jobs: apt-get -y update apt-get install unzip - - name: Set up JDK 11 + - name: Set up JDK 17 uses: actions/setup-java@v3 with: - distribution: 'zulu' - java-version: 11 - + distribution: "zulu" + java-version: 17 + - name: Cache SonarCloud packages uses: actions/cache@v3 with: @@ -67,7 +66,9 @@ jobs: - name: Install DOLFINx uses: ./.github/actions/install-dolfinx - + with: + petsc_arch: ${PETSC_ARCH} + - name: Run build-wrapper run: | cmake -S ./cpp -B build-mpc @@ -87,4 +88,4 @@ jobs: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} SONAR_TOKEN: ${{ secrets.SONAR_TOKEN }} run: | - sonar-scanner --define sonar.host.url="${{ env.SONAR_SERVER_URL }}" --define sonar.cfamily.build-wrapper-output="${{ env.BUILD_WRAPPER_OUT_DIR }}" \ No newline at end of file + sonar-scanner --define sonar.host.url="${{ env.SONAR_SERVER_URL }}" --define sonar.cfamily.build-wrapper-output="${{ env.BUILD_WRAPPER_OUT_DIR }}" diff --git a/.github/workflows/test_mpc.yml b/.github/workflows/test_mpc.yml index 6b334b49..407e5b4d 100644 --- a/.github/workflows/test_mpc.yml +++ b/.github/workflows/test_mpc.yml @@ -14,12 +14,18 @@ on: jobs: build: runs-on: ubuntu-latest - container: ghcr.io/fenics/test-env:nightly-mpich + container: ghcr.io/fenics/dolfinx/dev-env:v0.7.0 strategy: matrix: build_mode: [Release, Debug] - petsc_arch: [real, complex] + petsc_arch: + [ + linux-gnu-complex128-32, + linux-gnu-complex64-32, + linux-gnu-real64-32, + linux-gnu-real32-32, + ] # Due to: https://gitlab.com/petsc/petsc/-/issues/1288 CXX: [c++] #, clang++] @@ -31,13 +37,13 @@ jobs: # - CC: clang # CXX: c++ env: - DOLFINX_BRANCH: "main" - BASIX_BRANCH: "main" - UFL_BRANCH: "main" - FFCX_BRANCH: "main" + DOLFINX_BRANCH: "v0.7.1" + BASIX_BRANCH: "v0.7.0" + UFL_BRANCH: "2023.2.0" + FFCX_BRANCH: "v0.7.0" CC: ${{ matrix.CC }} CXX: ${{ matrix.CXX }} - PETSC_ARCH: "linux-gnu-${{ matrix.petsc_arch }}-32" + PETSC_ARCH: "${{ matrix.petsc_arch }}" OMPI_ALLOW_RUN_AS_ROOT: 1 OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1 OMPI_MCA_rmaps_base_oversubscribe: 1 @@ -73,10 +79,9 @@ jobs: cd python python3 -m mypy . --exclude=build - - name: Install h5py (https://github.com/h5py/h5py/issues/2222) + - name: Install h5py run: | - python3 -m pip install mpi4py cython numpy - python3 -m pip install --no-cache-dir --no-binary=h5py h5py --no-build-isolation + python3 -m pip install --no-cache-dir --no-binary=h5py h5py - name: Install DOLFINx uses: ./.github/actions/install-dolfinx @@ -94,10 +99,11 @@ jobs: cmake --install build-dir - name: Install DOLFINx-MPC (Python) - run: CXX_FLAGS="${MPC_CMAKE_CXX_FLAGS}" python3 -m pip -v install python/[test] + run: CXX_FLAGS="${MPC_CMAKE_CXX_FLAGS}" python3 -m pip -v install -e python/[test] + - name: Run tests (serial) - run: coverage run --rcfile=.coveragerc -m mpi4py -m pytest python/tests/ + run: coverage run --rcfile=.coveragerc -m mpi4py -m pytest python/tests/ -vs - name: Run tests (2 processes) run: mpirun -n 2 coverage run --rcfile=.coveragerc -m mpi4py -m pytest python/tests/ -vs diff --git a/Changelog.md b/Changelog.md index 8c70a7fa..2e2ea90e 100644 --- a/Changelog.md +++ b/Changelog.md @@ -1,37 +1,45 @@ # Changelog -## main +## v0.7.0 + - **API**: - Change input of `dolfinx_mpc.MultiPointConstraint.homogenize` and `dolfinx_mpc.backsubstitution` to `dolfinx.fem.Function` instead of `PETSc.Vec`. - - Add support for more floating types (`dtype` can now be specified as input for coefficients). In theory this means one could support single precision. Not thoroughly tested. + - **New feature**: Add support for more floating types (float32, float64, complex64, complex128). The floating type of a MPC is related to the mesh geometry. - This resulted in a minor refactoring of the pybindings, meaning that tte class `dolfinx_mpc.cpp.mpc.MultiPointConstraint` is replaced by `dolfinx_mpc.cpp.mpc.MultiPointConstraint_{dtype}` + - Casting scalar-type with `dolfinx.default_scalar_type` instead of `PETSc.ScalarType` + - Remove usage of `VectorFunctionSpace`. Use blocked basix element instead. - **DOLFINX API-changes**: - Use `dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))` instead of `dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 1))` as the latter is being deprecated. - Use `basix.ufl.element` in favor of `ufl.FiniteElement` as the latter is deprecated in DOLFINx. ## v0.6.1 (30.01.2023) - - Fixes for CI - - Add auto-publishing CI - - Fixes for `h5py` installation + +- Fixes for CI +- Add auto-publishing CI +- Fixes for `h5py` installation ## v0.6.0 (27.01.2023) - - Remove `dolfinx::common::impl::copy_N` in favor of `std::copy_n` by @jorgensd in #24 - - Improving and fixing `demo_periodic_gep.py` by @fmonteghetti in #22 and @conpierce8 in #30 - - Remove xtensor by @jorgensd in #25 - - Complex valued periodic constraint (scale) by @jorgensd in #34 - - Implement Hermitian pre-multiplication by @conpierce8 in #38 - - Fixes for packaging by @mirk in #41, #42, #4 - - Various updates to dependencies + +- Remove `dolfinx::common::impl::copy_N` in favor of `std::copy_n` by @jorgensd in #24 +- Improving and fixing `demo_periodic_gep.py` by @fmonteghetti in #22 and @conpierce8 in #30 +- Remove xtensor by @jorgensd in #25 +- Complex valued periodic constraint (scale) by @jorgensd in #34 +- Implement Hermitian pre-multiplication by @conpierce8 in #38 +- Fixes for packaging by @mirk in #41, #42, #4 +- Various updates to dependencies ## v0.5.0 (12.08.2022) - - Minimal C++ standard is now [C++20](https://en.cppreference.com/w/cpp/20) - - Deprecating GMSH IO functions from `dolfinx_mpc.utils`, see: [DOLFINx PR: 2261](https://github.com/FEniCS/dolfinx/pull/2261) for details. - - Various API changes in DOLFINx relating to `dolfinx.common.IndexMap`. - - Made code [mypy](https://mypy.readthedocs.io/en/stable/)-compatible (tests added to CI). - - Made code [PEP-561](https://peps.python.org/pep-0561/) compatible. + +- Minimal C++ standard is now [C++20](https://en.cppreference.com/w/cpp/20) +- Deprecating GMSH IO functions from `dolfinx_mpc.utils`, see: [DOLFINx PR: 2261](https://github.com/FEniCS/dolfinx/pull/2261) for details. +- Various API changes in DOLFINx relating to `dolfinx.common.IndexMap`. +- Made code [mypy](https://mypy.readthedocs.io/en/stable/)-compatible (tests added to CI). +- Made code [PEP-561](https://peps.python.org/pep-0561/) compatible. ## v0.4.0 (30.04.2022) + - **API**: + - **New feature**: Support for nonlinear problems (by @nate-sime) for mpc, see `test_nonlinear_assembly.py` for usage - Updated user interface for `dolfinx_mpc.create_slip_constraint`. See documentation for details. - **New feature**: Support for periodic constraints on sub-spaces. See `dolfinx_mpc.create_periodic_constraint` for details. @@ -50,11 +58,13 @@ - `slave_cells` does now longer exist as it can be gotten implicitly from `cell_to_slaves`. - **Performance**: + - Major rewrite of periodic boundary conditions. On average at least a 5 x performance speed-up. - The C++ assembler has been fully rewritten. - Various improvements to `ContactConstraint`. - **Bugs** + - Resolved issue where `create_facet_normal_approximation` would give you a 0 normal for a surface dof it was not owned by any of the cells with facets on the surface. - **DOLFINX API-changes**: @@ -67,19 +77,22 @@ - Various internal changes to handle new way of JIT-compliation of `dolfinx::fem::Form_{scalar_type}` ## 0.3.0 (25.08.2021) + - Minor internal changes ## 0.2.0 (06.08.2021) + - Add new MPC constraint: Periodic boundary condition constrained geometrically. See `demo_periodic_geometrical.py` for use-case. -- New: `demo_periodic_gep.py` proposed and initally implemented by [fmonteghetti](https://github.com/fmonteghetti) using SLEPc for eigen-value problems. - This demo illustrates the usage of the new `diagval` keyword argument in the `assemble_matrix` class. +- New: `demo_periodic_gep.py` proposed and initally implemented by [fmonteghetti](https://github.com/fmonteghetti) using SLEPc for eigen-value problems. + This demo illustrates the usage of the new `diagval` keyword argument in the `assemble_matrix` class. - **API**: + - Renaming and clean-up of `assemble_matrix` in C++ - Renaming of Periodic constraint due to additional geometrical constraint, `mpc.create_periodic_constraint` -> `mpc.create_periodic_constraint_geometrical/topological`. - Introduce new class `dolfinx_mpc.LinearProblem` mimicking the DOLFINx class (Usage illustrated in `demo_periodic_geometrical.py`) - Additional `kwarg` `b: PETSc.Vec` for `assemble_vector` to be able to re-use Vector. - - Additional `kwargs`: `form_compiler_parameters` and `jit_parameters` to `assemble_matrix`, `assemble_vector`, to allow usage of fast math etc. + - Additional `kwargs`: `form_compiler_parameters` and `jit_parameters` to `assemble_matrix`, `assemble_vector`, to allow usage of fast math etc. - **Performance**: - Slip condition constructor moved to C++ (Speedup for large problems) @@ -90,6 +103,6 @@ - `dolfinx.cpp.la.scatter_forward(x)` is replaced by `x.scatter_forward()` - Various interal updates to match DOLFINx API (including dof transformations moved outside of ffcx kernel) - ## 0.1.0 (11.05.2021) + - First tagged release of dolfinx_mpc, compatible with [DOLFINx 0.1.0](https://github.com/FEniCS/dolfinx/releases/tag/0.1.0). diff --git a/docs/_config.yml b/_config.yml similarity index 85% rename from docs/_config.yml rename to _config.yml index eae59601..4ce16a23 100644 --- a/docs/_config.yml +++ b/_config.yml @@ -5,11 +5,10 @@ title: "DOLFINx-MPC: An extension to DOLFINx for multi point constraints" author: Jørgen S. Dokken copyright: "2022" -# Force re-execution of notebooks on each build. -# See https://jupyterbook.org/content/execute.html execute: - execute_notebooks: force + execute_notebooks: cache +only_build_toc_files: true # Information about where the book exists on the web repository: @@ -24,6 +23,7 @@ launch_buttons: html: use_issues_button: true use_repository_button: true + parse: myst_enable_extensions: @@ -37,5 +37,10 @@ sphinx: - 'sphinx.ext.napoleon' - 'sphinx.ext.viewcode' + config: + nb_custom_formats: + .py: + - jupytext.reads + - fmt: py diff --git a/_toc.yml b/_toc.yml new file mode 100644 index 00000000..6fabfc58 --- /dev/null +++ b/_toc.yml @@ -0,0 +1,12 @@ +format: jb-book +root: index + +parts: + - caption: "Demos" + chapters: + - file: "python/demos/demo_stokes.py" + - caption: "Python API" + chapters: + - file: "docs/api" + - file: "docs/utils" + - file: "docs/numba" diff --git a/cpp/PeriodicConstraint.h b/cpp/PeriodicConstraint.h index 32c866f4..62c38ed7 100644 --- a/cpp/PeriodicConstraint.h +++ b/cpp/PeriodicConstraint.h @@ -73,7 +73,7 @@ dolfinx_mpc::mpc_data _create_periodic_condition( // Tolerance for adding scaled basis values to MPC. Any scaled basis // value with lower absolute value than the tolerance is ignored - const U tol = 1e-13; + const U tol = 500 * std::numeric_limits::epsilon(); auto mesh = V.mesh(); auto dofmap = V.dofmap(); @@ -355,7 +355,7 @@ dolfinx_mpc::mpc_data _create_periodic_condition( num_masters_per_slave_remote.reserve(bs * coords_recvb.size() / 3); std::vector remote_cell_collisions - = dolfinx_mpc::find_local_collisions(*mesh, tree, coords_recvb, 1e-20); + = dolfinx_mpc::find_local_collisions(*mesh, tree, coords_recvb, tol); auto [remote_basis_valuesb, r_basis_shape] = dolfinx_mpc::evaluate_basis_functions(V, coords_recvb, remote_cell_collisions); diff --git a/cpp/utils.h b/cpp/utils.h index b5101624..a8cb0ab7 100644 --- a/cpp/utils.h +++ b/cpp/utils.h @@ -1112,7 +1112,8 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, else { // Pull-back physical point xp to reference coordinate Xp - cmaps[0].pull_back_nonaffine(Xp, xp, coord_dofs); + cmaps[0].pull_back_nonaffine(Xp, xp, coord_dofs, + 500 * std::numeric_limits::epsilon(), 15); cmaps[0].tabulate(1, std::span(Xpb.data(), tdim), {1, tdim}, phi_b); dolfinx::fem::CoordinateElement::compute_jacobian(dphi, coord_dofs, diff --git a/docker/Dockerfile b/docker/Dockerfile index 422e509b..97f5f8f6 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -1,4 +1,4 @@ -FROM ghcr.io/fenics/dolfinx/dolfinx:v0.6.0-r1 +FROM ghcr.io/fenics/dolfinx/dolfinx:v0.7.0 WORKDIR /tmp diff --git a/docs/_toc.yml b/docs/_toc.yml deleted file mode 100644 index 5497cd6e..00000000 --- a/docs/_toc.yml +++ /dev/null @@ -1,9 +0,0 @@ -format: jb-book -root: index - -parts: - - caption: "Python API" - chapters: - - file: "api" - - file: "utils" - - file: "numba" diff --git a/docs/index.md b/index.md similarity index 100% rename from docs/index.md rename to index.md diff --git a/python/benchmarks/bench_contact_3D.py b/python/benchmarks/bench_contact_3D.py index d30377eb..71a2a501 100644 --- a/python/benchmarks/bench_contact_3D.py +++ b/python/benchmarks/bench_contact_3D.py @@ -8,27 +8,36 @@ # between two cubes. from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser +from pathlib import Path +import basix.ufl import numpy as np +from dolfinx import default_scalar_type, default_real_type from dolfinx.common import Timer, TimingType, list_timings, timing from dolfinx.cpp.mesh import entities_to_geometry -from dolfinx.fem import (Constant, Function, VectorFunctionSpace, dirichletbc, - form, locate_dofs_topological, set_bc) +from dolfinx.fem import (Constant, Function, dirichletbc, form, functionspace, + locate_dofs_topological, set_bc) from dolfinx.io import XDMFFile from dolfinx.mesh import (CellType, compute_midpoints, create_mesh, create_unit_cube, locate_entities_boundary, meshtags, refine) -from dolfinx_mpc import (MultiPointConstraint, apply_lifting, assemble_matrix, - assemble_vector) -from dolfinx_mpc.utils import (create_normal_approximation, log_info, - rigid_motions_nullspace, rotation_matrix) from mpi4py import MPI from petsc4py import PETSc from ufl import (Cell, Identity, Mesh, TestFunction, TrialFunction, VectorElement, dx, grad, inner, sym, tr) -from pathlib import Path + +from dolfinx_mpc import (MultiPointConstraint, apply_lifting, assemble_matrix, + assemble_vector) +from dolfinx_mpc.utils import (create_normal_approximation, log_info, + rigid_motions_nullspace, rotation_matrix) +import warnings comm = MPI.COMM_WORLD +if default_real_type == np.float32: + warnings.warn( + "Demo not supported in single precision as reading from XDMF only support double precision meshes") + exit(0) + def mesh_3D_dolfin(theta=0, ct=CellType.tetrahedron, ext="tetrahedron", num_refinements=0, N0=5): timer = Timer("Create mesh") @@ -82,7 +91,7 @@ def over_plane(p0, p1, p2): cell = Cell(ext, geometric_dimension=points.shape[1]) domain = Mesh(VectorElement("Lagrange", cell, 1)) # Rotate mesh - points = np.dot(r_matrix, points.T).T + points = np.dot(r_matrix, points.T).T.astype(default_real_type) mesh = create_mesh(MPI.COMM_SELF, cells, points, domain) with XDMFFile(MPI.COMM_SELF, tmp_mesh_name, "w") as xdmf: @@ -91,7 +100,6 @@ def over_plane(p0, p1, p2): MPI.COMM_WORLD.barrier() with XDMFFile(MPI.COMM_WORLD, tmp_mesh_name, "r") as xdmf: mesh = xdmf.read_mesh() - # Refine coarse mesh for i in range(num_refinements): mesh.topology.create_entities(mesh.topology.dim - 2) @@ -156,7 +164,7 @@ def over_plane(p0, p1, p2): mesh_dir.mkdir(exist_ok=True) fname = mesh_dir / f"mesh_{ext}_{theta:.2f}.xdmf" - with XDMFFile(MPI.COMM_WORLD, fname, "w") as o_f: + with XDMFFile(mesh.comm, fname, "w") as o_f: o_f.write_mesh(mesh) o_f.write_meshtags(ct, x=mesh.geometry) o_f.write_meshtags(mt, x=mesh.geometry) @@ -182,7 +190,8 @@ def demo_stacked_cubes(theta, ct, noslip, num_refinements, N0, timings=False): mesh.name = f"mesh_{celltype}_{theta:.2f}{type_ext:s}" # Create functionspaces - V = VectorFunctionSpace(mesh, ("Lagrange", 1)) + el = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 1, shape=(mesh.geometry.dim, )) + V = functionspace(mesh, el) # Define boundary conditions # Bottom boundary is fixed in all directions @@ -210,7 +219,7 @@ def top_v(x): return values u_top = Function(V) u_top.interpolate(top_v) - u_top.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD) + u_top.x.scatter_forward() top_dofs = locate_dofs_topological(V, fdim, mt.find(3)) bc_top = dirichletbc(u_top, top_dofs) @@ -218,7 +227,7 @@ def top_v(x): bcs = [bc_bottom, bc_top] # Elasticity parameters - E = PETSc.ScalarType(1.0e3) + E = default_scalar_type(1.0e3) nu = 0 mu = Constant(mesh, E / (2.0 * (1.0 + nu))) lmbda = Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) @@ -230,7 +239,7 @@ def sigma(v): u = TrialFunction(V) v = TestFunction(V) a = inner(sigma(u), grad(v)) * dx - rhs = inner(Constant(mesh, PETSc.ScalarType((0, 0, 0))), v) * dx + rhs = inner(Constant(mesh, default_scalar_type((0, 0, 0))), v) * dx log_info("Create constraints") @@ -258,12 +267,12 @@ def sigma(v): with Timer(f"{num_dofs}: Assemble-vector (C++)"): b = assemble_vector(linear_form, mpc) apply_lifting(b, [bilinear_form], [bcs], mpc) - b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore set_bc(b, bcs) list_timings(MPI.COMM_WORLD, [TimingType.wall]) # Solve Linear problem - opts = PETSc.Options() + opts = PETSc.Options() # type: ignore # opts["ksp_rtol"] = 1.0e-8 opts["pc_type"] = "gamg" # opts["pc_gamg_type"] = "agg" @@ -281,7 +290,7 @@ def sigma(v): # Create functionspace and build near nullspace A.setNearNullSpace(null_space) - solver = PETSc.KSP().create(comm) + solver = PETSc.KSP().create(comm) # type: ignore solver.setOperators(A) solver.setFromOptions() uh = Function(mpc.function_space) diff --git a/python/benchmarks/bench_elasticity.py b/python/benchmarks/bench_elasticity.py index c1db6ff6..b15ff399 100644 --- a/python/benchmarks/bench_elasticity.py +++ b/python/benchmarks/bench_elasticity.py @@ -9,11 +9,12 @@ from pathlib import Path from typing import Optional +import basix.ufl import h5py import numpy as np from dolfinx.common import Timer, TimingType, list_timings -from dolfinx.fem import (Constant, Function, VectorFunctionSpace, dirichletbc, - form, locate_dofs_topological, set_bc) +from dolfinx.fem import (Constant, Function, functionspace, dirichletbc, form, + locate_dofs_topological, set_bc) from dolfinx.io import XDMFFile from dolfinx.mesh import (create_unit_cube, locate_entities_boundary, meshtags, refine) @@ -36,7 +37,8 @@ def bench_elasticity_one(r_lvl: int = 0, out_hdf5: Optional[h5py.File] = None, mesh = refine(mesh, redistribute=True) fdim = mesh.topology.dim - 1 - V = VectorFunctionSpace(mesh, ("Lagrange", 1)) + el = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 1, shape=(mesh.geometry.dim,)) + V = functionspace(mesh, el) # Generate Dirichlet BC on lower boundary (Fixed) u_bc = Function(V) @@ -83,7 +85,7 @@ def sigma(v): # Create MPC with Timer("~Elasticity: Init constraint"): def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([1, 0, 0]): {l2b([1, 0, 1]): 0.5}} mpc = MultiPointConstraint(V) mpc.create_general_constraint(s_m_c, 2, 2) @@ -97,14 +99,14 @@ def l2b(li): b = assemble_vector(linear_form, mpc) # Apply boundary conditions apply_lifting(b, [bilinear_form], [bcs], mpc) - b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore set_bc(b, bcs) # Create functionspace and function for mpc vector # Solve Linear problem - solver = PETSc.KSP().create(MPI.COMM_WORLD) - opts = PETSc.Options() + solver = PETSc.KSP().create(mesh.comm) # type: ignore + opts = PETSc.Options() # type: ignore if boomeramg: opts["ksp_type"] = "cg" opts["ksp_rtol"] = 1.0e-5 @@ -135,7 +137,7 @@ def l2b(li): uh = b.copy() uh.set(0) solver.solve(b, uh) - uh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) + uh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) # type: ignore mpc.backsubstitution(uh) solver_time = timer.elapsed() @@ -165,7 +167,7 @@ def l2b(li): outdir = Path("results") outdir.mkdir(exist_ok=True, parents=True) fname = outdir / f"bench_elasticity_{r_lvl}.xdmf" - with XDMFFile(MPI.COMM_WORLD, fname, "w") as outfile: + with XDMFFile(mesh.comm, fname, "w") as outfile: outfile.write_mesh(mesh) outfile.write_function(u_h) diff --git a/python/benchmarks/bench_elasticity_edge.py b/python/benchmarks/bench_elasticity_edge.py index c617c022..df79e88a 100644 --- a/python/benchmarks/bench_elasticity_edge.py +++ b/python/benchmarks/bench_elasticity_edge.py @@ -6,24 +6,26 @@ import resource from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser +from pathlib import Path +import basix.ufl import h5py import numpy as np +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings -from dolfinx.fem import (Constant, Function, VectorFunctionSpace, dirichletbc, - form, locate_dofs_topological, set_bc) +from dolfinx.fem import (Constant, Function, functionspace, dirichletbc, form, + locate_dofs_topological, set_bc) from dolfinx.io import XDMFFile from dolfinx.mesh import (CellType, create_unit_cube, locate_entities_boundary, meshtags) -from dolfinx_mpc import (MultiPointConstraint, apply_lifting, assemble_matrix, - assemble_vector) -from dolfinx_mpc.utils import log_info, rigid_motions_nullspace from mpi4py import MPI from petsc4py import PETSc from ufl import (Identity, SpatialCoordinate, TestFunction, TrialFunction, as_vector, ds, dx, grad, inner, sym, tr) -from pathlib import Path +from dolfinx_mpc import (MultiPointConstraint, apply_lifting, assemble_matrix, + assemble_vector) +from dolfinx_mpc.utils import log_info, rigid_motions_nullspace def bench_elasticity_edge(tetra: bool = True, r_lvl: int = 0, out_hdf5=None, xdmf: bool = False, @@ -33,18 +35,16 @@ def bench_elasticity_edge(tetra: bool = True, r_lvl: int = 0, out_hdf5=None, xdm N *= 2 ct = CellType.tetrahedron if tetra else CellType.hexahedron mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, ct) - # Get number of unknowns on each edge - V = VectorFunctionSpace(mesh, ("Lagrange", int(degree))) + el = basix.ufl.element("Lagrange", mesh.topology.cell_name(), int(degree), shape=(mesh.geometry.dim,)) + V = functionspace(mesh, el) # Generate Dirichlet BC (Fixed) u_bc = Function(V) - with u_bc.vector.localForm() as u_local: - u_local.set(0.0) - u_bc.vector.destroy() + u_bc.x.array[:] = 0 def boundaries(x): - return np.isclose(x[0], np.finfo(float).eps) + return np.isclose(x[0], 0, 500 * np.finfo(x.dtype).resolution) fdim = mesh.topology.dim - 1 facets = locate_entities_boundary(mesh, fdim, boundaries) topological_dofs = locate_dofs_topological(V, fdim, facets) @@ -68,7 +68,7 @@ def periodic_relation(x): mpc = MultiPointConstraint(V) mpc.create_periodic_constraint_topological(V, periodic_mt, 2, periodic_relation, bcs, - scale=PETSc.ScalarType(0.5)) + scale=default_scalar_type(0.5)) mpc.finalize() # Create traction meshtag @@ -81,13 +81,13 @@ def traction_boundary(x): mt = meshtags(mesh, fdim, t_facets[arg_sort], facet_values) # Elasticity parameters - E = PETSc.ScalarType(1.0e4) + E = 1.0e4 nu = 0.1 - mu = Constant(mesh, E / (2.0 * (1.0 + nu))) - lmbda = Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) - g = Constant(mesh, PETSc.ScalarType((0, 0, -1e2))) + mu = Constant(mesh, default_scalar_type(E / (2.0 * (1.0 + nu)))) + lmbda = Constant(mesh, default_scalar_type(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))) + g = Constant(mesh, default_scalar_type((0, 0, -1e2))) x = SpatialCoordinate(mesh) - f = Constant(mesh, PETSc.ScalarType(1e3)) * as_vector((0, -(x[2] - 0.5)**2, (x[1] - 0.5)**2)) + f = Constant(mesh, default_scalar_type(1e3)) * as_vector((0, -(x[2] - 0.5)**2, (x[1] - 0.5)**2)) # Stress computation def epsilon(v): @@ -117,10 +117,10 @@ def sigma(v): # Apply boundary conditions apply_lifting(b, [bilinear_form], [bcs], mpc) - b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore set_bc(b, bcs) - opts = PETSc.Options() + opts = PETSc.Options() # type: ignore if boomeramg: opts["ksp_type"] = "cg" opts["ksp_rtol"] = 1.0e-5 @@ -145,8 +145,8 @@ def sigma(v): # opts["ksp_view"] = None # List progress of solver # Setup PETSc solver - solver = PETSc.KSP().create(MPI.COMM_WORLD) - solver.setFromOptions() + solver = PETSc.KSP().create(mesh.comm) # type: ignore + solver.setFromOptions() # type: ignore if info: log_info(f"Run {r_lvl}: Solving") @@ -186,7 +186,7 @@ def sigma(v): results = Path("results").absolute() results.mkdir(exist_ok=True) fname = results / f"bench_elasticity_edge_{r_lvl}.xdmf" - with XDMFFile(MPI.COMM_WORLD, fname, "w") as outfile: + with XDMFFile(mesh.comm, fname, "w") as outfile: outfile.write_mesh(mesh) outfile.write_function(u_h) diff --git a/python/benchmarks/bench_periodic.py b/python/benchmarks/bench_periodic.py index 73f42d54..e26dc1d3 100644 --- a/python/benchmarks/bench_periodic.py +++ b/python/benchmarks/bench_periodic.py @@ -13,24 +13,26 @@ from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser +from pathlib import Path import h5py import numpy as np +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings -from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form, +from dolfinx.fem import (Function, dirichletbc, form, functionspace, locate_dofs_geometrical, set_bc) from dolfinx.io import XDMFFile from dolfinx.mesh import (CellType, create_unit_cube, locate_entities_boundary, meshtags) -from dolfinx_mpc import (MultiPointConstraint, apply_lifting, assemble_matrix, - assemble_vector) -from dolfinx_mpc.utils import log_info from mpi4py import MPI from petsc4py import PETSc -from pathlib import Path from ufl import (SpatialCoordinate, TestFunction, TrialFunction, dx, exp, grad, inner, pi, sin) +from dolfinx_mpc import (MultiPointConstraint, apply_lifting, assemble_matrix, + assemble_vector) +from dolfinx_mpc.utils import log_info + def demo_periodic3D(tetra, r_lvl=0, out_hdf5=None, xdmf=False, boomeramg=False, kspview=False, degree=1): @@ -43,7 +45,7 @@ def demo_periodic3D(tetra, r_lvl=0, out_hdf5=None, N *= 2 mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, ct) - V = FunctionSpace(mesh, ("CG", degree)) + V = functionspace(mesh, ("CG", degree)) # Create Dirichlet boundary condition @@ -53,7 +55,7 @@ def dirichletboundary(x): mesh.topology.create_connectivity(2, 1) geometrical_dofs = locate_dofs_geometrical(V, dirichletboundary) - bc = dirichletbc(PETSc.ScalarType(0), geometrical_dofs, V) + bc = dirichletbc(default_scalar_type(0), geometrical_dofs, V) bcs = [bc] def PeriodicBoundary(x): @@ -134,7 +136,7 @@ def periodic_relation(x): # Solve linear problem log_info(f"Run {r_lvl}: Solving") - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) with Timer("~Periodic: Solve") as timer: # Create solver, set operator and options PETSc.Mat.setNearNullSpace(A, nullspace) @@ -179,7 +181,7 @@ def periodic_relation(x): results = Path("results").absolute() results.mkdir(exist_ok=True) fname = results / f"bench_periodic3d_{r_lvl}_{ext}.xdmf" - with XDMFFile(MPI.COMM_WORLD, fname, "w") as out_xdmf: + with XDMFFile(mesh.comm, fname, "w") as out_xdmf: out_xdmf.write_mesh(mesh) out_xdmf.write_function(u_h, 0.0, f"Xdmf/Domain/Grid[@Name='{mesh.name}'][1]") diff --git a/python/benchmarks/ref_elasticity.py b/python/benchmarks/ref_elasticity.py index 259e207f..f0c3a6bd 100644 --- a/python/benchmarks/ref_elasticity.py +++ b/python/benchmarks/ref_elasticity.py @@ -10,11 +10,13 @@ from time import perf_counter from typing import Optional +import basix.ufl import h5py import numpy as np +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings -from dolfinx.fem import (Constant, Function, VectorFunctionSpace, dirichletbc, - form, locate_dofs_topological) +from dolfinx.fem import (Constant, Function, dirichletbc, form, functionspace, + locate_dofs_topological) from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, set_bc) from dolfinx.io import XDMFFile @@ -47,7 +49,8 @@ def ref_elasticity(tetra: bool = True, r_lvl: int = 0, out_hdf5: Optional[h5py.F # set_log_level(LogLevel.ERROR) N = degree * N fdim = mesh.topology.dim - 1 - V = VectorFunctionSpace(mesh, ("Lagrange", int(degree))) + el = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 1, shape=(mesh.geometry.dim,)) + V = functionspace(mesh, el) # Generate Dirichlet BC on lower boundary (Fixed) u_bc = Function(V) @@ -71,13 +74,13 @@ def traction_boundary(x): mt = meshtags(mesh, fdim, t_facets[arg_sort], facet_values[arg_sort]) # Elasticity parameters - E = PETSc.ScalarType(1.0e4) + E = default_scalar_type(1.0e4) nu = 0.1 mu = Constant(mesh, E / (2.0 * (1.0 + nu))) lmbda = Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) - g = Constant(mesh, PETSc.ScalarType((0, 0, -1e2))) + g = Constant(mesh, default_scalar_type((0, 0, -1e2))) x = SpatialCoordinate(mesh) - f = Constant(mesh, PETSc.ScalarType(1e4)) * \ + f = Constant(mesh, default_scalar_type(1e4)) * \ as_vector((0, -(x[2] - 0.5)**2, (x[1] - 0.5)**2)) # Stress computation @@ -104,9 +107,9 @@ def sigma(v): linear_form = form(rhs) L_org = assemble_vector(linear_form) apply_lifting(L_org, [bilinear_form], [bcs]) - L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore set_bc(L_org, bcs) - opts = PETSc.Options() + opts = PETSc.Options() # type: ignore if boomeramg: opts["ksp_type"] = "cg" opts["ksp_rtol"] = 1.0e-5 @@ -132,7 +135,7 @@ def sigma(v): # opts["ksp_view"] = None # List progress of solver # Create solver, set operator and options - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) # type: ignore solver.setFromOptions() solver.setOperators(A_org) @@ -174,7 +177,7 @@ def sigma(v): outdir.mkdir(exist_ok=True, parents=True) fname = outdir / "ref_elasticity_{0:d}.xdmf".format(r_lvl) - with XDMFFile(MPI.COMM_WORLD, fname, "w") as out_xdmf: + with XDMFFile(mesh.comm, fname, "w") as out_xdmf: out_xdmf.write_mesh(mesh) out_xdmf.write_function(u_, 0.0, "Xdmf/Domain/Grid[@Name='{0:s}'][1]".format(mesh.name)) diff --git a/python/benchmarks/ref_periodic.py b/python/benchmarks/ref_periodic.py index bcc4fc4c..bf689248 100644 --- a/python/benchmarks/ref_periodic.py +++ b/python/benchmarks/ref_periodic.py @@ -24,8 +24,9 @@ import h5py import numpy as np +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings -from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form, +from dolfinx.fem import (Function, functionspace, dirichletbc, form, locate_dofs_geometrical) from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, set_bc) @@ -57,7 +58,7 @@ def reference_periodic(tetra: bool, r_lvl: int = 0, out_hdf5: Optional[h5py.File N *= 2 mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, CellType.hexahedron) - V = FunctionSpace(mesh, ("CG", degree)) + V = functionspace(mesh, ("CG", degree)) # Create Dirichlet boundary condition @@ -69,7 +70,7 @@ def dirichletboundary(x): mesh.topology.create_connectivity(2, 1) geometrical_dofs = locate_dofs_geometrical(V, dirichletboundary) - bc = dirichletbc(PETSc.ScalarType(0), geometrical_dofs, V) + bc = dirichletbc(default_scalar_type(0), geometrical_dofs, V) bcs = [bc] # Define variational problem @@ -90,15 +91,15 @@ def dirichletboundary(x): A_org.assemble() L_org = assemble_vector(linear_form) apply_lifting(L_org, [bilinear_form], [bcs]) - L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore set_bc(L_org, bcs) # Create PETSc nullspace - nullspace = PETSc.NullSpace().create(constant=True) - PETSc.Mat.setNearNullSpace(A_org, nullspace) + nullspace = PETSc.NullSpace().create(constant=True) # type: ignore + PETSc.Mat.setNearNullSpace(A_org, nullspace) # type: ignore # Set PETSc options - opts = PETSc.Options() + opts = PETSc.Options() # type: ignore if boomeramg: opts["ksp_type"] = "cg" opts["ksp_rtol"] = 1.0e-5 @@ -121,7 +122,7 @@ def dirichletboundary(x): # opts["ksp_view"] = None # List progress of solver # Initialize PETSc solver, set options and operator - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) # type: ignore solver.setFromOptions() solver.setOperators(A_org) @@ -131,8 +132,8 @@ def dirichletboundary(x): with Timer("Solve"): solver.solve(L_org, u_.vector) end = perf_counter() - u_.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, - mode=PETSc.ScatterMode.FORWARD) + u_.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, # type: ignore + mode=PETSc.ScatterMode.FORWARD) # type: ignore if kspview: solver.view() @@ -157,7 +158,7 @@ def dirichletboundary(x): fname = outdir / "reference_periodic_{0:d}_{1:s}.xdmf".format( r_lvl, ext) u_.name = "u_" + ext + "_unconstrained" - with XDMFFile(MPI.COMM_WORLD, fname, "w") as out_periodic: + with XDMFFile(mesh.comm, fname, "w") as out_periodic: out_periodic.write_mesh(mesh) out_periodic.write_function(u_, 0.0, "Xdmf/Domain/" diff --git a/python/demos/demo_contact_2D.py b/python/demos/demo_contact_2D.py index 0891e87a..d1cfa4da 100644 --- a/python/demos/demo_contact_2D.py +++ b/python/demos/demo_contact_2D.py @@ -13,14 +13,16 @@ # added to the to left corner of the top cube. +import warnings from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser from pathlib import Path import numpy as np import scipy.sparse.linalg from create_and_export_mesh import gmsh_2D_stacked, mesh_2D_dolfin +from dolfinx import default_scalar_type, default_real_type from dolfinx.common import Timer, TimingType, list_timings -from dolfinx.fem import (Constant, VectorFunctionSpace, dirichletbc, form, +from dolfinx.fem import (Constant, dirichletbc, form, functionspace, locate_dofs_geometrical) from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, set_bc) @@ -55,6 +57,9 @@ def demo_stacked_cubes(outfile: XDMFFile, theta: float, gmsh: bool = True, quad: mesh.name = f"mesh_{celltype}_{theta:.2f}_gmsh" else: + if default_real_type == np.float32: + warnings.warn("Demo does not run for single float precision due to limited xdmf support") + exit(0) mesh_name = "mesh" filename = meshdir / f"mesh_{celltype}_{theta:.2f}.xdmf" @@ -69,26 +74,25 @@ def demo_stacked_cubes(outfile: XDMFFile, theta: float, gmsh: bool = True, quad: mt = xdmf.read_meshtags(mesh, name="facet_tags") # Helper until meshtags can be read in from xdmf - V = VectorFunctionSpace(mesh, ("Lagrange", 1)) - + V = functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) r_matrix = rotation_matrix([0, 0, 1], theta) g_vec = np.dot(r_matrix, [0, -1.25e2, 0]) - g = Constant(mesh, PETSc.ScalarType(g_vec[:2])) + g = Constant(mesh, default_scalar_type(g_vec[:2])) def bottom_corner(x): - return np.isclose(x, [[0], [0], [0]]).all(axis=0) + return np.isclose(x, [[0], [0], [0]], atol=5e2 * np.finfo(default_scalar_type).resolution).all(axis=0) # Fix bottom corner - bc_value = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType) + bc_value = np.array((0,) * mesh.geometry.dim, dtype=default_scalar_type) # type: ignore bottom_dofs = locate_dofs_geometrical(V, bottom_corner) bc_bottom = dirichletbc(bc_value, bottom_dofs, V) bcs = [bc_bottom] # Elasticity parameters - E = PETSc.ScalarType(1.0e3) + E = 1.0e3 nu = 0 - mu = Constant(mesh, E / (2.0 * (1.0 + nu))) - lmbda = Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) + mu = Constant(mesh, default_scalar_type(E / (2.0 * (1.0 + nu)))) + lmbda = Constant(mesh, default_scalar_type(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))) # Stress computation def sigma(v): @@ -99,17 +103,18 @@ def sigma(v): v = TestFunction(V) a = inner(sigma(u), grad(v)) * dx ds = Measure("ds", domain=mesh, subdomain_data=mt, subdomain_id=3) - rhs = inner(Constant(mesh, PETSc.ScalarType((0, 0))), v) * dx + inner(g, v) * ds + rhs = inner(Constant(mesh, default_scalar_type((0, 0))), v) * dx + inner(g, v) * ds # type: ignore + tol = float(5e2 * np.finfo(default_scalar_type).resolution) def left_corner(x): - return np.isclose(x.T, np.dot(r_matrix, [0, 2, 0])).all(axis=1) + return np.isclose(x.T, np.dot(r_matrix, [0, 2, 0]), atol=tol).all(axis=1) # Create multi point constraint mpc = MultiPointConstraint(V) with Timer("~Contact: Create contact constraint"): nh = create_normal_approximation(V, mt, 4) - mpc.create_contact_slip_condition(mt, 4, 9, nh) + mpc.create_contact_slip_condition(mt, 4, 9, nh, eps2=tol) with Timer("~Contact: Add non-slip condition at bottom interface"): bottom_normal = facet_normal_approximation(V, mt, 5) @@ -123,8 +128,8 @@ def left_corner(x): mpc.create_slip_constraint(V, (mtv, 6), tangent, bcs=bcs) mpc.finalize() - rtol = 3e-8 - petsc_options = {"ksp_rtol": 1e-9, "ksp_atol": 1e-9, "pc_type": "gamg", + tol = float(5e2 * np.finfo(default_scalar_type).resolution) + petsc_options = {"ksp_rtol": tol, "ksp_atol": tol, "pc_type": "gamg", "pc_gamg_type": "agg", "pc_gamg_square_graph": 2, "pc_gamg_threshold": 0.02, "pc_gamg_coarse_eq_limit": 1000, "pc_gamg_sym_graph": True, "mg_levels_ksp_type": "chebyshev", "mg_levels_pc_type": "jacobi", @@ -166,7 +171,7 @@ def left_corner(x): A_org.assemble() L_org = assemble_vector(form(rhs)) apply_lifting(L_org, [form(a)], [bcs]) - L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore set_bc(L_org, bcs) root = 0 @@ -186,7 +191,7 @@ def left_corner(x): d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc, rtol=rtol, atol=rtol) + assert np.allclose(uh_numpy, u_mpc, rtol=tol, atol=tol) L_org.destroy() diff --git a/python/demos/demo_contact_3D.py b/python/demos/demo_contact_3D.py index eb5f883f..9aa1bc58 100644 --- a/python/demos/demo_contact_3D.py +++ b/python/demos/demo_contact_3D.py @@ -8,6 +8,7 @@ # between two cubes. +import warnings from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser from pathlib import Path @@ -15,6 +16,7 @@ import numpy as np import scipy.sparse.linalg from create_and_export_mesh import gmsh_3D_stacked, mesh_3D_dolfin +from dolfinx import default_real_type, default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.io import XDMFFile from dolfinx.mesh import CellType @@ -46,6 +48,9 @@ def demo_stacked_cubes(outfile: XDMFFile, theta: float, gmsh: bool = False, ct: mesh.topology.create_connectivity(tdim, tdim) mesh.topology.create_connectivity(fdim, tdim) else: + if default_real_type == np.float32: + warnings.warn("Demo does not run for single float precision due to limited xdmf support") + exit(0) mesh_3D_dolfin(theta, ct, celltype, res) MPI.COMM_WORLD.barrier() with XDMFFile(MPI.COMM_WORLD, f"meshes/mesh_{celltype}_{theta:.2f}.xdmf", "r") as xdmf: @@ -58,22 +63,22 @@ def demo_stacked_cubes(outfile: XDMFFile, theta: float, gmsh: bool = False, ct: mesh.name = f"mesh_{celltype}_{theta:.2f}{type_ext}{mesh_ext}" # Create functionspaces - V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) # Define boundary conditions # Bottom boundary is fixed in all directions bottom_dofs = fem.locate_dofs_topological(V, fdim, mt.find(5)) - u_bc = np.array((0, ) * mesh.geometry.dim, dtype=PETSc.ScalarType) + u_bc = np.array((0, ) * mesh.geometry.dim, dtype=default_scalar_type) bc_bottom = fem.dirichletbc(u_bc, bottom_dofs, V) - g_vec = np.array([0, 0, -4.25e-1], dtype=PETSc.ScalarType) + g_vec = np.array([0, 0, -4.25e-1], dtype=default_scalar_type) if not noslip: # Helper for orienting traction r_matrix = rotation_matrix([1 / np.sqrt(2), 1 / np.sqrt(2), 0], -theta) # Top boundary has a given deformation normal to the interface - g_vec = np.dot(r_matrix, g_vec) + g_vec = np.dot(r_matrix, g_vec).astype(default_scalar_type) top_dofs = fem.locate_dofs_topological(V, fdim, mt.find(3)) bc_top = fem.dirichletbc(g_vec, top_dofs, V) @@ -81,10 +86,10 @@ def demo_stacked_cubes(outfile: XDMFFile, theta: float, gmsh: bool = False, ct: bcs = [bc_bottom, bc_top] # Elasticity parameters - E = PETSc.ScalarType(1.0e3) + E = 1.0e3 nu = 0 - mu = fem.Constant(mesh, E / (2.0 * (1.0 + nu))) - lmbda = fem.Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) + mu = fem.Constant(mesh, default_scalar_type(E / (2.0 * (1.0 + nu)))) + lmbda = fem.Constant(mesh, default_scalar_type(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))) # Stress computation def sigma(v): @@ -95,21 +100,22 @@ def sigma(v): v = TestFunction(V) a = inner(sigma(u), grad(v)) * dx # NOTE: Traction deactivated until we have a way of fixing nullspace - # g = fem.Constant(mesh, PETSc.ScalarType(g_vec)) + # g = fem.Constant(mesh, default_scalar_type(g_vec)) # ds = Measure("ds", domain=mesh, subdomain_data=mt, subdomain_id=3) - rhs = inner(fem.Constant(mesh, PETSc.ScalarType((0, 0, 0))), v) * dx + rhs = inner(fem.Constant(mesh, default_scalar_type((0, 0, 0))), v) * dx # + inner(g, v) * ds bilinear_form = fem.form(a) linear_form = fem.form(rhs) mpc = MultiPointConstraint(V) + tol = float(5e2 * np.finfo(default_scalar_type).resolution) if noslip: with Timer("~~Contact: Create non-elastic constraint"): - mpc.create_contact_inelastic_condition(mt, 4, 9) + mpc.create_contact_inelastic_condition(mt, 4, 9, eps2=tol) else: with Timer("~Contact: Create contact constraint"): nh = create_normal_approximation(V, mt, 4) - mpc.create_contact_slip_condition(mt, 4, 9, nh) + mpc.create_contact_slip_condition(mt, 4, 9, nh, eps2=tol) with Timer("~~Contact: Add data and finialize MPC"): mpc.finalize() @@ -123,11 +129,11 @@ def sigma(v): b = assemble_vector(linear_form, mpc) apply_lifting(b, [bilinear_form], [bcs], mpc) - b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore fem.petsc.set_bc(b, bcs) # Solve Linear problem - opts = PETSc.Options() + opts = PETSc.Options() # type: ignore opts["ksp_rtol"] = 1.0e-8 opts["pc_type"] = "gamg" opts["pc_gamg_type"] = "agg" @@ -144,7 +150,7 @@ def sigma(v): # Create functionspace and build near nullspace A.setNearNullSpace(null_space) - solver = PETSc.KSP().create(mesh.comm) + solver = PETSc.KSP().create(mesh.comm) # type: ignore solver.setOperators(A) solver.setFromOptions() @@ -190,7 +196,7 @@ def sigma(v): A_org.assemble() L_org = fem.petsc.assemble_vector(linear_form) fem.petsc.apply_lifting(L_org, [bilinear_form], [bcs]) - L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore fem.petsc.set_bc(L_org, bcs) root = 0 diff --git a/python/demos/demo_elasticity.py b/python/demos/demo_elasticity.py index 722b81f7..c64f0454 100644 --- a/python/demos/demo_elasticity.py +++ b/python/demos/demo_elasticity.py @@ -8,6 +8,7 @@ import dolfinx.fem as fem import numpy as np import scipy.sparse.linalg +from dolfinx import default_scalar_type from dolfinx.common import Timer from dolfinx.io import XDMFFile from dolfinx.mesh import create_unit_square, locate_entities_boundary @@ -23,7 +24,7 @@ def demo_elasticity(): mesh = create_unit_square(MPI.COMM_WORLD, 10, 10) - V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) # Generate Dirichlet BC on lower boundary (Fixed) @@ -31,7 +32,7 @@ def boundaries(x): return np.isclose(x[0], np.finfo(float).eps) facets = locate_entities_boundary(mesh, 1, boundaries) topological_dofs = fem.locate_dofs_topological(V, 1, facets) - bc = fem.dirichletbc(np.array([0, 0], dtype=PETSc.ScalarType), topological_dofs, V) + bc = fem.dirichletbc(np.array([0, 0], dtype=default_scalar_type), topological_dofs, V) bcs = [bc] # Define variational problem @@ -39,10 +40,10 @@ def boundaries(x): v = TestFunction(V) # Elasticity parameters - E = PETSc.ScalarType(1.0e4) + E = 1.0e4 nu = 0.0 - mu = fem.Constant(mesh, E / (2.0 * (1.0 + nu))) - lmbda = fem.Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) + mu = fem.Constant(mesh, default_scalar_type(E / (2.0 * (1.0 + nu)))) + lmbda = fem.Constant(mesh, default_scalar_type(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))) # Stress computation def sigma(v): @@ -58,7 +59,7 @@ def sigma(v): # Create MPC def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([1, 0]): {l2b([1, 1]): 0.9}} mpc = MultiPointConstraint(V) mpc.create_general_constraint(s_m_c, 1, 1) @@ -72,7 +73,7 @@ def l2b(li): outdir = Path("results") outdir.mkdir(exist_ok=True, parents=True) - with XDMFFile(MPI.COMM_WORLD, outdir / "demo_elasticity.xdmf", "w") as outfile: + with XDMFFile(mesh.comm, outdir / "demo_elasticity.xdmf", "w") as outfile: outfile.write_mesh(mesh) outfile.write_function(u_h) @@ -86,7 +87,7 @@ def l2b(li): fem.petsc.apply_lifting(L_org, [bilinear_form], [bcs]) L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) fem.petsc.set_bc(L_org, bcs) - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) solver.setType(PETSc.KSP.Type.PREONLY) solver.getPC().setType(PETSc.PC.Type.LU) solver.setOperators(A_org) @@ -95,7 +96,7 @@ def l2b(li): u_.x.scatter_forward() u_.name = "u_unconstrained" - with XDMFFile(MPI.COMM_WORLD, outdir / "demo_elasticity.xdmf", "a") as outfile: + with XDMFFile(mesh.comm, outdir / "demo_elasticity.xdmf", "a") as outfile: outfile.write_function(u_) outfile.close() @@ -117,7 +118,7 @@ def l2b(li): d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + assert np.allclose(uh_numpy, u_mpc, atol=500 * np.finfo(u_mpc.dtype).resolution) # Print out master-slave connectivity for the first slave master_owner = None diff --git a/python/demos/demo_elasticity_disconnect.py b/python/demos/demo_elasticity_disconnect.py index d4f6333b..ffc5f229 100644 --- a/python/demos/demo_elasticity_disconnect.py +++ b/python/demos/demo_elasticity_disconnect.py @@ -12,12 +12,11 @@ import basix.ufl import gmsh import numpy as np -from dolfinx.fem import (Constant, Function, FunctionSpace, - VectorFunctionSpace, dirichletbc, +from dolfinx import default_scalar_type +from dolfinx.fem import (Constant, Function, dirichletbc, functionspace, locate_dofs_topological) from dolfinx.io import XDMFFile, gmshio from mpi4py import MPI -from petsc4py import PETSc from ufl import (Identity, Measure, SpatialCoordinate, TestFunction, TrialFunction, as_vector, grad, inner, sym, tr) @@ -107,12 +106,12 @@ gmsh.finalize() MPI.COMM_WORLD.barrier() -V = VectorFunctionSpace(mesh, ("Lagrange", 1)) +V = functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) tdim = mesh.topology.dim fdim = tdim - 1 -DG0 = FunctionSpace(mesh, ("DG", 0)) +DG0 = functionspace(mesh, ("DG", 0)) outer_dofs = locate_dofs_topological(DG0, tdim, ct.find(outer_tag)) inner_dofs = locate_dofs_topological(DG0, tdim, ct.find(inner_tag)) @@ -146,16 +145,16 @@ def sigma(v): dx = Measure("dx", domain=mesh, subdomain_data=ct) a = inner(sigma(u), grad(v)) * dx x = SpatialCoordinate(mesh) -rhs = inner(Constant(mesh, PETSc.ScalarType((0, 0, 0))), v) * dx -rhs += inner(Constant(mesh, PETSc.ScalarType((0.01, 0.02, 0))), v) * dx(outer_tag) -rhs += inner(as_vector(PETSc.ScalarType((0, 0, -9.81e-2))), v) * dx(inner_tag) +rhs = inner(Constant(mesh, default_scalar_type((0, 0, 0))), v) * dx +rhs += inner(Constant(mesh, default_scalar_type((0.01, 0.02, 0))), v) * dx(outer_tag) +rhs += inner(as_vector((0, 0, -9.81e-2)), v) * dx(inner_tag) # Create dirichletbc owning_processor, bc_dofs = determine_closest_block(V, -np.array([-r2, 0, 0])) bc_dofs = [] if bc_dofs is None else bc_dofs -u_fixed = np.array([0, 0, 0], dtype=PETSc.ScalarType) +u_fixed = np.array([0, 0, 0], dtype=default_scalar_type) bc_fixed = dirichletbc(u_fixed, np.asarray(bc_dofs, dtype=np.int32), V) bcs = [bc_fixed] @@ -176,7 +175,8 @@ def sigma(v): # Create nullspace null_space = rigid_motions_nullspace(mpc.function_space) -petsc_options = {"ksp_rtol": 1.0e-8, "pc_type": "gamg", "pc_gamg_type": "agg", +ksp_rtol = 5e2 * np.finfo(default_scalar_type).resolution +petsc_options = {"ksp_rtol": ksp_rtol, "pc_type": "gamg", "pc_gamg_type": "agg", "pc_gamg_coarse_eq_limit": 1000, "pc_gamg_sym_graph": True, "mg_levels_ksp_type": "chebyshev", "mg_levels_pc_type": "jacobi", "mg_levels_esteig_ksp_type": "cg", "matptap_via": "scalable", @@ -197,14 +197,15 @@ def sigma(v): print("Number of iterations: {0:d}".format(it)) # Write solution to file -V_out = FunctionSpace(mesh, basix.ufl.element("Lagrange", mesh.topology.cell_name(), - mesh.geometry.cmaps[0].degree, mesh.geometry.cmaps[0].variant, +V_out = functionspace(mesh, basix.ufl.element("Lagrange", mesh.topology.cell_name(), + mesh.geometry.cmaps[0].degree, + lagrange_variant=basix.LagrangeVariant(mesh.geometry.cmaps[0].variant), gdim=mesh.geometry.dim, shape=(V.dofmap.bs,))) u_out = Function(V_out) u_out.interpolate(u_h) u_out.name = "uh" out_path = Path("results") out_path.mkdir(exist_ok=True, parents=True) -with XDMFFile(MPI.COMM_WORLD, out_path / "demo_elasticity_disconnect.xdmf", "w") as xdmf: +with XDMFFile(mesh.comm, out_path / "demo_elasticity_disconnect.xdmf", "w") as xdmf: xdmf.write_mesh(mesh) xdmf.write_function(u_out) diff --git a/python/demos/demo_elasticity_disconnect_2D.py b/python/demos/demo_elasticity_disconnect_2D.py index 8b8732fd..7c3261a8 100644 --- a/python/demos/demo_elasticity_disconnect_2D.py +++ b/python/demos/demo_elasticity_disconnect_2D.py @@ -10,19 +10,20 @@ import gmsh import numpy as np -from dolfinx.fem import (Constant, Function, VectorFunctionSpace, dirichletbc, +from dolfinx import default_scalar_type +from dolfinx.fem import (Constant, Function, FunctionSpaceBase, dirichletbc, functionspace, locate_dofs_geometrical, - locate_dofs_topological, FunctionSpaceBase) + locate_dofs_topological) from dolfinx.io import XDMFFile, gmshio from dolfinx.mesh import locate_entities_boundary -from dolfinx_mpc import LinearProblem, MultiPointConstraint -from dolfinx_mpc.utils import (create_point_to_point_constraint, - rigid_motions_nullspace) from mpi4py import MPI -from petsc4py import PETSc from ufl import (Identity, Measure, SpatialCoordinate, TestFunction, TrialFunction, grad, inner, sym, tr) +from dolfinx_mpc import LinearProblem, MultiPointConstraint +from dolfinx_mpc.utils import (create_point_to_point_constraint, + rigid_motions_nullspace) + # Mesh parameters for creating a mesh consisting of two disjoint rectangles right_tag = 1 left_tag = 2 @@ -57,7 +58,7 @@ with XDMFFile(mesh.comm, "test.xdmf", "w") as xdmf: xdmf.write_mesh(mesh) -V = VectorFunctionSpace(mesh, ("Lagrange", 1)) +V = functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) tdim = mesh.topology.dim fdim = tdim - 1 @@ -96,13 +97,13 @@ def sigma(v): dx = Measure("dx", domain=mesh, subdomain_data=ct) a = inner(sigma(u), grad(v)) * dx x = SpatialCoordinate(mesh) -rhs = inner(Constant(mesh, PETSc.ScalarType((0, 0))), v) * dx +rhs = inner(Constant(mesh, default_scalar_type((0, 0))), v) * dx # Set boundary conditions -u_push = np.array([0.1, 0], dtype=PETSc.ScalarType) +u_push = np.array([0.1, 0], dtype=default_scalar_type) dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0)) bc_push = dirichletbc(u_push, dofs, V) -u_fix = np.array([0, 0], dtype=PETSc.ScalarType) +u_fix = np.array([0, 0], dtype=default_scalar_type) bc_fix = dirichletbc(u_fix, locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 2.1)), V) bcs = [bc_push, bc_fix] @@ -118,7 +119,7 @@ def gather_dof_coordinates(V: FunctionSpaceBase, dofs: np.ndarray): glob_num_nodes = MPI.COMM_WORLD.allreduce(num_nodes, op=MPI.SUM) recvbuf = None if MPI.COMM_WORLD.rank == 0: - recvbuf = np.zeros(3 * glob_num_nodes, dtype=np.float64) + recvbuf = np.zeros(3 * glob_num_nodes, dtype=V.mesh.geometry.x.dtype) sendbuf = coords.reshape(-1) sendcounts = np.array(MPI.COMM_WORLD.gather(len(sendbuf), 0)) MPI.COMM_WORLD.Gatherv(sendbuf, (recvbuf, sendcounts), root=0) @@ -187,6 +188,6 @@ def gather_dof_coordinates(V: FunctionSpaceBase, dofs: np.ndarray): u_h.name = "u" outdir = Path("results") outdir.mkdir(exist_ok=True, parents=True) -with XDMFFile(MPI.COMM_WORLD, outdir / "demo_elasticity_disconnect_2D.xdmf", "w") as xdmf: +with XDMFFile(mesh.comm, outdir / "demo_elasticity_disconnect_2D.xdmf", "w") as xdmf: xdmf.write_mesh(mesh) xdmf.write_function(u_h) diff --git a/python/demos/demo_periodic3d_topological.py b/python/demos/demo_periodic3d_topological.py index 298c4c4b..80e64dfc 100644 --- a/python/demos/demo_periodic3d_topological.py +++ b/python/demos/demo_periodic3d_topological.py @@ -22,13 +22,13 @@ import dolfinx.fem as fem import numpy as np import scipy.sparse.linalg +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.io import VTXWriter from dolfinx.mesh import (CellType, create_unit_cube, locate_entities_boundary, meshtags) from mpi4py import MPI from numpy.typing import NDArray -from petsc4py import PETSc from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx, exp, grad, inner, pi, sin) @@ -36,7 +36,7 @@ from dolfinx_mpc import LinearProblem # Get PETSc int and scalar types -complex_mode = True if np.dtype(PETSc.ScalarType).kind == 'c' else False +complex_mode = True if np.dtype(default_scalar_type).kind == 'c' else False def demo_periodic3D(celltype: CellType): @@ -45,19 +45,19 @@ def demo_periodic3D(celltype: CellType): # Tet setup N = 10 mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N) - V = fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) else: # Hex setup N = 10 mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, CellType.hexahedron) - V = fem.VectorFunctionSpace(mesh, ("CG", 2)) + V = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, ))) - def dirichletboundary(x: NDArray[np.float64]) -> NDArray[np.bool_]: + def dirichletboundary(x: NDArray[Union[np.float32, np.float64]]) -> NDArray[np.bool_]: return np.logical_or(np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)), np.logical_or(np.isclose(x[2], 0), np.isclose(x[2], 1))) # Create Dirichlet boundary condition - zero = PETSc.ScalarType([0, 0, 0]) + zero = default_scalar_type([0, 0, 0]) geometrical_dofs = fem.locate_dofs_geometrical(V, dirichletboundary) bc = fem.dirichletbc(zero, geometrical_dofs, V) bcs = [bc] @@ -77,7 +77,7 @@ def periodic_relation(x): return out_x with Timer("~~Periodic: Compute mpc condition"): mpc = dolfinx_mpc.MultiPointConstraint(V) - mpc.create_periodic_constraint_topological(V.sub(0), mt, 2, periodic_relation, bcs, PETSc.ScalarType(1)) + mpc.create_periodic_constraint_topological(V.sub(0), mt, 2, periodic_relation, bcs, default_scalar_type(1)) mpc.finalize() # Define variational problem u = TrialFunction(V) @@ -94,14 +94,14 @@ def periodic_relation(x): rhs = inner(f, v) * dx petsc_options: Dict[str, Union[str, float, int]] - if complex_mode: - rtol = 1e-16 + tol = float(5e2 * np.finfo(default_scalar_type).resolution) + if complex_mode or default_scalar_type == np.float32: petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} else: - rtol = 1e-8 - petsc_options = {"ksp_type": "cg", "ksp_rtol": rtol, "pc_type": "hypre", "pc_hypre_type": "boomeramg", + petsc_options = {"ksp_type": "cg", "ksp_rtol": str(tol), "pc_type": "hypre", "pc_hypre_type": "boomeramg", "pc_hypre_boomeramg_max_iter": 1, "pc_hypre_boomeramg_cycle_type": "v", "pc_hypre_boomeramg_print_statistics": 1} + problem = LinearProblem(a, rhs, mpc, bcs, petsc_options=petsc_options) u_h = problem.solve() @@ -130,7 +130,7 @@ def periodic_relation(x): outdir = Path("results") outdir.mkdir(exist_ok=True, parents=True) fname = outdir / f"demo_periodic3d_{ext}.bp" - out_periodic = VTXWriter(MPI.COMM_WORLD, fname, u_out) + out_periodic = VTXWriter(mesh.comm, fname, u_out, engine="BP4") out_periodic.write(0) out_periodic.close() @@ -152,7 +152,7 @@ def periodic_relation(x): d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc, rtol=rtol) + assert np.allclose(uh_numpy, u_mpc, rtol=tol, atol=tol) if __name__ == "__main__": diff --git a/python/demos/demo_periodic_geometrical.py b/python/demos/demo_periodic_geometrical.py index 0439f6d7..d9cb3d64 100644 --- a/python/demos/demo_periodic_geometrical.py +++ b/python/demos/demo_periodic_geometrical.py @@ -14,10 +14,10 @@ from pathlib import Path from typing import Union - import dolfinx.fem as fem import numpy as np import scipy.sparse.linalg +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.io import XDMFFile from dolfinx.mesh import create_unit_square, locate_entities_boundary @@ -30,13 +30,13 @@ from dolfinx_mpc import LinearProblem, MultiPointConstraint # Get PETSc int and scalar types -complex_mode = True if np.dtype(PETSc.ScalarType).kind == 'c' else False +complex_mode = True if np.dtype(default_scalar_type).kind == 'c' else False # Create mesh and finite element NX = 50 NY = 100 mesh = create_unit_square(MPI.COMM_WORLD, NX, NY) -V = fem.VectorFunctionSpace(mesh, ("CG", 1)) +V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) def dirichletboundary(x): @@ -46,7 +46,7 @@ def dirichletboundary(x): # Create Dirichlet boundary condition facets = locate_entities_boundary(mesh, 1, dirichletboundary) topological_dofs = fem.locate_dofs_topological(V, 1, facets) -zero = np.array([0, 0], dtype=PETSc.ScalarType) +zero = np.array([0, 0], dtype=default_scalar_type) bc = fem.dirichletbc(zero, topological_dofs, V) bcs = [bc] @@ -92,7 +92,7 @@ def periodic_relation(x): solver.setOptionsPrefix(solver_prefix) petsc_options: dict[str, Union[str, int, float]] -if complex_mode: +if complex_mode or default_scalar_type == np.float32: petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} else: petsc_options = {"ksp_type": "cg", "ksp_rtol": 1e-6, "pc_type": "hypre", "pc_hypre_type": "boomeramg", @@ -101,7 +101,7 @@ def periodic_relation(x): } # Set PETSc options -opts = PETSc.Options() +opts = PETSc.Options() # type: ignore opts.prefixPush(solver_prefix) if petsc_options is not None: for k, v in petsc_options.items(): @@ -121,7 +121,7 @@ def periodic_relation(x): outdir.mkdir(exist_ok=True, parents=True) uh.name = "u_mpc" -outfile = XDMFFile(MPI.COMM_WORLD, outdir / "demo_periodic_geometrical.xdmf", "w") +outfile = XDMFFile(mesh.comm, outdir / "demo_periodic_geometrical.xdmf", "w") outfile.write_mesh(mesh) outfile.write_function(uh) @@ -134,7 +134,7 @@ def periodic_relation(x): linear_form = fem.form(rhs) L_org = fem.petsc.assemble_vector(linear_form) fem.petsc.apply_lifting(L_org, [bilinear_form], [bcs]) -L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) +L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore fem.petsc.set_bc(L_org, bcs) solver.setOperators(A_org) u_ = fem.Function(V) @@ -165,6 +165,6 @@ def periodic_relation(x): d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + assert np.allclose(uh_numpy, u_mpc, atol=float(np.finfo(u_mpc.dtype).resolution)) list_timings(MPI.COMM_WORLD, [TimingType.wall]) L_org.destroy() diff --git a/python/demos/demo_periodic_gep.py b/python/demos/demo_periodic_gep.py index 2babfe8a..a8562900 100755 --- a/python/demos/demo_periodic_gep.py +++ b/python/demos/demo_periodic_gep.py @@ -32,6 +32,7 @@ import dolfinx.fem as fem import numpy as np +from dolfinx import default_scalar_type from dolfinx.io import XDMFFile from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags from mpi4py import MPI @@ -96,15 +97,16 @@ def EPS_print_results(EPS: SLEPc.EPS): for i in range(nconv): k = EPS.getEigenpair(i, vr, vi) error = EPS.computeError(i) - if k.imag != 0.0: + if not np.isclose(k.imag, 0.0): print0(f" {k.real:2.2e} + {k.imag:2.2e}j {error:1.1e}") else: pad = " " * 11 print0(f" {k.real:2.2e} {pad} {error:1.1e}") -def EPS_get_spectrum(EPS: SLEPc.EPS, - mpc: MultiPointConstraint) -> Tuple[List[complex], List[PETSc.Vec], List[PETSc.Vec]]: +def EPS_get_spectrum( + EPS: SLEPc.EPS, + mpc: MultiPointConstraint) -> Tuple[List[complex], List[PETSc.Vec], List[PETSc.Vec]]: # type: ignore """ Retrieve eigenvalues and eigenfunctions from SLEPc EPS object. Parameters ---------- @@ -137,7 +139,7 @@ def EPS_get_spectrum(EPS: SLEPc.EPS, return (eigval, eigvec_r, eigvec_i) -def solve_GEP_shiftinvert(A: PETSc.Mat, B: PETSc.Mat, +def solve_GEP_shiftinvert(A: PETSc.Mat, B: PETSc.Mat, # type: ignore #type: ignore problem_type: SLEPc.EPS.ProblemType = SLEPc.EPS.ProblemType.GNHEP, solver: SLEPc.EPS.Type = SLEPc.EPS.Type.KRYLOVSCHUR, nev: int = 10, tol: float = 1e-7, max_it: int = 10, @@ -225,7 +227,7 @@ def assemble_and_solve(boundary_condition: List[str] = ["dirichlet", "periodic"] # Create mesh and finite element N = 50 mesh = create_unit_square(comm, N, N) - V = fem.FunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1)) fdim = mesh.topology.dim - 1 bcs = [] @@ -321,11 +323,13 @@ def pbc_slave_to_master_map(x): # which is far from the region of interest. diagval_A = 1e2 diagval_B = 1e-2 + tol = float(5e2 * np.finfo(default_scalar_type).resolution) + A = assemble_matrix(mass_form, mpc, bcs=bcs, diagval=diagval_A) B = assemble_matrix(stiffness_form, mpc, bcs=bcs, diagval=diagval_B) EPS = solve_GEP_shiftinvert(A, B, problem_type=SLEPc.EPS.ProblemType.GHEP, solver=SLEPc.EPS.Type.KRYLOVSCHUR, - nev=Nev, tol=1e-7, max_it=10, + nev=Nev, tol=tol, max_it=10, target=1.5, shift=1.5, comm=comm) (eigval, eigvec_r, eigvec_i) = EPS_get_spectrum(EPS, mpc) # update slave DoF diff --git a/python/demos/demo_stokes.py b/python/demos/demo_stokes.py index 68c2d300..4dc5c691 100644 --- a/python/demos/demo_stokes.py +++ b/python/demos/demo_stokes.py @@ -1,31 +1,39 @@ -# Copyright (C) 2022 Jørgen S. Dokken +# # Stokes flow with slip conditions +# **Author** Jørgen S. Dokken # -# This file is part of DOLFINX_MPC -# -# SPDX-License-Identifier: MIT +# **License** MIT # # This demo illustrates how to apply a slip condition on an # interface not aligned with the coordiante axis. -# The demos solves the Stokes problem +# We start by the various modules required for this demo + +# + from pathlib import Path +from typing import Union +import basix.ufl +import dolfinx_mpc.utils import gmsh import numpy as np import scipy.sparse.linalg -from dolfinx import fem, io -from dolfinx.common import Timer, TimingType, list_timings +from dolfinx import common, default_scalar_type, fem, io from dolfinx.io import gmshio from mpi4py import MPI from numpy.typing import NDArray from petsc4py import PETSc -from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunctions, - TrialFunctions, VectorElement, div, dot, dx, grad, inner, - outer, sym) +from ufl import (FacetNormal, Identity, Measure, TestFunctions, TrialFunctions, + div, dot, dx, grad, inner, outer, sym) from ufl.core.expr import Expr -import dolfinx_mpc.utils from dolfinx_mpc import LinearProblem, MultiPointConstraint +# - +# ## Mesh generation +# +# Next we create the computational domain, a channel titled with respect to the coordinate axis. +# We use GMSH and the DOLFINx GMSH-IO to convert the GMSH model into a DOLFINx mesh + +# + def create_mesh_gmsh(L: int = 2, H: int = 1, res: float = 0.1, theta: float = np.pi / 5, @@ -105,41 +113,80 @@ def create_mesh_gmsh(L: int = 2, H: int = 1, res: float = 0.1, theta: float = np return mesh, ft -# ------------------- Mesh and function space creation ------------------------ mesh, mt = create_mesh_gmsh(res=0.1) - fdim = mesh.topology.dim - 1 -# Create the function space -P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2) -P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1) -TH = P2 * P1 -W = fem.FunctionSpace(mesh, TH) -V, V_to_W = W.sub(0).collapse() + +# The next step is the create the function spaces for the fluid velocit and pressure. +# We will use a mixed-formulation, and we use `basix.ufl` to create the Taylor-Hood finite element pair + +P2 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 2, gdim=mesh.geometry.dim, shape=(mesh.geometry.dim, )) +P1 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 1) +TH = basix.ufl.mixed_element([P2, P1]) +W = fem.functionspace(mesh, TH) +# - + +# ## Boundary conditions +# Next we want to prescribe an inflow condition on a set of facets, those marked by 3 in GMSH +# To do so, we start by creating a function in the **collapsed** subspace of V, and create the +# expression of choice for the boundary condition. +# We use a Python-function for this, that takes in an array of shape `(num_points, 3)` and +# returns the function values as an `(num_points, 2)` array. +# Note that we therefore can use vectorized Numpy operations to compute the inlet values for a sequence of points +# with a single function call + +V, _ = W.sub(0).collapse() Q, _ = W.sub(1).collapse() +inlet_velocity = fem.Function(V) -def inlet_velocity_expression(x: NDArray[np.float64]) -> NDArray[np.bool_]: +def inlet_velocity_expression(x: NDArray[Union[np.float32, np.float64]]) -> NDArray[Union[np.float32, + np.float64, + np.complex64, + np.complex128]]: return np.stack((np.sin(np.pi * np.sqrt(x[0]**2 + x[1]**2)), - 5 * x[1] * np.sin(np.pi * np.sqrt(x[0]**2 + x[1]**2)))) + 5 * x[1] * np.sin(np.pi * np.sqrt(x[0]**2 + x[1]**2)))).astype(default_scalar_type) -# ----------------------Defining boundary conditions---------------------- -# Inlet velocity Dirichlet BC -inlet_velocity = fem.Function(V) inlet_velocity.interpolate(inlet_velocity_expression) inlet_velocity.x.scatter_forward() + +# Next, we have to create the Dirichlet boundary condition. As the function is in the collapsed space, +# we send in the un-collapsed space and the collapsed space to locate dofs topological to find the +# appropriate dofs and the mapping from `V` to `W0`. + +# + W0 = W.sub(0) dofs = fem.locate_dofs_topological((W0, V), 1, mt.find(3)) bc1 = fem.dirichletbc(inlet_velocity, dofs, W0) - -# Collect Dirichlet boundary conditions bcs = [bc1] -# Slip conditions for walls +# - + +# Next, we want to create the slip conditions at the side walls of the channel. +# To do so, we compute an approximation of the normal at all the degrees of freedom that +# are associated with the facets marked with `1`, either by being on one of the vertices of a +# facet marked with `1`, or on the facet itself (for instance at a midpoint). +# If a dof is associated with a vertex that is connected to mutliple facets marked with `1`, +# we define the normal as the average between the normals at the vertex viewed from each facet. + +# + n = dolfinx_mpc.utils.create_normal_approximation(V, mt, 1) -with Timer("~Stokes: Create slip constraint"): +# - + +# Next, we can create the multipoint-constraint, enforcing that $\mathbf{u}\cdot\mathbf{n}=0$, +# except where we have already enforced a Dirichlet boundary condition + +# + +with common.Timer("~Stokes: Create slip constraint"): mpc = MultiPointConstraint(W) mpc.create_slip_constraint(W.sub(0), (mt, 1), n, bcs=bcs) mpc.finalize() +# - + +# ## Variational formulation +# We start by creating some convenience functions, including the tangential projection of a vector, +# the symmetric gradient and traction. + +# + def tangential_proj(u: Expr, n: Expr): @@ -157,50 +204,70 @@ def sym_grad(u: Expr): def T(u: Expr, p: Expr, mu: Expr): return 2 * mu * sym_grad(u) - p * Identity(u.ufl_shape[0]) +# - + + +# Next, we define the classical terms of the bilinear form `a` and linear form `L`. +# We note that we use the symmetric formulation after integration by parts. -# --------------------------Variational problem--------------------------- -# Traditional terms -mu = 1 -f = fem.Constant(mesh, PETSc.ScalarType((0, 0))) +# + +mu = fem.Constant(mesh, default_scalar_type(1.)) +f = fem.Constant(mesh, default_scalar_type((0, 0))) (u, p) = TrialFunctions(W) (v, q) = TestFunctions(W) a = (2 * mu * inner(sym_grad(u), sym_grad(v)) - inner(p, div(v)) - inner(div(u), q)) * dx L = inner(f, v) * dx +# - -# No prescribed shear stress +# We could prescibe some shear stress at the slip boundaries. However, in this demo, we set it to `0`, +# but include it in the variational formulation. We add the appropriate terms due to the slip condition, +# as explained in https://arxiv.org/pdf/2001.10639.pdf + +# + n = FacetNormal(mesh) -g_tau = tangential_proj(fem.Constant(mesh, PETSc.ScalarType(((0, 0), (0, 0)))) * n, n) +g_tau = tangential_proj(fem.Constant(mesh, default_scalar_type(((0, 0), (0, 0)))) * n, n) ds = Measure("ds", domain=mesh, subdomain_data=mt, subdomain_id=1) -# Terms due to slip condition -# Explained in for instance: https://arxiv.org/pdf/2001.10639.pdf a -= inner(outer(n, n) * dot(T(u, p, mu), n), v) * ds L += inner(g_tau, v) * ds +# - + +# ## Solve the variational problem +# We use the MUMPS solver provided through the DOLFINX_MPC PETSc interface to solve the variational problem -# Solve linear problem +# + petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_solver_type": "mumps"} problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options) U = problem.solve() +# - + +# ## Visualization +# We store the solution to the `VTX` file format, which can be opend with the +# `ADIOS2VTXReader` in Paraview -# ------------------------------ Output ---------------------------------- +# + u = U.sub(0).collapse() p = U.sub(1).collapse() u.name = "u" p.name = "p" -outdir = Path("results") +outdir = Path("results").absolute() outdir.mkdir(exist_ok=True, parents=True) -with io.VTXWriter(mesh.comm, outdir / "demo_stokes_u.bp", u) as vtx: +with io.VTXWriter(mesh.comm, outdir / "demo_stokes_u.bp", u, engine="BP4") as vtx: vtx.write(0.0) -with io.VTXWriter(mesh.comm, outdir / "demo_stokes_p.bp", p) as vtx: +with io.VTXWriter(mesh.comm, outdir / "demo_stokes_p.bp", p, engine="BP4") as vtx: vtx.write(0.0) +# - -# -------------------- Verification -------------------------------- -# Transfer data from the MPC problem to numpy arrays for comparison -with Timer("~Stokes: Verification of problem by global matrix reduction"): +# ## Verification +# We verify that the MPC implementation is correct by creating the global reduction matrix +# and performing global matrix-matrix-matrix multiplication using numpy + +# + +with common.Timer("~Stokes: Verification of problem by global matrix reduction"): # Solve the MPC problem using a global transformation matrix # and numpy solvers to get reference values @@ -212,7 +279,7 @@ def T(u: Expr, p: Expr, mu: Expr): L_org = fem.petsc.assemble_vector(linear_form) fem.petsc.apply_lifting(L_org, [bilinear_form], [bcs]) - L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore fem.petsc.set_bc(L_org, bcs) root = 0 dolfinx_mpc.utils.compare_mpc_lhs(A_org, problem.A, mpc, root=root) @@ -223,15 +290,20 @@ def T(u: Expr, p: Expr, mu: Expr): K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(L_org, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(U.vector, root=root) - + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system - d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) + d = scipy.sparse.linalg.spsolve(KTAK.astype(scipy_dtype), reduced_L.astype(scipy_dtype)) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + uh_numpy = K.astype(scipy_dtype) @ d.astype(scipy_dtype) + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, atol=float(5e2 + * np.finfo(u_mpc.dtype).resolution)) +# - + +# ## Timings +# Finally, we list the execution time of various operations used in the demo -# -------------------- List timings -------------------------- -list_timings(MPI.COMM_WORLD, [TimingType.wall]) +common.list_timings(MPI.COMM_WORLD, [common.TimingType.wall]) diff --git a/python/demos/demo_stokes_nest.py b/python/demos/demo_stokes_nest.py index 74dd76a0..5de4c4fd 100644 --- a/python/demos/demo_stokes_nest.py +++ b/python/demos/demo_stokes_nest.py @@ -18,6 +18,7 @@ import numpy as np import scipy.sparse.linalg import ufl +from dolfinx import default_scalar_type from dolfinx.io import gmshio from mpi4py import MPI from petsc4py import PETSc @@ -114,8 +115,8 @@ def create_mesh_gmsh(L: int = 2, H: int = 1, res: float = 0.1, theta: float = np Ve = basix.ufl.element(basix.ElementFamily.P, cellname, 2, shape=(mesh.geometry.dim,)) Qe = basix.ufl.element(basix.ElementFamily.P, cellname, 1) -V = dolfinx.fem.FunctionSpace(mesh, Ve) -Q = dolfinx.fem.FunctionSpace(mesh, Qe) +V = dolfinx.fem.functionspace(mesh, Ve) +Q = dolfinx.fem.functionspace(mesh, Qe) def inlet_velocity_expression(x): @@ -164,7 +165,7 @@ def T(u: Expr, p: Expr, mu: Expr): # --------------------------Variational problem--------------------------- # Traditional terms mu = 1 -f = dolfinx.fem.Constant(mesh, PETSc.ScalarType((0, 0))) +f = dolfinx.fem.Constant(mesh, default_scalar_type((0, 0))) (u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q) (v, q) = ufl.TestFunction(V), ufl.TestFunction(Q) a00 = 2 * mu * ufl.inner(sym_grad(u), sym_grad(v)) * ufl.dx @@ -173,12 +174,12 @@ def T(u: Expr, p: Expr, mu: Expr): a11 = None L0 = ufl.inner(f, v) * ufl.dx -L1 = ufl.inner(dolfinx.fem.Constant(mesh, PETSc.ScalarType(0.0)), q) * ufl.dx +L1 = ufl.inner(dolfinx.fem.Constant(mesh, default_scalar_type(0.0)), q) * ufl.dx # No prescribed shear stress n = ufl.FacetNormal(mesh) g_tau = tangential_proj(dolfinx.fem.Constant( - mesh, PETSc.ScalarType(((0, 0), (0, 0)))) * n, n) + mesh, default_scalar_type(((0, 0), (0, 0)))) * n, n) ds = ufl.Measure("ds", domain=mesh, subdomain_data=mt, subdomain_id=1) # Terms due to slip condition @@ -204,7 +205,7 @@ def T(u: Expr, p: Expr, mu: Expr): # Set Dirichlet boundary condition values in the RHS dolfinx.fem.petsc.apply_lifting_nest(b, a, bcs) for b_sub in b.getNestSubVecs(): - b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) + b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # type: ignore # bcs0 = dolfinx.cpp.fem.bcs_rows( # dolfinx.fem.assemble._create_cpp_form(L), bcs) @@ -213,19 +214,19 @@ def T(u: Expr, p: Expr, mu: Expr): # Preconditioner P11 = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(p * q * ufl.dx)) -P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]]) +P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]]) # type: ignore P.assemble() # ---------------------- Solve variational problem ----------------------- -ksp = PETSc.KSP().create(mesh.comm) +ksp = PETSc.KSP().create(mesh.comm) # type: ignore ksp.setOperators(A, P) ksp.setMonitor( - lambda ctx, it, r: PETSc.Sys.Print( + lambda ctx, it, r: PETSc.Sys.Print( # type: ignore f"Iteration: {it:>4d}, |r| = {r:.3e}")) ksp.setType("minres") ksp.setTolerances(rtol=1e-8) ksp.getPC().setType("fieldsplit") -ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE) +ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE) # type: ignore nested_IS = P.getNestISs() ksp.getPC().setFieldSplitIS( @@ -244,8 +245,8 @@ def T(u: Expr, p: Expr, mu: Expr): ksp.solve(b, Uh) for Uh_sub in Uh.getNestSubVecs(): - Uh_sub.ghostUpdate(addv=PETSc.InsertMode.INSERT, - mode=PETSc.ScatterMode.FORWARD) + Uh_sub.ghostUpdate(addv=PETSc.InsertMode.INSERT, # type: ignore + mode=PETSc.ScatterMode.FORWARD) # type: ignore # ----------------------------- Put NestVec into DOLFINx Function - --------- uh = dolfinx.fem.Function(mpc.function_space) uh.vector.setArray(Uh.getNestSubVecs()[0].array) @@ -274,12 +275,12 @@ def T(u: Expr, p: Expr, mu: Expr): outfile.write_function(uh) outfile.write_function(ph) -with dolfinx.io.VTXWriter(mesh.comm, outdir / "stokes_nest_uh.bp", uh) as vtx: +with dolfinx.io.VTXWriter(mesh.comm, outdir / "stokes_nest_uh.bp", uh, engine="BP4") as vtx: vtx.write(0.0) # -------------------- Verification -------------------------------- # Transfer data from the MPC problem to numpy arrays for comparison with dolfinx.common.Timer("~Stokes: Verification of problem by global matrix reduction"): - W = dolfinx.fem.FunctionSpace(mesh, ufl.MixedElement([Ve, Qe])) + W = dolfinx.fem.functionspace(mesh, basix.ufl.mixed_element([Ve, Qe])) V, V_to_W = W.sub(0).collapse() _, Q_to_W = W.sub(1).collapse() @@ -324,7 +325,7 @@ def T(u: Expr, p: Expr, mu: Expr): L_org = dolfinx.fem.petsc.assemble_vector(Lf) dolfinx.fem.petsc.apply_lifting(L_org, [af], [bcs]) - L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore dolfinx.fem.petsc.set_bc(L_org, bcs) root = 0 diff --git a/python/dolfinx_mpc/assemble_matrix.py b/python/dolfinx_mpc/assemble_matrix.py index 4a89c82c..7f1fefca 100644 --- a/python/dolfinx_mpc/assemble_matrix.py +++ b/python/dolfinx_mpc/assemble_matrix.py @@ -19,8 +19,8 @@ def assemble_matrix(form: _fem.Form, constraint: Union[MultiPointConstraint, Sequence[MultiPointConstraint]], bcs: Optional[Sequence[_fem.DirichletBC]] = None, - diagval: _PETSc.ScalarType = 1, - A: Optional[_PETSc.Mat] = None) -> _PETSc.Mat: + diagval: _PETSc.ScalarType = 1, # type: ignore + A: Optional[_PETSc.Mat] = None) -> _PETSc.Mat: # type: ignore """ Assemble a compiled DOLFINx bilinear form into a PETSc matrix with corresponding multi point constraints and Dirichlet boundary conditions. @@ -33,7 +33,7 @@ def assemble_matrix(form: _fem.Form, A: PETSc matrix to assemble into Returns: - _PETSc.Mat: The matrix with the assembled bi-linear form + _PETSc.Mat: The matrix with the assembled bi-linear form #type: ignore """ bcs = [] if bcs is None else [bc._cpp_object for bc in bcs] if not isinstance(constraint, Sequence): @@ -52,8 +52,8 @@ def assemble_matrix(form: _fem.Form, # Add one on diagonal for Dirichlet boundary conditions if form.function_spaces[0] is form.function_spaces[1]: - A.assemblyBegin(_PETSc.Mat.AssemblyType.FLUSH) - A.assemblyEnd(_PETSc.Mat.AssemblyType.FLUSH) + A.assemblyBegin(_PETSc.Mat.AssemblyType.FLUSH) # type: ignore + A.assemblyEnd(_PETSc.Mat.AssemblyType.FLUSH) # type: ignore _cpp.fem.petsc.insert_diagonal(A, form.function_spaces[0], bcs, diagval) A.assemble() @@ -107,17 +107,17 @@ def create_matrix_nest( A_[i][j] = cpp.mpc.create_matrix( a[i][j]._cpp_object, constraints[i]._cpp_object, constraints[j]._cpp_object) - A = _PETSc.Mat().createNest( + A = _PETSc.Mat().createNest( # type: ignore A_, comm=constraints[0].function_space.mesh.comm) return A def assemble_matrix_nest( - A: _PETSc.Mat, + A: _PETSc.Mat, # type: ignore a: Sequence[Sequence[_fem.Form]], constraints: Sequence[MultiPointConstraint], bcs: Sequence[_fem.DirichletBC] = [], - diagval: _PETSc.ScalarType = 1): + diagval: _PETSc.ScalarType = 1): # type: ignore """ Assemble a compiled DOLFINx bilinear form into a PETSc matrix of type "nest" with corresponding multi point constraints and Dirichlet boundary diff --git a/python/dolfinx_mpc/assemble_vector.py b/python/dolfinx_mpc/assemble_vector.py index bbf7b63f..9f10afd9 100644 --- a/python/dolfinx_mpc/assemble_vector.py +++ b/python/dolfinx_mpc/assemble_vector.py @@ -19,8 +19,8 @@ from .multipointconstraint import MultiPointConstraint -def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.DirichletBC]], - constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], scale: float = 1.0): +def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.DirichletBC]], # type: ignore + constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], scale: float = 1.0): # type: ignore """ Apply lifting to vector b, i.e. :math:`b = b - scale \\cdot K^T (A_j (g_j - x0_j))` @@ -46,7 +46,7 @@ def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.Diri def assemble_vector(form: ufl.form.Form, constraint: MultiPointConstraint, - b: Optional[_PETSc.Vec] = None) -> _PETSc.Vec: + b: Optional[_PETSc.Vec] = None) -> _PETSc.Vec: # type: ignore """ Assemble a linear form into vector `b` with corresponding multi point constraint @@ -73,7 +73,7 @@ def assemble_vector(form: ufl.form.Form, constraint: MultiPointConstraint, def create_vector_nest( L: Sequence[_fem.Form], - constraints: Sequence[MultiPointConstraint]) -> _PETSc.Vec: + constraints: Sequence[MultiPointConstraint]) -> _PETSc.Vec: # type: ignore """ Create a PETSc vector of type "nest" appropriate for the provided multi point constraints @@ -83,7 +83,7 @@ def create_vector_nest( constraints: An ordered list of multi point constraints Returns: - PETSc.Vec: A PETSc vector of type "nest" + PETSc.Vec: A PETSc vector of type "nest" #type: ignore """ assert len(constraints) == len(L) @@ -94,7 +94,7 @@ def create_vector_nest( def assemble_vector_nest( - b: _PETSc.Vec, + b: _PETSc.Vec, # type: ignore L: Sequence[_fem.Form], constraints: Sequence[MultiPointConstraint]): """ diff --git a/python/dolfinx_mpc/dictcondition.py b/python/dolfinx_mpc/dictcondition.py index 7e81ab0b..bb4f7562 100644 --- a/python/dolfinx_mpc/dictcondition.py +++ b/python/dolfinx_mpc/dictcondition.py @@ -6,13 +6,14 @@ import typing +import dolfinx import dolfinx.fem as fem import numpy as np -import numpy.typing -from petsc4py import PETSc +from dolfinx import default_scalar_type -def close_to(point: np.typing.NDArray[np.float64]): +def close_to(point: np.typing.NDArray[typing.Union[np.float64, np.float32]], + atol=1000 * np.finfo(dolfinx.default_real_type).resolution): """ Convenience function for locating a point [x,y,z] within an array x [[x0,...,xN],[y0,...,yN], [z0,...,zN]]. @@ -20,11 +21,11 @@ def close_to(point: np.typing.NDArray[np.float64]): Args: point: The point should be padded to 3D """ - return lambda x: np.isclose(x, point).all(axis=0) + return lambda x: np.isclose(x, point, atol=atol).all(axis=0) @typing.no_type_check -def create_dictionary_constraint(V: fem.FunctionSpace, slave_master_dict: +def create_dictionary_constraint(V: fem.functionspace, slave_master_dict: typing.Dict[bytes, typing.Dict[bytes, float]], subspace_slave: typing.Optional[int] = None, subspace_master: typing.Optional[int] = None): @@ -46,11 +47,10 @@ def create_dictionary_constraint(V: fem.FunctionSpace, slave_master_dict: .. highlight:: python .. code-block:: python - {np.array([d0, d1], dtype=numpy.float64).tobytes(): - {numpy.array([e0, e1], dtype=numpy.float64).tobytes(): alpha, - numpy.array([f0, f1], dtype=numpy.float64).tobytes(): beta}} + {np.array([d0, d1], dtype=mesh.geometry.x.dtype).tobytes(): + {numpy.array([e0, e1], dtype=mesh.geometry.x.dtype).tobytes(): alpha, + numpy.array([f0, f1], dtype=mesh.geometry.x.dtype).tobytes(): beta}} """ - dfloat = np.float64 comm = V.mesh.comm bs = V.dofmap.index_map_bs local_size = V.dofmap.index_map.size_local * bs @@ -60,6 +60,7 @@ def create_dictionary_constraint(V: fem.FunctionSpace, slave_master_dict: non_local_entities = {} slaves_local = {} slaves_ghost = {} + dfloat = V.mesh.geometry.x.dtype slave_point_nd = np.zeros((3, 1), dtype=dfloat) for i, slave_point in enumerate(slave_master_dict.keys()): num_masters = len(list(slave_master_dict[slave_point].keys())) @@ -82,14 +83,14 @@ def create_dictionary_constraint(V: fem.FunctionSpace, slave_master_dict: if slave_dofs[0] < local_size: slaves_local[i] = slave_dofs[0] owned_entities[i] = {"masters": np.full(num_masters, -1, dtype=np.int64), - "coeffs": np.full(num_masters, -1, dtype=np.float64), + "coeffs": np.full(num_masters, -1, dtype=dolfinx.default_scalar_type), "owners": np.full(num_masters, -1, dtype=np.int32), "master_count": 0, "local_index": []} slave_status = 1 else: slaves_ghost[i] = slave_dofs[0] ghosted_entities[i] = {"masters": np.full(num_masters, -1, dtype=np.int64), - "coeffs": np.full(num_masters, -1, dtype=np.float64), + "coeffs": np.full(num_masters, -1, dtype=dolfinx.default_scalar_type), "owners": np.full(num_masters, -1, dtype=np.int32), "master_count": 0, "local_index": []} slave_status = 0 @@ -210,5 +211,5 @@ def create_dictionary_constraint(V: fem.FunctionSpace, slave_master_dict: coeffs.extend(ghosted_slaves[slave_index]["coeffs"]) # type: ignore offsets.append(len(masters)) return (np.asarray(slaves, dtype=np.int32), np.asarray(masters, dtype=np.int64), - np.asarray(coeffs, dtype=PETSc.ScalarType), np.asarray(owners, dtype=np.int32), + np.asarray(coeffs, dtype=default_scalar_type), np.asarray(owners, dtype=np.int32), np.asarray(offsets, dtype=np.int32)) diff --git a/python/dolfinx_mpc/multipointconstraint.py b/python/dolfinx_mpc/multipointconstraint.py index 58145f4c..7ce797fc 100644 --- a/python/dolfinx_mpc/multipointconstraint.py +++ b/python/dolfinx_mpc/multipointconstraint.py @@ -230,7 +230,7 @@ def create_slip_constraint(self, space: _fem.FunctionSpaceBase, facet_marker: Tu .. highlight:: python .. code-block:: python - V = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = dolfinx.fem.functionspace(mesh, ("CG", 1)) mpc = MultiPointConstraint(V) n = dolfinx.fem.Function(V) mpc.create_slip_constaint(V, (mt, i), n) @@ -299,9 +299,9 @@ def create_general_constraint(self, slave_master_dict: Dict[bytes, Dict[bytes, f .. highlight:: python .. code-block:: python - {numpy.array([d0, d1], dtype=numpy.float64).tobytes(): - {numpy.array([e0, e1], dtype=numpy.float64).tobytes(): alpha, - numpy.array([f0, f1], dtype=numpy.float64).tobytes(): beta}} + {numpy.array([d0, d1], dtype=mesh.geometry.x.dtype).tobytes(): + {numpy.array([e0, e1], dtype=mesh.geometry.x.dtype).tobytes(): alpha, + numpy.array([f0, f1], dtype=mesh.geometry.x.dtype).tobytes(): beta}} """ slaves, masters, coeffs, owners, offsets = create_dictionary_constraint( self.V, slave_master_dict, subspace_slave, subspace_master) @@ -427,7 +427,7 @@ def function_space(self): self._not_finalized() return self.V - def backsubstitution(self, u: Union[_fem.Function, _PETSc.Vec]) -> None: + def backsubstitution(self, u: Union[_fem.Function, _PETSc.Vec]) -> None: # type: ignore """ For a Function, impose the multi-point constraint by backsubstiution. This function is used after solving the reduced problem to obtain the values @@ -445,7 +445,7 @@ def backsubstitution(self, u: Union[_fem.Function, _PETSc.Vec]) -> None: except AttributeError: with u.localForm() as vector_local: self._cpp_object.backsubstitution(vector_local.array_w) - u.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD) + u.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD) # type: ignore def homogenize(self, u: _fem.Function) -> None: """ diff --git a/python/dolfinx_mpc/numba/assemble_matrix.py b/python/dolfinx_mpc/numba/assemble_matrix.py index 53a0db3e..9ed04b22 100644 --- a/python/dolfinx_mpc/numba/assemble_matrix.py +++ b/python/dolfinx_mpc/numba/assemble_matrix.py @@ -12,6 +12,7 @@ import numba import numpy import numpy.typing as npt +import dolfinx from dolfinx.common import Timer from dolfinx_mpc.assemble_matrix import create_sparsity_pattern from dolfinx_mpc.multipointconstraint import MultiPointConstraint @@ -20,14 +21,14 @@ from .helpers import _bcs, _forms, extract_slave_cells, pack_slave_facet_info from .numba_setup import initialize_petsc, sink -mode = _PETSc.InsertMode.ADD_VALUES -insert = _PETSc.InsertMode.INSERT_VALUES +mode = _PETSc.InsertMode.ADD_VALUES # type: ignore +insert = _PETSc.InsertMode.INSERT_VALUES # type: ignore ffi, set_values_local = initialize_petsc() def assemble_matrix(form: _forms, constraint: MultiPointConstraint, - bcs: Optional[List[_bcs]] = None, diagval: _PETSc.ScalarType = 1., - A: Optional[_PETSc.Mat] = None): + bcs: Optional[List[_bcs]] = None, diagval: _PETSc.ScalarType = 1., # type: ignore + A: Optional[_PETSc.Mat] = None): # type: ignore """ Assembles a compiled DOLFINx form with given a multi point constraint and possible Dirichlet boundary conditions. @@ -106,8 +107,17 @@ def assemble_matrix(form: _forms, constraint: MultiPointConstraint, if e0.needs_dof_transformations or e1.needs_dof_transformations: raise NotImplementedError("Dof transformations not implemented") - is_complex = numpy.issubdtype(_PETSc.ScalarType, numpy.complexfloating) - nptype = "complex128" if is_complex else "float64" + if _PETSc.ScalarType == numpy.float32: # type: ignore + nptype = "float32" + elif _PETSc.ScalarType == numpy.float64: # type: ignore + nptype = "float64" + elif _PETSc.ScalarType == numpy.complex64: # type: ignore + nptype = "complex64" + elif _PETSc.ScalarType == numpy.complex128: # type: ignore + nptype = "complex128" + else: + raise RuntimeError(f"Unsupported scalar type {_PETSc.ScalarType}.") # type: ignore + ufcx_form = form.ufcx_form if num_cell_integrals > 0: # NOTE: This depends on enum ordering in ufcx.h @@ -157,8 +167,8 @@ def assemble_matrix(form: _forms, constraint: MultiPointConstraint, # Add one on diagonal for diriclet bc and slave dofs # NOTE: In the future one could use a constant in the dirichletbc if form.function_spaces[0] is form.function_spaces[1]: - A.assemblyBegin(_PETSc.Mat.AssemblyType.FLUSH) - A.assemblyEnd(_PETSc.Mat.AssemblyType.FLUSH) + A.assemblyBegin(_PETSc.Mat.AssemblyType.FLUSH) # type: ignore + A.assemblyEnd(_PETSc.Mat.AssemblyType.FLUSH) # type: ignore _cpp.fem.petsc.insert_diagonal(A, form.function_spaces[0], bcs, diagval) A.assemble() @@ -167,13 +177,13 @@ def assemble_matrix(form: _forms, constraint: MultiPointConstraint, @numba.njit -def add_diagonal(A: int, dofs: npt.NDArray[numpy.int32], diagval: _PETSc.ScalarType = 1): +def add_diagonal(A: int, dofs: npt.NDArray[numpy.int32], diagval: _PETSc.ScalarType = 1): # type: ignore """ Insert value on diagonal of matrix for given dofs. """ ffi_fb = ffi.from_buffer dof_list = numpy.zeros(1, dtype=numpy.int32) - dof_value = numpy.full(1, diagval, dtype=_PETSc.ScalarType) + dof_value = numpy.full(1, diagval, dtype=_PETSc.ScalarType) # type: ignore for dof in dofs: dof_list[0] = dof ierr_loc = set_values_local(A, 1, ffi_fb(dof_list), 1, ffi_fb(dof_list), ffi_fb(dof_value), mode) @@ -186,14 +196,14 @@ def assemble_slave_cells(A: int, kernel: cffi.FFI, active_cells: npt.NDArray[numpy.int32], mesh: Tuple[npt.NDArray[numpy.int32], - npt.NDArray[numpy.float64]], - coeffs: npt.NDArray[_PETSc.ScalarType], - constants: npt.NDArray[_PETSc.ScalarType], + npt.NDArray[dolfinx.default_real_type]], + coeffs: npt.NDArray[_PETSc.ScalarType], # type: ignore + constants: npt.NDArray[_PETSc.ScalarType], # type: ignore permutation_info: npt.NDArray[numpy.uint32], dofmap: npt.NDArray[numpy.int32], block_size: int, num_dofs_per_element: int, - mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], + mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], # type: ignore npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32]], is_bc: npt.NDArray[numpy.bool_]): @@ -211,9 +221,9 @@ def assemble_slave_cells(A: int, # NOTE: All cells are assumed to be of the same type num_xdofs_per_cell = x_dofmap.shape[1] - geometry = numpy.zeros((num_xdofs_per_cell, 3)) + geometry = numpy.zeros((num_xdofs_per_cell, 3), dtype=dolfinx.default_real_type) A_local = numpy.zeros((block_size * num_dofs_per_element, block_size - * num_dofs_per_element), dtype=_PETSc.ScalarType) + * num_dofs_per_element), dtype=_PETSc.ScalarType) # type: ignore masters, coefficients, offsets, c_to_s, c_to_s_off, is_slave = mpc # Loop over all cells @@ -238,7 +248,7 @@ def assemble_slave_cells(A: int, A_local[j * block_size + k, :] = 0 A_local[:, j * block_size + k] = 0 - A_local_copy: numpy.typing.NDArray[_PETSc.ScalarType] = A_local.copy() + A_local_copy: numpy.typing.NDArray[_PETSc.ScalarType] = A_local.copy() # type: ignore # Find local position of slaves slaves = c_to_s[c_to_s_off[cell]: c_to_s_off[cell + 1]] @@ -263,9 +273,9 @@ def assemble_slave_cells(A: int, @numba.njit def modify_mpc_cell(A: int, num_dofs: int, block_size: int, - Ae: npt.NDArray[_PETSc.ScalarType], + Ae: npt.NDArray[_PETSc.ScalarType], # type: ignore local_blocks: npt.NDArray[numpy.int32], - mpc_cell: Tuple[npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], + mpc_cell: Tuple[npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], # type: ignore npt.NDArray[_PETSc.ScalarType], npt.NDArray[numpy.int32], npt.NDArray[numpy.int8]]): """ @@ -287,7 +297,7 @@ def modify_mpc_cell(A: int, num_dofs: int, block_size: int, num_flattened_masters += offsets[slave + 1] - offsets[slave] # Strip a copy of Ae of all columns and rows belonging to a slave Ae_original = numpy.copy(Ae) - Ae_stripped = numpy.zeros((block_size * num_dofs, block_size * num_dofs), dtype=_PETSc.ScalarType) + Ae_stripped = numpy.zeros((block_size * num_dofs, block_size * num_dofs), dtype=_PETSc.ScalarType) # type: ignore for i in range(num_dofs): for b in range(block_size): is_slave0 = is_slave[local_blocks[i] * block_size + b] @@ -299,7 +309,7 @@ def modify_mpc_cell(A: int, num_dofs: int, block_size: int, j * block_size + c] flattened_masters = numpy.zeros(num_flattened_masters, dtype=numpy.int32) flattened_slaves = numpy.zeros(num_flattened_masters, dtype=numpy.int32) - flattened_coeffs = numpy.zeros(num_flattened_masters, dtype=_PETSc.ScalarType) + flattened_coeffs = numpy.zeros(num_flattened_masters, dtype=_PETSc.ScalarType) # type: ignore c = 0 for i, slave in enumerate(slaves): local_masters = masters[offsets[slave]: offsets[slave + 1]] @@ -312,9 +322,9 @@ def modify_mpc_cell(A: int, num_dofs: int, block_size: int, c += num_masters m0 = numpy.zeros(1, dtype=numpy.int32) m1 = numpy.zeros(1, dtype=numpy.int32) - Am0m1 = numpy.zeros((1, 1), dtype=_PETSc.ScalarType) - Arow = numpy.zeros(block_size * num_dofs, dtype=_PETSc.ScalarType) - Acol = numpy.zeros(block_size * num_dofs, dtype=_PETSc.ScalarType) + Am0m1 = numpy.zeros((1, 1), dtype=_PETSc.ScalarType) # type: ignore + Arow = numpy.zeros(block_size * num_dofs, dtype=_PETSc.ScalarType) # type: ignore + Acol = numpy.zeros(block_size * num_dofs, dtype=_PETSc.ScalarType) # type: ignore mpc_dofs = numpy.zeros(block_size * num_dofs, dtype=numpy.int32) ffi_fb = ffi.from_buffer for i in range(num_flattened_masters): @@ -360,14 +370,14 @@ def modify_mpc_cell(A: int, num_dofs: int, block_size: int, def assemble_exterior_slave_facets(A: int, kernel: cffi.FFI, mesh: Tuple[npt.NDArray[numpy.int32], npt.NDArray[numpy.float64]], - coeffs: npt.NDArray[_PETSc.ScalarType], - consts: npt.NDArray[_PETSc.ScalarType], + coeffs: npt.NDArray[_PETSc.ScalarType], # type: ignore + consts: npt.NDArray[_PETSc.ScalarType], # type: ignore perm: npt.NDArray[numpy.uint32], dofmap: npt.NDArray[numpy.int32], block_size: int, num_dofs_per_element: int, facet_info: npt.NDArray[numpy.int32], - mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], + mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], # type: ignore npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32]], is_bc: npt.NDArray[numpy.bool_], @@ -384,11 +394,11 @@ def assemble_exterior_slave_facets(A: int, kernel: cffi.FFI, facet_perm = numpy.zeros(1, dtype=numpy.uint8) # NOTE: All cells are assumed to be of the same type - geometry = numpy.zeros((x_dofmap.shape[1], 3)) + geometry = numpy.zeros((x_dofmap.shape[1], 3), dtype=x.dtype) # Numpy data used in facet loop A_local = numpy.zeros((num_dofs_per_element * block_size, - num_dofs_per_element * block_size), dtype=_PETSc.ScalarType) + num_dofs_per_element * block_size), dtype=_PETSc.ScalarType) # type: ignore local_dofs = numpy.zeros(block_size * num_dofs_per_element, dtype=numpy.int32) # Permutation info @@ -420,7 +430,7 @@ def assemble_exterior_slave_facets(A: int, kernel: cffi.FFI, A_local[j * block_size + k, :] = 0 A_local[:, j * block_size + k] = 0 - A_local_copy: numpy.typing.NDArray[_PETSc.ScalarType] = A_local.copy() + A_local_copy: numpy.typing.NDArray[_PETSc.ScalarType] = A_local.copy() # type: ignore slaves = c_to_s[c_to_s_off[cell_index]: c_to_s_off[cell_index + 1]] mpc_cell = (slaves, masters, coefficients, offsets, is_slave) modify_mpc_cell(A, num_dofs_per_element, block_size, A_local, local_blocks, mpc_cell) diff --git a/python/dolfinx_mpc/numba/assemble_vector.py b/python/dolfinx_mpc/numba/assemble_vector.py index 70d2a08a..d73dcee6 100644 --- a/python/dolfinx_mpc/numba/assemble_vector.py +++ b/python/dolfinx_mpc/numba/assemble_vector.py @@ -7,6 +7,7 @@ from typing import Optional, Tuple import cffi +import dolfinx import dolfinx.cpp as _cpp import dolfinx.fem as _fem import dolfinx.la as _la @@ -25,7 +26,8 @@ ffi, _ = initialize_petsc() -def assemble_vector(form: _forms, constraint: MultiPointConstraint, b: Optional[_PETSc.Vec] = None) -> _PETSc.Vec: +def assemble_vector(form: _forms, + constraint: MultiPointConstraint, b: Optional[_PETSc.Vec] = None) -> _PETSc.Vec: # type: ignore """ Assemble a compiled DOLFINx form into vector b. @@ -88,8 +90,16 @@ def assemble_vector(form: _forms, constraint: MultiPointConstraint, b: Optional[ subdomain_ids = form._cpp_object.integral_ids(_fem.IntegralType.cell) num_cell_integrals = len(subdomain_ids) - is_complex = numpy.issubdtype(_PETSc.ScalarType, numpy.complexfloating) - nptype = "complex128" if is_complex else "float64" + if _PETSc.ScalarType == numpy.float32: # type: ignore + nptype = "float32" + elif _PETSc.ScalarType == numpy.float64: # type: ignore + nptype = "float64" + elif _PETSc.ScalarType == numpy.complex64: # type: ignore + nptype = "complex64" + elif _PETSc.ScalarType == numpy.complex128: # type: ignore + nptype = "complex128" + else: + raise RuntimeError(f"Unsupported scalar type {_PETSc.ScalarType}.") # type: ignore ufcx_form = form.ufcx_form if num_cell_integrals > 0: V.mesh.topology.create_entity_permutations() @@ -134,17 +144,17 @@ def assemble_vector(form: _forms, constraint: MultiPointConstraint, b: Optional[ @numba.njit -def assemble_cells(b: npt.NDArray[_PETSc.ScalarType], +def assemble_cells(b: npt.NDArray[_PETSc.ScalarType], # type: ignore kernel: cffi.FFI, active_cells: npt.NDArray[numpy.int32], mesh: Tuple[npt.NDArray[numpy.int32], - npt.NDArray[numpy.float64]], - coeffs: npt.NDArray[_PETSc.ScalarType], - constants: npt.NDArray[_PETSc.ScalarType], + npt.NDArray[dolfinx.default_real_type]], + coeffs: npt.NDArray[_PETSc.ScalarType], # type: ignore + constants: npt.NDArray[_PETSc.ScalarType], # type: ignore permutation_info: npt.NDArray[numpy.uint32], dofmap: npt.NDArray[numpy.int32], block_size: int, num_dofs_per_element: int, - mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], + mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], # type: ignore npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32]]): """Assemble additional MPC contributions for cell integrals""" @@ -158,8 +168,8 @@ def assemble_cells(b: npt.NDArray[_PETSc.ScalarType], x_dofmap, x = mesh # NOTE: All cells are assumed to be of the same type - geometry = numpy.zeros((x_dofmap.shape[1], 3)) - b_local = numpy.zeros(block_size * num_dofs_per_element, dtype=_PETSc.ScalarType) + geometry = numpy.zeros((x_dofmap.shape[1], 3), dtype=dolfinx.default_real_type) + b_local = numpy.zeros(block_size * num_dofs_per_element, dtype=_PETSc.ScalarType) # type: ignore for cell_index in active_cells: @@ -184,18 +194,18 @@ def assemble_cells(b: npt.NDArray[_PETSc.ScalarType], @numba.njit -def assemble_exterior_slave_facets(b: npt.NDArray[_PETSc.ScalarType], +def assemble_exterior_slave_facets(b: npt.NDArray[_PETSc.ScalarType], # type: ignore kernel: cffi.FFI, facet_info: npt.NDArray[numpy.int32], mesh: Tuple[npt.NDArray[numpy.int32], npt.NDArray[numpy.float64]], - coeffs: npt.NDArray[_PETSc.ScalarType], - constants: npt.NDArray[_PETSc.ScalarType], + coeffs: npt.NDArray[_PETSc.ScalarType], # type: ignore + constants: npt.NDArray[_PETSc.ScalarType], # type: ignore permutation_info: npt.NDArray[numpy.uint32], dofmap: npt.NDArray[numpy.int32], block_size: int, num_dofs_per_element: int, - mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], + mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], # type: ignore npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32]], num_facets_per_cell: int): @@ -210,8 +220,8 @@ def assemble_exterior_slave_facets(b: npt.NDArray[_PETSc.ScalarType], # Unpack mesh data x_dofmap, x = mesh - geometry = numpy.zeros((x_dofmap.shape[1], 3)) - b_local = numpy.zeros(block_size * num_dofs_per_element, dtype=_PETSc.ScalarType) + geometry = numpy.zeros((x_dofmap.shape[1], 3), dtype=x.dtype) + b_local = numpy.zeros(block_size * num_dofs_per_element, dtype=_PETSc.ScalarType) # type: ignore for i in range(facet_info.shape[0]): # Extract cell index (local to process) and facet index (local to cell) for kernel cell_index, local_facet = facet_info[i] @@ -240,10 +250,10 @@ def assemble_exterior_slave_facets(b: npt.NDArray[_PETSc.ScalarType], @numba.njit(cache=True) -def modify_mpc_contributions(b: npt.NDArray[_PETSc.ScalarType], cell_index: int, - b_local: npt.NDArray[_PETSc.ScalarType], - b_copy: npt.NDArray[_PETSc.ScalarType], - mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], +def modify_mpc_contributions(b: npt.NDArray[_PETSc.ScalarType], cell_index: int, # type: ignore + b_local: npt.NDArray[_PETSc.ScalarType], # type: ignore + b_copy: npt.NDArray[_PETSc.ScalarType], # type: ignore + mpc: Tuple[npt.NDArray[numpy.int32], npt.NDArray[_PETSc.ScalarType], # type: ignore npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32], npt.NDArray[numpy.int32]], dofmap: npt.NDArray[numpy.int32], diff --git a/python/dolfinx_mpc/numba/numba_setup.py b/python/dolfinx_mpc/numba/numba_setup.py index 94887b1e..173ef709 100644 --- a/python/dolfinx_mpc/numba/numba_setup.py +++ b/python/dolfinx_mpc/numba/numba_setup.py @@ -29,10 +29,10 @@ def initialize_petsc() -> typing.Tuple[cffi.FFI, typing.Any]: petsc_arch = petsc4py.lib.getPathArchPETSc()[1] # Get PETSc int and scalar types - cmplx = True if np.dtype(PETSc.ScalarType).kind == 'c' else False + cmplx = True if np.dtype(PETSc.ScalarType).kind == 'c' else False # type: ignore - scalar_size = np.dtype(PETSc.ScalarType).itemsize - index_size = np.dtype(PETSc.IntType).itemsize + scalar_size = np.dtype(PETSc.ScalarType).itemsize # type: ignore + index_size = np.dtype(PETSc.IntType).itemsize # type: ignore if index_size == 8: c_int_t = "int64_t" ctypes_index = ctypes.c_int64 # type: ignore @@ -60,7 +60,6 @@ def initialize_petsc() -> typing.Tuple[cffi.FFI, typing.Any]: raise RuntimeError( "Cannot translate PETSc scalar type to a C type, complex: {} size: {}." .format(complex, scalar_size)) - # Load PETSc library via ctypes petsc_lib_name = ctypes.util.find_library("petsc") if petsc_lib_name is not None: diff --git a/python/dolfinx_mpc/problem.py b/python/dolfinx_mpc/problem.py index 27f5a54b..ef5306e0 100644 --- a/python/dolfinx_mpc/problem.py +++ b/python/dolfinx_mpc/problem.py @@ -36,7 +36,7 @@ class LinearProblem(dolfinx.fem.petsc.LinearProblem): .. code-block:: python u = dolfinx.fem.Function(mpc.function_space) - petsc_options: Parameters that is passed to the linear algebra backend PETSc. + petsc_options: Parameters that is passed to the linear algebra backend PETSc. #type: ignore For available choices for the 'petsc_options' kwarg, see the PETSc-documentation https://www.mcs.anl.gov/petsc/documentation/index.html. form_compiler_options: Parameters used in FFCx compilation of this form. Run `ffcx --help` at @@ -59,10 +59,10 @@ class LinearProblem(dolfinx.fem.petsc.LinearProblem): _a: _fem.Form _L: _fem.Form _mpc: MultiPointConstraint - _A: PETSc.Mat - _b: PETSc.Vec - _solver: PETSc.KSP - _x: PETSc.Vec + _A: PETSc.Mat # type: ignore + _b: PETSc.Vec # type: ignore + _solver: PETSc.KSP # type: ignore + _x: PETSc.Vec # type: ignore bcs: typing.List[_fem.DirichletBC] __slots__ = tuple(__annotations__) @@ -101,7 +101,7 @@ def __init__(self, a: ufl.Form, L: ufl.Form, mpc: MultiPointConstraint, self._mpc.function_space.dofmap.index_map_bs) self.bcs = [] if bcs is None else bcs - self._solver = PETSc.KSP().create(self.u.function_space.mesh.comm) + self._solver = PETSc.KSP().create(self.u.function_space.mesh.comm) # type: ignore self._solver.setOperators(self._A) # Give PETSc solver options a unique prefix @@ -109,7 +109,7 @@ def __init__(self, a: ufl.Form, L: ufl.Form, mpc: MultiPointConstraint, self._solver.setOptionsPrefix(solver_prefix) # Set PETSc options - opts = PETSc.Options() + opts = PETSc.Options() # type: ignore opts.prefixPush(solver_prefix) if petsc_options is not None: for k, v in petsc_options.items(): @@ -136,7 +136,7 @@ def solve(self) -> _fem.Function: # Apply boundary conditions to the rhs apply_lifting(self._b, [self._a], [self.bcs], self._mpc) - self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) + self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # type: ignore _fem.petsc.set_bc(self._b, self.bcs) # Solve linear system and update ghost values in the solution diff --git a/python/dolfinx_mpc/utils/mpc_utils.py b/python/dolfinx_mpc/utils/mpc_utils.py index 3d727543..c7b9a652 100644 --- a/python/dolfinx_mpc/utils/mpc_utils.py +++ b/python/dolfinx_mpc/utils/mpc_utils.py @@ -13,12 +13,14 @@ import dolfinx.la as _la import dolfinx.log as _log import dolfinx.mesh as _mesh -import dolfinx_mpc.cpp import numpy as np import ufl +from dolfinx import default_scalar_type as _dt from mpi4py import MPI from petsc4py import PETSc +import dolfinx_mpc.cpp + __all__ = ["rotation_matrix", "facet_normal_approximation", "log_info", "rigid_motions_nullspace", "determine_closest_block", "create_normal_approximation", "create_point_to_point_constraint"] @@ -111,8 +113,8 @@ def tangential_proj(u, n): _cpp.fem.petsc.assemble_matrix(A, bilinear_form._cpp_object, form_consts, form_coeffs, [bc_deac._cpp_object]) if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]: - A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH) - A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH) + A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH) # type: ignore + A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH) # type: ignore _cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [bc_deac._cpp_object], 1.0) A.assemble() linear_form = _fem.form(L, jit_options=jit_options, @@ -120,16 +122,16 @@ def tangential_proj(u, n): b = _fem.petsc.assemble_vector(linear_form) _fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]]) - b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) # type: ignore _fem.petsc.set_bc(b, [bc_deac]) # Solve Linear problem - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(V.mesh.comm) # type: ignore solver.setType("cg") solver.rtol = 1e-8 solver.setOperators(A) solver.solve(b, nh.vector) - nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) + nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) # type: ignore timer.stop() return nh @@ -191,8 +193,8 @@ def rigid_motions_nullspace(V: _fem.FunctionSpaceBase): basis[5][dofs[1]] = -x2 _la.orthonormalize(nullspace_basis) - assert _la.is_orthonormal(nullspace_basis) - return PETSc.NullSpace().create(vectors=nullspace_basis) + assert _la.is_orthonormal(nullspace_basis, float(500 * np.finfo(_x.x.array.dtype).resolution)) + return PETSc.NullSpace().create(vectors=nullspace_basis) # type: ignore def determine_closest_block(V, point): @@ -214,9 +216,9 @@ def determine_closest_block(V, point): midpoint_tree = _geometry.create_midpoint_tree(V.mesh, tdim, boundary_cells) # Find facet closest - + point = np.reshape(point, (1, 3)).astype(V.mesh.geometry.x.dtype) closest_cell = _geometry.compute_closest_entity( - bb_tree, midpoint_tree, V.mesh, np.reshape(point, (1, 3)))[0] + bb_tree, midpoint_tree, V.mesh, point)[0] # Set distance high if cell is not owned if cell_imap.size_local < closest_cell or closest_cell == -1: @@ -330,7 +332,7 @@ def create_point_to_point_constraint(V, slave_point, master_point, vector=None): if vector is None: masters = masters_as_glob owners = np.full(len(masters), master_proc, dtype=np.int32) - coeffs = np.ones(len(masters), dtype=PETSc.ScalarType) + coeffs = np.ones(len(masters), dtype=_dt) offsets = np.arange(0, len(masters) + 1, dtype=np.int32) else: for i in range(len(masters_as_glob)): @@ -391,7 +393,7 @@ def create_point_to_point_constraint(V, slave_point, master_point, vector=None): ghost_slaves[i] = local_size * block_size + idx slaves = np.asarray(np.append(slaves, ghost_slaves), dtype=np.int32) masters = np.asarray(np.append(masters, ghost_masters), dtype=np.int64) - coeffs = np.asarray(np.append(coeffs, ghost_coeffs), dtype=np.int32) + coeffs = np.asarray(np.append(coeffs, ghost_coeffs), dtype=_dt) owners = np.asarray(np.append(owners, ghost_owners), dtype=np.int32) offsets = np.asarray(np.append(offsets, ghost_offsets), dtype=np.int32) return slaves, masters, coeffs, owners, offsets diff --git a/python/dolfinx_mpc/utils/test.py b/python/dolfinx_mpc/utils/test.py index f3d59280..ae493997 100644 --- a/python/dolfinx_mpc/utils/test.py +++ b/python/dolfinx_mpc/utils/test.py @@ -14,6 +14,7 @@ import scipy.sparse import dolfinx_mpc import dolfinx.common +from typing import Any @pytest.fixture @@ -66,7 +67,7 @@ def gather_constants(constraint, root=0): g_consts = MPI.COMM_WORLD.gather(constants[:l_range[1] - l_range[0]], root=root) if MPI.COMM_WORLD.rank == root: block_size = constraint.function_space().dofmap.index_map_bs - global_consts = np.zeros(imap.size_global * block_size, dtype=PETSc.ScalarType) + global_consts = np.zeros(imap.size_global * block_size, dtype=constraint.coefficients()[0].dtype) for (r, vals) in zip(ranges, g_consts): global_consts[r[0]:r[1]] = vals return global_consts @@ -144,7 +145,7 @@ def gather_transformation_matrix(constraint, root=0): cols.append(dof - sum(dof > all_slaves)) # Gather K to root - K_vals = MPI.COMM_WORLD.gather(np.asarray(K_val, dtype=PETSc.ScalarType), root=root) + K_vals = MPI.COMM_WORLD.gather(np.asarray(K_val, dtype=coeffs.dtype), root=root) rows_g = MPI.COMM_WORLD.gather(np.asarray(rows, dtype=np.int64), root=root) cols_g = MPI.COMM_WORLD.gather(np.asarray(cols, dtype=np.int64), root=root) @@ -153,20 +154,21 @@ def gather_transformation_matrix(constraint, root=0): return K_sparse -def petsc_to_local_CSR(A: PETSc.Mat, mpc: dolfinx_mpc.MultiPointConstraint): +def petsc_to_local_CSR(A: PETSc.Mat, mpc: dolfinx_mpc.MultiPointConstraint): # type: ignore """ Convert a PETSc matrix to a local CSR matrix (scipy) including ghost entries """ - global_indices = np.asarray(mpc.function_space.dofmap.index_map.global_indices(), dtype=PETSc.IntType) + global_indices = np.asarray(mpc.function_space.dofmap.index_map.global_indices(), + dtype=PETSc.IntType) # type: ignore sort_index = np.argsort(global_indices) - is_A = PETSc.IS().createGeneral(global_indices[sort_index]) + is_A = PETSc.IS().createGeneral(global_indices[sort_index]) # type: ignore A_loc = A.createSubMatrices(is_A)[0] ai, aj, av = A_loc.getValuesCSR() A_csr = scipy.sparse.csr_matrix((av, aj, ai)) return A_csr[global_indices[:, None], global_indices] -def gather_PETScMatrix(A: PETSc.Mat, root=0) -> scipy.sparse.csr_matrix: +def gather_PETScMatrix(A: PETSc.Mat, root=0) -> scipy.sparse.csr_matrix: # type: ignore """ Given a distributed PETSc matrix, gather in on process 'root' in a scipy CSR matrix @@ -184,7 +186,7 @@ def gather_PETScMatrix(A: PETSc.Mat, root=0) -> scipy.sparse.csr_matrix: (np.hstack(av_all), np.hstack(aj_all), ai_cum), shape=A.getSize()) # type: ignore -def gather_PETScVector(vector: PETSc.Vec, root=0) -> np.ndarray: +def gather_PETScVector(vector: PETSc.Vec, root=0) -> np.ndarray: # type: ignore """ Gather a PETScVector from different processors on process 'root' as an numpy array @@ -204,8 +206,9 @@ def compare_CSR(A: scipy.sparse.csr_matrix, B: scipy.sparse.csr_matrix, atol=1e- assert diff.max() < atol -def compare_mpc_lhs(A_org: PETSc.Mat, A_mpc: PETSc.Mat, - mpc: dolfinx_mpc.MultiPointConstraint, root: int = 0): +def compare_mpc_lhs(A_org: PETSc.Mat, A_mpc: PETSc.Mat, # type: ignore + mpc: dolfinx_mpc.MultiPointConstraint, root: int = 0, + atol: np.floating[Any] = 5e3 * np.finfo(dolfinx.default_scalar_type).resolution): """ Compare an unmodified matrix for the problem with the one assembled with a multi point constraint. @@ -216,14 +219,16 @@ def compare_mpc_lhs(A_org: PETSc.Mat, A_mpc: PETSc.Mat, comm = mpc.V.mesh.comm V = mpc.V assert root < comm.size - + is_complex = np.issubdtype(mpc.coefficients()[0].dtype, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 K = gather_transformation_matrix(mpc, root=root) A_csr = gather_PETScMatrix(A_org, root=root) - # Get global slaves glob_slaves = _gather_slaves_global(mpc) A_mpc_csr = gather_PETScMatrix(A_mpc, root=root) if MPI.COMM_WORLD.rank == root: + K = K.astype(scipy_dtype) + A_csr = A_csr.astype(scipy_dtype) KTAK = np.conj(K.T) * A_csr * K # Remove identity rows of MPC matrix @@ -233,12 +238,13 @@ def compare_mpc_lhs(A_org: PETSc.Mat, A_mpc: PETSc.Mat, mpc_without_slaves = A_mpc_csr[cols_except_slaves[:, None], cols_except_slaves] # Compute difference - compare_CSR(KTAK, mpc_without_slaves) + compare_CSR(KTAK, mpc_without_slaves, atol=atol) timer.stop() -def compare_mpc_rhs(b_org: PETSc.Vec, b: PETSc.Vec, constraint: dolfinx_mpc.MultiPointConstraint, root: int = 0): +def compare_mpc_rhs(b_org: PETSc.Vec, b: PETSc.Vec, # type: ignore + constraint: dolfinx_mpc.MultiPointConstraint, root: int = 0): """ Compare an unconstrained RHS with an MPC rhs. """ diff --git a/python/setup.py b/python/setup.py index b34e863c..f8fdf170 100644 --- a/python/setup.py +++ b/python/setup.py @@ -16,10 +16,10 @@ VERSION = "0.7.0" -REQUIREMENTS = ["numpy>=1.21", "fenics-dolfinx>0.6.0.dev0"] +REQUIREMENTS = ["numpy>=1.21", "fenics-dolfinx>=0.7.0"] extras = { - 'docs': ['jupyter-book'], + 'docs': ['jupyter-book', 'jupytext'], "test": ["pytest", "coverage"] } diff --git a/python/tests/mwe123.py b/python/tests/mwe123.py new file mode 100644 index 00000000..f9037882 --- /dev/null +++ b/python/tests/mwe123.py @@ -0,0 +1,29 @@ +from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags +from dolfinx import default_scalar_type +from dolfinx import fem +import dolfinx_mpc +from mpi4py import MPI +import numpy as np +mesh = create_unit_square(MPI.COMM_WORLD, 1, 1) +V = fem.functionspace(mesh, ("Lagrange", 1)) + +dolfinx_mpc.cpp.mpc.test_tabulate(V._cpp_object) + + +# Create multipoint constraint +def periodic_relation(x): + out_x = np.copy(x) + out_x[0] = 1 - x[0] + return out_x + + +def PeriodicBoundary(x): + return np.isclose(x[0], 1) + + +facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, PeriodicBoundary) +arg_sort = np.argsort(facets) +mt = meshtags(mesh, mesh.topology.dim - 1, facets[arg_sort], np.full(len(facets), 2, dtype=np.int32)) + +mpc = dolfinx_mpc.MultiPointConstraint(V) +mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, [], default_scalar_type(1.)) diff --git a/python/tests/nitsche_ufl.py b/python/tests/nitsche_ufl.py index 9ddd450a..42890e18 100644 --- a/python/tests/nitsche_ufl.py +++ b/python/tests/nitsche_ufl.py @@ -101,7 +101,7 @@ def nitsche_ufl(mesh: dmesh.Mesh, mesh_data: Tuple[dmesh.MeshTags, int, int], g_vec = [i for i in range(mesh.geometry.dim)] g_vec[mesh.geometry.dim - 1] = gap - V = _fem.VectorFunctionSpace(mesh, ("CG", 1)) + V = _fem.functionspace(mesh, ("CG", 1)) u = _fem.Function(V) v = ufl.TestFunction(V) @@ -110,7 +110,7 @@ def nitsche_ufl(mesh: dmesh.Mesh, mesh_data: Tuple[dmesh.MeshTags, int, int], ds = ufl.Measure("ds", domain=mesh, metadata=metadata, subdomain_data=facet_marker) a = ufl.inner(sigma(u), epsilon(v)) * dx - zero = np.asarray([0, ] * mesh.geometry.dim, dtype=_PETSc.ScalarType) + zero = np.asarray([0, ] * mesh.geometry.dim, dtype=_PETSc.ScalarType) # type: ignore L = ufl.inner(_fem.Constant(mesh, zero), v) * dx # Derivation of one sided Nitsche with gap function @@ -179,7 +179,7 @@ def _u_D(x): # xdmf.write_mesh(mesh) # xdmf.close() - solver = _nls.petsc.NewtonSolver(mesh.comm, problem) + solver = _nls.petsc.NewtonSolver(mesh.comm, problem) # type: ignore null_space = rigid_motions_nullspace(V) solver.A.setNearNullSpace(null_space) @@ -201,11 +201,11 @@ def _u_initial(x): # Define solver and options ksp = solver.krylov_solver - opts = _PETSc.Options() - option_prefix = ksp.getOptionsPrefix() + opts = _PETSc.Options() # type: ignore + option_prefix = ksp.getOptionsPrefix() # type: ignore # Set PETSc options - opts = _PETSc.Options() + opts = _PETSc.Options() # type: ignore opts.prefixPush(option_prefix) for k, v in petsc_options.items(): opts[k] = v diff --git a/python/tests/test_cube_contact.py b/python/tests/test_cube_contact.py index e5e7cec2..7ef5bea6 100644 --- a/python/tests/test_cube_contact.py +++ b/python/tests/test_cube_contact.py @@ -8,19 +8,20 @@ # between two cubes. import dolfinx.fem as fem -import dolfinx_mpc.utils import gmsh import numpy as np import pytest import scipy.sparse.linalg import ufl +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.io import gmshio -from dolfinx_mpc.utils import get_assemblers # noqa: F401 from mpi4py import MPI from petsc4py import PETSc import dolfinx_mpc +import dolfinx_mpc.utils +from dolfinx_mpc.utils import get_assemblers # noqa: F401 theta = np.pi / 5 @@ -212,10 +213,10 @@ def top_v(x): bcs = [bc_bottom, bc_top] # Elasticity parameters - E = PETSc.ScalarType(1.0e3) + E = 1.0e3 nu = 0 - mu = fem.Constant(mesh, E / (2.0 * (1.0 + nu))) - lmbda = fem.Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) + mu = fem.Constant(mesh, default_scalar_type(E / (2.0 * (1.0 + nu)))) + lmbda = fem.Constant(mesh, default_scalar_type(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))) # Stress computation def sigma(v): @@ -226,7 +227,7 @@ def sigma(v): u = ufl.TrialFunction(V) v = ufl.TestFunction(V) a = ufl.inner(sigma(u), ufl.grad(v)) * ufl.dx - rhs = ufl.inner(fem.Constant(mesh, PETSc.ScalarType((0, 0, 0))), v) * ufl.dx + rhs = ufl.inner(fem.Constant(mesh, default_scalar_type((0, 0, 0))), v) * ufl.dx bilinear_form = fem.form(a) linear_form = fem.form(rhs) @@ -241,11 +242,11 @@ def sigma(v): mpc = dolfinx_mpc.MultiPointConstraint(V) if nonslip: with Timer("~Contact: Create non-elastic constraint"): - mpc.create_contact_inelastic_condition(mt, 4, 9) + mpc.create_contact_inelastic_condition(mt, 4, 9, eps2=500 * np.finfo(default_scalar_type).resolution) else: with Timer("~Contact: Create contact constraint"): nh = dolfinx_mpc.utils.create_normal_approximation(V, mt, 4) - mpc.create_contact_slip_condition(mt, 4, 9, nh) + mpc.create_contact_slip_condition(mt, 4, 9, nh, eps2=500 * np.finfo(default_scalar_type).resolution) mpc.finalize() @@ -297,7 +298,8 @@ def sigma(v): d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + atol = 1000 * np.finfo(default_scalar_type).resolution + assert np.allclose(uh_numpy, u_mpc, atol=atol) L_org.destroy() list_timings(comm, [TimingType.wall]) diff --git a/python/tests/test_integration_domains.py b/python/tests/test_integration_domains.py index 14f58e8a..815bfb22 100644 --- a/python/tests/test_integration_domains.py +++ b/python/tests/test_integration_domains.py @@ -5,18 +5,19 @@ # SPDX-License-Identifier: MIT -from dolfinx import fem -import dolfinx_mpc import numpy as np import pytest import scipy.sparse.linalg import ufl +from dolfinx import default_scalar_type, fem from dolfinx.common import Timer, TimingType, list_timings from dolfinx.mesh import compute_midpoints, create_unit_square, meshtags -from dolfinx_mpc.utils import get_assemblers # noqa: F401 from mpi4py import MPI from petsc4py import PETSc +import dolfinx_mpc +from dolfinx_mpc.utils import get_assemblers # noqa: F401 + @pytest.mark.parametrize("get_assemblers", ["C++", "numba"], indirect=True) def test_cell_domains(get_assemblers): # noqa: F811 @@ -27,7 +28,7 @@ def test_cell_domains(get_assemblers): # noqa: F811 N = 5 # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 15, N) - V = fem.FunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1)) def left_side(x): return x[0] < 0.5 @@ -44,16 +45,16 @@ def left_side(x): u = ufl.TrialFunction(V) v = ufl.TestFunction(V) x = ufl.SpatialCoordinate(mesh) - c1 = fem.Constant(mesh, PETSc.ScalarType(2)) - c2 = fem.Constant(mesh, PETSc.ScalarType(10)) + c1 = fem.Constant(mesh, default_scalar_type(2)) + c2 = fem.Constant(mesh, default_scalar_type(10)) dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct) a = c1 * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(1) +\ c2 * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(2)\ - + 0.01 * ufl.inner(u, v) * dx(1) + + 0.87 * ufl.inner(u, v) * dx(1) rhs = ufl.inner(x[1], v) * dx(1) + \ - ufl.inner(fem.Constant(mesh, PETSc.ScalarType(1)), v) * dx(2) + ufl.inner(fem.Constant(mesh, default_scalar_type(1)), v) * dx(2) bilinear_form = fem.form(a) linear_form = fem.form(rhs) @@ -64,7 +65,7 @@ def left_side(x): L_org.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {} for i in range(0, N + 1): @@ -80,9 +81,10 @@ def l2b(li): b = assemble_vector(linear_form, mpc) b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) solver.setType(PETSc.KSP.Type.PREONLY) - solver.getPC().setType(PETSc.PC.Type.LU) + pc = solver.getPC() + pc.setType(PETSc.PC.Type.LU) solver.setOperators(A) # Solve @@ -94,23 +96,32 @@ def l2b(li): root = 0 comm = mesh.comm + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore with Timer("~TEST: Compare"): dolfinx_mpc.utils.compare_mpc_lhs(A_org, A, mpc, root=root) dolfinx_mpc.utils.compare_mpc_rhs(L_org, b, mpc, root=root) - # Gather LHS, RHS and solution on one process + # Gather LHS, RHS and solution on one process (work with high precision ) + scipy_dtype = np.complex128 if is_complex else np.float64 A_csr = dolfinx_mpc.utils.gather_PETScMatrix(A_org, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(L_org, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(uh.vector, root=root) if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system - d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) + d = scipy.sparse.linalg.spsolve(KTAK.astype(scipy_dtype), reduced_L.astype(scipy_dtype)) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) - uh.vector.destroy() + uh_numpy = K.astype(scipy_dtype) @ d.astype(scipy_dtype) + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, rtol=500 + * np.finfo(default_scalar_type).resolution) + solver.destroy() + b.destroy() + del uh + A.destroy() + pc.destroy() + L_org.destroy() + A_org.destroy() list_timings(comm, [TimingType.wall]) diff --git a/python/tests/test_lifting.py b/python/tests/test_lifting.py index 885015f0..22123fa6 100644 --- a/python/tests/test_lifting.py +++ b/python/tests/test_lifting.py @@ -29,7 +29,7 @@ def test_lifting(get_assemblers): # noqa: F811 # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 1, 1, CellType.quadrilateral) - V = fem.FunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1)) # Solve Problem without MPC for reference u = ufl.TrialFunction(V) @@ -67,7 +67,7 @@ def dirichletboundary(x): # Create multipoint constraint def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([0, 0]): {l2b([0, 1]): 1}} mpc = dolfinx_mpc.MultiPointConstraint(V) @@ -81,7 +81,7 @@ def l2b(li): fem.petsc.set_bc(b, bcs) - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) solver.setType(PETSc.KSP.Type.PREONLY) solver.getPC().setType(PETSc.PC.Type.LU) solver.setOperators(A) diff --git a/python/tests/test_linear_problem.py b/python/tests/test_linear_problem.py index 9d606102..4d6ccd48 100644 --- a/python/tests/test_linear_problem.py +++ b/python/tests/test_linear_problem.py @@ -5,29 +5,30 @@ # SPDX-License-Identifier: MIT -import dolfinx_mpc -import dolfinx_mpc.utils import numpy as np import pytest import scipy.sparse.linalg import ufl -from dolfinx import fem +from dolfinx import default_scalar_type, fem from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags from mpi4py import MPI from petsc4py import PETSc +import dolfinx_mpc +import dolfinx_mpc.utils + @pytest.mark.parametrize("u_from_mpc", [True, False]) def test_pipeline(u_from_mpc): # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 5, 5) - V = fem.FunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1)) # Solve Problem without MPC for reference u = ufl.TrialFunction(V) v = ufl.TestFunction(V) - d = fem.Constant(mesh, PETSc.ScalarType(0.01)) + d = fem.Constant(mesh, default_scalar_type(0.08)) x = ufl.SpatialCoordinate(mesh) f = ufl.sin(2 * ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - d * ufl.inner(u, v) * ufl.dx @@ -69,19 +70,25 @@ def PeriodicBoundary(x): dolfinx_mpc.utils.compare_mpc_rhs(L_org, problem.b, mpc, root=root) # Gather LHS, RHS and solution on one process + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 A_csr = dolfinx_mpc.utils.gather_PETScMatrix(A_org, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(L_org, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(uh.vector, root=root) if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + uh_numpy = K.astype(scipy_dtype) @ d + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, + rtol=500 + * np.finfo(default_scalar_type).resolution, + atol=500 + * np.finfo(default_scalar_type).resolution) else: uh = fem.Function(V) diff --git a/python/tests/test_matrix_assembly.py b/python/tests/test_matrix_assembly.py index 5c1b136a..d5759349 100644 --- a/python/tests/test_matrix_assembly.py +++ b/python/tests/test_matrix_assembly.py @@ -28,7 +28,7 @@ def test_mpc_assembly(master_point, degree, celltype, get_assemblers): # noqa: # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 5, 3, celltype) - V = fem.FunctionSpace(mesh, ("Lagrange", degree)) + V = fem.functionspace(mesh, ("Lagrange", degree)) # Test against generated code and general assembler u = ufl.TrialFunction(V) @@ -37,7 +37,7 @@ def test_mpc_assembly(master_point, degree, celltype, get_assemblers): # noqa: bilinear_form = fem.form(a) def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([1, 0]): {l2b([0, 1]): 0.43, l2b([1, 1]): 0.11}, l2b([0, 0]): {l2b(master_point): 0.69}} mpc = dolfinx_mpc.MultiPointConstraint(V) @@ -63,14 +63,14 @@ def test_slave_on_same_cell(master_point, degree, celltype, get_assemblers): # # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 1, 8, celltype) - V = fem.FunctionSpace(mesh, ("Lagrange", degree)) + V = fem.functionspace(mesh, ("Lagrange", degree)) # Build master slave map - s_m_c = {np.array([1, 0], dtype=np.float64).tobytes(): - {np.array([0, 1], dtype=np.float64).tobytes(): 0.43, - np.array([1, 1], dtype=np.float64).tobytes(): 0.11}, - np.array([0, 0], dtype=np.float64).tobytes(): { - np.array(master_point, dtype=np.float64).tobytes(): 0.69}} + s_m_c = {np.array([1, 0], dtype=mesh.geometry.x.dtype).tobytes(): + {np.array([0, 1], dtype=mesh.geometry.x.dtype).tobytes(): 0.43, + np.array([1, 1], dtype=mesh.geometry.x.dtype).tobytes(): 0.11}, + np.array([0, 0], dtype=mesh.geometry.x.dtype).tobytes(): { + np.array(master_point, dtype=mesh.geometry.x.dtype).tobytes(): 0.69}} with Timer("~TEST: MPC INIT"): mpc = dolfinx_mpc.MultiPointConstraint(V) mpc.create_general_constraint(s_m_c) diff --git a/python/tests/test_mpc_pipeline.py b/python/tests/test_mpc_pipeline.py index ad56f0be..5bb3dbef 100644 --- a/python/tests/test_mpc_pipeline.py +++ b/python/tests/test_mpc_pipeline.py @@ -5,19 +5,20 @@ # SPDX-License-Identifier: MIT -import dolfinx_mpc -import dolfinx_mpc.utils import numpy as np import pytest import scipy.sparse.linalg import ufl -from dolfinx import fem +from dolfinx import default_scalar_type, fem from dolfinx.common import Timer, TimingType, list_timings from dolfinx.mesh import create_unit_square -from dolfinx_mpc.utils import get_assemblers # noqa: F401 from mpi4py import MPI from petsc4py import PETSc +import dolfinx_mpc +import dolfinx_mpc.utils +from dolfinx_mpc.utils import get_assemblers # noqa: F401 + @pytest.mark.parametrize("get_assemblers", ["C++", "numba"], indirect=True) @pytest.mark.parametrize("master_point", [[1, 1], [0, 1]]) @@ -26,13 +27,13 @@ def test_pipeline(master_point, get_assemblers): # noqa: F811 # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 3, 5) - V = fem.FunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1)) # Solve Problem without MPC for reference u = ufl.TrialFunction(V) v = ufl.TestFunction(V) - d = fem.Constant(mesh, PETSc.ScalarType(1.5)) - c = fem.Constant(mesh, PETSc.ScalarType(2)) + d = fem.Constant(mesh, default_scalar_type(1.5)) + c = fem.Constant(mesh, default_scalar_type(2)) x = ufl.SpatialCoordinate(mesh) f = c * ufl.sin(2 * ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) g = fem.Function(V) @@ -53,7 +54,7 @@ def test_pipeline(master_point, get_assemblers): # noqa: F811 # Create multipoint constraint def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([1, 0]): {l2b([0, 1]): 0.43, l2b([1, 1]): 0.11}, l2b([0, 0]): {l2b(master_point): 0.69}} @@ -65,7 +66,7 @@ def l2b(li): b = assemble_vector(linear_form, mpc) b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) solver.setType(PETSc.KSP.Type.PREONLY) solver.getPC().setType(PETSc.PC.Type.LU) solver.setOperators(A) @@ -85,19 +86,22 @@ def l2b(li): dolfinx_mpc.utils.compare_mpc_rhs(L_org, b, mpc, root=root) # Gather LHS, RHS and solution on one process + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 A_csr = dolfinx_mpc.utils.gather_PETScMatrix(A_org, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(L_org, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(uh.vector, root=root) if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + uh_numpy = K.astype(scipy_dtype) @ d + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, rtol=500 + * np.finfo(default_scalar_type).resolution) list_timings(comm, [TimingType.wall]) @@ -107,13 +111,13 @@ def test_linearproblem(master_point): # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 3, 5) - V = fem.FunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1)) # Solve Problem without MPC for reference u = ufl.TrialFunction(V) v = ufl.TestFunction(V) - d = fem.Constant(mesh, PETSc.ScalarType(1.5)) - c = fem.Constant(mesh, PETSc.ScalarType(2)) + d = fem.Constant(mesh, default_scalar_type(1.5)) + c = fem.Constant(mesh, default_scalar_type(2)) x = ufl.SpatialCoordinate(mesh) f = c * ufl.sin(2 * ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) g = fem.Function(V) @@ -131,7 +135,7 @@ def test_linearproblem(master_point): # Create multipoint constraint def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([1, 0]): {l2b([0, 1]): 0.43, l2b([1, 1]): 0.11}, @@ -149,22 +153,22 @@ def l2b(li): root = 0 comm = mesh.comm with Timer("~TEST: Compare"): - dolfinx_mpc.utils.compare_mpc_lhs(A_org, problem._A, mpc, root=root) - dolfinx_mpc.utils.compare_mpc_rhs(L_org, problem._b, mpc, root=root) - # Gather LHS, RHS and solution on one process + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 A_csr = dolfinx_mpc.utils.gather_PETScMatrix(A_org, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(L_org, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(uh.vector, root=root) if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + uh_numpy = K.astype(scipy_dtype) @ d + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, rtol=500 + * np.finfo(default_scalar_type).resolution) list_timings(comm, [TimingType.wall]) diff --git a/python/tests/test_nonlinear_assembly.py b/python/tests/test_nonlinear_assembly.py index e804df42..c1340ec7 100644 --- a/python/tests/test_nonlinear_assembly.py +++ b/python/tests/test_nonlinear_assembly.py @@ -26,7 +26,7 @@ def __init__(self, F, u, mpc, bcs=[], J=None, form_compiler_options={}, form_compiler_options=form_compiler_options, jit_options=jit_options) - def F(self, x: PETSc.Vec, F: PETSc.Vec): + def F(self, x: PETSc.Vec, F: PETSc.Vec): # type: ignore with F.localForm() as F_local: F_local.set(0.0) dolfinx_mpc.assemble_vector(self._L, self.mpc, b=F) @@ -34,10 +34,10 @@ def F(self, x: PETSc.Vec, F: PETSc.Vec): # Apply boundary condition dolfinx_mpc.apply_lifting(F, [self._a], bcs=[self.bcs], constraint=self.mpc, x0=[x], scale=-1.0) - F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) + F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # type: ignore dolfinx.fem.petsc.set_bc(F, self.bcs, x, -1.0) - def J(self, x: PETSc.Vec, A: PETSc.Mat): + def J(self, x: PETSc.Vec, A: PETSc.Mat): # type: ignore A.zeroEntries() dolfinx_mpc.assemble_matrix(self._a, self.mpc, bcs=self.bcs, A=A) A.assemble() @@ -65,17 +65,17 @@ def __init__(self, comm: MPI.Intracomm, problem: NonlinearMPCProblem, self.set_update(self.update) def update(self, solver: dolfinx.nls.petsc.NewtonSolver, - dx: PETSc.Vec, x: PETSc.Vec): + dx: PETSc.Vec, x: PETSc.Vec): # type: ignore # We need to use a vector created on the MPC's space to update ghosts self.u_mpc.vector.array = x.array_r self.u_mpc.vector.axpy(-1.0, dx) - self.u_mpc.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, - mode=PETSc.ScatterMode.FORWARD) + self.u_mpc.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, # type: ignore + mode=PETSc.ScatterMode.FORWARD) # type: ignore self.mpc.homogenize(self.u_mpc) self.mpc.backsubstitution(self.u_mpc) x.array = self.u_mpc.vector.array_r - x.ghostUpdate(addv=PETSc.InsertMode.INSERT, - mode=PETSc.ScatterMode.FORWARD) + x.ghostUpdate(addv=PETSc.InsertMode.INSERT, # type: ignore + mode=PETSc.ScatterMode.FORWARD) # type: ignore def solve(self, u: dolfinx.fem.Function): """Solve non-linear problem into function u. Returns the number @@ -85,17 +85,17 @@ def solve(self, u: dolfinx.fem.Function): return n, converged @property - def A(self) -> PETSc.Mat: + def A(self) -> PETSc.Mat: # type: ignore """Jacobian matrix""" return self._A @property - def b(self) -> PETSc.Vec: + def b(self) -> PETSc.Vec: # type: ignore """Residual vector""" return self._b -@pytest.mark.skipif(np.issubdtype(PETSc.ScalarType, np.complexfloating), +@pytest.mark.skipif(np.issubdtype(dolfinx.default_scalar_type, np.complexfloating), reason="This test does not work in complex mode.") @pytest.mark.parametrize("poly_order", [1, 2, 3]) def test_nonlinear_poisson(poly_order): @@ -110,19 +110,17 @@ def test_nonlinear_poisson(poly_order): for run_no, N in enumerate(N_vals): mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) - V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", poly_order)) + V = dolfinx.fem.functionspace(mesh, ("Lagrange", poly_order)) def boundary(x): return np.ones_like(x[0], dtype=np.int8) u_bc = dolfinx.fem.Function(V) - with u_bc.vector.localForm() as u_local: - u_local.set(0.0) - u_bc.vector.destroy() + u_bc.x.array[:] = 0.0 facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, boundary) topological_dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets) - zero = np.array(0, dtype=PETSc.ScalarType) + zero = np.array(0, dtype=dolfinx.default_scalar_type) bc = dolfinx.fem.dirichletbc(zero, topological_dofs, V) bcs = [bc] @@ -141,12 +139,12 @@ def boundary(x): # -- Impose the pi/2 rotational symmetry of the solution as a constraint, # -- except at the centre DoF def periodic_boundary(x): - eps = 1e-10 - return np.isclose(x[0], 0.5) & ( + eps = 1000 * np.finfo(x.dtype).resolution + return np.isclose(x[0], 0.5, atol=eps) & ( (x[1] < 0.5 - eps) | (x[1] > 0.5 + eps)) def periodic_relation(x): - out_x = np.zeros(x.shape) + out_x = np.zeros_like(x) out_x[0] = x[1] out_x[1] = x[0] out_x[2] = x[2] @@ -167,6 +165,8 @@ def periodic_relation(x): problem = NonlinearMPCProblem(F, u, mpc, bcs=bcs, J=J) solver = NewtonSolverMPC(mesh.comm, problem, mpc) + solver.atol = 1e1 * np.finfo(u.x.array.dtype).resolution + solver.rtol = 1e1 * np.finfo(u.x.array.dtype).resolution # Ensure the solver works with nonzero initial guess u.interpolate(lambda x: x[0]**2 * x[1]**2) @@ -199,7 +199,7 @@ def test_homogenize(tensor_order, poly_order): cellname = mesh.ufl_cell().cellname() el = basix.ufl.element(basix.ElementFamily.P, cellname, poly_order, shape=shape) - V = dolfinx.fem.FunctionSpace(mesh, el) + V = dolfinx.fem.functionspace(mesh, el) def periodic_boundary(x): return np.isclose(x[0], 0.0) diff --git a/python/tests/test_rectangular_assembly.py b/python/tests/test_rectangular_assembly.py index 57a8fc70..382ae8ce 100644 --- a/python/tests/test_rectangular_assembly.py +++ b/python/tests/test_rectangular_assembly.py @@ -31,9 +31,9 @@ def test_mixed_element(cell_type, ghost_mode): # Inlet velocity Dirichlet BC bc_facets = dolfinx.mesh.locate_entities_boundary( - mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 0.0)) + mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 0.0, atol=500 * np.finfo(x.dtype).resolution)) other_facets = dolfinx.mesh.locate_entities_boundary( - mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 1.0)) + mesh, mesh.topology.dim - 1, lambda x: np.isclose(x[0], 1.0, atol=500 * np.finfo(x.dtype).resolution)) arg_sort = np.argsort(other_facets) mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, other_facets[arg_sort], np.full_like(other_facets, 1)) @@ -52,7 +52,7 @@ def test_mixed_element(cell_type, ghost_mode): V = dolfinx.fem.functionspace(mesh, Ve) Q = dolfinx.fem.functionspace(mesh, Qe) - W = dolfinx.fem.functionspace(mesh, Ve * Qe) + W = dolfinx.fem.functionspace(mesh, basix.ufl.mixed_element([Ve, Qe])) inlet_velocity = dolfinx.fem.Function(V) inlet_velocity.interpolate( @@ -73,7 +73,7 @@ def test_mixed_element(cell_type, ghost_mode): mpc_q = dolfinx_mpc.MultiPointConstraint(Q) mpc_q.finalize() - f = dolfinx.fem.Constant(mesh, PETSc.ScalarType((0, 0))) + f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0))) (u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q) (v, q) = ufl.TestFunction(V), ufl.TestFunction(Q) a00 = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx @@ -83,7 +83,7 @@ def test_mixed_element(cell_type, ghost_mode): L0 = ufl.inner(f, v) * ufl.dx L1 = ufl.inner( - dolfinx.fem.Constant(mesh, PETSc.ScalarType(0.0)), q) * ufl.dx + dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0)), q) * ufl.dx n = ufl.FacetNormal(mesh) g_tau = ufl.as_vector((0.0, 0.0)) @@ -141,7 +141,7 @@ def test_mixed_element(cell_type, ghost_mode): mpc_vq.create_slip_constraint(W.sub(0), (mt, 1), n_approx, bcs=bcs) mpc_vq.finalize() - f = dolfinx.fem.Constant(mesh, PETSc.ScalarType((0, 0))) + f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0, 0))) (u, p) = ufl.TrialFunctions(W) (v, q) = ufl.TestFunctions(W) a = ( @@ -151,7 +151,7 @@ def test_mixed_element(cell_type, ghost_mode): ) L = ufl.inner(f, v) * ufl.dx + ufl.inner( - dolfinx.fem.Constant(mesh, PETSc.ScalarType(0.0)), q) * ufl.dx + dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0)), q) * ufl.dx # No prescribed shear stress n = ufl.FacetNormal(mesh) diff --git a/python/tests/test_surface_integral.py b/python/tests/test_surface_integral.py index f338e602..f1f2ef10 100644 --- a/python/tests/test_surface_integral.py +++ b/python/tests/test_surface_integral.py @@ -6,18 +6,20 @@ import dolfinx.fem as fem -import dolfinx_mpc -import dolfinx_mpc.utils import numpy as np import pytest import scipy.sparse.linalg import ufl +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags -from dolfinx_mpc.utils import get_assemblers # noqa: F401 from mpi4py import MPI from petsc4py import PETSc +import dolfinx_mpc +import dolfinx_mpc.utils +from dolfinx_mpc.utils import get_assemblers # noqa: F401 + @pytest.mark.parametrize("get_assemblers", ["C++", "numba"], indirect=True) def test_surface_integrals(get_assemblers): # noqa: F811 @@ -30,34 +32,31 @@ def test_surface_integrals(get_assemblers): # noqa: F811 # Fixed Dirichlet BC on the left wall def left_wall(x): - return np.isclose(x[0], np.finfo(float).eps) + return np.isclose(x[0], 0, atol=100 * np.finfo(x.dtype).resolution) fdim = mesh.topology.dim - 1 left_facets = locate_entities_boundary(mesh, fdim, left_wall) bc_dofs = fem.locate_dofs_topological(V, 1, left_facets) u_bc = fem.Function(V) - with u_bc.vector.localForm() as u_local: - u_local.set(0.0) - u_bc.vector.destroy() - + u_bc.x.array[:] = 0 bc = fem.dirichletbc(u_bc, bc_dofs) bcs = [bc] # Traction on top of domain def top(x): - return np.isclose(x[1], 1) + return np.isclose(x[1], 1, atol=100 * np.finfo(x.dtype).resolution) top_facets = locate_entities_boundary(mesh, 1, top) arg_sort = np.argsort(top_facets) mt = meshtags(mesh, fdim, top_facets[arg_sort], np.full(len(top_facets), 3, dtype=np.int32)) ds = ufl.Measure("ds", domain=mesh, subdomain_data=mt, subdomain_id=3) - g = fem.Constant(mesh, PETSc.ScalarType((0, -9.81e2))) + g = fem.Constant(mesh, default_scalar_type((0, -9.81e2))) # Elasticity parameters - E = PETSc.ScalarType(1.0e4) + E = 1.0e2 nu = 0.0 - mu = fem.Constant(mesh, E / (2.0 * (1.0 + nu))) - lmbda = fem.Constant(mesh, E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) + mu = fem.Constant(mesh, default_scalar_type(E / (2.0 * (1.0 + nu)))) + lmbda = fem.Constant(mesh, default_scalar_type(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))) # Stress computation def sigma(v): @@ -68,25 +67,26 @@ def sigma(v): u = ufl.TrialFunction(V) v = ufl.TestFunction(V) a = ufl.inner(sigma(u), ufl.grad(v)) * ufl.dx - rhs = ufl.inner(fem.Constant(mesh, PETSc.ScalarType((0, 0))), v) * ufl.dx\ + rhs = ufl.inner(fem.Constant(mesh, default_scalar_type((0, 0))), v) * ufl.dx\ + ufl.inner(g, v) * ds bilinear_form = fem.form(a) linear_form = fem.form(rhs) # Setup LU solver - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) solver.setType(PETSc.KSP.Type.PREONLY) solver.getPC().setType(PETSc.PC.Type.LU) # Setup multipointconstraint def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {} for i in range(1, N): s_m_c[l2b([1, i / N])] = {l2b([1, 1]): 0.8} mpc = dolfinx_mpc.MultiPointConstraint(V) mpc.create_general_constraint(s_m_c, 1, 1) mpc.finalize() + with Timer("~TEST: Assemble matrix old"): A = assemble_matrix(bilinear_form, mpc, bcs=bcs) with Timer("~TEST: Assemble vector"): @@ -116,23 +116,24 @@ def l2b(li): root = 0 comm = mesh.comm with Timer("~TEST: Compare"): - dolfinx_mpc.utils.compare_mpc_lhs(A_org, A, mpc, root=root) - dolfinx_mpc.utils.compare_mpc_rhs(L_org, b, mpc, root=root) - # Gather LHS, RHS and solution on one process + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 A_csr = dolfinx_mpc.utils.gather_PETScMatrix(A_org, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(L_org, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(uh.vector, root=root) if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + uh_numpy = K.astype(scipy_dtype) @ d + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, rtol=500 + * np.finfo(default_scalar_type).resolution) + L_org.destroy() b.destroy() list_timings(comm, [TimingType.wall]) @@ -161,8 +162,8 @@ def top(x): mt = meshtags(mesh, mesh.topology.dim - 1, np.array(indices[sort], dtype=np.intc), np.array(values[sort], dtype=np.intc)) ds = ufl.Measure("ds", domain=mesh, subdomain_data=mt) - g = fem.Constant(mesh, PETSc.ScalarType((2, 1))) - h = fem.Constant(mesh, PETSc.ScalarType((3, 2))) + g = fem.Constant(mesh, default_scalar_type((2, 1))) + h = fem.Constant(mesh, default_scalar_type((3, 2))) # Define variational problem u = ufl.TrialFunction(V) v = ufl.TestFunction(V) @@ -173,7 +174,7 @@ def top(x): # Create multipoint constraint and assemble system def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {} for i in range(1, N): s_m_c[l2b([1, i / N])] = {l2b([1, 1]): 0.3} diff --git a/python/tests/test_vector_assembly.py b/python/tests/test_vector_assembly.py index 3b85e954..b8a18fda 100644 --- a/python/tests/test_vector_assembly.py +++ b/python/tests/test_vector_assembly.py @@ -28,7 +28,7 @@ def test_mpc_assembly(master_point, degree, celltype, get_assemblers): # noqa: # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, 3, 5, celltype) - V = fem.FunctionSpace(mesh, ("Lagrange", degree)) + V = fem.functionspace(mesh, ("Lagrange", degree)) # Generate reference vector v = ufl.TestFunction(V) @@ -38,7 +38,7 @@ def test_mpc_assembly(master_point, degree, celltype, get_assemblers): # noqa: linear_form = fem.form(rhs) def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([1, 0]): {l2b([0, 1]): 0.43, l2b([1, 1]): 0.11}, l2b([0, 0]): {l2b(master_point): 0.69}} diff --git a/python/tests/test_vector_poisson.py b/python/tests/test_vector_poisson.py index 0d9827c7..ca2c03af 100644 --- a/python/tests/test_vector_poisson.py +++ b/python/tests/test_vector_poisson.py @@ -5,18 +5,20 @@ # SPDX-License-Identifier: MIT import dolfinx.fem as fem -import dolfinx_mpc -import dolfinx_mpc.utils import numpy as np import pytest import scipy.sparse.linalg import ufl +from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.mesh import create_unit_square -from dolfinx_mpc.utils import get_assemblers # noqa: F401 from mpi4py import MPI from petsc4py import PETSc +import dolfinx_mpc +import dolfinx_mpc.utils +from dolfinx_mpc.utils import get_assemblers # noqa: F401 + @pytest.mark.parametrize("get_assemblers", ["C++", "numba"], indirect=True) @pytest.mark.parametrize("Nx", [4]) @@ -32,7 +34,7 @@ def test_vector_possion(Nx, Ny, slave_space, master_space, get_assemblers): # n V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,))) def boundary(x): - return np.isclose(x.T, [0, 0, 0]).all(axis=1) + return np.isclose(x.T, [0, 0, 0], atol=500 * np.finfo(x.dtype).resolution).all(axis=1) # Define boundary conditions (HAS TO BE NON-MASTER NODES) u_bc = fem.Function(V) @@ -56,13 +58,15 @@ def boundary(x): linear_form = fem.form(rhs) # Setup LU solver - solver = PETSc.KSP().create(MPI.COMM_WORLD) + solver = PETSc.KSP().create(mesh.comm) solver.setType(PETSc.KSP.Type.PREONLY) - solver.getPC().setType(PETSc.PC.Type.LU) + pc = solver.getPC() + pc.setType(PETSc.PC.Type.LU) + pc.setFactorSolverType("mumps") # Create multipoint constraint def l2b(li): - return np.array(li, dtype=np.float64).tobytes() + return np.array(li, dtype=mesh.geometry.x.dtype).tobytes() s_m_c = {l2b([1, 0]): {l2b([1, 1]): 0.1, l2b([0.5, 1]): 0.3}} mpc = dolfinx_mpc.MultiPointConstraint(V) mpc.create_general_constraint(s_m_c, slave_space, master_space) @@ -101,19 +105,23 @@ def l2b(li): dolfinx_mpc.utils.compare_mpc_rhs(L_org, b, mpc, root=root) # Gather LHS, RHS and solution on one process + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 A_csr = dolfinx_mpc.utils.gather_PETScMatrix(A_org, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(L_org, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(uh.vector, root=root) if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc) + uh_numpy = K.astype(scipy_dtype) @ d + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, rtol=500 + * np.finfo(default_scalar_type).resolution) + b.destroy() L_org.destroy() list_timings(comm, [TimingType.wall])