Scaling of condition number in 1D, 2D, and 3D

Linear Algebra
, ,
5

To further supply on top of @Stein’s code and comments, you can also extract a sub-matrix from your PETSc matrix, and use that as input to your eigenvalue problem.

I here illustrate this with standard PETSc KSP, but that is just to show the procedure:

# # Extract a sub matrix without boundary condition dofs and use it with PETSc
# Author: Jørgen S. Dokken
# SPDX License identifier: MIT

from mpi4py import MPI
from petsc4py import PETSc
import dolfinx.fem.petsc
import ufl
import numpy as np
M = 30
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, M, M)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u, v= ufl.TrialFunction(V), ufl.TestFunction(V)

a = dolfinx.fem.form(ufl.inner(u, v) * ufl.dx)

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
exterior_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim-1, exterior_facets)

bcs = [dolfinx.fem.dirichletbc(0., exterior_dofs, V)]
A = dolfinx.fem.petsc.assemble_matrix(a, bcs, diag=0.0)
A.assemble()

num_dofs_local = V.dofmap.index_map.size_local*V.dofmap.index_map_bs
num_ghosts = V.dofmap.index_map.num_ghosts*V.dofmap.index_map_bs
not_dirichlet = np.full(num_dofs_local+num_ghosts, 1, dtype=np.int32)
not_dirichlet[exterior_dofs] = 0

idx_set = np.flatnonzero(not_dirichlet).astype(np.int32)
idx_set_row = idx_set[idx_set<num_dofs_local]
idx_set_col = idx_set[idx_set<num_dofs_local]
isx_row = PETSc.IS(A.getComm()).createGeneral(idx_set_row)
isx_col = PETSc.IS(A.getComm()).createGeneral(idx_set_col)
lgm_row = A.getLGMap()[0].applyIS(isx_row)
lgm_col = A.getLGMap()[1].applyIS(isx_col)



A_sub = A.createSubMatrix(lgm_row, lgm_col)
A_sub.assemble()

assert np.isclose(A.norm(0), A_sub.norm(0))
assert np.isclose(A.norm(2), A_sub.norm(2))
assert np.isclose(A.norm(3), A_sub.norm(3))

x = ufl.SpatialCoordinate(mesh)
L = dolfinx.fem.form(ufl.inner(x[0]+2*x[1], v)*ufl.dx)
b = dolfinx.fem.petsc.assemble_vector(L)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

b_sub = b.getSubVector(lgm_row)

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A_sub)
ksp.setType("preonly")
ksp.setErrorIfNotConverged(True)
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.getPC().getFactorMatrix().setMumpsIcntl(14, 1000)

x= b.duplicate()
x_sub = x.getSubVector(lgm_row)

ksp.solve(b_sub, x_sub)

x.restoreSubVector(lgm_row, x_sub)


uh = dolfinx.fem.Function(V)
x.copy(uh.x.petsc_vec)
uh.x.scatter_forward()
with dolfinx.io.VTXWriter(uh.function_space.mesh.comm, "u_sub.bp", [uh]) as bp:
    bp.write(0.0)
2 Likes
Structural Modal Analysis