diff --git a/cpp/dolfinx/la/MatrixCSR.h b/cpp/dolfinx/la/MatrixCSR.h index 3a5a8a26068..09296819d99 100644 --- a/cpp/dolfinx/la/MatrixCSR.h +++ b/cpp/dolfinx/la/MatrixCSR.h @@ -357,7 +357,7 @@ class MatrixCSR /// Block size /// @return block sizes for rows and columns - const std::array& block_size() const { return _bs; } + std::array block_size() const { return _bs; } private: // Maps for the distribution of the ows and columns diff --git a/cpp/dolfinx/la/Vector.h b/cpp/dolfinx/la/Vector.h index 97de3cd3d98..0201349af78 100644 --- a/cpp/dolfinx/la/Vector.h +++ b/cpp/dolfinx/la/Vector.h @@ -7,6 +7,7 @@ #pragma once #include "utils.h" +#include #include #include #include @@ -305,53 +306,60 @@ auto norm(const V& x, Norm type = Norm::l2) /// @param[in,out] basis The set of vectors to orthonormalise. The /// vectors must have identical parallel layouts. The vectors are /// modified in-place. -/// @param[in] tol The tolerance used to detect a linear dependency +/// @tparam V dolfinx::la::Vector template -void orthonormalize(std::span basis, double tol = 1.0e-10) +void orthonormalize(std::vector> basis) { using T = typename V::value_type; + using U = typename dolfinx::scalar_value_type_t; + // Loop over each vector in basis for (std::size_t i = 0; i < basis.size(); ++i) { // Orthogonalize vector i with respect to previously orthonormalized // vectors + V& bi = basis[i].get(); for (std::size_t j = 0; j < i; ++j) { + const V& bj = basis[j].get(); + // basis_i <- basis_i - dot_ij basis_j - T dot_ij = inner_product(basis[i], basis[j]); - std::transform(basis[j].array().begin(), basis[j].array().end(), - basis[i].array().begin(), basis[i].mutable_array().begin(), + auto dot_ij = inner_product(bi, bj); + std::transform(bj.array().begin(), bj.array().end(), bi.array().begin(), + bi.mutable_array().begin(), [dot_ij](auto xj, auto xi) { return xi - dot_ij * xj; }); } // Normalise basis function - double norm = la::norm(basis[i], la::Norm::l2); - if (norm < tol) + auto norm = la::norm(bi, la::Norm::l2); + if (norm * norm < std::numeric_limits::epsilon()) { throw std::runtime_error( "Linear dependency detected. Cannot orthogonalize."); } - std::transform(basis[i].array().begin(), basis[i].array().end(), - basis[i].mutable_array().begin(), + std::transform(bi.array().begin(), bi.array().end(), + bi.mutable_array().begin(), [norm](auto x) { return x / norm; }); } } -/// Test if basis is orthonormal -/// @param[in] basis The set of vectors to check -/// @param[in] tol The tolerance used to test for orthonormality -/// @return True is basis is orthonormal, otherwise false +/// @brief Test if basis is orthonormal. +/// @param[in] basis The set of vectors to check. +/// @return True is basis is orthonormal, otherwise false. template -bool is_orthonormal(std::span basis, double tol = 1.0e-10) +bool is_orthonormal(std::vector> basis) { using T = typename V::value_type; + using U = typename dolfinx::scalar_value_type_t; + + // auto tol = std::sqrt(T(std::numeric_limits::epsilon())); for (std::size_t i = 0; i < basis.size(); i++) { for (std::size_t j = i; j < basis.size(); j++) { - const double delta_ij = (i == j) ? 1.0 : 0.0; - T dot_ij = inner_product(basis[i], basis[j]); - if (std::abs(delta_ij - dot_ij) > tol) + T delta_ij = (i == j) ? T(1) : T(0); + auto dot_ij = inner_product(basis[i].get(), basis[j].get()); + if (std::norm(delta_ij - dot_ij) > std::numeric_limits::epsilon()) return false; } } diff --git a/cpp/dolfinx/la/petsc.cpp b/cpp/dolfinx/la/petsc.cpp index 1712e0aaaa5..d3d938877a1 100644 --- a/cpp/dolfinx/la/petsc.cpp +++ b/cpp/dolfinx/la/petsc.cpp @@ -75,15 +75,24 @@ Vec la::petsc::create_vector(MPI_Comm comm, std::array range, // Get local size assert(range[1] >= range[0]); - const std::int32_t local_size = range[1] - range[0]; + std::int32_t local_size = range[1] - range[0]; Vec x; - const std::vector _ghosts(ghosts.begin(), ghosts.end()); - ierr = VecCreateGhostBlock(comm, bs, bs * local_size, PETSC_DETERMINE, - _ghosts.size(), _ghosts.data(), &x); - CHECK_ERROR("VecCreateGhostBlock"); - assert(x); + std::vector _ghosts(ghosts.begin(), ghosts.end()); + if (bs == 1) + { + ierr = VecCreateGhost(comm, local_size, PETSC_DETERMINE, _ghosts.size(), + _ghosts.data(), &x); + CHECK_ERROR("VecCreateGhost"); + } + else + { + ierr = VecCreateGhostBlock(comm, bs, bs * local_size, PETSC_DETERMINE, + _ghosts.size(), _ghosts.data(), &x); + CHECK_ERROR("VecCreateGhostBlock"); + } + assert(x); return x; } //----------------------------------------------------------------------------- @@ -94,8 +103,23 @@ Vec la::petsc::create_vector_wrap(const common::IndexMap& map, int bs, const std::int64_t size_global = bs * map.size_global(); const std::vector ghosts(map.ghosts().begin(), map.ghosts().end()); Vec vec; - VecCreateGhostBlockWithArray(map.comm(), bs, size_local, size_global, - ghosts.size(), ghosts.data(), x.data(), &vec); + PetscErrorCode ierr; + if (bs == 1) + { + ierr + = VecCreateGhostWithArray(map.comm(), size_local, size_global, + ghosts.size(), ghosts.data(), x.data(), &vec); + CHECK_ERROR("VecCreateGhostWithArray"); + } + else + { + ierr = VecCreateGhostBlockWithArray(map.comm(), bs, size_local, size_global, + ghosts.size(), ghosts.data(), x.data(), + &vec); + CHECK_ERROR("VecCreateGhostBlockWithArray"); + } + + assert(vec); return vec; } //----------------------------------------------------------------------------- @@ -368,7 +392,6 @@ void petsc::options::clear() } //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -//----------------------------------------------------------------------------- petsc::Vector::Vector(const common::IndexMap& map, int bs) : _x(la::petsc::create_vector(map, bs)) { @@ -688,18 +711,6 @@ int petsc::KrylovSolver::solve(Vec x, const Vec b, bool transpose) const petsc::error(ierr, __FILE__, "KSPSolve"); } - // FIXME: Remove ghost updating? - // Update ghost values in solution vector - Vec xg; - VecGhostGetLocalForm(x, &xg); - const bool is_ghosted = xg ? true : false; - VecGhostRestoreLocalForm(x, &xg); - if (is_ghosted) - { - VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD); - VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD); - } - // Get the number of iterations PetscInt num_iterations = 0; ierr = KSPGetIterationNumber(_ksp, &num_iterations); diff --git a/python/demo/demo_axis/demo_axis.py b/python/demo/demo_axis/demo_axis.py index 3b082a6001b..6c24d2d1c08 100644 --- a/python/demo/demo_axis/demo_axis.py +++ b/python/demo/demo_axis/demo_axis.py @@ -18,15 +18,14 @@ from functools import partial import numpy as np -from mesh_sphere_axis import generate_mesh_sphere_axis -from scipy.special import jv, jvp - import ufl from basix.ufl import element, mixed_element -from dolfinx import fem, io, mesh, plot - +from dolfinx.fem.petsc import LinearProblem +from mesh_sphere_axis import generate_mesh_sphere_axis from mpi4py import MPI -from petsc4py import PETSc +from scipy.special import jv, jvp + +from dolfinx import default_scalar_type, fem, io, mesh, plot try: from dolfinx.io import VTXWriter @@ -52,7 +51,7 @@ # The time-harmonic Maxwell equation is complex-valued PDE. PETSc must # therefore have compiled with complex scalars. -if not np.issubdtype(PETSc.ScalarType, np.complexfloating): +if not np.issubdtype(default_scalar_type, np.complexfloating): print("Demo should only be executed with DOLFINx complex mode") exit(0) @@ -491,7 +490,7 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg): phi = np.pi / 4 # Initialize phase term -phase = fem.Constant(msh, PETSc.ScalarType(np.exp(1j * 0 * phi))) +phase = fem.Constant(msh, default_scalar_type(np.exp(1j * 0 * phi))) # - # We now solve the problem: @@ -518,7 +517,7 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg): + k0 ** 2 * ufl.inner(eps_pml * Es_m, v_m) * rho * dPml a, L = ufl.lhs(F), ufl.rhs(F) - problem = fem.petsc.LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) + problem = LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) Esh_m = problem.solve() # Scattered magnetic field diff --git a/python/demo/demo_biharmonic.py b/python/demo/demo_biharmonic.py index 5fe5f9d152e..86af5bbe649 100644 --- a/python/demo/demo_biharmonic.py +++ b/python/demo/demo_biharmonic.py @@ -115,6 +115,7 @@ import ufl from dolfinx import fem, io, mesh, plot +from dolfinx.fem.petsc import LinearProblem from dolfinx.mesh import CellType, GhostMode from ufl import (CellDiameter, FacetNormal, avg, div, dS, dx, grad, inner, jump, pi, sin) @@ -208,14 +209,13 @@ L = inner(f, v) * dx # - -# We create a {py:class}`LinearProblem ` +# We create a {py:class}`LinearProblem ` # object that brings together the variational problem, the Dirichlet # boundary condition, and which specifies the linear solver. In this # case we use a direct (LU) solver. The {py:func}`solve -# ` will compute a solution. +# ` will compute a solution. -problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", - "pc_type": "lu"}) +problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) uh = problem.solve() # The solution can be written to a {py:class}`XDMFFile diff --git a/python/demo/demo_elasticity.py b/python/demo/demo_elasticity.py index c25c73f502f..a0d5bfb18fd 100644 --- a/python/demo/demo_elasticity.py +++ b/python/demo/demo_elasticity.py @@ -33,11 +33,12 @@ from dolfinx.io import XDMFFile from dolfinx.mesh import (CellType, GhostMode, create_box, locate_entities_boundary) +from mpi4py import MPI from petsc4py import PETSc from ufl import dx, grad, inner -from dolfinx import default_real_type, la -from mpi4py import MPI +import dolfinx +from dolfinx import la dtype = PETSc.ScalarType # - @@ -58,33 +59,32 @@ def build_nullspace(V): # Create vectors that will span the nullspace bs = V.dofmap.index_map_bs length0 = V.dofmap.index_map.size_local - length1 = length0 + V.dofmap.index_map.num_ghosts - basis = [np.zeros(bs * length1, dtype=dtype) for i in range(6)] + basis = [la.vector(V.dofmap.index_map, bs=bs, dtype=dtype) for i in range(6)] + b = [b.array for b in basis] # Get dof indices for each subspace (x, y and z dofs) dofs = [V.sub(i).dofmap.list.flatten() for i in range(3)] # Set the three translational rigid body modes for i in range(3): - basis[i][dofs[i]] = 1.0 + b[i][dofs[i]] = 1.0 # Set the three rotational rigid body modes x = V.tabulate_dof_coordinates() dofs_block = V.dofmap.list.flatten() x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2] - basis[3][dofs[0]] = -x1 - basis[3][dofs[1]] = x0 - basis[4][dofs[0]] = x2 - basis[4][dofs[2]] = -x0 - basis[5][dofs[2]] = x1 - basis[5][dofs[1]] = -x2 - - # Create PETSc Vec objects (excluding ghosts) and normalise - basis_petsc = [PETSc.Vec().createWithArray(x[:bs * length0], bsize=3, comm=V.mesh.comm) for x in basis] - la.orthonormalize(basis_petsc) - assert la.is_orthonormal(basis_petsc, eps=1.0e3 * np.finfo(default_real_type).eps) - - # Create and return a PETSc nullspace + b[3][dofs[0]] = -x1 + b[3][dofs[1]] = x0 + b[4][dofs[0]] = x2 + b[4][dofs[2]] = -x0 + b[5][dofs[2]] = x1 + b[5][dofs[1]] = -x2 + + _basis = [x._cpp_object for x in basis] + dolfinx.cpp.la.orthonormalize(_basis) + assert dolfinx.cpp.la.is_orthonormal(_basis) + + basis_petsc = [PETSc.Vec().createWithArray(x[:bs * length0], bsize=3, comm=V.mesh.comm) for x in b] return PETSc.NullSpace().create(vectors=basis_petsc) @@ -167,6 +167,7 @@ def σ(v): ns = build_nullspace(V) A.setNearNullSpace(ns) +A.setOption(PETSc.Mat.Option.SPD, True) # Set PETSc solver options, create a PETSc Krylov solver, and attach the # matrix `A` to the solver: @@ -175,7 +176,7 @@ def σ(v): # Set solver options opts = PETSc.Options() opts["ksp_type"] = "cg" -opts["ksp_rtol"] = 1.0e-10 +opts["ksp_rtol"] = 1.0e-8 opts["pc_type"] = "gamg" # Use Chebyshev smoothing for multigrid @@ -183,8 +184,7 @@ def σ(v): opts["mg_levels_pc_type"] = "jacobi" # Improve estimate of eigenvalues for Chebyshev smoothing -opts["mg_levels_esteig_ksp_type"] = "cg" -opts["mg_levels_ksp_chebyshev_esteig_steps"] = 20 +opts["mg_levels_ksp_chebyshev_esteig_steps"] = 10 # Create PETSc Krylov solver and turn convergence monitoring on solver = PETSc.KSP().create(msh.comm) diff --git a/python/demo/demo_lagrange_variants.py b/python/demo/demo_lagrange_variants.py index 99f29902ef9..b14bb06acb3 100644 --- a/python/demo/demo_lagrange_variants.py +++ b/python/demo/demo_lagrange_variants.py @@ -17,21 +17,17 @@ # # We begin this demo by importing the required modules. -import matplotlib.pylab as plt -import numpy as np - # + import basix import basix.ufl +import matplotlib.pylab as plt +import numpy as np import ufl -from dolfinx import default_scalar_type, fem, mesh -from ufl import ds, dx, grad, inner - +from dolfinx.fem.petsc import LinearProblem from mpi4py import MPI +from ufl import ds, dx, grad, inner -if np.issubdtype(default_scalar_type, np.complexfloating): - print("Demo should only be executed with DOLFINx real mode") - exit(0) +from dolfinx import default_scalar_type, fem, mesh # - # Note that Basix and the Basix UFL wrapper are imported directly. Basix @@ -132,7 +128,7 @@ a = inner(grad(u), grad(v)) * dx L = inner(f, v) * dx + inner(g, v) * ds -problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) +problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) uh = problem.solve() # - diff --git a/python/demo/demo_mixed-poisson.py b/python/demo/demo_mixed-poisson.py index 4ff02275b00..e80c1ec26ec 100644 --- a/python/demo/demo_mixed-poisson.py +++ b/python/demo/demo_mixed-poisson.py @@ -85,14 +85,14 @@ # + import numpy as np - from basix.ufl import element, mixed_element -from dolfinx import fem, io, mesh +from dolfinx.fem.petsc import LinearProblem +from mpi4py import MPI +from petsc4py import PETSc from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner) -from mpi4py import MPI -from petsc4py import PETSc +from dolfinx import fem, io, mesh domain = mesh.create_unit_square( MPI.COMM_WORLD, @@ -159,8 +159,8 @@ def f2(x): bcs = [bc_top, bc_bottom] -problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={ - "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}) +problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps"}) try: w_h = problem.solve() except PETSc.Error as e: diff --git a/python/demo/demo_navier_stokes.py b/python/demo/demo_navier_stokes.py index 09b35f0bf41..1ee284132a1 100644 --- a/python/demo/demo_navier_stokes.py +++ b/python/demo/demo_navier_stokes.py @@ -158,7 +158,8 @@ # + import numpy as np -from dolfinx import fem, io, mesh, default_real_type +from dolfinx import default_real_type, fem, io, mesh +from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block from ufl import (CellDiameter, FacetNormal, TestFunction, TrialFunction, avg, conditional, div, dot, dS, ds, dx, grad, gt, inner, outer) @@ -282,9 +283,9 @@ def jump(phi, n): bcs = [bc_u] # Assemble Stokes problem -A = fem.petsc.assemble_matrix_block(a, bcs=bcs) +A = assemble_matrix_block(a, bcs=bcs) A.assemble() -b = fem.petsc.assemble_vector_block(L, a, bcs=bcs) +b = assemble_vector_block(L, a, bcs=bcs) # Create and configure solver ksp = PETSc.KSP().create(msh.comm) diff --git a/python/demo/demo_pml/demo_pml.py b/python/demo/demo_pml/demo_pml.py index cc81085f6b2..bd1bc87974f 100644 --- a/python/demo/demo_pml/demo_pml.py +++ b/python/demo/demo_pml/demo_pml.py @@ -33,23 +33,22 @@ from functools import partial from typing import Tuple, Union -from efficiencies_pml_demo import calculate_analytical_efficiencies -from mesh_wire_pml import generate_mesh_wire - import ufl from basix.ufl import element -from dolfinx import fem, mesh, plot +from dolfinx.fem.petsc import LinearProblem from dolfinx.io import VTXWriter, gmshio - +from efficiencies_pml_demo import calculate_analytical_efficiencies +from mesh_wire_pml import generate_mesh_wire from mpi4py import MPI -from petsc4py import PETSc + +from dolfinx import default_scalar_type, fem, mesh, plot # - # Since we want to solve time-harmonic Maxwell's equation, we require # that the demo is executed with DOLFINx (PETSc) complex mode. -if not np.issubdtype(PETSc.ScalarType, np.complexfloating): +if not np.issubdtype(default_scalar_type, np.complexfloating): print("Demo should only be executed with DOLFINx complex mode") exit(0) @@ -412,7 +411,7 @@ def create_eps_mu(pml: ufl.tensors.ListTensor, a, L = ufl.lhs(F), ufl.rhs(F) -problem = fem.petsc.LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) +problem = LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) Esh = problem.solve() # - diff --git a/python/demo/demo_poisson.py b/python/demo/demo_poisson.py index d6fb650864e..05b6a559f57 100644 --- a/python/demo/demo_poisson.py +++ b/python/demo/demo_poisson.py @@ -67,12 +67,14 @@ # + import numpy as np + import ufl -from mpi4py import MPI -from petsc4py.PETSc import ScalarType +from dolfinx import fem, io, mesh, plot +from dolfinx.fem.petsc import LinearProblem from ufl import ds, dx, grad, inner -from dolfinx import fem, io, mesh, plot +from mpi4py import MPI +from petsc4py.PETSc import ScalarType # - @@ -129,14 +131,14 @@ L = inner(f, v) * dx + inner(g, v) * ds # - -# A {py:class}`LinearProblem ` object is +# A {py:class}`LinearProblem ` object is # created that brings together the variational problem, the Dirichlet # boundary condition, and which specifies the linear solver. In this # case an LU solver us sued. The {py:func}`solve -# ` computes the solution. +# ` computes the solution. # + -problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) +problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) uh = problem.solve() # - diff --git a/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py b/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py index c8731f03ac1..42015813e9a 100644 --- a/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py +++ b/python/demo/demo_scattering_boundary_conditions/demo_scattering_boundary_conditions.py @@ -47,15 +47,14 @@ have_pyvista = False from typing import Tuple -from analytical_efficiencies_wire import calculate_analytical_efficiencies -from mesh_wire import generate_mesh_wire - import ufl +from analytical_efficiencies_wire import calculate_analytical_efficiencies from basix.ufl import element -from dolfinx import fem, io, plot - +from dolfinx.fem.petsc import LinearProblem +from mesh_wire import generate_mesh_wire from mpi4py import MPI -from petsc4py import PETSc + +from dolfinx import default_scalar_type, fem, io, plot # - @@ -63,7 +62,7 @@ # solve a complex-valued PDE, and therefore need to use PETSc compiled # with complex numbers. -if not np.issubdtype(PETSc.ScalarType, np.complexfloating): +if not np.issubdtype(default_scalar_type, np.complexfloating): print("Demo should only be executed with DOLFINx complex mode") exit(0) @@ -407,7 +406,7 @@ def curl_2d(f: fem.Function): # `Esh`: a, L = ufl.lhs(F), ufl.rhs(F) -problem = fem.petsc.LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) +problem = LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) Esh = problem.solve() # We save the solution as an [ADIOS2 diff --git a/python/demo/demo_stokes.py b/python/demo/demo_stokes.py index 7eba180cc47..5718beda308 100644 --- a/python/demo/demo_stokes.py +++ b/python/demo/demo_stokes.py @@ -91,6 +91,7 @@ from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc, extract_function_spaces, form, locate_dofs_topological) +from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block from dolfinx.io import XDMFFile from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary from ufl import div, dx, grad, inner @@ -196,6 +197,13 @@ def nested_iterative_solver(): P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]]) P.assemble() + A00 = A.getNestSubMatrix(0, 0) + A00.setOption(PETSc.Mat.Option.SPD, True) + + P00, P11 = P.getNestSubMatrix(0, 0), P.getNestSubMatrix(1, 1) + P00.setOption(PETSc.Mat.Option.SPD, True) + P11.setOption(PETSc.Mat.Option.SPD, True) + # Assemble right-hand side vector b = fem.petsc.assemble_vector_nest(L) @@ -279,8 +287,8 @@ def nested_iterative_solver(): norm_u = u.x.norm() norm_p = p.x.norm() if MPI.COMM_WORLD.rank == 0: - print(f"(A) Norm of velocity coefficient vector (blocked, iterative): {norm_u}") - print(f"(A) Norm of pressure coefficient vector (blocked, iterative): {norm_p}") + print(f"(A) Norm of velocity coefficient vector (nested, iterative): {norm_u}") + print(f"(A) Norm of pressure coefficient vector (nested, iterative): {norm_p}") return norm_u, norm_p @@ -297,11 +305,11 @@ def block_operators(): # Assembler matrix operator, preconditioner and RHS vector into # single objects but preserving block structure - A = fem.petsc.assemble_matrix_block(a, bcs=bcs) + A = assemble_matrix_block(a, bcs=bcs) A.assemble() - P = fem.petsc.assemble_matrix_block(a_p, bcs=bcs) + P = assemble_matrix_block(a_p, bcs=bcs) P.assemble() - b = fem.petsc.assemble_vector_block(L, a, bcs=bcs) + b = assemble_vector_block(L, a, bcs=bcs) # Set the nullspace for pressure (since pressure is determined only # up to a constant) @@ -429,8 +437,8 @@ def block_direct_solver(): # Compute the $L^2$ norms of the u and p vectors norm_u, norm_p = u.x.norm(), p.x.norm() if MPI.COMM_WORLD.rank == 0: - print(f"(C) Norm of velocity coefficient vector (blocked, iterative): {norm_u}") - print(f"(C) Norm of pressure coefficient vector (blocked, iterative): {norm_p}") + print(f"(C) Norm of velocity coefficient vector (blocked, direct): {norm_u}") + print(f"(C) Norm of pressure coefficient vector (blocked, direct): {norm_p}") return norm_u, norm_p @@ -515,8 +523,8 @@ def mixed_direct(): # Compute norms norm_u, norm_p = u.x.norm(), p.x.norm() if MPI.COMM_WORLD.rank == 0: - print(f"(D) Norm of velocity coefficient vector (blocked, iterative): {norm_u}") - print(f"(D) Norm of pressure coefficient vector (blocked, iterative): {norm_p}") + print(f"(D) Norm of velocity coefficient vector (monolithic, direct): {norm_u}") + print(f"(D) Norm of pressure coefficient vector (monolithic, direct): {norm_p}") return norm_u, norm_u diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py index c5d4343d363..ac9e484c536 100644 --- a/python/demo/demo_tnt-elements.py +++ b/python/demo/demo_tnt-elements.py @@ -21,20 +21,19 @@ # We begin this demo by importing the required modules. # + +import basix +import basix.ufl import matplotlib import matplotlib.pylab as plt import numpy as np - -import basix -import basix.ufl -from dolfinx import fem, mesh +from dolfinx.fem.petsc import LinearProblem +from mpi4py import MPI from ufl import (SpatialCoordinate, TestFunction, TrialFunction, cos, div, dx, grad, inner, sin) -from mpi4py import MPI +from dolfinx import fem, mesh matplotlib.use('agg') - # - # ## Defining a degree 1 TNT element @@ -209,7 +208,7 @@ def poisson_error(V): bc = fem.dirichletbc(u_bc, bdofs) # Solve - problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_rtol": 1e-12}) + problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_rtol": 1e-12}) uh = problem.solve() M = (u_exact - uh)**2 * dx diff --git a/python/demo/demo_types.py b/python/demo/demo_types.py index 688fb02f2fa..832976b39fe 100644 --- a/python/demo/demo_types.py +++ b/python/demo/demo_types.py @@ -208,8 +208,8 @@ def σ(v): display_scalar(uh, "poisson", np.imag) -# uh = elasticity(dtype=np.float32) -# uh = elasticity(dtype=np.float64) -# uh = elasticity(dtype=np.complex64) -# uh = elasticity(dtype=np.complex128) -# display_vector(uh, "elasticity", np.real) +uh = elasticity(dtype=np.float32) +uh = elasticity(dtype=np.float64) +uh = elasticity(dtype=np.complex64) +uh = elasticity(dtype=np.complex128) +display_vector(uh, "elasticity", np.real) diff --git a/python/demo/demo_waveguide/demo_half_loaded_waveguide.py b/python/demo/demo_waveguide/demo_half_loaded_waveguide.py index efef4591d80..35b02374f3c 100644 --- a/python/demo/demo_waveguide/demo_half_loaded_waveguide.py +++ b/python/demo/demo_waveguide/demo_half_loaded_waveguide.py @@ -45,12 +45,12 @@ import ufl from analytical_modes import verify_mode from basix.ufl import element, mixed_element +from dolfinx.fem.petsc import assemble_matrix from dolfinx.mesh import (CellType, create_rectangle, exterior_facet_indices, locate_entities) from mpi4py import MPI -from petsc4py.PETSc import ScalarType -from dolfinx import fem, io, plot +from dolfinx import default_scalar_type, fem, io, plot try: import pyvista @@ -104,8 +104,8 @@ def Omega_v(x): cells_v = locate_entities(msh, msh.topology.dim, Omega_v) cells_d = locate_entities(msh, msh.topology.dim, Omega_d) -eps.x.array[cells_d] = np.full_like(cells_d, eps_d, dtype=ScalarType) -eps.x.array[cells_v] = np.full_like(cells_v, eps_v, dtype=ScalarType) +eps.x.array[cells_d] = np.full_like(cells_d, eps_d, dtype=default_scalar_type) +eps.x.array[cells_v] = np.full_like(cells_v, eps_v, dtype=default_scalar_type) # - # In order to find the weak form of our problem, the starting point are @@ -230,9 +230,9 @@ def Omega_v(x): # Now we can solve the problem with SLEPc. First of all, we need to # assemble our $A$ and $B$ matrices with PETSc in this way: -A = fem.petsc.assemble_matrix(a, bcs=[bc]) +A = assemble_matrix(a, bcs=[bc]) A.assemble() -B = fem.petsc.assemble_matrix(b, bcs=[bc]) +B = assemble_matrix(b, bcs=[bc]) B.assemble() # Now, we need to create the eigenvalue problem in SLEPc. Our problem is diff --git a/python/dolfinx/__init__.py b/python/dolfinx/__init__.py index abab765c5d4..48eb8cb912d 100644 --- a/python/dolfinx/__init__.py +++ b/python/dolfinx/__init__.py @@ -26,10 +26,14 @@ import sys -from petsc4py import PETSc as _PETSc - -default_scalar_type = _PETSc.ScalarType -default_real_type = _PETSc.RealType +try: + from petsc4py import PETSc as _PETSc + default_scalar_type = _PETSc.ScalarType + default_real_type = _PETSc.RealType +except ImportError: + import numpy as _np + default_scalar_type = _np.float64 + default_real_type = _np.float64 from dolfinx import common from dolfinx import cpp as _cpp diff --git a/python/dolfinx/common.py b/python/dolfinx/common.py index ac7db4bfbdd..2572d239640 100644 --- a/python/dolfinx/common.py +++ b/python/dolfinx/common.py @@ -3,7 +3,7 @@ # This file is part of DOLFINx (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -"""General tools for timing and configuration""" +"""General tools for timing and configuration.""" import functools import typing diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py index 532019091dc..e92e7ba81a2 100644 --- a/python/dolfinx/fem/__init__.py +++ b/python/dolfinx/fem/__init__.py @@ -9,7 +9,6 @@ from dolfinx.cpp.fem import (IntegralType, create_nonmatching_meshes_interpolation_data) from dolfinx.cpp.fem import create_sparsity_pattern as _create_sparsity_pattern -from dolfinx.fem import petsc from dolfinx.fem.assemble import (apply_lifting, assemble_matrix, assemble_scalar, assemble_vector, create_matrix, set_bc, create_vector) @@ -43,4 +42,4 @@ def create_sparsity_pattern(a: Form): "DirichletBC", "dirichletbc", "bcs_by_block", "DofMap", "Form", "form", "IntegralType", "create_vector", "locate_dofs_geometrical", "locate_dofs_topological", - "extract_function_spaces", "petsc", "create_nonmatching_meshes_interpolation_data"] + "extract_function_spaces", "create_nonmatching_meshes_interpolation_data"] diff --git a/python/dolfinx/fem/assemble.py b/python/dolfinx/fem/assemble.py index fa7b937746e..6793c552fb5 100644 --- a/python/dolfinx/fem/assemble.py +++ b/python/dolfinx/fem/assemble.py @@ -12,15 +12,15 @@ import typing import numpy as np - -import dolfinx -from dolfinx import cpp as _cpp -from dolfinx import la from dolfinx.cpp.fem import pack_coefficients as _pack_coefficients from dolfinx.cpp.fem import pack_constants as _pack_constants from dolfinx.fem.bcs import DirichletBC from dolfinx.fem.forms import Form +import dolfinx +from dolfinx import cpp as _cpp +from dolfinx import la + def pack_constants(form: typing.Union[Form, typing.Sequence[Form]]) -> typing.Union[np.ndarray, typing.Sequence[np.ndarray]]: diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 6dc0eb25050..b33d6641efc 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -394,7 +394,8 @@ def copy(self) -> Function: def vector(self): """PETSc vector holding the degrees-of-freedom.""" if self._petsc_x is None: - self._petsc_x = _cpp.la.petsc.create_vector_wrap(self._cpp_object.x) + from dolfinx.la import create_petsc_vector_wrap + self._petsc_x = create_petsc_vector_wrap(self.x) return self._petsc_x @property diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py index 24d1e7e3ca4..f96ea3496e8 100644 --- a/python/dolfinx/fem/petsc.py +++ b/python/dolfinx/fem/petsc.py @@ -17,9 +17,9 @@ import os import typing +import petsc4py +import petsc4py.lib import ufl -from dolfinx import cpp as _cpp -from dolfinx import la from dolfinx.cpp.fem import pack_coefficients as _pack_coefficients from dolfinx.cpp.fem import pack_constants as _pack_constants from dolfinx.fem import assemble @@ -29,11 +29,19 @@ from dolfinx.fem.forms import extract_function_spaces as _extract_spaces from dolfinx.fem.forms import form as _create_form from dolfinx.fem.function import Function as _Function - -import petsc4py -import petsc4py.lib +from dolfinx.la import create_petsc_vector from petsc4py import PETSc +from dolfinx import cpp as _cpp +from dolfinx import la + +__all__ = ["create_vector", "create_vector_block", "create_vector_nest", + "create_matrix", "create_matrix_block", "create_matrix_nest", + "assemble_vector", "assemble_vector_nest", "assemble_vector_block", + "assemble_matrix", "assemble_matrix_nest", "assemble_matrix_block", + "apply_lifting", "apply_lifting_nest", "set_bc", "set_bc_nest", + "LinearProblem", "NonlinearProblem"] + def _extract_function_spaces(a: typing.List[typing.List[Form]]): """From a rectangular array of bilinear forms, extract the function @@ -75,11 +83,11 @@ def create_vector(L: Form) -> PETSc.Vec: L: A linear form. Returns: - A PETSc vector with a layout that is compatible with `L`. + A PETSc vector with a layout that is compatible with ``L``. """ dofmap = L.function_spaces[0].dofmap - return la.create_petsc_vector(dofmap.index_map, dofmap.index_map_bs) + return create_petsc_vector(dofmap.index_map, dofmap.index_map_bs) def create_vector_block(L: typing.List[Form]) -> PETSc.Vec: @@ -89,7 +97,7 @@ def create_vector_block(L: typing.List[Form]) -> PETSc.Vec: L: List of linear forms. Returns: - A PETSc vector with a layout that is compatible with `L`. + A PETSc vector with a layout that is compatible with ``L``. """ maps = [(form.function_spaces[0].dofmap.index_map, @@ -98,14 +106,14 @@ def create_vector_block(L: typing.List[Form]) -> PETSc.Vec: def create_vector_nest(L: typing.List[Form]) -> PETSc.Vec: - """Create a PETSc netsted vector (``VecNest``) that is compaible with a list of linear forms. + """Create a PETSc nested vector (``VecNest``) that is compatible with a list of linear forms. Args: L: List of linear forms. Returns: A PETSc nested vector (``VecNest``) with a layout that is - compatible with `L`. + compatible with ``L``. """ maps = [(form.function_spaces[0].dofmap.index_map, @@ -116,14 +124,14 @@ def create_vector_nest(L: typing.List[Form]) -> PETSc.Vec: # -- Matrix instantiation ---------------------------------------------------- def create_matrix(a: Form, mat_type=None) -> PETSc.Mat: - """Create a PETSc matrix that is compaible with a bilinear form. + """Create a PETSc matrix that is compatible with a bilinear form. Args: a: A bilinear form. mat_type: The PETSc matrix type (``MatType``). Returns: - A PETSc matrix with a layout that is compatible with `a`. + A PETSc matrix with a layout that is compatible with ``a``. """ if mat_type is None: @@ -133,13 +141,14 @@ def create_matrix(a: Form, mat_type=None) -> PETSc.Mat: def create_matrix_block(a: typing.List[typing.List[Form]]) -> PETSc.Mat: - """Create a PETSc matrix that is compaible with a rectangular array of bilinear forms. + """Create a PETSc matrix that is compatible with a rectangular array of bilinear forms. Args: - a: A rectangular array of bilinear forms. + a: Rectangular array of bilinear forms. Returns: - A PETSc matrix with a blocked layout that is compatible with `a`. + A PETSc matrix with a blocked layout that is compatible with + ``a``. """ _a = [[None if form is None else form._cpp_object for form in arow] for arow in a] @@ -147,13 +156,13 @@ def create_matrix_block(a: typing.List[typing.List[Form]]) -> PETSc.Mat: def create_matrix_nest(a: typing.List[typing.List[Form]]) -> PETSc.Mat: - """Create a PETSc matrix (``MatNest``) that is compaible with a rectangular array of bilinear forms. + """Create a PETSc matrix (``MatNest``) that is compatible with a rectangular array of bilinear forms. Args: - a: A rectangular array of bilinear forms. + a: Rectangular array of bilinear forms. Returns: - A PETSc matrix ('MatNest``) that is compatible with `a`. + A PETSc matrix (``MatNest``) that is compatible with ``a``. """ _a = [[None if form is None else form._cpp_object for form in arow] for arow in a] @@ -164,11 +173,6 @@ def create_matrix_nest(a: typing.List[typing.List[Form]]) -> PETSc.Mat: @functools.singledispatch def assemble_vector(L: typing.Any, constants=None, coeffs=None) -> PETSc.Vec: - return _assemble_vector_form(L, constants, coeffs) - - -@assemble_vector.register(Form) -def _assemble_vector_form(L: Form, constants=None, coeffs=None) -> PETSc.Vec: """Assemble linear form into a new PETSc vector. Note: @@ -182,8 +186,8 @@ def _assemble_vector_form(L: Form, constants=None, coeffs=None) -> PETSc.Vec: An assembled vector. """ - b = la.create_petsc_vector(L.function_spaces[0].dofmap.index_map, - L.function_spaces[0].dofmap.index_map_bs) + b = create_petsc_vector(L.function_spaces[0].dofmap.index_map, + L.function_spaces[0].dofmap.index_map_bs) with b.localForm() as b_local: assemble._assemble_vector_array(b_local.array_w, L, constants, coeffs) return b @@ -200,7 +204,7 @@ def _assemble_vector_vec(b: PETSc.Vec, L: Form, constants=None, coeffs=None) -> Args: b: Vector to assemble the contribution of the linear form into. - L: A linear form to assemble into `b`. + L: A linear form to assemble into ``b``. Returns: An assembled vector. @@ -213,12 +217,7 @@ def _assemble_vector_vec(b: PETSc.Vec, L: Form, constants=None, coeffs=None) -> @functools.singledispatch def assemble_vector_nest(L: typing.Any, constants=None, coeffs=None) -> PETSc.Vec: - return _assemble_vector_nest_forms(L, constants, coeffs) - - -@assemble_vector_nest.register(list) -def _assemble_vector_nest_forms(L: typing.List[Form], constants=None, coeffs=None) -> PETSc.Vec: - """Assemble linear forms into a new nested PETSc (VecNest) vector. + """Assemble linear forms into a new nested PETSc (``VecNest``) vector. The returned vector is not finalised, i.e. ghost values are not accumulated on the owning processes. @@ -234,7 +233,7 @@ def _assemble_vector_nest_forms(L: typing.List[Form], constants=None, coeffs=Non @assemble_vector_nest.register def _assemble_vector_nest_vec(b: PETSc.Vec, L: typing.List[Form], constants=None, coeffs=None) -> PETSc.Vec: - """Assemble linear forms into a nested PETSc (VecNest) vector. The + """Assemble linear forms into a nested PETSc (``VecNest``) vector. The vector is not zeroed before assembly and it is not finalised, i.e. ghost values are not accumulated on the owning processes. @@ -327,27 +326,34 @@ def _assemble_vector_block_vec(b: PETSc.Vec, # -- Matrix assembly --------------------------------------------------------- @functools.singledispatch def assemble_matrix(a: typing.Any, bcs: typing.List[DirichletBC] = [], - diagonal: float = 1.0, - constants=None, coeffs=None) -> PETSc.Mat: - return _assemble_matrix_form(a, bcs, diagonal, constants, coeffs) - - -@assemble_matrix.register(Form) -def _assemble_matrix_form(a: Form, bcs: typing.List[DirichletBC] = [], - diagonal: float = 1.0, - constants=None, coeffs=None) -> PETSc.Mat: + diagonal: float = 1.0, constants=None, coeffs=None): """Assemble bilinear form into a matrix. The returned matrix is not finalised, i.e. ghost values are not accumulated. + Note: + The returned matrix is not 'assembled', i.e. ghost contributions + have not been communicated. + + Args: + a: Bilinear form to assembled into a matrix. + bc: Dirichlet boundary conditions applied to the system. + diagonal: Value to set on matrix diagonal for Dirichlet boundary + condition constrained degrees-of-freedom. + constants: Constants appearing the in the form. + coeffs: Coefficients appearing the in the form. + + Returns: + Matrix representing the bilinear form. + """ A = _cpp.fem.petsc.create_matrix(a._cpp_object) - _assemble_matrix_mat(A, a, bcs, diagonal, constants, coeffs) + assemble_matrix_mat(A, a, bcs, diagonal, constants, coeffs) return A @assemble_matrix.register -def _assemble_matrix_mat(A: PETSc.Mat, a: Form, bcs: typing.List[DirichletBC] = [], - diagonal: float = 1.0, constants=None, coeffs=None) -> PETSc.Mat: +def assemble_matrix_mat(A: PETSc.Mat, a: Form, bcs: typing.List[DirichletBC] = [], + diagonal: float = 1.0, constants=None, coeffs=None) -> PETSc.Mat: """Assemble bilinear form into a matrix. The returned matrix is not finalised, i.e. ghost values are not accumulated. @@ -367,9 +373,23 @@ def _assemble_matrix_mat(A: PETSc.Mat, a: Form, bcs: typing.List[DirichletBC] = @functools.singledispatch def assemble_matrix_nest(a: typing.List[typing.List[Form]], bcs: typing.List[DirichletBC] = [], mat_types=[], - diagonal: float = 1.0, - constants=None, coeffs=None) -> PETSc.Mat: - """Assemble bilinear forms into matrix""" + diagonal: float = 1.0, constants=None, coeffs=None) -> PETSc.Mat: + """Create a nested matrix and assembled bilinear forms into the matrix. + + Args: + a: Rectangular (list-of-lists) array for bilinear forms. + bcs: Dirichlet boundary conditions. + mat_types: PETSc matrix type for each matrix block. + diagonal: Value to set on matrix diagonal for Dirichlet boundary + condition constrained degrees-of-freedom. + constants: Constants appearing the in the form. + coeffs: Coefficients appearing the in the form. + + Returns: + PETSc matrix (``MatNest``) representing the block of bilinear + forms. + + """ _a = [[None if form is None else form._cpp_object for form in arow] for arow in a] A = _cpp.fem.petsc.create_matrix_nest(_a, mat_types) _assemble_matrix_nest_mat(A, a, bcs, diagonal, constants, coeffs) @@ -380,7 +400,24 @@ def assemble_matrix_nest(a: typing.List[typing.List[Form]], def _assemble_matrix_nest_mat(A: PETSc.Mat, a: typing.List[typing.List[Form]], bcs: typing.List[DirichletBC] = [], diagonal: float = 1.0, constants=None, coeffs=None) -> PETSc.Mat: - """Assemble bilinear forms into matrix""" + """Assemble bilinear forms into a nested matrix + + Args: + A: PETSc ``MatNest`` matrix. Matrix must have been correctly + initialized for the bilinear forms. + a: Rectangular (list-of-lists) array for bilinear forms. + bcs: Dirichlet boundary conditions. + mat_types: PETSc matrix type for each matrix block. + diagonal: Value to set on matrix diagonal for Dirichlet boundary + condition constrained degrees-of-freedom. + constants: Constants appearing the in the form. + coeffs: Coefficients appearing the in the form. + + Returns: + PETSc matrix (``MatNest``) representing the block of bilinear + forms. + + """ constants = [[form and _pack_constants(form._cpp_object) for form in forms] for forms in a] if constants is None else constants coeffs = [[{} if form is None else _pack_coefficients( @@ -389,7 +426,7 @@ def _assemble_matrix_nest_mat(A: PETSc.Mat, a: typing.List[typing.List[Form]], 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_mat(Asub, a_block, bcs, diagonal, const, coeff) + assemble_matrix_mat(Asub, a_block, bcs, diagonal, const, coeff) elif i == j: for bc in bcs: row_forms = [row_form for row_form in a_row if row_form is not None] @@ -407,7 +444,7 @@ def assemble_matrix_block(a: typing.List[typing.List[Form]], bcs: typing.List[DirichletBC] = [], diagonal: float = 1.0, constants=None, coeffs=None) -> PETSc.Mat: # type: ignore - """Assemble bilinear forms into matrix""" + """Assemble bilinear forms into a blocked matrix.""" _a = [[None if form is None else form._cpp_object for form in arow] for arow in a] A = _cpp.fem.petsc.create_matrix_block(_a) return _assemble_matrix_block_mat(A, a, bcs, diagonal, constants, coeffs) @@ -417,7 +454,7 @@ def assemble_matrix_block(a: typing.List[typing.List[Form]], def _assemble_matrix_block_mat(A: PETSc.Mat, a: typing.List[typing.List[Form]], bcs: typing.List[DirichletBC] = [], diagonal: float = 1.0, constants=None, coeffs=None) -> PETSc.Mat: - """Assemble bilinear forms into matrix""" + """Assemble bilinear forms into a blocked matrix.""" constants = [[form and _pack_constants(form._cpp_object) for form in forms] for forms in a] if constants is None else constants @@ -520,8 +557,10 @@ def __init__(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[DirichletBC] = [], """Initialize solver for a linear variational problem. Args: - a: A bilinear UFL form, the left hand side of the variational problem. - L: A linear UFL form, the right hand side of the variational problem. + a: A bilinear UFL form, the left hand side of the + variational problem. + L: A linear UFL form, the right hand side of the variational + problem. bcs: A list of Dirichlet boundary conditions. u: The solution function. It will be created if not provided. petsc_options: Options that are passed to the linear @@ -538,14 +577,12 @@ def __init__(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[DirichletBC] = [], Example:: - problem = LinearProblem(a, L, [bc0, bc1], - petsc_options={"ksp_type": "preonly", - "pc_type": "lu", - "pc_factor_mat_solver_type": "mumps"}) + problem = LinearProblem(a, L, [bc0, bc1], petsc_options={"ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps"}) """ self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options) self._A = create_matrix(self._a) - self._L = _create_form(L, form_compiler_options=form_compiler_options, jit_options=jit_options) self._b = create_vector(self._L) @@ -557,7 +594,7 @@ def __init__(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[DirichletBC] = [], else: self.u = u - self._x = _cpp.la.petsc.create_vector_wrap(self.u.x._cpp_object) + self._x = la.create_petsc_vector_wrap(self.u.x) self.bcs = bcs self._solver = PETSc.KSP().create(self.u.function_space.mesh.comm) @@ -592,7 +629,7 @@ def solve(self) -> _Function: # Assemble lhs self._A.zeroEntries() - _assemble_matrix_mat(self._A, self._a, bcs=self.bcs) + assemble_matrix_mat(self._A, self._a, bcs=self.bcs) self._A.assemble() # Assemble rhs @@ -646,7 +683,7 @@ class NonlinearProblem: def __init__(self, F: ufl.form.Form, u: _Function, bcs: typing.List[DirichletBC] = [], J: ufl.form.Form = None, form_compiler_options={}, jit_options={}): - """Initialize solver for solving a non-linear problem using Newton's method, :math:`dF/du(u) du = -F(u)`. + """Initialize solver for solving a non-linear problem using Newton's method, :math:`(dF/du)(u) du = -F(u)`. Args: F: The PDE residual F(u, v) @@ -689,7 +726,7 @@ def a(self) -> Form: """Compiled bilinear form (the Jacobian form)""" return self._a - def form(self, x: PETSc.Vec): + def form(self, x: PETSc.Vec) -> None: """This function is called before the residual or Jacobian is computed. This is usually used to update ghost values. @@ -699,7 +736,7 @@ def form(self, x: PETSc.Vec): """ x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD) - def F(self, x: PETSc.Vec, b: PETSc.Vec): + def F(self, x: PETSc.Vec, b: PETSc.Vec) -> None: """Assemble the residual F into the vector b. Args: @@ -717,7 +754,7 @@ def F(self, x: PETSc.Vec, b: PETSc.Vec): b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) set_bc(b, self.bcs, x, -1.0) - def J(self, x: PETSc.Vec, A: PETSc.Mat): + def J(self, x: PETSc.Vec, A: PETSc.Mat) -> None: """Assemble the Jacobian matrix. Args: @@ -725,19 +762,19 @@ def J(self, x: PETSc.Vec, A: PETSc.Mat): """ A.zeroEntries() - _assemble_matrix_mat(A, self._a, self.bcs) + assemble_matrix_mat(A, self._a, self.bcs) A.assemble() def load_petsc_lib(loader: typing.Callable[[str], typing.Any]) -> typing.Any: - """ - Load PETSc shared library using loader callable, e.g. ctypes.CDLL. + """Load PETSc shared library using loader callable, e.g. ctypes.CDLL. Args: loader: A callable that accepts a library path and returns a wrapped library. Returns: A wrapped library of the type returned by the callable. + """ petsc_lib = None diff --git a/python/dolfinx/geometry.py b/python/dolfinx/geometry.py index 45de5f96b9a..0263f7be6f1 100644 --- a/python/dolfinx/geometry.py +++ b/python/dolfinx/geometry.py @@ -28,7 +28,7 @@ class BoundingBoxTree: """Bounding box trees used in collision detection.""" def __init__(self, tree): - """Wrap a C++ BoundingBoxTree + """Wrap a C++ BoundingBoxTree. Note: This initializer should not be used in user code. Use @@ -39,14 +39,14 @@ def __init__(self, tree): @property def num_bboxes(self) -> int: - """Number of bounding boxes""" + """Number of bounding boxes.""" return self._cpp_object.num_bboxes def get_bbox(self, i) -> npt.NDArray[np.floating]: """Get lower and upper corners of the ith bounding box. Args: - i: Index of the box + i: Index of the box. Returns: The 'lower' and 'upper' points of the bounding box. @@ -64,11 +64,11 @@ def bb_tree(mesh: Mesh, dim: int, entities: typing.Optional[npt.NDArray[np.int32 """Create a bounding box tree for use in collision detection. Args: - mesh: The mesh - dim: The dimension of the mesh entities - entities: List of entity indices (local to process). If not supplied, - all owned and ghosted entities are used. - padding: Padding for each bounding box + mesh: The mesh. + dim: Dimension of the mesh entities to build bounding box for. + entities: List of entity indices (local to process). If not + supplied, all owned and ghosted entities are used. + padding: Padding for each bounding box. Returns: Bounding box tree. @@ -93,8 +93,8 @@ def compute_collisions_trees(tree0: BoundingBoxTree, tree1: BoundingBoxTree) -> """Compute all collisions between two bounding box trees. Args: - tree0: First bounding box tree - tree1: Second bounding box tree + tree0: First bounding box tree. + tree1: Second bounding box tree. Returns: List of pairs of intersecting box indices from each tree. Shape @@ -111,8 +111,8 @@ def compute_collisions_points(tree: BoundingBoxTree, x: npt.NDArray[np.floating] than one box. Args: - tree: Bounding box tree - x: Points (`shape=(num_points, 3)`) + tree: Bounding box tree. + x: Points (``shape=(num_points, 3)``). Returns: For each point, the bounding box leaves that collide with the @@ -127,14 +127,15 @@ def compute_closest_entity(tree: BoundingBoxTree, midpoint_tree: BoundingBoxTree """Compute closest mesh entity to a point. Args: - tree: bounding box tree for the entities + tree: bounding box tree for the entities. midpoint_tree: A bounding box tree with the midpoints of all the mesh entities. This is used to accelerate the search. - mesh: The mesh - points: The points to check for collision, shape=(num_points, 3) + mesh: The mesh. + points: The points to check for collision, ``shape=(num_points, + 3)``. Returns: - Mesh entity index for each point in `points`. Returns -1 for a + Mesh entity index for each point in ``points``. Returns -1 for a point if the bounding box tree is empty. """ @@ -160,14 +161,16 @@ def compute_colliding_cells(mesh: Mesh, candidates: AdjacencyList_int32, x: npt. """From a mesh, find which cells collide with a set of points. Args: - mesh: The mesh + mesh: The mesh. candidate_cells: Adjacency list of candidate colliding cells for - the ith point in `x` - points: The points to check for collision shape=(num_points, 3) + the ith point in ``x``. + points: The points to check for collision ``shape=(num_points, + 3)``, Returns: Adjacency list where the ith node is the list of entities that - collide with the ith point + collide with the ith point. + """ return _cpp.geometry.compute_colliding_cells(mesh._cpp_object, candidates, x) @@ -179,14 +182,14 @@ def squared_distance(mesh: Mesh, dim: int, entities: typing.List[int], points: n input entity. Args: - mesh: Mesh containing the entities - dim: Topological dimension of the mesh entities - entities: Indices of the mesh entities (local to process) + mesh: Mesh containing the entities. + dim: Topological dimension of the mesh entities. + entities: Indices of the mesh entities (local to process). points: Points to compute the shortest distance from - (``shape=(num_points, 3)``) + (``shape=(num_points, 3)``). Returns: - Squared shortest distance from points[i] to entities[i] + Squared shortest distance from ``points[i]`` to ``entities[i]``. """ return _cpp.geometry.squared_distance(mesh._cpp_object, dim, entities, points) @@ -198,8 +201,8 @@ def compute_distance_gjk(p: npt.NDArray[np.floating], q: npt.NDArray[np.floating Uses the Gilbert–Johnson–Keerthi (GJK) distance algorithm. Args: - p: Body 1 list of points (``shape=(num_points, gdim)``) - q: Body 2 list of points (``shape=(num_points, gdim)``) + p: Body 1 list of points (``shape=(num_points, gdim)``). + q: Body 2 list of points (``shape=(num_points, gdim)``). Returns: Shortest vector between the two bodies. diff --git a/python/dolfinx/graph.py b/python/dolfinx/graph.py index 4ebd8ba6c29..9827a23179c 100644 --- a/python/dolfinx/graph.py +++ b/python/dolfinx/graph.py @@ -3,7 +3,7 @@ # This file is part of DOLFINx (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -"""Graph representations and operations on graphs""" +"""Graph representations and operations on graphs.""" from __future__ import annotations @@ -39,12 +39,12 @@ def create_adjacencylist(data: np.ndarray, offsets=None): Args: data: The adjacency array. If the array is one-dimensional, - offsets should be supplied. If the array is two-dimensional the - number of edges per node is the second dimension. + offsets should be supplied. If the array is two-dimensional + the number of edges per node is the second dimension. offsets: The offsets array with the number of edges per node. Returns: - An adjacency list + An adjacency list. """ if offsets is None: diff --git a/python/dolfinx/io/utils.py b/python/dolfinx/io/utils.py index 364e468ae0c..0ab8239b2e7 100644 --- a/python/dolfinx/io/utils.py +++ b/python/dolfinx/io/utils.py @@ -4,24 +4,23 @@ # This file is part of DOLFINx (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -"""IO module for input data and post-processing file output""" +"""IO module for input data and post-processing file output.""" import typing -import numpy as np -import numpy.typing as npt - import basix import basix.ufl +import numpy as np +import numpy.typing as npt import ufl -from dolfinx import cpp as _cpp from dolfinx.cpp.io import perm_gmsh as cell_perm_gmsh # noqa F401 from dolfinx.cpp.io import perm_vtk as cell_perm_vtk # noqa F401 from dolfinx.fem import Function from dolfinx.mesh import GhostMode, Mesh, MeshTags - from mpi4py import MPI as _MPI +from dolfinx import cpp as _cpp + __all__ = ["VTKFile", "XDMFFile", "cell_perm_gmsh", "cell_perm_vtk", "distribute_entity_data"] @@ -39,14 +38,13 @@ def _extract_cpp_functions(functions: typing.Union[typing.List[Function], Functi __all__ = __all__ + ["FidesWriter", "VTXWriter"] class VTXWriter: - """Interface to VTK files for ADIOS2 + """Writer for VTX files, using ADIOS2 to create the files. VTX supports arbitrary order Lagrange finite elements for the geometry description and arbitrary order (discontinuous) Lagrange finite elements for Functions. - The files can be displayed by Paraview. The storage backend uses - ADIOS2. + The files can be displayed by Paraview. """ @@ -105,15 +103,14 @@ def close(self): self._cpp_object.close() class FidesWriter: - """Interface to Fides file format. + """Writer for Fides files, using ADIOS2 to create the files. - Fides supports first order Lagrange finite elements for the - geometry description and first order Lagrange finite elements - for functions. All functions has to be of the same element - family and same order. + Fides (https://fides.readthedocs.io/) supports first order + Lagrange finite elements for the geometry description and first + order Lagrange finite elements for functions. All functions have + to be of the same element family and same order. - The files can be displayed by Paraview. The storage backend uses - ADIOS2. + The files can be displayed by Paraview. """ @@ -124,9 +121,9 @@ def __init__(self, comm: _MPI.Comm, filename: str, output: typing.Union[Mesh, ty element family and degree Args: - comm: The MPI communicator - filename: The output filename - output: The data to output. Either a mesh, a single + comm: MPI communicator. + filename: Output filename. + output: Data to output. Either a mesh, a single first order Lagrange function or list of first order Lagrange functions. engine: ADIOS2 engine to use for output. See @@ -168,7 +165,7 @@ def close(self): class VTKFile(_cpp.io.VTKFile): - """Interface to VTK files + """Interface to VTK files. VTK supports arbitrary order Lagrange finite elements for the geometry description. XDMF is the preferred format for geometry @@ -212,14 +209,17 @@ def write_function(self, u: Function, t: float = 0.0, mesh_xpath="/Xdmf/Domain/G """Write function to file for a given time. Note: - Function is interpolated onto the mesh nodes, as a Nth order Lagrange function, - where N is the order of the coordinate map. - If the Function is a cell-wise constant, it is saved as a cell-wise constant. + Function is interpolated onto the mesh nodes, as a Nth order + Lagrange function, where N is the order of the coordinate + map. If the Function is a cell-wise constant, it is saved as + a cell-wise constant. Args: - u: The Function to write to file. + u: Function to write to file. t: Time associated with Function output . - mesh_xpath: Path to mesh associated with the Function in the XDMFFile. + mesh_xpath: Path to mesh associated with the Function in the + XDMFFile. + """ super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath) diff --git a/python/dolfinx/jit.py b/python/dolfinx/jit.py index 033000537d3..0ac28a13100 100644 --- a/python/dolfinx/jit.py +++ b/python/dolfinx/jit.py @@ -17,7 +17,7 @@ from mpi4py import MPI -__all__ = ["ffcx_jit", "get_options"] +__all__ = ["ffcx_jit", "get_options", "mpi_jit_decorator"] DOLFINX_DEFAULT_JIT_OPTIONS = { "cache_dir": @@ -121,13 +121,14 @@ def get_options(priority_options: Optional[dict] = None) -> dict: """Return a copy of the merged JIT option values for DOLFINx. Args: - priority_options: take priority over all other option values (see notes) + priority_options: Take priority over all other option values + (see notes). Returns: - dict: merged option values + dict: Merged option values. Notes: - See ffcx_jit for user facing documentation. + See :func:`ffcx_jit` for user facing documentation. """ options = {} @@ -153,15 +154,15 @@ def ffcx_jit(ufl_object, form_compiler_options={}, jit_options={}): """Compile UFL object with FFCx and CFFI. Args: - ufl_object: Object to compile, e.g. ufl.Form + ufl_object: Object to compile, e.g. ``ufl.Form``. form_compiler_options: Options used in FFCx compilation of - this form. Run `ffcx --help` at the commandline to see all - available options. Takes priority over all other option + this form. Run ``ffcx --help`` at the command line to see + all available options. Takes priority over all other option values. jit_options: Options used in CFFI JIT compilation of C code - generated by FFCx. See `python/dolfinx/jit.py` for all - available options. Takes priority over all other - option values. + generated by FFCx. See ``python/dolfinx/jit.py`` for all + available options. Takes priority over all other option + values. Returns: (compiled object, module, (header code, implementation code)) @@ -195,7 +196,6 @@ def ffcx_jit(ufl_object, form_compiler_options={}, jit_options={}): **{ "cffi_extra_compile_args": ["-O2", "-march=native" ], "cffi_verbose": True }** """ - # Prepare form compiler options with priority options p_ffcx = ffcx.get_options(form_compiler_options) p_jit = get_options(jit_options) diff --git a/python/dolfinx/la.py b/python/dolfinx/la.py index c7642b42f95..8cef960cc13 100644 --- a/python/dolfinx/la.py +++ b/python/dolfinx/la.py @@ -10,12 +10,11 @@ import numpy as np from dolfinx.cpp.common import IndexMap from dolfinx.cpp.la import BlockMode, InsertMode, Norm -from dolfinx.cpp.la.petsc import create_vector as create_petsc_vector from dolfinx import cpp as _cpp -__all__ = ["orthonormalize", "is_orthonormal", "create_petsc_vector", "matrix_csr", "vector", - "MatrixCSR", "Norm", "InsertMode", "Vector", ] +__all__ = ["orthonormalize", "is_orthonormal", "matrix_csr", "vector", + "MatrixCSR", "Norm", "InsertMode", "Vector", "create_petsc_vector"] class MatrixCSR: @@ -23,11 +22,7 @@ def __init__(self, A): """A distributed sparse matrix that uses compressed sparse row storage. Args: - sp: The sparsity pattern that defines the nonzero structure - of the matrix the parallel distribution of the matrix - bm: The block mode (compact or expanded). Relevant only if - block size is greater than one. - + A: The C++/pybind11 matrix object. Note: Objects of this type should be created using :func:`matrix_csr` and not created using the class @@ -39,7 +34,7 @@ def __init__(self, A): def index_map(self, i: int) -> IndexMap: """Index map for row/column. - Arg: + Args: i: 0 for row map, 1 for column map. """ @@ -48,12 +43,14 @@ def index_map(self, i: int) -> IndexMap: @property def block_size(self): """Block sizes for the matrix.""" - return self._cpp_object.block_size + return self._cpp_object.bs def add(self, x, rows, cols, bs=1) -> None: + """Add a block of values in the matrix.""" self._cpp_object.add(x, rows, cols, bs) def set(self, x, rows, cols, bs=1) -> None: + """Set a block of values in the matrix.""" self._cpp_object.set(x, rows, cols, bs) def set_value(self, x: np.floating) -> None: @@ -63,9 +60,10 @@ def set_value(self, x: np.floating) -> None: x: The value to set all non-zero entries to. """ - self._cpp_object.set(x) + self._cpp_object.set_value(x) def finalize(self) -> None: + """Scatter and accumulate ghost values.""" self._cpp_object.finalize() def squared_norm(self) -> np.floating: @@ -118,7 +116,7 @@ def to_scipy(self, ghosted=False): """ from scipy.sparse import csr_matrix as _csr, bsr_matrix as _bsr - bs0, bs1 = self._cpp_object.block_size + bs0, bs1 = self._cpp_object.bs ncols = self.index_map(1).size_local + self.index_map(1).num_ghosts if ghosted: nrows = self.index_map(0).size_local + self.index_map(0).num_ghosts @@ -179,22 +177,42 @@ def __init__(self, x): @property def index_map(self) -> IndexMap: + """Index map that describes size and parallel distribution.""" return self._cpp_object.index_map + @property + def block_size(self) -> int: + """Block size for the vector.""" + return self._cpp_object.bs + @property def array(self) -> np.ndarray: + """Local representation of the vector.""" return self._cpp_object.array - def set(self, value: np.floating) -> None: - self._cpp_object.set(value) - def scatter_forward(self) -> None: + """Update ghost entries.""" self._cpp_object.scatter_forward() - def scatter_reverse(self, mode: InsertMode): + def scatter_reverse(self, mode: InsertMode) -> None: + """Scatter ghost entries to owner. + + Args: + mode: Control how scattered values are set/accumulated by owner. + + """ self._cpp_object.scatter_reverse(mode) def norm(self, type=_cpp.la.Norm.l2) -> np.floating: + """Compute a norm of the vector. + + Args: + type: Norm type to computed. + + Returns: + Computed norm. + + """ return self._cpp_object.norm(type) @@ -232,14 +250,34 @@ def create_petsc_vector_wrap(x: Vector): x: The vector to wrap as a PETSc vector. Returns: - A PETSc vector that shares data with `x`. + A PETSc vector that shares data with ``x``. Note: - The vector `x` must not be destroyed before the returned PETSc + The vector ``x`` must not be destroyed before the returned PETSc object. """ - return _cpp.la.petsc.create_vector_wrap(x._cpp_object) + from petsc4py import PETSc + map = x.index_map + ghosts = map.ghosts.astype(PETSc.IntType) + bs = x.block_size + size = (map.size_local * bs, map.size_global * bs) + return PETSc.Vec().createGhostWithArray(ghosts, x.array, size=size, bsize=bs, comm=map.comm) + + +def create_petsc_vector(map, bs: int): + """Create a distributed PETSc vector. + + Args: + map: Index map that describes the size and parallel layout of + the vector to create. + bs: Block size of the vector. + + """ + from petsc4py import PETSc + ghosts = map.ghosts.astype(PETSc.IntType) + size = (map.size_local * bs, map.size_global * bs) + return PETSc.Vec().createGhost(ghosts, size=size, bsize=bs, comm=map.comm) def orthonormalize(basis): diff --git a/python/dolfinx/log.py b/python/dolfinx/log.py index 5398b9726f6..d927867c704 100644 --- a/python/dolfinx/log.py +++ b/python/dolfinx/log.py @@ -4,7 +4,7 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later -"""Logging module""" +"""Logging module.""" # Import pybind11 wrapped code intp dolfinx.log from dolfinx.cpp.log import (LogLevel, get_log_level, log, # noqa diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 4a6aa59f9c3..58597819243 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -29,7 +29,8 @@ "GhostMode", "build_dual_graph", "cell_dim", "compute_midpoints", "exterior_facet_indices", "compute_incident_entities", "create_cell_partitioner", "create_interval", "create_unit_interval", "create_rectangle", "create_unit_square", - "create_box", "create_unit_cube", "to_type", "to_string"] + "create_box", "create_unit_cube", "to_type", "to_string", "refine_plaza", + "transfer_meshtag"] def compute_incident_entities(topology, entities: npt.NDArray[np.int32], d0: int, d1: int): @@ -49,7 +50,8 @@ def __init__(self, mesh: _cpp.mesh.Mesh, domain: ufl.Mesh): domain: The UFL domain Note: - Mesh objects should not usually be created using this class directly. + Mesh objects should not usually be created using this class + directly. """ self._cpp_object = mesh @@ -69,7 +71,7 @@ def name(self, value): self._cpp_object.name = value def ufl_cell(self) -> ufl.Cell: - """Return the UFL cell type""" + """Return the UFL cell type.""" return ufl.Cell(self.topology.cell_name(), geometric_dimension=self.geometry.dim) def ufl_domain(self) -> ufl.Mesh: @@ -94,15 +96,15 @@ def geometry(self): def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: - """Compute mesh entities satisfying a geometric marking function + """Compute mesh entities satisfying a geometric marking function. Args: - mesh: Mesh to locate entities on - dim: Topological dimension of the mesh entities to consider - marker: A function that takes an array of points `x` with shape - `(gdim, num_points)` and returns an array of booleans of - length `num_points`, evaluating to `True` for entities to be - located. + mesh: Mesh to locate entities on. + dim: Topological dimension of the mesh entities to consider. + marker: A function that takes an array of points ``x`` with + shape ``(gdim, num_points)`` and returns an array of + booleans of length ``num_points``, evaluating to `True` for + entities to be located. Returns: Indices (local to the process) of marked mesh entities. @@ -113,7 +115,7 @@ def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: """Compute mesh entities that are connected to an owned boundary - facet and satisfy a geometric marking function + facet and satisfy a geometric marking function. For vertices and edges, in parallel this function will not necessarily mark all entities that are on the exterior boundary. For @@ -125,12 +127,12 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n communication. Args: - mesh: Mesh to locate boundary entities on + mesh: Mesh to locate boundary entities on. dim: Topological dimension of the mesh entities to consider - marker: Function that takes an array of points `x` with shape - `(gdim, num_points)` and returns an array of booleans of - length `num_points`, evaluating to `True` for entities to be - located. + marker: Function that takes an array of points ``x`` with shape + ``(gdim, num_points)`` and returns an array of booleans of + length ``num_points``, evaluating to ``True`` for entities + to be located. Returns: Indices (local to the process) of marked mesh entities. @@ -157,11 +159,13 @@ def transfer_meshtag(meshtag: MeshTags, mesh1: Mesh, parent_cell: npt.NDArray[np """Generate cell mesh tags on a refined mesh from the mesh tags on the coarse parent mesh. Args: - meshtag: Mesh tags on the coarse, parent mesh - mesh1: The refined mesh - parent_cell: Index of the parent cell for each cell in the refined mesh + meshtag: Mesh tags on the coarse, parent mesh. + mesh1: The refined mesh. + parent_cell: Index of the parent cell for each cell in the + refined mesh. parent_facet: Index of the local parent facet for each cell - in the refined mesh. Only required for transfer tags on facets. + in the refined mesh. Only required for transfer tags on + facets. Returns: Mesh tags on the refined mesh. @@ -183,13 +187,13 @@ def refine(mesh: Mesh, edges: typing.Optional[np.ndarray] = None, redistribute: Args: mesh: Mesh from which to create the refined mesh. - edges: Indices of edges to split during refinement. If `None`, + edges: Indices of edges to split during refinement. If ``None``, uniform refinement is uses. redistribute: - Refined mesh is re-partitioned if `True` + Refined mesh is re-partitioned if ``True``. Returns: - Refined mesh + Refined mesh. """ if edges is None: @@ -208,10 +212,10 @@ def refine_plaza(mesh: Mesh, edges: typing.Optional[np.ndarray] = None, redistri Args: mesh: Mesh from which to create the refined mesh. - edges: Indices of edges to split during refinement. If `None`, + edges: Indices of edges to split during refinement. If ``None``, uniform refinement is uses. redistribute: - Refined mesh is re-partitioned if `True` + Refined mesh is re-partitioned if ``True``. option: Control computation of the parent-refined mesh data. @@ -236,8 +240,10 @@ def create_mesh(comm: _MPI.Comm, cells: typing.Union[np.ndarray, _cpp.graph.Adja Args: comm: MPI communicator to define the mesh on. - cells: Cells of the mesh. `cells[i]` is the 'nodes' of cell `i`. - x: Mesh geometry ('node' coordinates), with shape ``(num_nodes, gdim)`` + cells: Cells of the mesh. ``cells[i]`` is the 'nodes' of cell + ``i``. + x: Mesh geometry ('node' coordinates), with shape ``(num_nodes, + gdim)`` domain: UFL mesh. ghost_mode: The ghost mode used in the mesh partitioning. partitioner: Function that computes the parallel distribution of @@ -294,8 +300,9 @@ def __init__(self, meshtags): A Python mesh is passed to the initializer as it may have UFL data attached that is not attached the C++ Mesh that is - associated with the C++ `meshtags` object. If `mesh` is - passed, `mesh` and `meshtags` must share the same C++ mesh. + associated with the C++ ``meshtags`` object. If `mesh` is + passed, ``mesh`` and ``meshtags`` must share the same C++ + mesh. """ self._cpp_object = meshtags @@ -335,10 +342,10 @@ def find(self, value) -> npt.NDArray[np.int32]: """Get a list of all entity indices with a given value. Args: - value: Mesh tag value to search for + value: Mesh tag value to search for. Return: - Indices of entities with tag `value` + Indices of entities with tag ``value``. """ return self._cpp_object.find(value) @@ -349,15 +356,15 @@ def meshtags(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32], """Create a MeshTags object that associates data with a subset of mesh entities. Args: - mesh: The mesh - dim: Topological dimension of the mesh entity + mesh: The mesh. + dim: Topological dimension of the mesh entity. entities: Indices (local to process) of entities to associate values with. The array must be sorted and must not contain duplicates. - values: The corresponding value for each entity + values: The corresponding value for each entity. Returns: - A MeshTags object + A mesh tags object. Note: The type of the returned MeshTags is inferred from the type of @@ -392,14 +399,14 @@ def meshtags_from_entities(mesh: Mesh, dim: int, entities: _cpp.graph.AdjacencyL mesh entities, where the entities are defined by their vertices. Args: - mesh: The mesh - dim: Topological dimension of the mesh entity + mesh: The mesh. + dim: Topological dimension of the mesh entity. entities: Entities to associated values with, with entities - defined by their vertices - values: The corresponding value for each entity + defined by their vertices. + values: The corresponding value for each entity. Returns: - A MeshTags object + A mesh tags object. Note: The type of the returned MeshTags is inferred from the type of @@ -425,15 +432,15 @@ def create_interval(comm: _MPI.Comm, nx: int, points: npt.ArrayLike, comm: MPI communicator nx: Number of cells points: Coordinates of the end points - dtype: Float type for the mesh geometry (`numpy.float32` or - `numpy.float64`) + dtype: Float type for the mesh geometry (``numpy.float32`` or + ``numpy.float64``) ghost_mode: Ghost mode used in the mesh partitioning. Options - are `GhostMode.none' and `GhostMode.shared_facet`. + are ``GhostMode.none`` and ``GhostMode.shared_facet``. partitioner: Partitioning function to use for determining the parallel distribution of cells across MPI ranks Returns: - An interval mesh + An interval mesh. """ if partitioner is None: @@ -453,18 +460,18 @@ def create_unit_interval(comm: _MPI.Comm, nx: int, dtype: typing.Optional[npt.DT """Create a mesh on the unit interval. Args: - comm: MPI communicator - nx: Number of cells - points: Coordinates of the end points - dtype: Float type for the mesh geometry (`numpy.float32` or - `numpy.float64`) + comm: MPI communicator. + nx: Number of cells. + points: Coordinates of the end points. + dtype: Float type for the mesh geometry (``numpy.float32`` or + ``numpy.float64``). ghost_mode: Ghost mode used in the mesh partitioning. Options - are `GhostMode.none' and `GhostMode.shared_facet`. + are ``GhostMode.none`` and ``GhostMode.shared_facet``. partitioner: Partitioning function to use for determining the - parallel distribution of cells across MPI ranks + parallel distribution of cells across MPI ranks. Returns: - A unit interval mesh with end points at 0 and 1 + A unit interval mesh with end points at 0 and 1. """ if partitioner is None: @@ -479,22 +486,22 @@ def create_rectangle(comm: _MPI.Comm, points: npt.ArrayLike, n: npt.ArrayLike, """Create a rectangle mesh. Args: - comm: MPI communicator - points: Coordinates of the lower-left and upper-right corners of the - rectangle - n: Number of cells in each direction - cell_type: Mesh cell type - dtype: Float type for the mesh geometry (`numpy.float32` or - `numpy.float64`) - ghost_mode: Ghost mode used in the mesh partitioning + comm: MPI communicator. + points: Coordinates of the lower-left and upper-right corners of + the rectangle. + n: Number of cells in each direction. + cell_type: Mesh cell type. + dtype: Float type for the mesh geometry (``numpy.float32`` or + ``numpy.float64``) + ghost_mode: Ghost mode used in the mesh partitioning. partitioner: Function that computes the parallel distribution of - cells across MPI ranks + cells across MPI ranks. diagonal: Direction of diagonal of triangular meshes. The options are ``left``, ``right``, ``crossed``, ``left/right``, ``right/left``. Returns: - A mesh of a rectangle + A mesh of a rectangle. """ if partitioner is None: @@ -516,20 +523,20 @@ def create_unit_square(comm: _MPI.Comm, nx: int, ny: int, cell_type=CellType.tri """Create a mesh of a unit square. Args: - comm: MPI communicator - nx: Number of cells in the "x" direction - ny: Number of cells in the "y" direction - cell_type: Mesh cell type - dtype: Float type for the mesh geometry (`numpy.float32` or - `numpy.float64`) - ghost_mode: Ghost mode used in the mesh partitioning - partitioner:Function that computes the parallel distribution of cells across - MPI ranks + comm: MPI communicator. + nx: Number of cells in the "x" direction. + ny: Number of cells in the "y" direction. + cell_type: Mesh cell type. + dtype: Float type for the mesh geometry (``numpy.float32`` or + ``numpy.float64``). + ghost_mode: Ghost mode used in the mesh partitioning. + partitioner:Function that computes the parallel distribution of + cells across MPI ranks. diagonal: - Direction of diagonal + Direction of diagonal. Returns: - A mesh of a square with corners at (0, 0) and (1, 1) + A mesh of a square with corners at (0, 0) and (1, 1). """ if partitioner is None: @@ -546,19 +553,19 @@ def create_box(comm: _MPI.Comm, points: typing.List[npt.ArrayLike], n: list, """Create a box mesh. Args: - comm: MPI communicator + comm: MPI communicator. points: Coordinates of the 'lower-left' and 'upper-right' - corners of the box + corners of the box. n: List of cells in each direction - cell_type: The cell type - dtype: Float type for the mesh geometry (`numpy.float32` or - `numpy.float64`) - ghost_mode: The ghost mode used in the mesh partitioning + cell_type: The cell type. + dtype: Float type for the mesh geometry (``numpy.float32`` or + ``numpy.float64``). + ghost_mode: The ghost mode used in the mesh partitioning. partitioner: Function that computes the parallel distribution of - cells across MPI ranks + cells across MPI ranks. Returns: - A mesh of a box domain + A mesh of a box domain. """ if partitioner is None: @@ -579,20 +586,20 @@ def create_unit_cube(comm: _MPI.Comm, nx: int, ny: int, nz: int, cell_type=CellT """Create a mesh of a unit cube. Args: - comm: MPI communicator - nx: Number of cells in "x" direction - ny: Number of cells in "y" direction - nz: Number of cells in "z" direction + comm: MPI communicator. + nx: Number of cells in "x" direction. + ny: Number of cells in "y" direction. + nz: Number of cells in "z" direction. cell_type: Mesh cell type - dtype: Float type for the mesh geometry (`numpy.float32` or - `numpy.float64`) - ghost_mode: Ghost mode used in the mesh partitioning + dtype: Float type for the mesh geometry (``numpy.float32`` or + ``numpy.float64``). + ghost_mode: Ghost mode used in the mesh partitioning. partitioner: Function that computes the parallel distribution of - cells across MPI ranks + cells across MPI ranks. Returns: A mesh of an axis-aligned unit cube with corners at (0, 0, 0) - and (1, 1, 1) + and (1, 1, 1). """ if partitioner is None: diff --git a/python/dolfinx/nls/__init__.py b/python/dolfinx/nls/__init__.py index b1a70a2b453..770e4c30191 100644 --- a/python/dolfinx/nls/__init__.py +++ b/python/dolfinx/nls/__init__.py @@ -4,7 +4,3 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later """Tools for solving nonlinear problems.""" - -from dolfinx.nls import petsc - -__all__ = ["petsc"] diff --git a/python/dolfinx/nls/petsc.py b/python/dolfinx/nls/petsc.py index d4e377a6a1b..8b07c32fa45 100644 --- a/python/dolfinx/nls/petsc.py +++ b/python/dolfinx/nls/petsc.py @@ -18,6 +18,7 @@ from dolfinx import cpp as _cpp from dolfinx import fem +from dolfinx.fem.petsc import create_matrix, create_vector __all__ = ["NewtonSolver"] @@ -29,9 +30,9 @@ def __init__(self, comm: MPI.Intracomm, problem: NonlinearProblem): # Create matrix and vector to be used for assembly # of the non-linear problem - self._A = fem.petsc.create_matrix(problem.a) + self._A = create_matrix(problem.a) self.setJ(problem.J, self._A) - self._b = fem.petsc.create_vector(problem.L) + self._b = create_vector(problem.L) self.setF(problem.F, self._b) self.set_form(problem.form) diff --git a/python/dolfinx/pkgconfig.py b/python/dolfinx/pkgconfig.py index ee09d57947f..7cca48b30d4 100644 --- a/python/dolfinx/pkgconfig.py +++ b/python/dolfinx/pkgconfig.py @@ -24,13 +24,13 @@ def _pkgconfig_query(s): return (rc, out.rstrip().decode('utf-8')) -def exists(package): - "Test for the existence of a pkg-config file for a named package" - return (_pkgconfig_query("--exists " + package)[0] == 0) +def exists(package) -> bool: + """Test for the existence of a pkg-config file for a named package.""" + return _pkgconfig_query("--exists " + package)[0] == 0 def parse(package): - "Return a dict containing compile-time definitions" + """Return a dict containing compile-time definitions.""" parse_map = { '-D': 'define_macros', '-I': 'include_dirs', @@ -40,11 +40,11 @@ def parse(package): result = {x: [] for x in parse_map.values()} - # Execute the query to pkg-config and clean the result. + # Execute the query to pkg-config and clean the result out = _pkgconfig_query(package + ' --cflags --libs')[1] out = out.replace('\\"', '') - # Iterate through each token in the output. + # Iterate through each token in the output for token in out.split(): key = parse_map.get(token[:2]) if key: diff --git a/python/dolfinx/plot.py b/python/dolfinx/plot.py index b9028ebe306..519ab83dfb4 100644 --- a/python/dolfinx/plot.py +++ b/python/dolfinx/plot.py @@ -35,6 +35,14 @@ def create_vtk_mesh(msh: mesh.Mesh, dim: typing.Optional[int] = None, entities=N dimension. The vertex indices in the returned topology array are the indices for the associated entry in the mesh geometry. + Args: + mesh: Mesh to extract data from. + dim: Topological dimension of entities to extract. + entities: Entities to extract. Extract all if ``None``. + + Returns: + Topology, type for each cell, and geometry in VTK-ready format. + """ if dim is None: dim = msh.topology.dim @@ -83,9 +91,21 @@ def create_vtk_mesh(msh: mesh.Mesh, dim: typing.Optional[int] = None, entities=N @create_vtk_mesh.register(fem.FunctionSpace) def _(V: fem.FunctionSpace, entities=None): """Creates a VTK mesh topology (topology array and array of cell - types) that is based on the degree-of-freedom coordinates. Note that - this function supports Lagrange elements (continuous and - discontinuous) only. + types) that is based on the degree-of-freedom coordinates. + + This function supports visualisation when the degree of the finite + element space is different from the geometric degree of the mesh. + + Note: + This function supports Lagrange elements (continuous and + discontinuous) only. + + Args: + V: Mesh to extract data from. + entities: Entities to extract. Extract all if ``None``. + + Returns: + Topology, type for each cell, and geometry in VTK-ready format. """ if not (V.ufl_element().family() in ['Discontinuous Lagrange', "Lagrange", "DQ", "Q", "DP", "P"]): diff --git a/python/dolfinx/wrappers/common.cpp b/python/dolfinx/wrappers/common.cpp index 643da6ac8c0..bcbdd22e32b 100644 --- a/python/dolfinx/wrappers/common.cpp +++ b/python/dolfinx/wrappers/common.cpp @@ -88,6 +88,8 @@ void common(py::module& m) }), py::arg("comm"), py::arg("local_size"), py::arg("dest_src"), py::arg("ghosts"), py::arg("ghost_owners")) + .def_property_readonly("comm", [](const dolfinx::common::IndexMap& self) + { return MPICommWrapper(self.comm()); }) .def_property_readonly("size_local", &dolfinx::common::IndexMap::size_local) .def_property_readonly("size_global", diff --git a/python/dolfinx/wrappers/la.cpp b/python/dolfinx/wrappers/la.cpp index b02e92aa5e1..033d066f465 100644 --- a/python/dolfinx/wrappers/la.cpp +++ b/python/dolfinx/wrappers/la.cpp @@ -46,7 +46,6 @@ void declare_objects(py::module& m, const std::string& type) py::arg("vec")) .def_property_readonly("dtype", [](const dolfinx::la::Vector& self) { return py::dtype::of(); }) - .def("set", &dolfinx::la::Vector::set, py::arg("v")) .def( "norm", [](dolfinx::la::Vector& self, dolfinx::la::Norm type) @@ -92,8 +91,7 @@ void declare_objects(py::module& m, const std::string& type) py::arg("p"), py::arg("block_mode")) .def_property_readonly("dtype", [](const dolfinx::la::MatrixCSR& self) { return py::dtype::of(); }) - .def_property_readonly("block_size", - &dolfinx::la::MatrixCSR::block_size) + .def_property_readonly("bs", &dolfinx::la::MatrixCSR::block_size) .def("squared_norm", &dolfinx::la::MatrixCSR::squared_norm) .def("index_map", &dolfinx::la::MatrixCSR::index_map) .def("add", @@ -126,7 +124,7 @@ void declare_objects(py::module& m, const std::string& type) throw std::runtime_error( "Block size not supported in this function"); }) - .def("set", + .def("set_value", static_cast::*)(T)>( &dolfinx::la::MatrixCSR::set), py::arg("x")) @@ -167,13 +165,33 @@ void declare_objects(py::module& m, const std::string& type) .def("finalize_end", &dolfinx::la::MatrixCSR::finalize_end); } +// Declare objects that have multiple scalar types +template +void declare_functions(py::module& m) +{ + m.def( + "inner_product", + [](const dolfinx::la::Vector& x, const dolfinx::la::Vector& y) + { return dolfinx::la::inner_product(x, y); }, + py::arg("x"), py::arg("y")); + m.def( + "orthonormalize", + [](std::vector>> basis) + { dolfinx::la::orthonormalize(basis); }, + py::arg("basis")); + m.def( + "is_orthonormal", + [](std::vector>> + basis) { return dolfinx::la::is_orthonormal(basis); }, + py::arg("basis")); +} + } // namespace namespace dolfinx_wrappers { void la(py::module& m) { - py::enum_(m, "InsertMode") .value("add", PyInsertMode::add) .value("insert", PyInsertMode::insert); @@ -253,5 +271,10 @@ void la(py::module& m) declare_objects(m, "float64"); declare_objects>(m, "complex64"); declare_objects>(m, "complex128"); + + declare_functions(m); + declare_functions(m); + declare_functions>(m); + declare_functions>(m); } } // namespace dolfinx_wrappers diff --git a/python/dolfinx/wrappers/petsc.cpp b/python/dolfinx/wrappers/petsc.cpp index fabd1084751..9a47e688f37 100644 --- a/python/dolfinx/wrappers/petsc.cpp +++ b/python/dolfinx/wrappers/petsc.cpp @@ -23,7 +23,6 @@ #include #include #include - #include #include #include @@ -33,7 +32,6 @@ namespace { - // Declare assembler function that have multiple scalar types template void declare_petsc_discrete_operators(py::module& m) @@ -79,7 +77,6 @@ void declare_petsc_discrete_operators(py::module& m) return A; }, py::return_value_policy::take_ownership, py::arg("V0"), py::arg("V1")); - m.def( "interpolation_matrix", [](const dolfinx::fem::FunctionSpace& V0, @@ -123,17 +120,6 @@ void declare_petsc_discrete_operators(py::module& m) void petsc_la_module(py::module& m) { - m.def("create_vector", - py::overload_cast( - &dolfinx::la::petsc::create_vector), - py::return_value_policy::take_ownership, py::arg("index_map"), - py::arg("bs"), "Create a ghosted PETSc Vec for index map."); - m.def( - "create_vector_wrap", - [](dolfinx::la::Vector& x) - { return dolfinx::la::petsc::create_vector_wrap(x); }, - py::return_value_policy::take_ownership, py::arg("x"), - "Create a ghosted PETSc Vec that wraps a DOLFINx Vector"); m.def( "create_matrix", [](dolfinx_wrappers::MPICommWrapper comm, @@ -142,11 +128,9 @@ void petsc_la_module(py::module& m) py::return_value_policy::take_ownership, py::arg("comm"), py::arg("p"), py::arg("type") = std::string(), "Create a PETSc Mat from sparsity pattern."); - // TODO: check reference counting for index sets m.def("create_index_sets", &dolfinx::la::petsc::create_index_sets, py::arg("maps"), py::return_value_policy::take_ownership); - m.def( "scatter_local_vectors", [](Vec x, @@ -163,6 +147,7 @@ void petsc_la_module(py::module& m) py::arg("x"), py::arg("x_b"), py::arg("maps"), "Scatter the (ordered) list of sub vectors into a block " "vector."); + m.def( "get_local_vectors", [](const Vec x, @@ -344,7 +329,6 @@ void petsc_nls_module(py::module& m) namespace dolfinx_wrappers { - void petsc(py::module& m_fem, py::module& m_la, py::module& m_nls) { py::module petsc_fem_mod @@ -356,7 +340,7 @@ void petsc(py::module& m_fem, py::module& m_la, py::module& m_nls) petsc_la_module(petsc_la_mod); py::module petsc_nls_mod - = m_nls.def_submodule("petsc", "PETSc-specific non-linear solvers"); + = m_nls.def_submodule("petsc", "PETSc-specific nonlinear solvers"); petsc_nls_module(petsc_nls_mod); } } // namespace dolfinx_wrappers diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 65933c18c63..2350f55f0ae 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -17,11 +17,20 @@ VectorFunctionSpace, assemble_scalar, bcs_by_block, dirichletbc, extract_function_spaces, form, locate_dofs_geometrical, locate_dofs_topological) -from dolfinx.fem.petsc import (apply_lifting, apply_lifting_nest, - assemble_matrix, assemble_matrix_block, - assemble_matrix_nest, assemble_vector, - assemble_vector_block, assemble_vector_nest, - set_bc, set_bc_nest) +from dolfinx.fem.petsc import apply_lifting as petsc_apply_lifting +from dolfinx.fem.petsc import apply_lifting_nest as petsc_apply_lifting_nest +from dolfinx.fem.petsc import assemble_matrix as petsc_assemble_matrix +from dolfinx.fem.petsc import \ + assemble_matrix_block as petsc_assemble_matrix_block +from dolfinx.fem.petsc import \ + assemble_matrix_nest as petsc_assemble_matrix_nest +from dolfinx.fem.petsc import assemble_vector as petsc_assemble_vector +from dolfinx.fem.petsc import \ + assemble_vector_block as petsc_assemble_vector_block +from dolfinx.fem.petsc import \ + assemble_vector_nest as petsc_assemble_vector_nest +from dolfinx.fem.petsc import set_bc as petsc_set_bc +from dolfinx.fem.petsc import set_bc_nest as petsc_set_bc_nest from dolfinx.mesh import (CellType, GhostMode, create_mesh, create_rectangle, create_unit_cube, create_unit_square, locate_entities_boundary) @@ -75,81 +84,78 @@ def test_assemble_functional_ds(mode, dtype): assert value == pytest.approx(4.0, 1e-6) -def test_assemble_derivatives(): +@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128]) +def test_assemble_derivatives(dtype): """This test checks the original_coefficient_positions, which may change under differentiation (some coefficients and constants are eliminated)""" - mesh = create_unit_square(MPI.COMM_WORLD, 12, 12) + mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, dtype=dtype(0).real.dtype) Q = FunctionSpace(mesh, ("Lagrange", 1)) - u = Function(Q) + u = Function(Q, dtype=dtype) v = ufl.TestFunction(Q) du = ufl.TrialFunction(Q) - b = Function(Q) - c1 = Constant(mesh, np.array([[1.0, 0.0], [3.0, 4.0]], PETSc.ScalarType)) - c2 = Constant(mesh, PETSc.ScalarType(2.0)) + b = Function(Q, dtype=dtype) + c1 = Constant(mesh, np.array([[1.0, 0.0], [3.0, 4.0]], dtype=dtype)) + c2 = Constant(mesh, dtype(2.0)) b.x.array[:] = 2.0 # derivative eliminates 'u' and 'c1' L = ufl.inner(c1, c1) * v * dx + c2 * b * inner(u, v) * dx - a = form(derivative(L, u, du)) + a = form(derivative(L, u, du), dtype=dtype) - A1 = assemble_matrix(a) - A1.assemble() - a = form(c2 * b * inner(du, v) * dx) - A2 = assemble_matrix(a) - A2.assemble() - assert (A1 - A2).norm() == pytest.approx(0.0, rel=1e-12, abs=1e-12) - A1.destroy(), A2.destroy() + A1 = fem.assemble_matrix(a) + A1.finalize() + a = form(c2 * b * inner(du, v) * dx, dtype=dtype) + A2 = fem.assemble_matrix(a) + A2.finalize() + assert np.allclose(A1.data, A2.data) @pytest.mark.parametrize("mode", [GhostMode.none, GhostMode.shared_facet]) -def test_basic_assembly(mode): - mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, ghost_mode=mode) +@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128]) +def test_basic_assembly(mode, dtype): + mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, ghost_mode=mode, dtype=dtype(0).real.dtype) V = FunctionSpace(mesh, ("Lagrange", 1)) u, v = ufl.TrialFunction(V), ufl.TestFunction(V) - f = Function(V) + f = Function(V, dtype=dtype) f.x.array[:] = 10.0 a = inner(f * u, v) * dx + inner(u, v) * ds L = inner(f, v) * dx + inner(2.0, v) * ds - a, L = form(a), form(L) + a, L = form(a, dtype=dtype), form(L, dtype=dtype) # Initial assembly - A = assemble_matrix(a) - A.assemble() - assert isinstance(A, PETSc.Mat) - b = assemble_vector(L) - b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - assert isinstance(b, PETSc.Vec) + A = fem.assemble_matrix(a) + A.finalize() + assert isinstance(A, la.MatrixCSR) + b = fem.assemble_vector(L) + b.scatter_reverse(la.InsertMode.add) + assert isinstance(b, la.Vector) # Second assembly - normA = A.norm() - A.zeroEntries() - A = assemble_matrix(A, a) - A.assemble() - assert isinstance(A, PETSc.Mat) - assert normA == pytest.approx(A.norm()) + normA = A.squared_norm() + A.set_value(0) + A = fem.assemble_matrix(A, a) + A.finalize() + assert isinstance(A, la.MatrixCSR) + assert normA == pytest.approx(A.squared_norm()) normb = b.norm() - with b.localForm() as b_local: - b_local.set(0.0) - b = assemble_vector(b, L) - b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - assert isinstance(b, PETSc.Vec) + b.array[:] = 0 + fem.assemble_vector(b.array, L) + b.scatter_reverse(la.InsertMode.add) assert normb == pytest.approx(b.norm()) # Vector re-assembly - no zeroing (but need to zero ghost entries) - with b.localForm() as b_local: - b_local.array[b.local_size:] = 0.0 - assemble_vector(b, L) - b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) + b.array[b.index_map.size_local * b.block_size:] = 0 + fem.assemble_vector(b.array, L) + b.scatter_reverse(la.InsertMode.add) assert 2.0 * normb == pytest.approx(b.norm()) # Matrix re-assembly (no zeroing) - assemble_matrix(A, a) - A.assemble() - assert 2.0 * normA == pytest.approx(A.norm()) - A.destroy(), b.destroy() + fem.assemble_matrix(A, a) + A.finalize() + assert 4 * normA == pytest.approx(A.squared_norm()) @pytest.mark.parametrize("mode", [GhostMode.none, GhostMode.shared_facet]) @@ -157,13 +163,12 @@ def test_basic_assembly_petsc_matrixcsr(mode): mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, ghost_mode=mode) V = FunctionSpace(mesh, ("Lagrange", 1)) u, v = ufl.TrialFunction(V), ufl.TestFunction(V) - a = inner(u, v) * dx + inner(u, v) * ds - a = form(a) + a = form(inner(u, v) * dx + inner(u, v) * ds) A0 = fem.assemble_matrix(a) A0.finalize() assert isinstance(A0, la.MatrixCSR) - A1 = fem.petsc.assemble_matrix(a) + A1 = petsc_assemble_matrix(a) A1.assemble() assert isinstance(A1, PETSc.Mat) assert np.sqrt(A0.squared_norm()) == pytest.approx(A1.norm(), 1.0e-5) @@ -175,7 +180,7 @@ def test_basic_assembly_petsc_matrixcsr(mode): A0 = fem.assemble_matrix(a) A0.finalize() assert isinstance(A0, la.MatrixCSR) - A1 = fem.petsc.assemble_matrix(a) + A1 = petsc_assemble_matrix(a) A1.assemble() assert isinstance(A1, PETSc.Mat) assert np.sqrt(A0.squared_norm()) == pytest.approx(A1.norm(), rel=1.0e-8, abs=1.0e-5) @@ -187,36 +192,35 @@ def test_assembly_bcs(mode): mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, ghost_mode=mode) V = FunctionSpace(mesh, ("Lagrange", 1)) u, v = ufl.TrialFunction(V), ufl.TestFunction(V) - a = inner(u, v) * dx + inner(u, v) * ds - L = inner(1.0, v) * dx - a, L = form(a), form(L) + a = form(inner(u, v) * dx + inner(u, v) * ds) + L = form(inner(1.0, v) * dx) bdofsV = locate_dofs_geometrical(V, lambda x: np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0))) bc = dirichletbc(PETSc.ScalarType(1), bdofsV, V) # Assemble and apply 'global' lifting of bcs - A = assemble_matrix(a) + A = petsc_assemble_matrix(a) A.assemble() - b = assemble_vector(L) + b = petsc_assemble_vector(L) b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) g = b.duplicate() with g.localForm() as g_local: g_local.set(0.0) - set_bc(g, [bc]) + petsc_set_bc(g, [bc]) f = b - A * g - set_bc(f, [bc]) + petsc_set_bc(f, [bc]) # Assemble vector and apply lifting of bcs during assembly - b_bc = assemble_vector(L) - apply_lifting(b_bc, [a], [[bc]]) + b_bc = petsc_assemble_vector(L) + petsc_apply_lifting(b_bc, [a], [[bc]]) b_bc.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - set_bc(b_bc, [bc]) + petsc_set_bc(b_bc, [bc]) assert (f - b_bc).norm() == pytest.approx(0.0, rel=1e-6, abs=1e-6) A.destroy(), b.destroy(), g.destroy() @pytest.mark.skip_in_parallel -def test_assemble_manifold(): +def test_petsc_assemble_manifold(): """Test assembly of poisson problem on a mesh with topological dimension 1 but embedded in 2D (gdim=2)""" points = np.array([[0.0, 0.0], [0.2, 0.0], [0.4, 0.0], @@ -237,12 +241,12 @@ def test_assemble_manifold(): bcdofs = locate_dofs_geometrical(U, lambda x: np.isclose(x[0], 0.0)) bcs = [dirichletbc(PETSc.ScalarType(0), bcdofs, U)] - A = assemble_matrix(a, bcs=bcs) + A = petsc_assemble_matrix(a, bcs=bcs) A.assemble() - b = assemble_vector(L) - apply_lifting(b, [a], bcs=[bcs]) - set_bc(b, bcs) + b = petsc_assemble_vector(L) + petsc_apply_lifting(b, [a], bcs=[bcs]) + petsc_set_bc(b, bcs) assert np.isclose(b.norm(), 0.41231) assert np.isclose(A.norm(), 25.0199) @@ -300,32 +304,32 @@ def test_matrix_assembly_block(mode): def blocked(): """Monolithic blocked""" - A = assemble_matrix_block(a_block, bcs=[bc]) + A = petsc_assemble_matrix_block(a_block, bcs=[bc]) A.assemble() - b = assemble_vector_block(L_block, a_block, bcs=[bc]) + b = petsc_assemble_vector_block(L_block, a_block, bcs=[bc]) assert A.getType() != "nest" Anorm = A.norm() bnorm = b.norm() A.destroy(), b.destroy() with pytest.raises(RuntimeError): - assemble_matrix_block(a_block_none, bcs=[bc]) + petsc_assemble_matrix_block(a_block_none, bcs=[bc]) return Anorm, bnorm def nest(): """Nested (MatNest)""" - A = assemble_matrix_nest(a_block, bcs=[bc], mat_types=[["baij", "aij", "aij"], - ["aij", "", "aij"], - ["aij", "aij", "aij"]]) + A = petsc_assemble_matrix_nest(a_block, bcs=[bc], mat_types=[["baij", "aij", "aij"], + ["aij", "", "aij"], + ["aij", "aij", "aij"]]) A.assemble() with pytest.raises(RuntimeError): - assemble_matrix_nest(a_block_none, bcs=[bc]) + petsc_assemble_matrix_nest(a_block_none, bcs=[bc]) - b = assemble_vector_nest(L_block) - apply_lifting_nest(b, a_block, bcs=[bc]) + b = petsc_assemble_vector_nest(L_block) + petsc_apply_lifting_nest(b, a_block, bcs=[bc]) for b_sub in b.getNestSubVecs(): b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) bcs0 = bcs_by_block([L.function_spaces[0] for L in L_block], [bc]) - set_bc_nest(b, bcs0) + petsc_set_bc_nest(b, bcs0) b.assemble() bnorm = math.sqrt(sum([x.norm()**2 for x in b.getNestSubVecs()])) Anorm = nest_matrix_norm(A) @@ -345,12 +349,12 @@ def monolithic(): bdofsW_V1 = locate_dofs_topological(W.sub(1), mesh.topology.dim - 1, bndry_facets) bc = dirichletbc(u_bc, bdofsW_V1, W.sub(1)) - A = assemble_matrix(a, bcs=[bc]) + A = petsc_assemble_matrix(a, bcs=[bc]) A.assemble() - b = assemble_vector(L) - apply_lifting(b, [a], bcs=[[bc]]) + b = petsc_assemble_vector(L) + petsc_apply_lifting(b, [a], bcs=[[bc]]) b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - set_bc(b, [bc]) + petsc_set_bc(b, [bc]) assert A.getType() != "nest" Anorm = A.norm() bnorm = b.norm() @@ -411,8 +415,8 @@ def monitor(ksp, its, rnorm): def blocked(): """Blocked""" - A = assemble_matrix_block([[a00, a01], [a10, a11]], bcs=bcs) - b = assemble_vector_block([L0, L1], [[a00, a01], [a10, a11]], bcs=bcs) + A = petsc_assemble_matrix_block([[a00, a01], [a10, a11]], bcs=bcs) + b = petsc_assemble_vector_block([L0, L1], [[a00, a01], [a10, a11]], bcs=bcs) A.assemble() x = A.createVecLeft() ksp = PETSc.KSP() @@ -432,14 +436,14 @@ def blocked(): def nested(): """Nested (MatNest)""" - A = assemble_matrix_nest([[a00, a01], [a10, a11]], bcs=bcs, diagonal=1.0) + A = petsc_assemble_matrix_nest([[a00, a01], [a10, a11]], bcs=bcs, diagonal=1.0) A.assemble() - b = assemble_vector_nest([L0, L1]) - apply_lifting_nest(b, [[a00, a01], [a10, a11]], bcs=bcs) + b = petsc_assemble_vector_nest([L0, L1]) + petsc_apply_lifting_nest(b, [[a00, a01], [a10, a11]], bcs=bcs) for b_sub in b.getNestSubVecs(): b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) bcs0 = bcs_by_block([L0.function_spaces[0], L1.function_spaces[0]], bcs) - set_bc_nest(b, bcs0) + petsc_set_bc_nest(b, bcs0) b.assemble() x = b.copy() @@ -477,12 +481,12 @@ def monolithic(): bdofsW1_V1 = locate_dofs_topological(W.sub(1), facetdim, bndry_facets) bcs = [dirichletbc(u0_bc, bdofsW0_V0, W.sub(0)), dirichletbc(u1_bc, bdofsW1_V1, W.sub(1))] - A = assemble_matrix(a, bcs=bcs) + A = petsc_assemble_matrix(a, bcs=bcs) A.assemble() - b = assemble_vector(L) - apply_lifting(b, [a], [bcs]) + b = petsc_assemble_vector(L) + petsc_apply_lifting(b, [a], [bcs]) b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - set_bc(b, bcs) + petsc_set_bc(b, bcs) x = b.copy() ksp = PETSc.KSP() @@ -563,18 +567,18 @@ def boundary1(x): def nested_solve(): """Nested solver""" - A = assemble_matrix_nest(form([[a00, a01], [a10, a11]]), bcs=[bc0, bc1], - mat_types=[["baij", "aij"], ["aij", ""]]) + A = petsc_assemble_matrix_nest(form([[a00, a01], [a10, a11]]), bcs=[bc0, bc1], + mat_types=[["baij", "aij"], ["aij", ""]]) A.assemble() - P = assemble_matrix_nest(form([[p00, p01], [p10, p11]]), bcs=[bc0, bc1], - mat_types=[["aij", "aij"], ["aij", ""]]) + P = petsc_assemble_matrix_nest(form([[p00, p01], [p10, p11]]), bcs=[bc0, bc1], + mat_types=[["aij", "aij"], ["aij", ""]]) P.assemble() - b = assemble_vector_nest(form([L0, L1])) - apply_lifting_nest(b, form([[a00, a01], [a10, a11]]), [bc0, bc1]) + b = petsc_assemble_vector_nest(form([L0, L1])) + petsc_apply_lifting_nest(b, form([[a00, a01], [a10, a11]]), [bc0, bc1]) for b_sub in b.getNestSubVecs(): b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) bcs = bcs_by_block(extract_function_spaces(form([L0, L1])), [bc0, bc1]) - set_bc_nest(b, bcs) + petsc_set_bc_nest(b, bcs) b.assemble() ksp = PETSc.KSP() @@ -605,11 +609,11 @@ def monitor(ksp, its, rnorm): def blocked_solve(): """Blocked (monolithic) solver""" - A = assemble_matrix_block(form([[a00, a01], [a10, a11]]), bcs=[bc0, bc1]) + A = petsc_assemble_matrix_block(form([[a00, a01], [a10, a11]]), bcs=[bc0, bc1]) A.assemble() - P = assemble_matrix_block(form([[p00, p01], [p10, p11]]), bcs=[bc0, bc1]) + P = petsc_assemble_matrix_block(form([[p00, p01], [p10, p11]]), bcs=[bc0, bc1]) P.assemble() - b = assemble_vector_block(form([L0, L1]), form([[a00, a01], [a10, a11]]), bcs=[bc0, bc1]) + b = petsc_assemble_vector_block(form([L0, L1]), form([[a00, a01], [a10, a11]]), bcs=[bc0, bc1]) ksp = PETSc.KSP() ksp.create(mesh.comm) @@ -656,15 +660,15 @@ def monolithic_solve(): bc0 = dirichletbc(u0, bdofsW0_P2_0, W.sub(0)) bc1 = dirichletbc(u0, bdofsW0_P2_1, W.sub(0)) - A = assemble_matrix(a, bcs=[bc0, bc1]) + A = petsc_assemble_matrix(a, bcs=[bc0, bc1]) A.assemble() - P = assemble_matrix(p_form, bcs=[bc0, bc1]) + P = petsc_assemble_matrix(p_form, bcs=[bc0, bc1]) P.assemble() - b = assemble_vector(L) - apply_lifting(b, [a], bcs=[[bc0, bc1]]) + b = petsc_assemble_vector(L) + petsc_apply_lifting(b, [a], bcs=[[bc0, bc1]]) b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - set_bc(b, [bc0, bc1]) + petsc_set_bc(b, [bc0, bc1]) ksp = PETSc.KSP() ksp.create(mesh.comm) @@ -708,12 +712,12 @@ def test_basic_interior_facet_assembly(): u, v = ufl.TrialFunction(V), ufl.TestFunction(V) a = ufl.inner(ufl.avg(u), ufl.avg(v)) * ufl.dS a = form(a) - A = assemble_matrix(a) + A = petsc_assemble_matrix(a) A.assemble() assert isinstance(A, PETSc.Mat) L = ufl.conj(ufl.avg(v)) * ufl.dS L = form(L) - b = assemble_vector(L) + b = petsc_assemble_vector(L) b.assemble() assert isinstance(b, PETSc.Vec) A.destroy() @@ -739,18 +743,18 @@ def test_basic_assembly_constant(mode): a, L = form(a), form(L) # Initial assembly - A1 = assemble_matrix(a) + A1 = petsc_assemble_matrix(a) A1.assemble() - b1 = assemble_vector(L) + b1 = petsc_assemble_vector(L) b1.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) c.value = [[1.0, 2.0], [3.0, 4.0]] - A2 = assemble_matrix(a) + A2 = petsc_assemble_matrix(a) A2.assemble() - b2 = assemble_vector(L) + b2 = petsc_assemble_vector(L) b2.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) assert (A1 * 3.0 - A2 * 5.0).norm() == pytest.approx(0.0, abs=1.0e-6) assert (b1 * 3.0 - b2 * 5.0).norm() == pytest.approx(0.0, abs=1.0e-5) @@ -805,13 +809,13 @@ def test_pack_coefficients(): _F = form(F) # -- Test vector - b0 = assemble_vector(_F) + b0 = petsc_assemble_vector(_F) b0.assemble() constants = _cpp.fem.pack_constants(_F._cpp_object) coeffs = _cpp.fem.pack_coefficients(_F._cpp_object) with b0.localForm() as _b0: for c in [(None, None), (None, coeffs), (constants, None), (constants, coeffs)]: - b = assemble_vector(_F, c[0], c[1]) + b = petsc_assemble_vector(_F, c[0], c[1]) b.assemble() with b.localForm() as _b: assert (_b0.array_r == _b.array_r).all() @@ -822,7 +826,7 @@ def test_pack_coefficients(): coeff *= 5.0 with b0.localForm() as _b0: for c in [(None, coeffs), (constants, None), (constants, coeffs)]: - b = assemble_vector(_F, c[0], c[1]) + b = petsc_assemble_vector(_F, c[0], c[1]) b.assemble() with b.localForm() as _b: assert (_b0 - _b).norm() > 1.0e-5 @@ -832,13 +836,13 @@ def test_pack_coefficients(): J = ufl.derivative(F, u, du) J = form(J) - A0 = assemble_matrix(J) + A0 = petsc_assemble_matrix(J) A0.assemble() constants = _cpp.fem.pack_constants(J._cpp_object) coeffs = _cpp.fem.pack_coefficients(J._cpp_object) for c in [(None, None), (None, coeffs), (constants, None), (constants, coeffs)]: - A = assemble_matrix(J, constants=c[0], coeffs=c[1]) + A = petsc_assemble_matrix(J, constants=c[0], coeffs=c[1]) A.assemble() assert pytest.approx((A - A0).norm(), 1.0e-12) == 0.0 @@ -847,7 +851,7 @@ def test_pack_coefficients(): for coeff in coeffs.values(): coeff *= 5.0 for c in [(None, coeffs), (constants, None), (constants, coeffs)]: - A = assemble_matrix(J, constants=c[0], coeffs=c[1]) + A = petsc_assemble_matrix(J, constants=c[0], coeffs=c[1]) A.assemble() assert (A - A0).norm() > 1.0e-5 @@ -867,13 +871,13 @@ def test_coefficents_non_constant(): # -- Volume integral vector F = form((ufl.inner(u, v) - ufl.inner(x[0] * x[1]**2, v)) * dx) - b0 = assemble_vector(F) + b0 = petsc_assemble_vector(F) b0.assemble() assert np.linalg.norm(b0.array) == pytest.approx(0.0, abs=1.0e-7) # -- Exterior facet integral vector F = form((ufl.inner(u, v) - ufl.inner(x[0] * x[1]**2, v)) * ds) - b0 = assemble_vector(F) + b0 = petsc_assemble_vector(F) b0.assemble() assert np.linalg.norm(b0.array) == pytest.approx(0.0, abs=1.0e-7) @@ -890,7 +894,7 @@ def test_coefficents_non_constant(): F = (ufl.inner(u1('+') * u0('-'), ufl.avg(v)) - ufl.inner(x[0] * x[1]**2, ufl.avg(v))) * ufl.dS F = form(F) - b0 = assemble_vector(F) + b0 = petsc_assemble_vector(F) b0.assemble() assert np.linalg.norm(b0.array) == pytest.approx(0.0, abs=1.0e-7) @@ -973,9 +977,9 @@ def partitioner(comm, nparts, local_graph, num_ghost_nodes): assert sum == pytest.approx(6.0) # Assemble - A = assemble_matrix(a) + A = petsc_assemble_matrix(a) A.assemble() - b = assemble_vector(L) + b = petsc_assemble_vector(L) b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # Solve @@ -1001,16 +1005,16 @@ def test_matrix_assembly_rectangular(mode): def single(): a = form(ufl.inner(u, v0) * ufl.dx) - A = assemble_matrix(a, bcs=[]) + A = petsc_assemble_matrix(a, bcs=[]) A.assemble() return A def block(): a = form([[ufl.inner(u, v0) * ufl.dx], [ufl.inner(u, v1) * ufl.dx]]) - A0 = assemble_matrix_block(a, bcs=[]) + A0 = petsc_assemble_matrix_block(a, bcs=[]) A0.assemble() - A1 = assemble_matrix_nest(a, bcs=[]) + A1 = petsc_assemble_matrix_nest(a, bcs=[]) A1.assemble() return A0, A1 diff --git a/python/test/unit/io/test_xdmf_function.py b/python/test/unit/io/test_xdmf_function.py index fe160dedfce..ec2b679f4a4 100644 --- a/python/test/unit/io/test_xdmf_function.py +++ b/python/test/unit/io/test_xdmf_function.py @@ -48,7 +48,7 @@ def test_save_1d_scalar(tempdir, encoding, dtype, use_pathlib): mesh = create_unit_interval(MPI.COMM_WORLD, 32, dtype=xtype) V = FunctionSpace(mesh, ("Lagrange", 2)) u = Function(V, dtype=dtype) - u.x.set(1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0) with XDMFFile(mesh.comm, filename2, "w", encoding=encoding) as file: file.write_mesh(mesh) file.write_function(u) @@ -63,7 +63,7 @@ def test_save_2d_scalar(tempdir, encoding, dtype, cell_type): mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, cell_type, dtype=xtype) V = FunctionSpace(mesh, ("Lagrange", 2)) u = Function(V, dtype=dtype) - u.x.set(1.0) + u.x.array[:] = 1.0 with XDMFFile(mesh.comm, filename, "w", encoding=encoding) as file: file.write_mesh(mesh) file.write_function(u) @@ -78,7 +78,7 @@ def test_save_3d_scalar(tempdir, encoding, dtype, cell_type): mesh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 4, cell_type, dtype=xtype) V = FunctionSpace(mesh, ("Lagrange", 2)) u = Function(V, dtype=dtype) - u.x.set(1.0) + u.x.array[:] = 1.0 with XDMFFile(mesh.comm, filename, "w", encoding=encoding) as file: file.write_mesh(mesh) file.write_function(u) @@ -93,7 +93,7 @@ def test_save_2d_vector(tempdir, encoding, dtype, cell_type): mesh = create_unit_square(MPI.COMM_WORLD, 12, 13, cell_type, dtype=xtype) V = VectorFunctionSpace(mesh, ("Lagrange", 2)) u = Function(V, dtype=dtype) - u.x.set(1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0) with XDMFFile(mesh.comm, filename, "w", encoding=encoding) as file: file.write_mesh(mesh) file.write_function(u) @@ -107,7 +107,7 @@ def test_save_3d_vector(tempdir, encoding, dtype, cell_type): filename = Path(tempdir, "u_3Dv.xdmf") mesh = create_unit_cube(MPI.COMM_WORLD, 2, 2, 2, cell_type, dtype=xtype) u = Function(VectorFunctionSpace(mesh, ("Lagrange", 1)), dtype=dtype) - u.x.set(1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0) with XDMFFile(mesh.comm, filename, "w", encoding=encoding) as file: file.write_mesh(mesh) file.write_function(u) @@ -121,7 +121,7 @@ def test_save_2d_tensor(tempdir, encoding, dtype, cell_type): filename = Path(tempdir, "tensor.xdmf") mesh = create_unit_square(MPI.COMM_WORLD, 16, 16, cell_type, dtype=xtype) u = Function(TensorFunctionSpace(mesh, ("Lagrange", 2)), dtype=dtype) - u.x.set(1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0) with XDMFFile(mesh.comm, filename, "w", encoding=encoding) as file: file.write_mesh(mesh) file.write_function(u) @@ -135,7 +135,7 @@ def test_save_3d_tensor(tempdir, encoding, dtype, cell_type): filename = Path(tempdir, "u3t.xdmf") mesh = create_unit_cube(MPI.COMM_WORLD, 4, 4, 4, cell_type, dtype=xtype) u = Function(TensorFunctionSpace(mesh, ("Lagrange", 2)), dtype=dtype) - u.x.set(1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0) with XDMFFile(mesh.comm, filename, "w", encoding=encoding) as file: file.write_mesh(mesh) file.write_function(u) @@ -151,11 +151,11 @@ def test_save_3d_vector_series(tempdir, encoding, dtype, cell_type): u = Function(VectorFunctionSpace(mesh, ("Lagrange", 2)), dtype=dtype) with XDMFFile(mesh.comm, filename, "w", encoding=encoding) as file: file.write_mesh(mesh) - u.x.set(1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 1.0 + (1j if np.issubdtype(dtype, np.complexfloating) else 0) file.write_function(u, 0.1) - u.x.set(2.0 + (2j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 2.0 + (2j if np.issubdtype(dtype, np.complexfloating) else 0) file.write_function(u, 0.2) with XDMFFile(mesh.comm, filename, "a", encoding=encoding) as file: - u.x.set(3.0 + (3j if np.issubdtype(dtype, np.complexfloating) else 0)) + u.x.array[:] = 3.0 + (3j if np.issubdtype(dtype, np.complexfloating) else 0) file.write_function(u, 0.3) diff --git a/python/test/unit/la/test_nullspace.py b/python/test/unit/la/test_nullspace.py index 9f1dce9a9c9..e35e861d704 100644 --- a/python/test/unit/la/test_nullspace.py +++ b/python/test/unit/la/test_nullspace.py @@ -9,17 +9,17 @@ import numpy as np import pytest - import ufl -from dolfinx import la +from dolfinx.la import create_petsc_vector from dolfinx.fem import VectorFunctionSpace, form from dolfinx.fem.petsc import assemble_matrix from dolfinx.mesh import (CellType, GhostMode, create_box, create_unit_cube, create_unit_square) -from ufl import TestFunction, TrialFunction, dx, grad, inner - from mpi4py import MPI from petsc4py import PETSc +from ufl import TestFunction, TrialFunction, dx, grad, inner + +from dolfinx import la def build_elastic_nullspace(V): @@ -67,7 +67,7 @@ def build_broken_elastic_nullspace(V): """Function to build incorrect null space for 2D elasticity""" # Create list of vectors for null space - ns = [la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs) for i in range(4)] + ns = [create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs) for i in range(4)] with ExitStack() as stack: vec_local = [stack.enter_context(x.localForm()) for x in ns] diff --git a/python/test/unit/nls/test_newton.py b/python/test/unit/nls/test_newton.py index 7b150cf4322..4b3b8c3700e 100644 --- a/python/test/unit/nls/test_newton.py +++ b/python/test/unit/nls/test_newton.py @@ -11,13 +11,14 @@ locate_dofs_geometrical) from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, create_matrix, create_vector, set_bc) +from dolfinx.la import create_petsc_vector from dolfinx.mesh import create_unit_square from mpi4py import MPI from petsc4py import PETSc from ufl import TestFunction, TrialFunction, derivative, dx, grad, inner from dolfinx import cpp as _cpp -from dolfinx import default_real_type, la +from dolfinx import default_real_type class NonlinearPDEProblem: @@ -171,7 +172,7 @@ def test_nonlinear_pde_snes(): problem = NonlinearPDE_SNESProblem(F, u, bc) u.x.array[:] = 0.9 - b = la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs) + b = create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs) J = create_matrix(problem.a) # Create Newton solver and solve