|
| 1 | +""" |
| 2 | +.. _ref_mapdl_math_mapdl_vs_scipy: |
| 3 | +
|
| 4 | +Compute Eigenvalues using MAPDL or SciPy |
| 5 | +---------------------------------------- |
| 6 | +
|
| 7 | +This example shows: |
| 8 | +
|
| 9 | +- How to extract the stiffness and mass matrices from a MAPDL model. |
| 10 | +- How to use the ``Math`` module of PyMapdl to compute the first |
| 11 | + eigenvalues. |
| 12 | +- How to can get these matrices in the SciPy world, to get the same |
| 13 | + solutions using Python resources. |
| 14 | +- How MAPDL is really faster than SciPy :) |
| 15 | +""" |
| 16 | + |
| 17 | +############################################################################### |
| 18 | +# First load python packages we need for this example |
| 19 | +import time |
| 20 | +import math |
| 21 | + |
| 22 | +import matplotlib.pylab as plt |
| 23 | +import numpy as np |
| 24 | +import scipy |
| 25 | +from scipy.sparse.linalg import eigsh |
| 26 | + |
| 27 | +from ansys.mapdl.core import launch_mapdl |
| 28 | +from ansys.mapdl.core import examples |
| 29 | + |
| 30 | +############################################################################### |
| 31 | +# Next: |
| 32 | +# |
| 33 | +# - Load the ansys.mapdl module |
| 34 | +# - Get the ``Math`` module of PyMapdl |
| 35 | +# |
| 36 | +mapdl = launch_mapdl() |
| 37 | +print(mapdl) |
| 38 | +mm = mapdl.math |
| 39 | + |
| 40 | +############################################################################### |
| 41 | +# APDLMath EigenSolve |
| 42 | +# First load the input file using MAPDL. |
| 43 | +# |
| 44 | +print(mapdl.input(examples.examples.wing_model)) |
| 45 | + |
| 46 | + |
| 47 | +############################################################################### |
| 48 | +# Plot and mesh using the ``eplot`` method. |
| 49 | +mapdl.eplot() |
| 50 | + |
| 51 | + |
| 52 | +############################################################################### |
| 53 | +# Next, setup a modal Analysis and request the :math:`K` and math:`M` |
| 54 | +# matrices to be formed. MAPDL stores these matrices in a ``.FULL`` |
| 55 | +# file. |
| 56 | + |
| 57 | +print(mapdl.slashsolu()) |
| 58 | +print(mapdl.antype(antype='MODAL')) |
| 59 | +print(mapdl.modopt(method='LANB', nmode='10', freqb='1.')) |
| 60 | +print(mapdl.wrfull(ldstep='1')) |
| 61 | + |
| 62 | +# store the output of the solve command |
| 63 | +output = mapdl.solve() |
| 64 | + |
| 65 | + |
| 66 | +############################################################################### |
| 67 | +# Read the sparse matrices using PyMapdl. |
| 68 | +# |
| 69 | +mapdl.finish() |
| 70 | +mm.free() |
| 71 | +k = mm.stiff(fname='file.full') |
| 72 | +M = mm.mass(fname='file.full') |
| 73 | + |
| 74 | + |
| 75 | +############################################################################### |
| 76 | +# Solve the eigenproblem using PyMapdl with APDLMath. |
| 77 | +# |
| 78 | +nev = 10 |
| 79 | +A = mm.mat(k.nrow, nev) |
| 80 | + |
| 81 | +t1 = time.time() |
| 82 | +ev = mm.eigs(nev, k, M, phi=A, fmin=1.) |
| 83 | +t2 = time.time() |
| 84 | +mapdl_elapsed_time = t2-t1 |
| 85 | +print('\nElapsed time to solve this problem : ', mapdl_elapsed_time) |
| 86 | + |
| 87 | +############################################################################### |
| 88 | +# Print eigenfrequencies and accuracy. |
| 89 | +# |
| 90 | +# Accuracy : :math:`\frac{||(K-\lambda.M).\phi||_2}{||K.\phi||_2}` |
| 91 | +# |
| 92 | +mapdl_acc = np.empty(nev) |
| 93 | + |
| 94 | +for i in range(nev): |
| 95 | + f = ev[i] # Eigenfrequency (Hz) |
| 96 | + omega = 2*np.pi*f # omega = 2.pi.Frequency |
| 97 | + lam = omega**2 # lambda = omega^2 |
| 98 | + |
| 99 | + phi = A[i] # i-th eigenshape |
| 100 | + kphi = k.dot(phi) # K.Phi |
| 101 | + mphi = M.dot(phi) # M.Phi |
| 102 | + |
| 103 | + kphi_nrm = kphi.norm() # Normalization scalar value |
| 104 | + |
| 105 | + mphi *= lam # (K-\lambda.M).Phi |
| 106 | + kphi -= mphi |
| 107 | + |
| 108 | + mapdl_acc[i] = kphi.norm()/kphi_nrm # compute the residual |
| 109 | + print(f"[{i}] : Freq = {f:8.2f} Hz\t Residual = {mapdl_acc[i]:.5}") |
| 110 | + |
| 111 | + |
| 112 | +############################################################################### |
| 113 | +# Use SciPy to Solve the same Eigenproblem |
| 114 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 115 | +# |
| 116 | +# First get MAPDL sparse matrices into the Python memory as SciPy |
| 117 | +# matrices. |
| 118 | +# |
| 119 | +pk = k.asarray() |
| 120 | +pm = M.asarray() |
| 121 | + |
| 122 | +# get_ipython().run_line_magic('matplotlib', 'inline') |
| 123 | + |
| 124 | +fig, (ax1, ax2) = plt.subplots(1, 2) |
| 125 | +fig.suptitle('K and M Matrix profiles') |
| 126 | +ax1.spy(pk, markersize=0.01) |
| 127 | +ax1.set_title('K Matrix') |
| 128 | +ax2.spy(pm, markersize=0.01) |
| 129 | +ax2.set_title('M Matrix') |
| 130 | +plt.show(block=True) |
| 131 | + |
| 132 | + |
| 133 | +############################################################################### |
| 134 | +# Make the sparse matrices for SciPy unsymmetric as symmetric matrices in SciPy |
| 135 | +# are memory inefficient. |
| 136 | +# |
| 137 | +# :math:`K = K + K^T - diag(K)` |
| 138 | +# |
| 139 | +pkd = scipy.sparse.diags(pk.diagonal()) |
| 140 | +pK = pk + pk.transpose() - pkd |
| 141 | +pmd = scipy.sparse.diags(pm.diagonal()) |
| 142 | +pm = pm + pm.transpose() - pmd |
| 143 | + |
| 144 | + |
| 145 | +############################################################################### |
| 146 | +# Plot Matrices |
| 147 | +# |
| 148 | +fig, (ax1, ax2) = plt.subplots(1, 2) |
| 149 | +fig.suptitle('K and M Matrix profiles') |
| 150 | +ax1.spy(pk, markersize=0.01) |
| 151 | +ax1.set_title('K Matrix') |
| 152 | +ax2.spy(pm, markersize=0.01) |
| 153 | +ax2.set_title('M Matrix') |
| 154 | +plt.show(block=True) |
| 155 | + |
| 156 | + |
| 157 | +############################################################################### |
| 158 | +# Solve the eigenproblem |
| 159 | +# |
| 160 | +t3 = time.time() |
| 161 | +vals, vecs = eigsh(A=pK, M=pm, k=10, sigma=1, which='LA') |
| 162 | +t4 = time.time() |
| 163 | +scipy_elapsed_time = t4 - t3 |
| 164 | +print('\nElapsed time to solve this problem : ', scipy_elapsed_time) |
| 165 | + |
| 166 | + |
| 167 | +############################################################################### |
| 168 | +# Convert Lambda values to Frequency values: |
| 169 | +# :math:`freq = \frac{\sqrt(\lambda)}{2.\pi}` |
| 170 | +# |
| 171 | +freqs = np.sqrt(vals) / (2*math.pi) |
| 172 | + |
| 173 | + |
| 174 | +############################################################################### |
| 175 | +# Compute the residual error for SciPy. |
| 176 | +# |
| 177 | +# :math:`Err=\frac{||(K-\lambda.M).\phi||_2}{||K.\phi||_2}` |
| 178 | +# |
| 179 | +scipy_acc = np.zeros(nev) |
| 180 | + |
| 181 | +for i in range(nev): |
| 182 | + lam = vals[i] # i-th eigenvalue |
| 183 | + phi = vecs.T[i] # i-th eigenshape |
| 184 | + |
| 185 | + kphi = pk*phi.T # K.Phi |
| 186 | + mphi = pm*phi.T # M.Phi |
| 187 | + |
| 188 | + kphi_nrm = np.linalg.norm(kphi, 2) # Normalization scalar value |
| 189 | + |
| 190 | + mphi *= lam # (K-\lambda.M).Phi |
| 191 | + kphi -= mphi |
| 192 | + |
| 193 | + scipy_acc[i] = 1 - np.linalg.norm(kphi, 2)/kphi_nrm # compute the residual |
| 194 | + print(f"[{i}] : Freq = {freqs[i]:8.2f} Hz\t Residual = {scipy_acc[i]:.5}") |
| 195 | + |
| 196 | + |
| 197 | +############################################################################### |
| 198 | +# MAPDL is more accurate than SciPy. |
| 199 | +# |
| 200 | +fig = plt.figure(figsize=(12, 10)) |
| 201 | +ax = plt.axes() |
| 202 | +x = np.linspace(1, 10, 10) |
| 203 | +plt.title('Residual Error') |
| 204 | +plt.yscale('log') |
| 205 | +plt.xlabel('Mode') |
| 206 | +plt.ylabel('% Error') |
| 207 | +ax.bar(x, scipy_acc, label='SciPy Results') |
| 208 | +ax.bar(x, mapdl_acc, label='MAPDL Results') |
| 209 | +plt.legend(loc='lower right') |
| 210 | +plt.show() |
| 211 | + |
| 212 | +############################################################################### |
| 213 | +# MAPDL is faster than SciPy. |
| 214 | +# |
| 215 | +ratio = scipy_elapsed_time/mapdl_elapsed_time |
| 216 | +print(f"Mapdl is {ratio:.3} times faster") |
0 commit comments