From 30fa1ef434f78bcc6daf6ead6495d0c75203d3ad Mon Sep 17 00:00:00 2001 From: Musawer Ahmad Saqif Date: Sun, 28 Jun 2026 13:37:19 -0600 Subject: [PATCH] Check PETSc errors and validate BC assembly Add PETSc/SLEPc error checking in solver and linear algebra wrappers. Reject unsupported DirichletBC use for bilinear forms with different trial and test spaces, with MatrixCSR and PETSc regression coverage. Move the ADIOS2 availability check before importing VTXWriter in the optional test. --- cpp/dolfinx/fem/petsc.cpp | 8 ++- cpp/dolfinx/la/petsc.cpp | 98 +++++++++++++++++--------- cpp/dolfinx/la/slepc.cpp | 78 ++++++++++++-------- cpp/dolfinx/nls/NewtonSolver.cpp | 31 +++++--- python/dolfinx/fem/assemble.py | 28 +++++++- python/dolfinx/fem/petsc.py | 23 ++++-- python/test/unit/fem/test_assembler.py | 28 ++++++++ python/test/unit/io/test_adios2.py | 2 +- 8 files changed, 211 insertions(+), 85 deletions(-) diff --git a/cpp/dolfinx/fem/petsc.cpp b/cpp/dolfinx/fem/petsc.cpp index 305a20e9c2f..323c4ac80ea 100644 --- a/cpp/dolfinx/fem/petsc.cpp +++ b/cpp/dolfinx/fem/petsc.cpp @@ -54,7 +54,7 @@ Vec fem::petsc::create_vector_block( = VecCreateGhost(maps[0].first.get().comm(), local_size, PETSC_DETERMINE, _ghosts.size(), _ghosts.data(), &x); if (ierr != 0) - throw std::runtime_error("Call to PETSc VecCreateGhost failed."); + la::petsc::error(ierr, __FILE__, "VecCreateGhost"); return x; } @@ -76,8 +76,10 @@ Vec fem::petsc::create_vector_nest( // Create nested (VecNest) vector Vec y; - VecCreateNest(vecs.front()->comm(), petsc_vecs.size(), nullptr, - petsc_vecs.data(), &y); + PetscErrorCode ierr = VecCreateNest(vecs.front()->comm(), petsc_vecs.size(), + nullptr, petsc_vecs.data(), &y); + if (ierr != 0) + la::petsc::error(ierr, __FILE__, "VecCreateNest"); return y; } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/la/petsc.cpp b/cpp/dolfinx/la/petsc.cpp index cb23cfafc4e..eb41a6a0317 100644 --- a/cpp/dolfinx/la/petsc.cpp +++ b/cpp/dolfinx/la/petsc.cpp @@ -57,11 +57,15 @@ la::petsc::create_vectors(MPI_Comm comm, std::vector v(x.size()); for (std::size_t i = 0; i < v.size(); ++i) { - VecCreateMPI(comm, x[i].size(), PETSC_DETERMINE, &v[i]); - PetscScalar* data; - VecGetArray(v[i], &data); + PetscErrorCode ierr + = VecCreateMPI(comm, x[i].size(), PETSC_DETERMINE, &v[i]); + CHECK_ERROR("VecCreateMPI"); + PetscScalar* data = nullptr; + ierr = VecGetArray(v[i], &data); + CHECK_ERROR("VecGetArray"); std::ranges::copy(x[i], data); - VecRestoreArray(v[i], &data); + ierr = VecRestoreArray(v[i], &data); + CHECK_ERROR("VecRestoreArray"); } return v; @@ -140,7 +144,9 @@ std::vector la::petsc::create_index_sets( std::int32_t size = map.first.get().size_local() + map.first.get().num_ghosts(); IS _is; - ISCreateStride(PETSC_COMM_SELF, bs * size, offset, 1, &_is); + PetscErrorCode ierr + = ISCreateStride(PETSC_COMM_SELF, bs * size, offset, 1, &_is); + CHECK_ERROR("ISCreateStride"); is.push_back(_is); offset += bs * size; } @@ -160,11 +166,14 @@ std::vector> la::petsc::get_local_vectors( // Unwrap PETSc vector Vec x_local; - VecGhostGetLocalForm(x, &x_local); + PetscErrorCode ierr = VecGhostGetLocalForm(x, &x_local); + CHECK_ERROR("VecGhostGetLocalForm"); PetscInt n = 0; - VecGetSize(x_local, &n); + ierr = VecGetSize(x_local, &n); + CHECK_ERROR("VecGetSize"); const PetscScalar* array = nullptr; - VecGetArrayRead(x_local, &array); + ierr = VecGetArrayRead(x_local, &array); + CHECK_ERROR("VecGetArrayRead"); std::span _x(array, n); // Copy PETSc Vec data in to local vectors @@ -185,8 +194,10 @@ std::vector> la::petsc::get_local_vectors( offset_ghost += size_ghost; } - VecRestoreArrayRead(x_local, &array); - VecGhostRestoreLocalForm(x, &x_local); + ierr = VecRestoreArrayRead(x_local, &array); + CHECK_ERROR("VecRestoreArrayRead"); + ierr = VecGhostRestoreLocalForm(x, &x_local); + CHECK_ERROR("VecGhostRestoreLocalForm"); return x_b; } @@ -205,11 +216,14 @@ void la::petsc::scatter_local_vectors( offset_owned += map.first.get().size_local() * map.second; Vec x_local; - VecGhostGetLocalForm(x, &x_local); + PetscErrorCode ierr = VecGhostGetLocalForm(x, &x_local); + CHECK_ERROR("VecGhostGetLocalForm"); PetscInt n = 0; - VecGetSize(x_local, &n); + ierr = VecGetSize(x_local, &n); + CHECK_ERROR("VecGetSize"); PetscScalar* array = nullptr; - VecGetArray(x_local, &array); + ierr = VecGetArray(x_local, &array); + CHECK_ERROR("VecGetArray"); std::span _x(array, n); // Copy local vectors into PETSc Vec @@ -228,8 +242,10 @@ void la::petsc::scatter_local_vectors( offset_ghost += size_ghost; } - VecRestoreArray(x_local, &array); - VecGhostRestoreLocalForm(x, &x_local); + ierr = VecRestoreArray(x_local, &array); + CHECK_ERROR("VecRestoreArray"); + ierr = VecGhostRestoreLocalForm(x, &x_local); + CHECK_ERROR("VecGhostRestoreLocalForm"); } //----------------------------------------------------------------------------- Mat la::petsc::create_matrix(MPI_Comm comm, const SparsityPattern& sp, @@ -245,8 +261,12 @@ Mat la::petsc::create_matrix(MPI_Comm comm, const SparsityPattern& sp, std::array maps = {sp.index_map(0), sp.index_map(1)}; const std::array bs = {sp.block_size(0), sp.block_size(1)}; - if (type) - MatSetType(A, type->c_str()); + if (type and !type->empty()) + { + ierr = MatSetType(A, type->c_str()); + if (ierr != 0) + petsc::error(ierr, __FILE__, "MatSetType"); + } // Get global and local dimensions const std::int64_t M = bs[0] * maps[0]->size_global(); @@ -294,7 +314,7 @@ Mat la::petsc::create_matrix(MPI_Comm comm, const SparsityPattern& sp, ierr = MatXAIJSetPreallocation(A, _bs, _nnz_diag.data(), _nnz_offdiag.data(), nullptr, nullptr); if (ierr != 0) - petsc::error(ierr, __FILE__, "MatXIJSetPreallocation"); + petsc::error(ierr, __FILE__, "MatXAIJSetPreallocation"); // Set block sizes ierr = MatSetBlockSizes(A, bs[0], bs[1]); @@ -405,7 +425,10 @@ petsc::Vector::Vector(Vec x, bool inc_ref_count) : _x(x) { assert(x); if (inc_ref_count) - PetscObjectReference((PetscObject)_x); + { + PetscErrorCode ierr = PetscObjectReference((PetscObject)_x); + CHECK_ERROR("PetscObjectReference"); + } } //----------------------------------------------------------------------------- petsc::Vector::Vector(Vector&& v) noexcept : _x(std::exchange(v._x, nullptr)) {} @@ -413,7 +436,7 @@ petsc::Vector::Vector(Vector&& v) noexcept : _x(std::exchange(v._x, nullptr)) {} petsc::Vector::~Vector() { if (_x) - VecDestroy(&_x); + (void)VecDestroy(&_x); } //----------------------------------------------------------------------------- petsc::Vector& petsc::Vector::operator=(Vector&& v) noexcept @@ -425,10 +448,13 @@ petsc::Vector& petsc::Vector::operator=(Vector&& v) noexcept petsc::Vector petsc::Vector::copy() const { Vec _y; - VecDuplicate(_x, &_y); - VecCopy(_x, _y); + PetscErrorCode ierr = VecDuplicate(_x, &_y); + CHECK_ERROR("VecDuplicate"); + ierr = VecCopy(_x, _y); + CHECK_ERROR("VecCopy"); Vector y(_y, true); - VecDestroy(&_y); + ierr = VecDestroy(&_y); + CHECK_ERROR("VecDestroy"); return y; } //----------------------------------------------------------------------------- @@ -499,7 +525,10 @@ petsc::Operator::Operator(Mat A, bool inc_ref_count) : _matA(A) { assert(A); if (inc_ref_count) - PetscObjectReference((PetscObject)_matA); + { + PetscErrorCode ierr = PetscObjectReference((PetscObject)_matA); + CHECK_ERROR("PetscObjectReference"); + } } //----------------------------------------------------------------------------- petsc::Operator::Operator(Operator&& A) noexcept @@ -512,7 +541,7 @@ petsc::Operator::~Operator() // Decrease reference count (PETSc will destroy object once reference // counts reached zero) if (_matA) - MatDestroy(&_matA); + (void)MatDestroy(&_matA); } //----------------------------------------------------------------------------- petsc::Operator& petsc::Operator::operator=(Operator&& A) noexcept @@ -527,7 +556,7 @@ std::array petsc::Operator::size() const PetscInt m(0), n(0); PetscErrorCode ierr = MatGetSize(_matA, &m, &n); if (ierr != 0) - petsc::error(ierr, __FILE__, "MetGetSize"); + petsc::error(ierr, __FILE__, "MatGetSize"); return {{m, n}}; } //----------------------------------------------------------------------------- @@ -620,21 +649,24 @@ void petsc::Matrix::apply(AssemblyType type) void petsc::Matrix::set_options_prefix(const std::string& options_prefix) { assert(_matA); - MatSetOptionsPrefix(_matA, options_prefix.c_str()); + PetscErrorCode ierr = MatSetOptionsPrefix(_matA, options_prefix.c_str()); + CHECK_ERROR("MatSetOptionsPrefix"); } //----------------------------------------------------------------------------- std::string petsc::Matrix::get_options_prefix() const { assert(_matA); const char* prefix = nullptr; - MatGetOptionsPrefix(_matA, &prefix); + PetscErrorCode ierr = MatGetOptionsPrefix(_matA, &prefix); + CHECK_ERROR("MatGetOptionsPrefix"); return std::string(prefix); } //----------------------------------------------------------------------------- void petsc::Matrix::set_from_options() { assert(_matA); - MatSetFromOptions(_matA); + PetscErrorCode ierr = MatSetFromOptions(_matA); + CHECK_ERROR("MatSetFromOptions"); } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- @@ -666,7 +698,7 @@ petsc::KrylovSolver::KrylovSolver(KrylovSolver&& solver) noexcept petsc::KrylovSolver::~KrylovSolver() { if (_ksp) - KSPDestroy(&_ksp); + (void)KSPDestroy(&_ksp); } //----------------------------------------------------------------------------- petsc::KrylovSolver& @@ -695,11 +727,11 @@ int petsc::KrylovSolver::solve(Vec x, const Vec b, bool transpose) const // Get PETSc operators Mat _A, _P; - KSPGetOperators(_ksp, &_A, &_P); + PetscErrorCode ierr = KSPGetOperators(_ksp, &_A, &_P); + if (ierr != 0) + petsc::error(ierr, __FILE__, "KSPGetOperators"); assert(_A); - PetscErrorCode ierr; - // Solve linear system spdlog::info("PETSc Krylov solver starting to solve system."); diff --git a/cpp/dolfinx/la/slepc.cpp b/cpp/dolfinx/la/slepc.cpp index a942d5ccb18..5ee14191a2f 100644 --- a/cpp/dolfinx/la/slepc.cpp +++ b/cpp/dolfinx/la/slepc.cpp @@ -17,8 +17,22 @@ using namespace dolfinx; using namespace dolfinx::la; +namespace +{ +//----------------------------------------------------------------------------- +void check_petsc_error(PetscErrorCode ierr, const char* petsc_function) +{ + if (ierr != 0) + petsc::error(ierr, __FILE__, petsc_function); +} +//----------------------------------------------------------------------------- +} // namespace + //----------------------------------------------------------------------------- -SLEPcEigenSolver::SLEPcEigenSolver(MPI_Comm comm) { EPSCreate(comm, &_eps); } +SLEPcEigenSolver::SLEPcEigenSolver(MPI_Comm comm) +{ + check_petsc_error(EPSCreate(comm, &_eps), "EPSCreate"); +} //----------------------------------------------------------------------------- SLEPcEigenSolver::SLEPcEigenSolver(EPS eps, bool inc_ref_count) : _eps(eps) { @@ -29,8 +43,7 @@ SLEPcEigenSolver::SLEPcEigenSolver(EPS eps, bool inc_ref_count) : _eps(eps) if (inc_ref_count) { ierr = PetscObjectReference((PetscObject)_eps); - if (ierr != 0) - petsc::error(ierr, __FILE__, "PetscObjectReference"); + check_petsc_error(ierr, "PetscObjectReference"); } } //----------------------------------------------------------------------------- @@ -43,7 +56,7 @@ SLEPcEigenSolver::SLEPcEigenSolver(SLEPcEigenSolver&& solver) noexcept SLEPcEigenSolver::~SLEPcEigenSolver() { if (_eps) - EPSDestroy(&_eps); + (void)EPSDestroy(&_eps); } //----------------------------------------------------------------------------- SLEPcEigenSolver& @@ -56,7 +69,7 @@ SLEPcEigenSolver::operator=(SLEPcEigenSolver&& solver) noexcept void SLEPcEigenSolver::set_operators(const Mat A, const Mat B) { assert(_eps); - EPSSetOperators(_eps, A, B); + check_petsc_error(EPSSetOperators(_eps, A, B), "EPSSetOperators"); } //----------------------------------------------------------------------------- void SLEPcEigenSolver::solve() @@ -64,10 +77,10 @@ void SLEPcEigenSolver::solve() // Get operators Mat A, B; assert(_eps); - EPSGetOperators(_eps, &A, &B); + check_petsc_error(EPSGetOperators(_eps, &A, &B), "EPSGetOperators"); PetscInt m(0), n(0); - MatGetSize(A, &m, &n); + check_petsc_error(MatGetSize(A, &m, &n), "MatGetSize"); solve(m); } //----------------------------------------------------------------------------- @@ -77,35 +90,38 @@ void SLEPcEigenSolver::solve(std::int64_t n) // Get operators Mat A, B; assert(_eps); - EPSGetOperators(_eps, &A, &B); + check_petsc_error(EPSGetOperators(_eps, &A, &B), "EPSGetOperators"); PetscInt _m(0), _n(0); - MatGetSize(A, &_m, &_n); + check_petsc_error(MatGetSize(A, &_m, &_n), "MatGetSize"); assert(n <= _n); #endif // Set number of eigenpairs to compute assert(_eps); - EPSSetDimensions(_eps, n, PETSC_DECIDE, PETSC_DECIDE); + check_petsc_error(EPSSetDimensions(_eps, n, PETSC_DECIDE, PETSC_DECIDE), + "EPSSetDimensions"); // Set any options from the PETSc database - EPSSetFromOptions(_eps); + check_petsc_error(EPSSetFromOptions(_eps), "EPSSetFromOptions"); // Solve eigenvalue problem - EPSSolve(_eps); + check_petsc_error(EPSSolve(_eps), "EPSSolve"); // Check for convergence EPSConvergedReason reason; - EPSGetConvergedReason(_eps, &reason); + check_petsc_error(EPSGetConvergedReason(_eps, &reason), + "EPSGetConvergedReason"); if (reason < 0) spdlog::warn("Eigenvalue solver did not converge"); // Report solver status PetscInt num_iterations = 0; - EPSGetIterationNumber(_eps, &num_iterations); + check_petsc_error(EPSGetIterationNumber(_eps, &num_iterations), + "EPSGetIterationNumber"); EPSType eps_type = nullptr; - EPSGetType(_eps, &eps_type); + check_petsc_error(EPSGetType(_eps, &eps_type), "EPSGetType"); spdlog::info("Eigenvalue solver ({}) converged in {} iterations.", eps_type, num_iterations); } @@ -116,17 +132,20 @@ std::complex SLEPcEigenSolver::get_eigenvalue(int i) const // Get number of computed values PetscInt num_computed_eigenvalues; - EPSGetConverged(_eps, &num_computed_eigenvalues); + check_petsc_error(EPSGetConverged(_eps, &num_computed_eigenvalues), + "EPSGetConverged"); if (i < num_computed_eigenvalues) { #ifdef PETSC_USE_COMPLEX PetscScalar l; - EPSGetEigenvalue(_eps, i, &l, nullptr); + check_petsc_error(EPSGetEigenvalue(_eps, i, &l, nullptr), + "EPSGetEigenvalue"); return l; #else PetscScalar lr, li; - EPSGetEigenvalue(_eps, i, &lr, &li); + check_petsc_error(EPSGetEigenvalue(_eps, i, &lr, &li), + "EPSGetEigenvalue"); return std::complex(lr, li); #endif } @@ -145,9 +164,11 @@ void SLEPcEigenSolver::get_eigenpair(PetscScalar& lr, PetscScalar& lc, Vec r, // Get number of computed eigenvectors/values PetscInt num_computed_eigenvalues; - EPSGetConverged(_eps, &num_computed_eigenvalues); + check_petsc_error(EPSGetConverged(_eps, &num_computed_eigenvalues), + "EPSGetConverged"); if (ii < num_computed_eigenvalues) - EPSGetEigenpair(_eps, ii, &lr, &lc, r, c); + check_petsc_error(EPSGetEigenpair(_eps, ii, &lr, &lc, r, c), + "EPSGetEigenpair"); else { throw std::runtime_error("Requested eigenpair (" + std::to_string(i) @@ -159,7 +180,7 @@ std::int64_t SLEPcEigenSolver::get_number_converged() const { PetscInt num_conv; assert(_eps); - EPSGetConverged(_eps, &num_conv); + check_petsc_error(EPSGetConverged(_eps, &num_conv), "EPSGetConverged"); return num_conv; } //----------------------------------------------------------------------------- @@ -167,8 +188,7 @@ void SLEPcEigenSolver::set_options_prefix(const std::string& options_prefix) { assert(_eps); PetscErrorCode ierr = EPSSetOptionsPrefix(_eps, options_prefix.c_str()); - if (ierr != 0) - petsc::error(ierr, __FILE__, "EPSSetOptionsPrefix"); + check_petsc_error(ierr, "EPSSetOptionsPrefix"); } //----------------------------------------------------------------------------- std::string SLEPcEigenSolver::get_options_prefix() const @@ -176,8 +196,7 @@ std::string SLEPcEigenSolver::get_options_prefix() const assert(_eps); const char* prefix = nullptr; PetscErrorCode ierr = EPSGetOptionsPrefix(_eps, &prefix); - if (ierr != 0) - petsc::error(ierr, __FILE__, "EPSGetOptionsPrefix"); + check_petsc_error(ierr, "EPSGetOptionsPrefix"); return std::string(prefix); } //----------------------------------------------------------------------------- @@ -185,15 +204,15 @@ void SLEPcEigenSolver::set_from_options() const { assert(_eps); PetscErrorCode ierr = EPSSetFromOptions(_eps); - if (ierr != 0) - petsc::error(ierr, __FILE__, "EPSSetFromOptions"); + check_petsc_error(ierr, "EPSSetFromOptions"); } //----------------------------------------------------------------------------- int SLEPcEigenSolver::get_iteration_number() const { assert(_eps); PetscInt num_iter; - EPSGetIterationNumber(_eps, &num_iter); + check_petsc_error(EPSGetIterationNumber(_eps, &num_iter), + "EPSGetIterationNumber"); return num_iter; } //----------------------------------------------------------------------------- @@ -203,7 +222,8 @@ MPI_Comm SLEPcEigenSolver::comm() const { assert(_eps); MPI_Comm mpi_comm = MPI_COMM_NULL; - PetscObjectGetComm((PetscObject)_eps, &mpi_comm); + check_petsc_error(PetscObjectGetComm((PetscObject)_eps, &mpi_comm), + "PetscObjectGetComm"); return mpi_comm; } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/nls/NewtonSolver.cpp b/cpp/dolfinx/nls/NewtonSolver.cpp index 1a2984abc9e..80f84751d4c 100644 --- a/cpp/dolfinx/nls/NewtonSolver.cpp +++ b/cpp/dolfinx/nls/NewtonSolver.cpp @@ -19,6 +19,12 @@ using namespace dolfinx; namespace { //----------------------------------------------------------------------------- +void check_petsc_error(PetscErrorCode ierr, const std::string& petsc_function) +{ + if (ierr != 0) + la::petsc::error(ierr, __FILE__, petsc_function); +} +//----------------------------------------------------------------------------- /// Convergence test /// @param solver The Newton solver @@ -29,7 +35,7 @@ std::pair converged(const nls::petsc::NewtonSolver& solver, const Vec r) { PetscReal residual = 0; - VecNorm(r, NORM_2, &residual); + check_petsc_error(VecNorm(r, NORM_2, &residual), "VecNorm"); // Relative residual const double relative_residual = residual / solver.residual0(); @@ -62,7 +68,7 @@ std::pair converged(const nls::petsc::NewtonSolver& solver, void update_solution(const nls::petsc::NewtonSolver& solver, const Vec dx, Vec x) { - VecAXPY(x, -solver.relaxation_parameter, dx); + check_petsc_error(VecAXPY(x, -solver.relaxation_parameter, dx), "VecAXPY"); } //----------------------------------------------------------------------------- } // namespace @@ -85,13 +91,13 @@ nls::petsc::NewtonSolver::NewtonSolver(MPI_Comm comm) nls::petsc::NewtonSolver::~NewtonSolver() { if (_b) - VecDestroy(&_b); + (void)VecDestroy(&_b); if (_dx) - VecDestroy(&_dx); + (void)VecDestroy(&_dx); if (_matJ) - MatDestroy(&_matJ); + (void)MatDestroy(&_matJ); if (_matP) - MatDestroy(&_matP); + (void)MatDestroy(&_matP); } //----------------------------------------------------------------------------- void nls::petsc::NewtonSolver::setF(std::function F, @@ -99,7 +105,8 @@ void nls::petsc::NewtonSolver::setF(std::function F, { _fnF = std::move(F); _b = b; - PetscObjectReference((PetscObject)_b); + check_petsc_error(PetscObjectReference((PetscObject)_b), + "PetscObjectReference"); } //----------------------------------------------------------------------------- void nls::petsc::NewtonSolver::setJ(std::function J, @@ -107,7 +114,8 @@ void nls::petsc::NewtonSolver::setJ(std::function J, { _fnJ = std::move(J); _matJ = Jmat; - PetscObjectReference((PetscObject)_matJ); + check_petsc_error(PetscObjectReference((PetscObject)_matJ), + "PetscObjectReference"); } //----------------------------------------------------------------------------- void nls::petsc::NewtonSolver::setP(std::function P, @@ -115,7 +123,8 @@ void nls::petsc::NewtonSolver::setP(std::function P, { _fnP = std::move(P); _matP = Pmat; - PetscObjectReference((PetscObject)_matP); + check_petsc_error(PetscObjectReference((PetscObject)_matP), + "PetscObjectReference"); } //----------------------------------------------------------------------------- const la::petsc::KrylovSolver& @@ -195,7 +204,7 @@ std::pair nls::petsc::NewtonSolver::solve(Vec x) _solver.set_operators(_matJ, _matJ); if (!_dx) - MatCreateVecs(_matJ, &_dx, nullptr); + check_petsc_error(MatCreateVecs(_matJ, &_dx, nullptr), "MatCreateVecs"); // Start iterations while (!newton_converged and _iteration < max_it) @@ -224,7 +233,7 @@ std::pair nls::petsc::NewtonSolver::solve(Vec x) if (_iteration == 1) { PetscReal _r = 0; - VecNorm(_dx, NORM_2, &_r); + check_petsc_error(VecNorm(_dx, NORM_2, &_r), "VecNorm"); _residual0 = _r; } diff --git a/python/dolfinx/fem/assemble.py b/python/dolfinx/fem/assemble.py index 4a1b89e017d..9fc2dca9ae5 100644 --- a/python/dolfinx/fem/assemble.py +++ b/python/dolfinx/fem/assemble.py @@ -26,6 +26,25 @@ from dolfinx.fem.function import FunctionSpace +def _check_bcs_for_different_trial_test_spaces(a: Form, bcs: Sequence[DirichletBC]) -> None: + """Reject DirichletBCs for unsupported Petrov-Galerkin assembly.""" + if not bcs or len(a.function_spaces) != 2: + return + + test_space, trial_space = a.function_spaces + if test_space is trial_space: + return + if test_space.mesh is not trial_space.mesh: + return + + for bc in bcs: + if test_space.contains(bc.function_space) or trial_space.contains(bc.function_space): + raise RuntimeError( + "Dirichlet boundary conditions for bilinear forms with different test " + "and trial spaces are not supported by this matrix assembly path." + ) + + def pack_constants( form: Form | Sequence[Form], ) -> npt.NDArray | Sequence[npt.NDArray]: @@ -288,6 +307,7 @@ def assemble_matrix( accumulated. """ bcs = [] if bcs is None else bcs + _check_bcs_for_different_trial_test_spaces(a, bcs) A: la.MatrixCSR = create_matrix(a, block_mode) _assemble_matrix_csr(A, a, bcs, diag, constants, coeffs) return A @@ -325,15 +345,17 @@ def _assemble_matrix_csr( The returned matrix is not finalised, i.e. ghost values are not accumulated. """ - bcs = [] if bcs is None else [bc._cpp_object for bc in bcs] + bcs = [] if bcs is None else bcs + _check_bcs_for_different_trial_test_spaces(a, bcs) + _bcs = [bc._cpp_object for bc in bcs] constants = pack_constants(a) if constants is None else constants # type: ignore[assignment] coeffs = pack_coefficients(a) if coeffs is None else coeffs # type: ignore[assignment] - _cpp.fem.assemble_matrix(A._cpp_object, a._cpp_object, constants, coeffs, bcs) + _cpp.fem.assemble_matrix(A._cpp_object, a._cpp_object, constants, coeffs, _bcs) # If matrix is a 'diagonal'block, set diagonal entry for constrained # dofs if a.function_spaces[0] is a.function_spaces[1]: - _cpp.fem.insert_diagonal(A._cpp_object, a.function_spaces[0], bcs, diag) + _cpp.fem.insert_diagonal(A._cpp_object, a.function_spaces[0], _bcs, diag) return A diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py index e5dee8c5cf1..21e4dc60c8f 100644 --- a/python/dolfinx/fem/petsc.py +++ b/python/dolfinx/fem/petsc.py @@ -49,7 +49,10 @@ from dolfinx.cpp.fem.petsc import discrete_gradient as _discrete_gradient from dolfinx.cpp.fem.petsc import interpolation_matrix as _interpolation_matrix from dolfinx.fem import IntegralType, pack_coefficients, pack_constants -from dolfinx.fem.assemble import _assemble_vector_array +from dolfinx.fem.assemble import ( + _assemble_vector_array, + _check_bcs_for_different_trial_test_spaces, +) from dolfinx.fem.assemble import apply_lifting as _apply_lifting from dolfinx.fem.bcs import DirichletBC from dolfinx.fem.bcs import bcs_by_block as _bcs_by_block @@ -395,8 +398,12 @@ def assemble_matrix( Returns: Matrix representing the bilinear form. """ # noqa: D301 + bcs = [] if bcs is None else bcs + if not isinstance(a, Sequence): + _check_bcs_for_different_trial_test_spaces(a, bcs) + A = create_matrix(a, kind) - assemble_matrix(A, a, bcs, diag, constants, coeffs) # type: ignore[arg-type] + assemble_matrix(A, a, bcs, diag, constants, coeffs, _check_bcs=False) # type: ignore[arg-type] return A @@ -412,6 +419,8 @@ def _( | Sequence[Sequence[dict[tuple[IntegralType, int], npt.NDArray]]] | None ) = None, + *, + _check_bcs: bool = True, ) -> PETSc.Mat: # type: ignore[name-defined] """Assemble bilinear form into a matrix. @@ -422,6 +431,10 @@ def _( The returned matrix is not finalised, i.e. ghost values are not accumulated. """ + bcs = [] if bcs is None else bcs + if _check_bcs and not isinstance(a, Sequence): + _check_bcs_for_different_trial_test_spaces(a, bcs) + if A.getType() == PETSc.Mat.Type.NEST: # type: ignore[attr-defined] if not isinstance(a, Sequence): raise ValueError("Must provide a sequence of forms when assembling a nest matrix") @@ -431,7 +444,7 @@ def _( for j, (a_block, const, coeff) in enumerate(zip(a_row, const_row, coeff_row)): if a_block is not None: Asub = A.getNestSubMatrix(i, j) - assemble_matrix(Asub, a_block, bcs, diag, const, coeff) # type: ignore[arg-type] + assemble_matrix(Asub, a_block, bcs, diag, const, coeff, _check_bcs=False) # type: ignore[arg-type] elif i == j: for bc in bcs: row_forms = [row_form for row_form in a_row if row_form is not None] @@ -462,7 +475,7 @@ def _( [(Vsub.dofmaps[0].index_map, Vsub.dofmaps[0].index_map_bs) for Vsub in V[1]] # type: ignore[union-attr] ) - _bcs = [bc._cpp_object for bc in bcs] if bcs is not None else [] + _bcs = [bc._cpp_object for bc in bcs] for i, a_row in enumerate(a): for j, a_sub in enumerate(a_row): if a_sub is not None: @@ -501,7 +514,7 @@ def _( else: # Non-blocked constants = pack_constants(a) if constants is None else constants # type: ignore[assignment] coeffs = pack_coefficients(a) if coeffs is None else coeffs # type: ignore[assignment] - _bcs = [bc._cpp_object for bc in bcs] if bcs is not None else [] + _bcs = [bc._cpp_object for bc in bcs] _cpp.fem.petsc.assemble_matrix(A, a._cpp_object, constants, coeffs, _bcs) if a.function_spaces[0] is a.function_spaces[1]: A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH) # type: ignore[attr-defined] diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 55d77b88a51..2f2afa2b03f 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -178,6 +178,26 @@ def nest_matrix_norm(A): return math.sqrt(norm) +def _different_trial_test_space_form_and_bc(): + msh = create_unit_square(MPI.COMM_WORLD, 4, 4) + U = functionspace(msh, ("Lagrange", 1)) + V = functionspace(msh, ("Lagrange", 2)) + u = ufl.TrialFunction(U) + v = ufl.TestFunction(V) + a = form(inner(u, v) * dx) + + facets = locate_entities_boundary(msh, msh.topology.dim - 1, lambda x: np.isclose(x[0], 0.0)) + dofs = locate_dofs_topological(V, msh.topology.dim - 1, facets) + bc = dirichletbc(default_scalar_type(0), dofs, V) + return a, bc + + +def test_assembly_rejects_bcs_for_different_trial_test_spaces(): + a, bc = _different_trial_test_space_form_and_bc() + with pytest.raises(RuntimeError, match="different test and trial spaces"): + assemble_matrix(a, bcs=[bc]) + + @pytest.mark.petsc4py def test_vector_single_space_as_block(): from dolfinx.fem.petsc import create_vector as petsc_create_vector @@ -193,6 +213,14 @@ def test_vector_single_space_as_block(): class TestPETScAssemblers: """Test PETSc-based assemblers for matrices and vectors.""" + def test_assembly_rejects_bcs_for_different_trial_test_spaces(self): + """Test that PETSc assembly rejects same-domain Petrov-Galerkin BCs.""" + from dolfinx.fem.petsc import assemble_matrix as petsc_assemble_matrix + + a, bc = _different_trial_test_space_form_and_bc() + with pytest.raises(RuntimeError, match="different test and trial spaces"): + petsc_assemble_matrix(a, bcs=[bc]) + @pytest.mark.parametrize("mode", [GhostMode.none, GhostMode.shared_facet]) def test_basic_assembly_petsc_matrixcsr(self, mode): """Test basic assembly of PETSc Mat and compare with MatrixCSR assembly.""" diff --git a/python/test/unit/io/test_adios2.py b/python/test/unit/io/test_adios2.py index aeeb526a0c6..70a49a6c597 100644 --- a/python/test/unit/io/test_adios2.py +++ b/python/test/unit/io/test_adios2.py @@ -251,9 +251,9 @@ def test_vtx_reuse_mesh(self, tempdir, dim, simplex, reuse): def test_dg_0_data(self, tempdir): """Test that we can mix DG-0 and other Lagrange functions.""" + adios2 = import_adios2() from dolfinx.io import VTXWriter - adios2 = import_adios2() mesh = generate_mesh(2, False) v = Function(functionspace(mesh, ("Lagrange", 2, (2,)))) filename = Path(tempdir, "v.bp")