.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/ansys-math_vs_scipy.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_ansys-math_vs_scipy.py: .. _ref_ansys_math_vs_scipy: Compute Eigenvalues using PyAnsys Math and SciPy ------------------------------------------------ This example shows how to perform these tasks: - Extract the stiffness and mass matrices from an MAPDL model. - Use PyAnsys Math to compute the first eigenvalues. - Get these matrices using SciPy to obtain the same solutions using Python resources. - See if PyAnsys Math is more accurate and faster than SciPy. .. GENERATED FROM PYTHON SOURCE LINES 39-42 Perform required imports and start PyAnsys Math ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Perform required imports. .. GENERATED FROM PYTHON SOURCE LINES 42-57 .. code-block:: Python import math import time from ansys.mapdl.core.examples import examples import matplotlib.pylab as plt import numpy as np import scipy from scipy.sparse.linalg import eigsh import ansys.math.core.math as pymath # Start PyAnsys Math as a service. mm = pymath.AnsMath() .. GENERATED FROM PYTHON SOURCE LINES 58-61 Load the input file ~~~~~~~~~~~~~~~~~~~ Load the input file using MAPDL. .. GENERATED FROM PYTHON SOURCE LINES 61-65 .. code-block:: Python print(mm._mapdl.input(examples.wing_model)) .. rst-class:: sphx-glr-script-out .. code-block:: none /INPUT FILE= wing.dat LINE= 0 *** WARNING *** CP = 0.966 TIME= 15:01:41 The /BATCH command must be the first line of input. The /BATCH command is ignored. *** MAPDL - ENGINEERING ANALYSIS SYSTEM RELEASE 24.1 *** Ansys Mechanical Enterprise Academic Student 01055371 VERSION=LINUX x64 15:01:41 APR 22, 2024 CP= 0.967 ***** MAPDL ANALYSIS DEFINITION (PREP7) ***** *** WARNING *** CP = 0.967 TIME= 15:01:41 Deactivation of element shape checking is not recommended. *** WARNING *** CP = 1.004 TIME= 15:01:41 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 = 1.015 TIME= 15:01:41 CLEAR, SELECT, and MESH boundary condition commands are not possible after MODMSH. ***** ROUTINE COMPLETED ***** CP = 1.019 .. GENERATED FROM PYTHON SOURCE LINES 66-69 Plot and mesh ~~~~~~~~~~~~~ Plot and mesh using the ``eplot`` method. .. GENERATED FROM PYTHON SOURCE LINES 69-73 .. code-block:: Python mm._mapdl.eplot() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /examples/images/sphx_glr_ansys-math_vs_scipy_001.png :alt: ansys math vs scipy :srcset: /examples/images/sphx_glr_ansys-math_vs_scipy_001.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/runner/work/pyansys-math/pyansys-math/doc/source/examples/images/sphx_glr_ansys-math_vs_scipy_001.vtksz .. GENERATED FROM PYTHON SOURCE LINES 74-78 Set up modal analysis ~~~~~~~~~~~~~~~~~~~~~ Set up a modal analysis and form the :math:`K` and :math:`M` matrices. MAPDL stores these matrices in a ``.FULL`` file. .. GENERATED FROM PYTHON SOURCE LINES 78-88 .. code-block:: Python print(mm._mapdl.slashsolu()) print(mm._mapdl.antype(antype="MODAL")) print(mm._mapdl.modopt(method="LANB", nmode="10", freqb="1.")) print(mm._mapdl.wrfull(ldstep="1")) # store the output of the solve command output = mm._mapdl.solve() .. rst-class:: sphx-glr-script-out .. code-block:: none ***** 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 .. GENERATED FROM PYTHON SOURCE LINES 89-92 Read sparse matrices ~~~~~~~~~~~~~~~~~~~~ Read the sparse matrices using PyAnsys Math. .. GENERATED FROM PYTHON SOURCE LINES 92-99 .. code-block:: Python mm._mapdl.finish() mm.free() k = mm.stiff(fname="file.full") M = mm.mass(fname="file.full") .. GENERATED FROM PYTHON SOURCE LINES 100-104 Solve eigenproblem ~~~~~~~~~~~~~~~~~~ Solve the eigenproblme using PyAnsys Math. .. GENERATED FROM PYTHON SOURCE LINES 104-113 .. code-block:: Python 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() pymath_elapsed_time = t2 - t1 print("\nElapsed time to solve this problem : ", pymath_elapsed_time) .. rst-class:: sphx-glr-script-out .. code-block:: none Elapsed time to solve this problem : 0.4790050983428955 .. GENERATED FROM PYTHON SOURCE LINES 114-120 Print eigenfrequencies and accuracy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Print the eigenfrequencies and the accuracy. Accuracy : :math:`\frac{||(K-\lambda.M).\phi||_2}{||K.\phi||_2}` .. GENERATED FROM PYTHON SOURCE LINES 120-140 .. code-block:: Python pymath_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 pymath_acc[i] = kphi.norm() / kphi_nrm # compute the residual print(f"[{i}] : Freq = {f:8.2f} Hz\t Residual = {pymath_acc[i]:.5}") .. rst-class:: sphx-glr-script-out .. code-block:: none [0] : Freq = 352.39 Hz Residual = 7.1398e-09 [1] : Freq = 385.20 Hz Residual = 4.4352e-09 [2] : Freq = 656.78 Hz Residual = 7.6434e-09 [3] : Freq = 764.73 Hz Residual = 1.1473e-08 [4] : Freq = 825.49 Hz Residual = 1.4281e-08 [5] : Freq = 1039.25 Hz Residual = 1.8512e-08 [6] : Freq = 1143.60 Hz Residual = 1.4537e-08 [7] : Freq = 1258.00 Hz Residual = 6.6247e-08 [8] : Freq = 1334.23 Hz Residual = 2.154e-07 [9] : Freq = 1352.10 Hz Residual = 6.4539e-08 .. GENERATED FROM PYTHON SOURCE LINES 141-145 Use SciPy to solve the same eigenproblem ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Get the MAPDL sparse matrices into Python memory as SciPy matrices. .. GENERATED FROM PYTHON SOURCE LINES 145-160 .. code-block:: Python 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) .. image-sg:: /examples/images/sphx_glr_ansys-math_vs_scipy_002.png :alt: K and M Matrix profiles, K Matrix, M Matrix :srcset: /examples/images/sphx_glr_ansys-math_vs_scipy_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 161-166 Make the sparse matrices for SciPy unsymmetric because symmetric matrices in SciPy are memory inefficient. :math:`K = K + K^T - diag(K)` .. GENERATED FROM PYTHON SOURCE LINES 166-172 .. code-block:: Python pkd = scipy.sparse.diags(pk.diagonal()) pK = pk + pk.transpose() - pkd pmd = scipy.sparse.diags(pm.diagonal()) pm = pm + pm.transpose() - pmd .. GENERATED FROM PYTHON SOURCE LINES 173-175 Plot the matrices. .. GENERATED FROM PYTHON SOURCE LINES 175-184 .. code-block:: Python 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) .. image-sg:: /examples/images/sphx_glr_ansys-math_vs_scipy_003.png :alt: K and M Matrix profiles, K Matrix, M Matrix :srcset: /examples/images/sphx_glr_ansys-math_vs_scipy_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 185-187 Solve the eigenproblem. .. GENERATED FROM PYTHON SOURCE LINES 187-194 .. code-block:: Python 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) .. rst-class:: sphx-glr-script-out .. code-block:: none Elapsed time to solve this problem : 2.037623405456543 .. GENERATED FROM PYTHON SOURCE LINES 195-198 Convert lambda values to frequency values: :math:`freq = \frac{\sqrt(\lambda)}{2.\pi}` .. GENERATED FROM PYTHON SOURCE LINES 198-201 .. code-block:: Python freqs = np.sqrt(vals) / (2 * math.pi) .. GENERATED FROM PYTHON SOURCE LINES 202-206 Compute the residual error for SciPy. :math:`Err=\frac{||(K-\lambda.M).\phi||_2}{||K.\phi||_2}` .. GENERATED FROM PYTHON SOURCE LINES 206-224 .. code-block:: Python 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}") .. rst-class:: sphx-glr-script-out .. code-block:: none [0] : Freq = 352.39 Hz Residual = 8.0111e-05 [1] : Freq = 385.20 Hz Residual = 0.00010356 [2] : Freq = 656.78 Hz Residual = 0.00024284 [3] : Freq = 764.73 Hz Residual = 0.00016261 [4] : Freq = 825.49 Hz Residual = 0.00038862 [5] : Freq = 1039.25 Hz Residual = 0.00057672 [6] : Freq = 1143.60 Hz Residual = 0.0025674 [7] : Freq = 1258.00 Hz Residual = 0.00033886 [8] : Freq = 1334.23 Hz Residual = 0.00046688 [9] : Freq = 1352.10 Hz Residual = 0.0011237 .. GENERATED FROM PYTHON SOURCE LINES 225-228 See if PyAnsys Math is more accurate than SciPy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Plot residual error to see if PyAnsys Math is more accurate than SciPy. .. GENERATED FROM PYTHON SOURCE LINES 228-241 .. code-block:: Python 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, pymath_acc, label="PyAnsys Math Results") plt.legend(loc="lower right") plt.show() .. image-sg:: /examples/images/sphx_glr_ansys-math_vs_scipy_004.png :alt: Residual Error :srcset: /examples/images/sphx_glr_ansys-math_vs_scipy_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 242-245 See if PyAnsys Math is faster than SciPy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Plot elapsed time to see if PyAnsys Math is more accurate than SciPy. .. GENERATED FROM PYTHON SOURCE LINES 245-249 .. code-block:: Python ratio = scipy_elapsed_time / pymath_elapsed_time print(f"PyAnsys Math is {ratio:.3} times faster.") .. rst-class:: sphx-glr-script-out .. code-block:: none PyAnsys Math is 4.25 times faster. .. GENERATED FROM PYTHON SOURCE LINES 250-253 Stop PyAnsys Math ~~~~~~~~~~~~~~~~~ Stop PyAnsys Math. .. GENERATED FROM PYTHON SOURCE LINES 253-255 .. code-block:: Python mm._mapdl.exit() .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/ansys/mapdl/core/launcher.py:811: UserWarning: The environment variable 'PYMAPDL_START_INSTANCE' is set, hence the argument 'start_instance' is overwritten. warnings.warn( .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 7.563 seconds) .. _sphx_glr_download_examples_ansys-math_vs_scipy.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ansys-math_vs_scipy.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ansys-math_vs_scipy.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_