# 1D Oscillator with External Eigensolver¶

Here, we will analyze the oscillations of a square plate by performing an eigenvalue analysis in different configurations using an external eigensolver (Python script):

Here you can download the python script, mesh input file, final mesh and material file.

## Setup¶

We consider a square mass of dimension 1mx1m, which is constrained in $x$ direction on the left and is connected to a spring and a damper in the lower left corner:

The system constitutes a 1D oscillator in the $y$-direction, described by the following ordinary differential equation:

$$m\ddot{y}+c\dot{y}+ky=f(t) \label{eq:momentumbalance}$$

From the normalized form of \eqref{eq:momentumbalance} $$\ddot{y} + 2\xi\omega_0 \dot{y}+\omega_0^2 y=\frac{f(t)}{m}$$

with $\omega_0 = \sqrt{\frac{k}{m}}$ and $\xi = \frac{c}{2\omega_0m}$, we can find the following analytical expression for the eigenvalues: $$\lambda_{1,2} = -\xi\omega_0 \pm \omega_0 \sqrt{\xi^2-1}$$ Note that the eigenvalue for the undamped case ($\xi=0$) is simply $\omega_0^2$.

The corresponding eigenvalue problem in FE-formulation takes the form: $$(\boldsymbol{K}+\lambda\boldsymbol{C}+\lambda^2\boldsymbol{M})\boldsymbol{x}=\boldsymbol{0} \label{eq:quadraticEVP}$$ For the undamped case, \eqref{eq:quadraticEVP} simplifies to a generalized eigenvalue problem: $$\boldsymbol{K}\boldsymbol{x}=\lambda\boldsymbol{M}\boldsymbol{x}$$ The quadratic eigenvalue problem is also converted to a generalized one, similar to the transformation of second order ODEs to first order ODEs. The resulting generalized eigenvalue problem with increased system size is solved with the external eigensolver.

## Meshing using Gmsh¶

We use Gmsh as a preprocessor. ASCII *.geo input files are used to create geometry and mesh. In the following, the workflow for creating the input file for the square mass is illustrated:

• Define geometrical constants and choose element type.

• Create the square plate by first defining the corner points, connecting them with lines and then generating a closed surface.

• Define a structured mesh; in this case 1 element for the entire domain is sufficient, since we are not interested in the mechanical deformation of the plate.

• Define physical groups for relevant nodes, edges and the surface to define boundary conditions and the computational domain in the simulation.

• Export the mesh.

• The .msh file can be directly generated from the UnitSquare.geo input file. Note that only Gmsh version 2 is supported by openCFS, so we need to export the mesh in the legacy format by the following command:

gmsh -3 -format msh2 UnitSquare.geo

• Final mesh:

## External Eigensolver¶

As already mentioned, we use a Python script to solve the generalized eigenvalue problem for our oscillator:

• Import the required libraries for linear algebra functions and the sys module for taking arguments from the terminal:
from scipy.io import mmread, mmwrite
import scipy.sparse
import numpy as np
import sys

• Read in the arguments from the terminal call:
AFileName = sys.argv[1]
BFileName = sys.argv[2]
arg_sigma = complex(sys.argv[3])
arg_k = int(sys.argv[4])
arg_tol = sys.argv[5]

• Import the system matrices from .mtx format:
A = mmread(AFileName)

• Solve the generalized eigenvalue problem with the scipy.linalg.sparse package for sparse linear algebra:
eigenvalues, eigenvectors = scipy.sparse.linalg.eigs(A,k=arg_k,M=B,sigma=arg_sigma,tol=arg_tol)

• Reshape the arrays holding eigenvalues and eigenvectors since the openCFS input routine only supports vector-like structures:
eigenvalues = eigenvalues[np.newaxis].T
eigenvectors = eigenvectors.reshape((-1,1), order='F')

• Export the resulting arrays to the Matrix Market coordinate format:
eigenvalues = scipy.sparse.coo_matrix(eigenvalues)
eigenvectors = scipy.sparse.coo_matrix(eigenvectors)
mmwrite(sys.argv[6], eigenvalues)
mmwrite(sys.argv[7], eigenvectors)


## Undamped Oscillator¶

We start by setting up the XML file for the undamped case with $m=1kg$ and $k=1 \frac{N}{m}$. Thus, we expect an eigenvalue of $\lambda=1 \frac{1}{s^2}$.

• Define the .msh file as input. The results are stored in a *.hdf5 file. The material parameters are taken from the mat.xml file.
  <fileFormats>
<input>
<gmsh fileName="UnitSquare.msh"/>
</input>
<output>
<hdf5/>
</output>
<materialData file="mat.xml" format="xml"/>
</fileFormats>

• Define the computational domain as plane and assign the material:
  <domain geometryType='plane'  printGridInfo="no">
<regionList>
<region name="V_square" material="UnitMaterial" />
</regionList>
</domain>

• Define the eigenvalue analysis step: we have a generalized eigenvalue problem and specify a shift point around which the solver should look for the eigenvalues. We only expect one eigenvalue, so we set the number of computed eigenvalues to 1.
 <sequenceStep>
<sequenceStep>
<analysis>
<eigenValue>
<valuesAround>
<shiftPoint>
<Real>0.95</Real>
<Imag>0</Imag>
</shiftPoint>
<number>1</number>
</valuesAround>
<eigenVectors normalization="norm" side="right"/>
<problemType>
<Generalized>
<aMatrix>stiffness</aMatrix>
<bMatrix>mass</bMatrix>
</Generalized>
</problemType>
</eigenValue>
</analysis>


• Define plane strain solid mechanics as the governing partial differential equation (PDE) on our computational domain:
 <pdeList>
<mechanic subType='planeStrain'>
<regionList>
<region name="V_square"/>
</regionList>

• Constrain the plate in $x$ direction on the left side. Define a spring element acting in $y$-direction at the lower left node with the stiffness of $k=1 \frac{N}{m}$.
         <bcsAndLoads>
<fix name="S_fixed">
<comp dof="x"/>
</fix>
<concentratedElem name="N_concentrated-element" dof="y" stiffnessValue="1.0"/>

• Define the stored results:
        <storeResults>
<nodeResult type="mechDisplacement">
<allRegions/>
</nodeResult>
</storeResults>
</mechanic>
</pdeList>

• Choose the external eigensolver for our problem and indicate what parameters should be included in the terminal call:
  <linearSystems>
<system>
<solutionStrategy>
<standard>
<eigenSolver id="ext"/>
</standard>
</solutionStrategy>
<eigenSolverList>
<external id="ext">
<logging>yes</logging>
<cmd>python3 EigenSolver.py</cmd>
<arguments>
<AFileName/>
<BFileName/>
<shiftPoint/>
<number/>
<tolerance>1e-9</tolerance>
<EigenValuesFileName/>
<EigenVectorsFileName/>
</arguments>
</external>
</eigenSolverList>
</system>
</linearSystems>
</sequenceStep>

</cfsSimulation>

• The material data is found in the mat.xml file:
<?xml version='1.0' encoding='utf-8'?>
<cfsMaterialDataBase xmlns="http://www.cfs++.org/material">
<material name="UnitMaterial">
<mechanical>
<density>
<linear>
<real>1.0</real>
</linear>
</density>
<elasticity>
<linear>
<isotropic>
<elasticityModulus>
<real>1.0e+9</real>
</elasticityModulus>
<poissonNumber>
<real>0.25</real>
</poissonNumber>
</isotropic>
</linear>
</elasticity>
</mechanical>
</material>
</cfsMaterialDataBase>

• Running the simulation by the command cfs UndampedOscillator, gives us the following results:
++ Reading mesh ... "UnitSquare.msh" -> regular
++ Creating PDE 'mechanic' for analysis 'eigenValue'
++ External Eigensolver - Executing the command 'python3 EigenSolver.py UndampedOscillator_AMatrix.mtx UndampedOscillator_BMatrix.mtx 0.950000+0.000000j 1 1.000000e-09 UndampedOscillator_EigenValues.mtx UndampedOscillator_EigenVectors.mtx'
Import Olivia
Import Olivia

Mode |               Real Part |            Error Bounds |          Imaginary Part |
1 |                       1 |                   1e-09 |             2.23709e-18 |


The call to the Python script is printed with all arguments and can be used for debugging.

## Damped Oscillator¶

We set up the XML file for the damped case with the same mass and stiffness as before ($m=1kg$ and $k=1 \frac{N}{m}$) and introduce a damping $c=\frac{6}{5}\frac{Ns}{m}$ . This gives us $\lambda_{1,2}= -0.6 \pm 0.8j \frac{1}{s}$ for the analytical eigenvalues. Only few changes are necessary to achieve this setup starting from the undamped case:

• We define a quadratic eigenvalue problem and compute 2 instead of one eigenvalues around a chosen shift point:
<analysis>
<eigenValue>
<valuesAround>
<shiftPoint>
<Real>-0.6</Real>
<Imag>0</Imag>
</shiftPoint>
<number>2</number>
</valuesAround>
<eigenVectors normalization="norm" side="right"/>
<problemType>
<linear>damping</linear>
<constant>stiffness</constant>
</problemType>
</eigenValue>
</analysis>

• In addition to the stiffness, we define a damping at the lower left node
 <concentratedElem name="N_concentrated-element" dof="y" stiffnessValue="1.0" dampingValue="3/5*2"/>

• Define the solvers for the quadratic eigenvalue problem: the system is converted to a generalized eigenvalue problem (quadratic eigensolver), which is then solved with the external eigensolver:
<linearSystems>
<system>
<solutionStrategy>
<standard>
</standard>
</solutionStrategy>
<eigenSolverList>
<generalisedEigenSolver id="ext"/>
<linearisation>firstCompanion</linearisation>
<external id="ext">
<logging>yes</logging>
<cmd>python3 EigenSolver.py</cmd>
<arguments>
<AFileName/>
<BFileName/>
<shiftPoint/>
<number/>
<tolerance>1e-9</tolerance>
<EigenValuesFileName/>
<EigenVectorsFileName/>
</arguments>
</external>
</eigenSolverList>
</system>
</linearSystems>

• Running the simulation by the command cfs DampedOscillator, gives us the following results:
++ Reading mesh ... "UnitSquare.msh" -> regular
++ Creating PDE 'mechanic' for analysis 'eigenValue'

Time needed for the linearisation: 0.00065 sec

++ External Eigensolver - Executing the command 'python3 EigenSolver.py DampedOscillator_AMatrix.mtx DampedOscillator_BMatrix.mtx -0.600000+0.000000j 2 1.000000e-09 DampedOscillator_EigenValues.mtx DampedOscillator_EigenVectors.mtx'
Import Olivia
Import Olivia

Mode |               Real Part |            Error Bounds |          Imaginary Part |
1 |                    -0.6 |                   1e-09 |                    -0.8 |
2 |                    -0.6 |                   1e-09 |                     0.8 |