Hey all,
I am solving a linear system (A*U=F) and need to gather information about the A matrix to select an appropriate iterative solver for my simulation. Specifically, I want to know about the sparsity pattern of the matrix, including:
- the number of zeros in the matrix,
- their locations, and
- the condition number of the matrix.
Below is a minimal code example for elasticity where Matrix A is the one I want to analyze. However, I’m not sure how to extract this information. What is the best way to do that?
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
import dolfinx
# -------------------------Material Properties
E = 1 * 130e9
nu = 0.28
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# ......................................... Reading mesh file
h = 0.5
mesh= create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
[15, 15, 15], cell_type=CellType.hexahedron)
# ......................................... Define Loading
Force = 200.0
Area2 = 0.5**2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,)) # Displacement Element
Space = functionspace(mesh, el)
# ......................................... Define Trail and Test Functions
u = TrialFunction(Space)
del_u = TestFunction(Space)
# ....................................... Mark boundary regions
tol = 1 * 1e-16
def Bottom(x):
return np.isclose(x[2], 0, atol=tol)
# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], h, atol=tol))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
# # ...................................... BCs
# Clamp BC for left side
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
clamp_dofs = locate_dofs_topological(Space, mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(np.array([0, 0, 0], dtype=default_scalar_type), clamp_dofs, Space)
bc = [bc0] # BC 4
# ......................................... Define tensor
delta = Identity(3)
i, j, k, l = indices(4)
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])
# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l)) # Elasticity tensor
def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])
# ......................................... Define Variational Formulation
VarForm = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
a = VarForm
L = tr_S2[i] * del_u[i] * ds(2)
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = assemble_vector(linear_form)
apply_lifting(b, [bilinear_form], [bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bc)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
problem = LinearProblem(a, L, bcs=bc,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()
with VTXWriter(mesh.comm, ("Pyramid/DIS_Pyramid.bp"), [Result]) as vtx:
vtx.write(0.0)
I found a previous post here that seemed helpful, but I encountered the follwing error when I ran it. I am using DOLFINx version 0.8.0.
