# 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 :)

import math


First load python packages we need for this example

import time

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

from ansys.mapdl.core import examples, launch_mapdl


Next:

• Get the Math module of PyMapdl

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

Product:             Ansys Mechanical Enterprise
MAPDL Version:       22.2
ansys.mapdl Version: 0.64.0


APDLMath EigenSolve First load the input file using MAPDL.

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

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

***** MAPDL 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() 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()

*****  MAPDL 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)

Elapsed time to solve this problem :  0.8896512985229492


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}")

 : Freq =   352.39 Hz         Residual = 3.4723e-08
 : Freq =   385.16 Hz         Residual = 8.1511e-09
 : Freq =   656.78 Hz         Residual = 1.303e-08
 : Freq =   764.72 Hz         Residual = 1.3847e-08
 : Freq =   825.33 Hz         Residual = 1.1408e-08
 : Freq =  1039.21 Hz         Residual = 4.2396e-09
 : Freq =  1143.66 Hz         Residual = 1.1628e-08
 : Freq =  1258.00 Hz         Residual = 5.3364e-08
 : Freq =  1334.23 Hz         Residual = 5.1172e-08
 : Freq =  1352.06 Hz         Residual = 3.7675e-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)

Elapsed time to solve this problem :  4.250144720077515


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}")

 : Freq =   352.39 Hz         Residual = 8.0098e-05
 : Freq =   385.16 Hz         Residual = 0.00010369
 : Freq =   656.78 Hz         Residual = 0.00024262
 : Freq =   764.72 Hz         Residual = 0.00016258
 : Freq =   825.33 Hz         Residual = 0.00039002
 : Freq =  1039.21 Hz         Residual = 0.00057628
 : Freq =  1143.66 Hz         Residual = 0.0025916
 : Freq =  1258.00 Hz         Residual = 0.00033872
 : Freq =  1334.23 Hz         Residual = 0.00046633
 : Freq =  1352.06 Hz         Residual = 0.0011262


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")

Mapdl is 4.78 times faster


stop mapdl

mapdl.exit()


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

Gallery generated by Sphinx-Gallery