Conversion from array to PETSc and from PETSc to array

1

I am using this code for converting array to PETSc

def array2petsc4py(g):
    Xpt = PETSc.Mat().createAIJ(g.shape)
    Xpt.setUp()
    Xpt.setValues(range(0, g.shape[0]), range(0, g.shape[1]), g) 
    Xpt.assemble()
    return Xpt

I am using this code for converting PETSc to array

def petsc2array(v):
    s=v.getValues(range(0, v.getSize()[0]), range(0,  v.getSize()[1]))
    return s

The entire code is as under:

from dolfin import *
import numpy as np

L, B, H = 20., 0.5, 1.

Nx = 200
Ny = int(B/L*Nx)+1
Nz = int(H/L*Nx)+1

mesh = BoxMesh(Point(0.,0.,0.),Point(L,B,H), Nx, Ny, Nz)
E, nu = Constant(1e5), Constant(0.)
rho = Constant(1e-3)

# Lame coefficient for constitutive relation
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
u_ = TrialFunction(V)
du = TestFunction(V)


def left(x, on_boundary):
    return near(x[0],0.)

bc = DirichletBC(V, Constant((0.,0.,0.)), left)
k_form = inner(sigma(du),eps(u_))*dx
l_form = Constant(1.)*u_[0]*dx
K = PETScMatrix()
b = PETScVector()
assemble_system(k_form, l_form, bc, A_tensor=K, b_tensor=b)

m_form = rho*dot(du,u_)*dx
M = PETScMatrix()
assemble(m_form, tensor=M)
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.
N_eig = 6   # number of eigenvalues
print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(N_eig)

# Exact solution computation
from scipy.optimize import root
from math import cos, cosh
falpha = lambda x: cos(x)*cosh(x)+1
alpha = lambda n: root(falpha, (2*n+1)*pi/2.)['x'][0]

# Set up file for exporting results
file_results = XDMFFile("modal_analysis.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

# Extraction
eig_v = []
eig_l = []
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)

    # 3D eigenfrequency
    freq_3D = sqrt(r)/2/pi

    # Beam eigenfrequency
    if i % 2 == 0: # exact solution should correspond to weak axis bending
        I_bend = H*B**3/12.
    else:          #exact solution should correspond to strong axis bending
        I_bend = B*H**3/12.
    freq_beam = alpha(i/2)**2*sqrt(E*I_bend/(rho*B*H*L**4))/2/pi

    #print("Solid FE: {0:8.5f} [Hz]   Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
    print(freq_3D)
    
    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector "+str(i))
    eigenmode.vector()[:] = rx
    file_results.write(eigenmode, 0)
    eig_l.append(freq_3D)
    eig_v.append(rx)

X=np.array(eig_v).T
print(X)
X=array2petsc4py(X)
print(petsc2array(X))

When I am printing (eig_v).T, I am getting this result

[[-5.74603728e-09 -5.56168924e-08  1.39166771e-06  1.34571155e-05
   2.99130317e-05 -5.24574139e-05]
 [-1.46828847e-09  4.44368075e-10  4.13287260e-07 -1.15319026e-07
   1.00016719e-05 -2.05460734e-04]
 [ 1.12387171e-10  6.62216770e-09 -2.67054929e-08 -2.22954680e-06
  -5.71791541e-07 -9.43435280e-05]
 ...
 [ 3.50790481e-01  6.82910454e-01  1.21546520e+00  2.34954088e+00
  -1.98597854e+00 -2.75347126e-01]
 [-2.00053800e+01  1.97297796e-01 -2.00011482e+01  1.87148446e-01
   1.99670139e+01  2.22417356e+01]
 [ 1.95350971e-01  1.99766043e+01  1.79867559e-01  1.98325977e+01
  -1.67631860e-01  9.82987713e+00]]

When I am printing petsc2array(X), I am getting this:

array([[-5.74603728e-09, -1.46828847e-09,  1.12387171e-10,
        -7.20729164e-09, -1.85249140e-09,  1.37110354e-10],
       [-2.93643085e-03, -8.65676446e-04,  5.92329939e-05,
        -4.28157489e-03, -8.63416199e-04,  8.01453964e-05],
       [-3.35274371e-09, -1.22287034e-09,  6.93487603e-11,
        -7.22509339e-09, -1.88445084e-09,  1.31201876e-10],
       ...,
       [-1.84894696e-01,  2.22416506e+01, -8.43399827e+00,
        -2.03052265e-01,  2.22414448e+01, -4.78124768e+00],
       [-2.21061076e-01,  2.22413487e+01, -1.12857016e+00,
        -2.39041436e-01,  2.22414310e+01,  2.52412793e+00],
       [-2.57113443e-01,  2.22416240e+01,  6.17693519e+00,
        -2.75347126e-01,  2.22417356e+01,  9.82987713e+00]])

But both of the results have to be matched, but here, they are not matching.

2

I would strongly suggest making a minimal code example, as described in Read before posting: How do I get my question answered?
as your question is only about converting a numpy array <-> PETScMatrix, which does not have to include an eigenvalue problem.

In general, you can convert it using the dense array:


from petsc4py import PETSc

def petsc2array(v):
    s=v.getValues(range(0, v.getSize()[0]), range(0,  v.getSize()[1]))
    return s


print(X)
Xpetsc = PETSc.Mat().createDense(X.shape, array=X)

X_back = petsc2array(Xpetsc)
assert(np.allclose(X_back, X))
2 Likes
3

Sir, for a matrix of around 100000*100000 order, while converting array to petsc I am facing this error.

image (53)
image (53)765×294 96.6 KB

4

You should never use dense matrices when you have systems of that size. Then you need to use CSR matrices, as explained in: Converting to SciPy sparse matrix without Eigen backend - #2 by MiroK

If you have a 100 000 by 100 000 means 10 000 000 000 nonzeros.