# Compute Eigenvalues using MAPDL or SciPy¶

This example shows:

• How to extract the stiffness and mass matrices from a MAPDL model.

• How to use the Math module of PyMapdl to compute the first eigenvalues.

• How to can get these matrices in the SciPy world, to get the same solutions using Python resources.

• How MAPDL is really faster than SciPy :)

First load python packages we need for this example

import time
import math

import matplotlib.pylab as plt
import numpy as np
import scipy
from scipy.sparse.linalg import eigsh

from ansys.mapdl.core import launch_mapdl
from ansys.mapdl.core import examples


Next:

• Get the Math module of PyMapdl

mapdl = launch_mapdl()
print(mapdl)
mm = mapdl.math


Out:

Product:             Ansys Mechanical Enterprise
MAPDL Version:       21.2
ansys.mapdl Version: 0.60.5


APDLMath EigenSolve First load the input file using MAPDL.

print(mapdl.input(examples.examples.wing_model))


Out:

/INPUT FILE=    LINE=       0
*****ANSYS VERIFICATION RUN ONLY*****
DO NOT USE RESULTS FOR PRODUCTION

***** ANSYS ANALYSIS DEFINITION (PREP7) *****

*** WARNING ***                         CP =       0.000   TIME= 00:00:00
Deactivation of element shape checking is not recommended.

*** WARNING ***                         CP =       0.000   TIME= 00:00:00
The mesh of area 1 contains PLANE42 triangles, which are much too stiff
in bending.  Use quadratic (6- or 8-node) elements if possible.

*** WARNING ***                         CP =       0.000   TIME= 00:00:00
CLEAR, SELECT, and MESH boundary condition commands are not possible
after MODMSH.

***** ROUTINE COMPLETED *****  CP =         0.000


Plot and mesh using the eplot method.

mapdl.eplot()


Out:

/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/pyvista/core/dataset.py:1401: PyvistaDeprecationWarning: Use of point_arrays is deprecated. Use point_data instead.
warnings.warn(
/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/pyvista/core/dataset.py:1541: PyvistaDeprecationWarning: Use of cell_arrays is deprecated. Use cell_data instead.
warnings.warn(

[(3.933687241492801, 3.0105415922331633, 3.21223669211588),
(1.1499194488984512, 0.22677379963881444, 0.4284688995215311),
(0.0, 0.0, 1.0)]


Next, setup a modal Analysis and request the $$K$$ and math:M matrices to be formed. MAPDL stores these matrices in a .FULL file.

print(mapdl.slashsolu())
print(mapdl.antype(antype="MODAL"))
print(mapdl.modopt(method="LANB", nmode="10", freqb="1."))
print(mapdl.wrfull(ldstep="1"))

# store the output of the solve command
output = mapdl.solve()


Out:

*****  ANSYS SOLUTION ROUTINE  *****
PERFORM A MODAL ANALYSIS
THIS WILL BE A NEW ANALYSIS
USE SYM. BLOCK LANCZOS MODE EXTRACTION METHOD
EXTRACT    10 MODES
SHIFT POINT FOR EIGENVALUE CALCULATION=  1.0000
NORMALIZE THE MODE SHAPES TO THE MASS MATRIX
STOP SOLUTION AFTER FULL FILE HAS BEEN WRITTEN
LOADSTEP =    1 SUBSTEP =    1 EQ. ITER =    1


Read the sparse matrices using PyMapdl.

mapdl.finish()
mm.free()
k = mm.stiff(fname="file.full")
M = mm.mass(fname="file.full")


Solve the eigenproblem using PyMapdl with APDLMath.

nev = 10
A = mm.mat(k.nrow, nev)

t1 = time.time()
ev = mm.eigs(nev, k, M, phi=A, fmin=1.0)
t2 = time.time()
mapdl_elapsed_time = t2 - t1
print("\nElapsed time to solve this problem : ", mapdl_elapsed_time)


Out:

Elapsed time to solve this problem :  0.48459291458129883


Print eigenfrequencies and accuracy.

Accuracy : $$\frac{||(K-\lambda.M).\phi||_2}{||K.\phi||_2}$$

mapdl_acc = np.empty(nev)

for i in range(nev):
f = ev[i]  # Eigenfrequency (Hz)
omega = 2 * np.pi * f  # omega = 2.pi.Frequency
lam = omega ** 2  # lambda = omega^2

phi = A[i]  # i-th eigenshape
kphi = k.dot(phi)  # K.Phi
mphi = M.dot(phi)  # M.Phi

kphi_nrm = kphi.norm()  # Normalization scalar value

mphi *= lam  # (K-\lambda.M).Phi
kphi -= mphi

mapdl_acc[i] = kphi.norm() / kphi_nrm  # compute the residual
print(f"[{i}] : Freq = {f:8.2f} Hz\t Residual = {mapdl_acc[i]:.5}")


Out:

[0] : Freq =   352.40 Hz         Residual = 6.4363e-09
[1] : Freq =   385.20 Hz         Residual = 5.6669e-09
[2] : Freq =   656.81 Hz         Residual = 6.2898e-09
[3] : Freq =   764.73 Hz         Residual = 6.6152e-09
[4] : Freq =   825.47 Hz         Residual = 5.4891e-09
[5] : Freq =  1039.28 Hz         Residual = 6.2336e-09
[6] : Freq =  1143.61 Hz         Residual = 2.6386e-08
[7] : Freq =  1258.00 Hz         Residual = 4.5162e-08
[8] : Freq =  1334.23 Hz         Residual = 5.2739e-08
[9] : Freq =  1352.10 Hz         Residual = 3.6164e-08


## Use SciPy to Solve the same Eigenproblem¶

First get MAPDL sparse matrices into the Python memory as SciPy matrices.

pk = k.asarray()
pm = M.asarray()

# get_ipython().run_line_magic('matplotlib', 'inline')

fig, (ax1, ax2) = plt.subplots(1, 2)
fig.suptitle("K and M Matrix profiles")
ax1.spy(pk, markersize=0.01)
ax1.set_title("K Matrix")
ax2.spy(pm, markersize=0.01)
ax2.set_title("M Matrix")
plt.show(block=True)


Make the sparse matrices for SciPy unsymmetric as symmetric matrices in SciPy are memory inefficient.

$$K = K + K^T - diag(K)$$

pkd = scipy.sparse.diags(pk.diagonal())
pK = pk + pk.transpose() - pkd
pmd = scipy.sparse.diags(pm.diagonal())
pm = pm + pm.transpose() - pmd


Plot Matrices

fig, (ax1, ax2) = plt.subplots(1, 2)
fig.suptitle("K and M Matrix profiles")
ax1.spy(pk, markersize=0.01)
ax1.set_title("K Matrix")
ax2.spy(pm, markersize=0.01)
ax2.set_title("M Matrix")
plt.show(block=True)


Solve the eigenproblem

t3 = time.time()
vals, vecs = eigsh(A=pK, M=pm, k=10, sigma=1, which="LA")
t4 = time.time()
scipy_elapsed_time = t4 - t3
print("\nElapsed time to solve this problem : ", scipy_elapsed_time)


Out:

Elapsed time to solve this problem :  2.5159525871276855


Convert Lambda values to Frequency values: $$freq = \frac{\sqrt(\lambda)}{2.\pi}$$

freqs = np.sqrt(vals) / (2 * math.pi)


Compute the residual error for SciPy.

$$Err=\frac{||(K-\lambda.M).\phi||_2}{||K.\phi||_2}$$

scipy_acc = np.zeros(nev)

for i in range(nev):
lam = vals[i]  # i-th eigenvalue
phi = vecs.T[i]  # i-th eigenshape

kphi = pk * phi.T  # K.Phi
mphi = pm * phi.T  # M.Phi

kphi_nrm = np.linalg.norm(kphi, 2)  # Normalization scalar value

mphi *= lam  # (K-\lambda.M).Phi
kphi -= mphi

scipy_acc[i] = 1 - np.linalg.norm(kphi, 2) / kphi_nrm  # compute the residual
print(f"[{i}] : Freq = {freqs[i]:8.2f} Hz\t Residual = {scipy_acc[i]:.5}")


Out:

[0] : Freq =   352.40 Hz         Residual = 8.0071e-05
[1] : Freq =   385.20 Hz         Residual = 0.00010352
[2] : Freq =   656.81 Hz         Residual = 0.00024255
[3] : Freq =   764.73 Hz         Residual = 0.00016257
[4] : Freq =   825.47 Hz         Residual = 0.0003896
[5] : Freq =  1039.28 Hz         Residual = 0.00057568
[6] : Freq =  1143.61 Hz         Residual = 0.0025878
[7] : Freq =  1258.00 Hz         Residual = 0.00033875
[8] : Freq =  1334.24 Hz         Residual = 0.00046618
[9] : Freq =  1352.10 Hz         Residual = 0.0011268


MAPDL is more accurate than SciPy.

fig = plt.figure(figsize=(12, 10))
ax = plt.axes()
x = np.linspace(1, 10, 10)
plt.title("Residual Error")
plt.yscale("log")
plt.xlabel("Mode")
plt.ylabel("% Error")
ax.bar(x, scipy_acc, label="SciPy Results")
ax.bar(x, mapdl_acc, label="MAPDL Results")
plt.legend(loc="lower right")
plt.show()


MAPDL is faster than SciPy.

ratio = scipy_elapsed_time / mapdl_elapsed_time
print(f"Mapdl is {ratio:.3} times faster")


Out:

Mapdl is 5.19 times faster


Total running time of the script: ( 0 minutes 5.352 seconds)

Gallery generated by Sphinx-Gallery