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.

5

Hi, dokken. I have meet the same problem that I convert the array to the Mat of seqaij and convert the Mat back to array, it becomes different from the original array. Exactly the order of its elements seems to have changed. Iā€™d like to know what this is all about.

6

Without a minimal reproducible example, I am not sure I can give any further guidance.