waveguide.py

class waveguicsx.waveguide.Waveguide(comm: _MPI.Comm, M: Mat, K0: Mat, K1: Mat, K2: Mat)

Bases: object

A class for solving complex waveguide problems based on SLEPc eigensolver.

The full documentation is entirely defined in the `waveguide.py’ module.

The following matrix problem is considered: (K0-omega**2*M + 1j*k*(K1-K1^T) + k**2*K2)*U=F. This kind of problem typically stems from the so-called SAFE (Semi-Analytical Finite Element) method.

The class enables to deal with complex waveguides, two-dimensional (e.g. plates) or three-dimensional (arbitrarily shaped cross-section), inhomogeneous in the transverse directions, anisotropic. Complex-valued problems can be handled including the effects of non-propagating modes (evanescent, inhomogeneous), viscoelastic loss (complex material properties) or perfectly matched layers (PML) to simulate buried waveguides.

The free response (F=0) is an eigenvalue problem, solved iteratively by varying the parameter which can be the angular frequency omega or the wavenumber k. In the former case, the eigenvalue is k, while in the latter case, the eigenvalue is omega^2. The loops over the parameter (angular frequency or wavenumber) can be parallelized, as shown in some tutorials (using mpi4py).

Various modal properties (energy velocity, group velocity, excitability…) can be post-processed as a function of the frequency and plotted as dispersion curves.

The forced reponse (F is not 0) is solved in the frequency domain by expanding the solution as a sum of eigenmodes using biorthogonality relationship, leading to very fast computations of excited wavefields.

Example:

# In this example, the matrices M, K0, K1, K2 and the excitation vector F are supposed to be dimensional for simplicity
# Yet, in practice, the problem would better be normalized to avoid ill-conditioning (see tutorials)

from waveguicsx.waveguide import Waveguide

# Definition of the excitation signal (here, a toneburst)
excitation = Signal()
excitation.toneburst(fs=400e3, T=2e-3, fc=100e3, n=8) #central frequency 100 kHz, 8 cycles, duration 2 ms, sampling frequency 400 kHz
excitation.plot() #plot time signal
excitation.plot_spectrum() #plot spectrum
omega = 2*np.pi*excitation.frequency #angular frequency range

# Initialization of waveguide
wg = Waveguide(MPI.COMM_WORLD, M, K0, K1, K2)
wg.set_parameters(omega=omega) #set the parameter range (here, angular frequency)

# Free response (dispersion curves)
wg.solve(nev=20, target=0) #solution of eigenvalue problem (iteration over the parameter omega), 20 eigenvalues requested at each frequency
wg.compute_energy_velocity() #post-process energy velocity
wg.plot() #plot k vs. omega
wg.plot_energy_velocity() #plot ve vs. omega

# Computation of modal coefficients due to an excitation vector F
wg.compute_response_coefficient(F=F) #F should be a PETSc vector
wg.plot_coefficient() #plot modal coefficients vs. omega

# Forced response at degree of freedom dof and axial coordinates z (dof should be an integer)
frequency, response = wg.compute_response(dof=dof, z=[0.5, 1., 1.5, 2.], spectrum=excitation.spectrum, plot=False) #response in the frequency domain
response = Signal(frequency=frequency, spectrum=response) #define response as a Signal object
response.plot_spectrum() #plot frequency response
response.ifft() #response in the time domain
response.plot() #plot time response
plt.show()

Attributes

commmpi4py.MPI.Intracomm

MPI communicator (parallel processing)

M, K0, K1, K2petsc4py.PETSc.Mat

SAFE matrices

Fpetsc4py.PETSc.Vec

SAFE excitation vector

problem_typestr

problem_type is “omega” if the varying parameter is omega, “wavenumber” if this is k

two_sidedbool

if True, left eigenvectors will be also computed (otherwise, only right eigenvectors are computed)

target: complex number or user-defined function of the parameter

target around which eigenpairs are looked for (see method solve)

omega or wavenumbernumpy.ndarray

the parameter range specified by the user (see method set_parameters)

evpPEP or EPS instance (SLEPc object)

eigensolver parameters

eigenvalueslist of numpy arrays

list of wavenumbers or angular frequencies, access to components with eigenvalues[ip][imode] (ip: parameter index, imode: mode index)

eigenvectorslist of PETSc matrices

list of mode shapes, access to components with eigenvectors[ik][idof,imode] (ik: parameter index, imode: mode index, idof: dof index) or eigenvectors[ik].getColumnVector(imode)

eigenforceslist of PETSc matrices

list of eigenforces (acces to components: see eigenvectors)

opposite_goinglist of numpy arrays

list of opposite-going mode (acces to components: see eigenvectors)

energy_velocitylist of numpy arrays

list of energy velocity (access to component: see eigenvalues)

group_velocitylist of numpy arrays

list of group velocity (access to component: see eigenvalues)

traveling_directionlist of numpy arrays

list of traveling_direction (access to component: see eigenvalues)

pml_ratiolist of numpy arrays

list of pml ratio, used for filtering out PML modes (access to component: see eigenvalues)

coefficientlist of numpy arrays

list of response coefficient to excitation vector F (access to component: see eigenvalues)

excitabilitylist of numpy arrays

list of excitability to excitation vector F (access to component: see eigenvalues)

complex_powerlist of numpy arrays

list of complex power flow of individual modes (access to component: see eigenvalues)

plot_scalerdictionnary

dictionnary containing the scaling factors of various modal properties, useful to plot results in a dimensional form

Methods

__init__(comm:’_MPI.Comm’, M:PETSc.Mat, K0:PETSc.Mat, K1:PETSc.Mat, K2:PETSc.Mat):

Constructor, initialization of waveguide

set_parameters(omega=None, wavenumber=None, two_sided=False):

Set problem type (problem_type), the parameter range (omega or wavenumber) as well as default parameters of SLEPc eigensolver (evp); set two_sided to True to compute left eigenvectors also (left eigenvectors are the opposite-going modes)

solve(nev=1, target=0):

Solve the eigenvalue problem repeatedly for the parameter range, solutions are stored as attributes (names: eigenvalues, eigenvectors)

compute_eigenforce():

Compute the eigenforces for the whole parameter range and store them as an attribute (name: eigenforces)

compute_poynting_normalization():

Normalization of eigenvectors and eigenforces, so that U’=U/sqrt(|P|), where P is the normal component of complex Poynting vector

compute_opposite_going(plot=False):

Compute opposite-going mode pairs based on on wavenumber and biorthogonality for the whole parameter range and store them as attributes (name: opposite_going), set plot to True to visualize the biorthogonality values of detected pairs

compute_energy_velocity():

Compute the energy velocities for the whole parameter range and store them as an attribute (name: energy_velocity)

compute_group_velocity():

Compute the group velocities for the whole parameter range and store them as an attribute (name: energy_velocity)

compute_traveling_direction():

Compute the traveling directions for the whole parameter range and store them as an attribute (name: traveling_direction)

compute_pml_ratio():

Compute the pml ratios for the whole parameter range and store them as an attribute (name: pml_ratio)

compute_response_coefficient(F, spectrum=None, wavenumber_function=None, dof=None):

Compute the response coefficients due to excitation vector F for the whole parameter range and store them as an attribute (name: coefficient)

compute_complex_power(self):

Compute the individual complex power flow of modes given by P=-1j*omega/2*U^H*F

compute_response(dof, z, spectrum=None, wavenumber_function=None, plot=False):

Compute the response at the degree of freedom dof and the axial coordinate z for the whole frequency range

track_mode(omega_index, mode_index, threshold=0.9, plot=False):

Track a mode over the whole frequency range thanks to eigenvector similarity. The mode is specified by its index, mode_index, at a given angular frequency index, omega_index.

plot(direction=None, pml_threshold=None, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot dispersion curves Re(omega) vs. Re(wavenumber) using matplotlib

plot_phase_velocity(direction=None, pml_threshold=None, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot phase velocity dispersion curves, vp=Re(omega)/Re(wavenumber) vs. Re(omega)

plot_attenuation(direction=None, pml_threshold=None, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot attenuation dispersion curves, Im(wavenumber) vs. Re(omega) if omega is the parameter, or Im(omega) vs. Re(omega) if wavenumber is the parameter

plot_energy_velocity(direction=None, pml_threshold=None, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot energy velocity dispersion curves, ve vs. Re(omega)

plot_group_velocity(direction=None, pml_threshold=None, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot group velocity dispersion curves, vg vs. Re(omega)

plot_coefficient(direction=None, pml_threshold=None, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot response coefficients as a function of frequency, |q| vs. Re(omega)

plot_excitability(direction=None, pml_threshold=None, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot excitability as a function of frequency, |e| vs. Re(omega)

plot_spectrum(index=0, ax=None, color=”k”, marker=”o”, markersize=2, linestyle=””, **kwargs):

Plot the spectrum, Im(eigenvalues) vs. Re(eigenvalues), for the parameter index specified by the user

set_plot_scaler(length=1, time=1, mass=1, dim=3):

Define the characteristic length, time, mass, as well as dim, and calculate the scaling factors of modal properties, which are stored in the attribute name plot_scaler (useful to visualize plots in a dimensional form)

set_parameters(omega: Optional[ndarray[complex]] = None, wavenumber: Optional[ndarray[complex]] = None, two_sided: bool = False)

Set the parameter range (omega or wavenumber) as well as default parameters of the SLEPc eigensolver (evp). The user must specify the parameter omega or wavenumber, but not both. This method generates the attributes omega (or wavenumber) and evp. After calling this method, various SLEPc parameters can be set by changing the attribute evp manually. Set two_sided=True for solving left eigenvectors also.

Parameters

omega or wavenumbernumpy.ndarray

the parameter range specified by the user

two_sidedbool

False if left eigenvectiors are not needed, True if they must be solved also

solve(nev=1, target=0)

Solve the dispersion problem, i.e. the eigenvalue problem repeatedly for the parameter range (omega or wavenumber). The solutions are stored in the attributes eigenvalues and eigenvectors. If two_sided is True, left eigensolutions are also solved.

Note: left eigensolutions correspond to opposite-going modes and are hence added to the right eigensolutions (i.e. in eigenvalues and eigenvectors) after removing any possible duplicates.

Parameters

nevint

number of eigenpairs requested

targetcomplex number or user-defined function of the parameter, optional (default: 0)

target around which eigenpairs are looked for a small shift might sometimes prevent errors (e.g. zero pivot with dirichlet bc)

compute_eigenforces()

Post-process the eigenforces F=(K1^T+1j*k*K2)*U for every mode in the whole parameter range

compute_poynting_normalization()

Post-process the normalization of eigenvectors and eigenforces, so that U’=U/sqrt(|P|), where P is the normal component of complex Poynting vector (P=-1j*omega/2*U^H*F). After normalization, every mode is such that |P|=1 and the attribute _poynting_normalization is set to True. Normalization is not mandatory but, when applied, has to be done before any response coefficient computation.

compute_energy_velocity()

Post-process the energy velocity ve=Re(P)/Re(E) for every mode in the whole parameter range, where P is the normal component of complex Poynting vector and E is the total energy (cross-section time-averaged). Warning in case of PML: the integration is currently applied over the whole cross-section (including PML), the so-defined energy velocity is questionable.

compute_opposite_going(plot=False)

Post-process the pairing of opposite-going modes, based on wavenumber and biorthogonality criteria, and store them as a an attribute (name: opposite_going, -1 value for unpaired modes). Compute their biorthogonality normalization factors, Um^T*F-m - U-m^T*Fm, where m and -m denote opposite-going modes, for the whole parameter range and store them as an attribute (name: _biorthogonality_factor). If plot is set to True, the biorthogonality criterion found by the algorithm is plotted as a function of frequency index, allowing visual check that there is no values close to zero (a factor close to zero probably means a lack of biorthogonality). The biorthogonality criterion is defined as |biorthogonality_factor*omega/4|.

Notes:

• when an unpaired mode is found, the value -1 is stored in opposite_going (and NaN value in _biorthogonality_factor), meaning that this mode will be discarded in the computation of group velocity, traveling direction, coefficient and excitability (NaN values stored)

• if modes with lack of biorthogonality or two many unpaired modes occur, try to recompute the eigenproblem by increasing the accuracy (e.g. reducing the tolerance)

• lack of biorthogonality may be also due to multiple modes (*); in this case, try to use an unstructured mesh instead

• if two_sided is True, lack of biorthogonolity may occur for specific target: try another target (e.g. add a small imaginary part)

(*) e.g. flexural modes in a cylinder with structured mesh

compute_group_velocity()

Post-process the group velocity, vg=1/Re(dk/domega) for every mode in the whole parameter range (opposite-going modes required). For unpaired modes, NaN values are set.

compute_traveling_direction(delta=0.01)

Post-process the traveling direction, +1 or -1, for every mode in the whole parameter range, using the sign of Im(k + 1j*delta/v) where delta is the imaginary shift used for analytical continuation of k, and v is the group velocity (or, if not available, the energy velocity). This criterion is based on the limiting absorption principle (theoretically, vg should be used instead of ve). For unpaired modes, NaN values are set.

compute_pml_ratio()

Post-process the pml ratio (useful to filter out PML mode), given by 1-Im(Ek)/|Ek| where Ek denotes the “complex” kinetic energy, for every mode in the whole parameter range. Reminder: the pml ratio tends to 1 for mode shapes vanishing inside the PML.

compute_response_coefficient(F, spectrum=None, wavenumber_function=None, dof=None)

Computation of modal coefficients due to the excitation vector F for every mode in the whole omega range (opposite-going eigenvectors are required). Modal coefficients qm are defined from: U(z,omega) = sum qm(omega)*Um(omega)*exp(i*km*z), m=1…M, omega denotes the angular frequency. For unpaired modes, NaN values are set. Assumption: the source is centred at z=0.

Note: spectrum and wavenumber_function can be specified in compute_response(…) instead of compute_response_coefficient(…), but not in both functions in the same time (otherwise the excitation will be modulated twice)

Parameters

FPETSc vector

SAFE excitation vector

spectrumnumpy.ndarray

when specified, spectrum is a vector of length omega used to modulate F in terms of frequency (default: 1 for all frequencies)

wavenumber_function: python function

when specified, wavenumber_function is a python function used to modulate F in terms of wavenumber (example: wavenumber_function = lambda x: np.sin(x), default: 1 for all wavenumbers, i.e. source localized at z=0)

dofint

when specified, it calculates the modal excitability (stored in the attribute excitability), i.e. qm*Um at the degree of freedom dof and for a unit excitation vector (i.e. such that the sum of the elements of F is equal to 1)

compute_complex_power()

Post-process the individual complex power flow of modes given by P=-1j*omega/2*U^H*F (normal component of complex Poynting vector), where U and F denote the eigenvector and eigenforce of a single mode. The ‘usual’ power is given by the real part.

Notes:

• For lossy media (e.g. with viscoelasticity or with PML), it should be reminded that the individual power should be carefully handled: the power of a sum is not equal to the sum of powers because Auld’s complex biorthogonality relation does no longer hold

• This inequality also generally applies between two multiple modes if any(*); in this case, try to use an unstructured mesh instead?

• Warning in case of PML: integration is currently done over the whole cross-section (see also compute_energy_velocity), the so-defined complex power flow is questionable

(*) e.g. flexural modes in a cylinder with structured mesh, whether the medium is lossy or lossless

compute_response(dof, z, omega_index=None, spectrum=None, wavenumber_function=None, plot=False)

Post-process the response (modal expansion) at the degree of freedom dof and the axial coordinate z, for the whole frequency range. The outputs are frequency, a numpy 1d array of size len(omega), and response, a numpy 2d array of size len(dof or z)*len(omega). dof and z cannot be both vectors, except if omega_index is specified or omega is scalar (single frequency computation): in that case, the array response is of size len(z)*len(dof), which can be useful to plot the whole field at a single frequency.

The response at each frequency omega is calculated from: U(z,omega) = sum qm(omega)*Um(omega)*exp(i*km*z), m=1…M, where z is the receiver position along the waveguide axis. M is the number of modes traveling in the proper direction, positive if z is positive, negative if z is negative. The pairing of opposite-going eigenvectors is required, unpaired modes are discarded from the expansion.

The outputs frequency and response are made dimensional when values in plot_scaler are not set to 1.

Assumption: the source is assumed to be centred at z=0.

Warning: the response calculation is only valid if z lies oustide the source region.

Note: spectrum and wavenumber_function can be specified in compute_response_coefficient(…) instead of compute_response(…), but not in both functions in the same time (otherwise the excitation will be modulated twice).

Parameters

dofnumpy array of integer

dof where the response is computed

znumpy array

axial coordinate where the response is computed

omega_indexint

omega index to compute the response at a single frequency, allowing the consideration of multiple dof and z

spectrumnumpy.ndarray

when specified, spectrum is a vector of length omega used to modulate F in terms of frequency (default: 1 for all frequencies)

wavenumber_function: python function

when specified, wavenumber_function is a python function used to modulate F in terms of wavenumber (example: wavenumber_function = lambda x: np.sin(x), default: 1 for all wavenumbers, i.e. source localized at z=0)

plotbool

if set to True, the magnitude and phase of response are plotted as a function of frequency

Returns

frequencynumpy 1d array

the frequency vector, i.e. omega/(2*pi)

responsenumpy array (1d or 2d)

the matrix response

ll_abs : matplotlib list of lines for magnitude plot when plot is set to True ll_angle : matplotlib list of lines for phase plot when plot is set to True

track_mode(omega_index, mode_index, threshold=0.9, plot=False)

Track a mode over the whole frequency range. The mode is specified by its index, mode_index, at a given angular frequency index, omega_index. Tracking is performed thanks to similarity between eigenvectors and eigenforces (value between 0 and 1). Tracking is stopped if similarity becomes lower than threshold. It returns mode, the index list identifying the mode position at each frequency (index is set to -1 for frequencies at which the mode has not been successfully tracked due to low similarity). If plot is set to True, the real and imaginary parts of eigenvalue are plotted w.r.t. frequency index, for visual check that the desired mode has been properly tracked.

plot_phase_velocity(**kwargs)

Plot phase velocity dispersion curves, vp vs. Re(omega), where omega is replaced with frequency for dimensional results. Parameters and Returns: see plot(…).

plot_attenuation(**kwargs)

Plot attenuation dispersion curves, Im(wavenumber) vs. Re(omega) if omega is the parameter, or Im(omega) vs. Re(wavenumber) if wavenumber is the parameter, where omega is replaced with frequency for dimensional results. Parameters and Returns: see plot(…).

plot_energy_velocity(**kwargs)

Plot energy velocity dispersion curves, ve vs. Re(omega), where omega is replaced with frequency for dimensional results. Parameters and Returns: see plot(…).

plot_group_velocity(**kwargs)

Plot group velocity dispersion curves, vg vs. Re(omega), where omega is replaced with frequency for dimensional results. Parameters and Returns: see plot(…).

plot_coefficient(**kwargs)

Plot response coefficients as a function of frequency, |q| vs. Re(omega), where omega is replaced with frequency for dimensional results. Parameters and Returns: see plot(…).

plot_excitability(**kwargs)

Plot excitability as a function of frequency, |e| vs. Re(omega), where omega is replaced with frequency for dimensional results. Parameters and Returns: see plot(…).

plot_complex_power(**kwargs)

Plot complex power flow of individual modes as a function of frequency, Re(P) and Im(P) vs. Re(omega), where omega is replaced with frequency for dimensional results. Gray color is used for the imaginary part. Parameters and Returns: see plot(…).

plot(x=None, y=None, c=None, direction=None, pml_threshold=None, mode=None, ax=None, color='k', marker='o', markersize=2, **kwargs)

Plot dispersion curves y[1](y[0]) vs. x[1](x[0]) as scatter plot. If the index list, mode, is specified by the user, a single mode is plotted as a continuous single colored curve (the index list mode can be obtained from method track_mode(…)).

Parameters

x, y, c: list

x[0], y[0], c[0] are strings corresponding to modal properties for the x-axis, y-axis and marker colors respectively (these strings can be: ‘omega’, ‘wavenumber’, ‘energy_velocity’, ‘group_velocity’, ‘pml_ratio’, ‘eigenvalues’, ‘excitability’, ‘eigenvectors’, ‘eigenforces’, ‘coefficient’, ‘power flow’, ‘frequency’, ‘attenuation’, ‘phase_velocity’), x[1], y[1], c[1] are the functions applied to x[0], y[0] and c[0] respectively (e.g. np.abs, np.angle, np.real, np.imag, etc.). If x is None but not y, x is set to [‘omega’, np.real] if results are normalized, or set to [‘frequency’, np.real] if they are dimensional. If both x and are None, plot dispersion curves Re(omega) or Re(frequency) vs. Re(wavenumber). If c is None, a single color is used for coloring markers (given by the input variable color).

direction: int

+1 for positive-going modes, -1 for negative-going modes, None for plotting all modes

pml_threshold: float

threshold to filter out PML modes (modes such that pml_ratio<pml_threshold)

mode: index list

index list identifying the mode position at each frequency

ax: matplotlib axis

the matplotlib axis on which to plot data (created if None)

color: str, marker: str, markersize: int, linestyle: str, **kwargs are passed to ax.plot

Returns

sc: the matplotlib collection

set_plot_scaler(length=1, time=1, mass=1, dim=3)

Define the characteristic length, time and mass in order to visualize plots in a dimensional form (by default, they are equal to 1). Set dim=3 for three-dimensional waveguides, dim=2 for two-dimensional waveguides (e.g. plates). Scaling factors for ‘omega’, ‘wavenumber’, ‘energy_velocity’, ‘group_velocity’, ‘pml_ratio’, ‘eigenvalues’, ‘excitability’, ‘eigenvectors’, ‘eigenforces’, ‘coefficient’, ‘complex_power’, ‘frequency’, ‘attenuation’, ‘phase_velocity’ are stored in the attribute name plot_scaler. If poynting normalization has already been applied, then the scalers for ‘eigenvectors’, ‘eigenforces’ and ‘coefficient’ are such that the dimensional cross-section power flow of eigenmodes is equal to 1 Watt (if no poynting normalization applied, these scalers are left to 1). Reminder: while the dimension of U (displacement) is in meter, the dimension of F (force) is in Newton for 3D waveguides and in Newton/meter for 2D waveguides (F is in mass*length**(dim-2)/time**2).

plot_spectrum(index=0, c=None, ax=None, color='k', marker='o', markersize=2, **kwargs)

Plot the spectrum, Im(k) vs. Re(k) computed for omega[index] (if the parameter is the frequency), or Im(omega) vs. Re(omega) for wavenumber[index] (if the parameter is the wavenumber).

Parameters

index: int

parameter index

c: list

c[0] is a string (must be an attribute of self) and c[1] is a function used for coloring markers, a single color (given by the input variable color) is used if c is None

ax: matplotlib axis

the matplotlib axis on which to plot data (created if None)

color: str, marker: str, markersize: int, linestyle: str, **kwargs are passed to ax.plot

Returns

sc: the matplotlib collection

class waveguicsx.waveguide.Signal(time=None, waveform=None, frequency=None, spectrum=None, alpha=0)

Bases: object

A class for handling signals in the time domain and in the frequency domain

Reminder:

• the sampling frequency fs must be at least twice the highest excited frequency (fs>=2fmax)

• the time duration T must be large enough to capture the slowest wave at z, the source-receiver distance

Fourier transform definition used: X(f) = 2/T * integral of x(t)*exp(+i*omega*t)*dt

Two remarks:

1. this is not the numpy fft function convention, which is in exp(-i*omega*t)

2. the true amplitude of the Fourier transform, when needed, has to be obtained by multiplying the output (spectrum) by the scalar T/2, where T is the duration of the time signal (with the above definition: the division by T simplifies dimensionless analyses, and the factor 2 is used because only the positive part of the spectrum is considered)

Complex Fourier transform:

A complex Fourier transform is applied if alpha is set to a nonzero value. The frequency vector has then an imaginary part, constant and equal to alpha/(2*pi). Complex frequency computations can be useful for the analysis of long time duration signals (avoids aliasing). A good choice is alpha = log(50)/T. Note that the first frequency component is kept in that case (the frequency has a zero real part but non-zero imaginary part).

Example:

mysignal = Signal(alpha=0*np.log(50)/5e-4)
mysignal.toneburst(fs=5000e3, T=5e-4, fc=100e3, n=5)
mysignal.plot()
mysignal.fft()
mysignal.plot_spectrum()
mysignal.ifft(coeff=1)
mysignal.plot()
plt.show()

Attributes

timenumpy 1d array

time vector

waveformnumpy nd array

waveform vectors stacked as rows (waveform is an array of size number_of_signals*len(time))

frequencynumpy 1d array

frequency vector

spectrumnumpy nd array

spectrum vectors stacked as rows (spectrum is an array of size number_of_signals*len(frequency))

alphafloat

decaying parameter to apply complex Fourier transform (useful for long time duration signal)

Methods

__init__(time=None, waveform=None, frequency=None, spectrum=None, alpha=0):

Constructor, initialization of signal (specify either waveform vs. time or spectrum vs. frequency)

fft():

Compute Fourier transform, results are stored as attributes (names: frequency, spectrum)

ifft(coeff=1):

Compute inverse Fourier transform, results are stored as attributes (names: time, waveform)

ricker(fs, T, fc):

Generate a Ricker signal

toneburst(fs, T, fc, n):

Generate a toneburst signal

chirp(fs, T, f0, f1, chirp_duration):

Generate a chirp signal

plot(ax=None, color=”k”, linewidth=1, linestyle=”-”, **kwargs):

Plot time waveform (waveform vs. time)

plot_spectrum(ax=None, color=”k”, linewidth=1, linestyle=”-”, **kwargs):

Plot the spectrum (spectrum vs. frequency), in magnitude and phase

fft()

Compute Fourier transform (positive frequency part only, time waveform are assumed to be real). If the number of time steps is odd, one point is added. The zero frequency, if any, is suppressed. Results are stored as attributes (names: frequency, spectrum). spectrum is an array of size number_of_signals*len(frequency)

ifft(coeff=1)

Compute inverse Fourier transform (only the positive frequency part is needed, time waveform are assumed to be real). Zero padding is applied in the low-frequency range (if missing) and in the high-frequency range (if coeff is greater than 1). Zero padding in the high frequency range is applied up to the frequency coeff*max(frequency). Results are stored as attributes (names: time, waveform). waveform is an array of size number_of_signals*len(time).

ricker(fs, T, fc)

Generate a Ricker wavelet signal of unit amplitude (fs: sampling frequency, T: time duration, fc: Ricker central frequency)

Note that for better accuracy:

• fs is rounded so that fs/fc is an integer

• T is adjusted so that the number of points is even

toneburst(fs, T, fc, n)

Generate a toneburst signal (fs: sampling frequency, T: time duration, fc: central frequency, n: number of cycles). This signal is a Hanning-modulated n cycles sinusoidal toneburst centred at fc Hz (with unit amplitude). For this kind of excitation, fmax can be considered as 2*fc roughly, hence one should choose fs>=4fc.

Note that for better accuracy:

• fs is rounded so that fs/fc is an integer

• T is adjusted so that the number of points is even

chirp(fs, T, f0, f1, chirp_duration)

Generate a chirp of unit amplitude (fs: sampling frequency, T: time duration, f0: first frequency, f1: last frequency, chirp_duration: time to sweep from f0 to f1). Note that for better accuracy, T is adjusted so that the number of points is even.

plot(ax=None, color='k', linewidth=1, linestyle='-', **kwargs)

Plot time waveform (waveform vs. time)

Parameters

ax: matplotlib axis

the matplotlib axis on which to plot data (created if None)

color: str, linewidth: int, linestyle: str, **kwargs are passed to ax.plot

Returns

ll: the matplotlib list of lines

plot_spectrum(ax=None, color='k', linewidth=1, linestyle='-', **kwargs)

Plot the spectrum (spectrum vs. frequency), in magnitude and phase

Parameters

ax: matplotlib axis

the matplotlib axis on which to plot data (created if None)

color: str, linewidth: int, linestyle: str, **kwargs are passed to ax.plot

Returns

ll_abs: the matplotlib list of lines for magnitude plot ll_angle: same but for phase plot

scattering.py

class waveguicsx.scattering.Scattering(comm: _MPI.Comm, M: Mat, K: Mat, C: Mat, tbcs: list)

Bases: object

A class for solving scattering problems by local inhomogeneities in complex waveguides based on PETSc.

The full documentation is entirely defined in the `scattering.py’ module.

The following matrix problem is considered: (K-omega**2*M-1j*omega*C)*U=F. This kind of problem typically typically stems from a finite element (FE) model of a small portion of waveguide including a local inhomogeneity (e.g. defects). The cross-section extremities of the truncated FE model are then handled as transparent boundary conditions (BCs) to reproduce semi-infinite waveguides. The so-obtained scattering problem is solved repeatedly for each frequency. The loops over the angular frequency can be parallelized, as shown in some tutorials (using mpi4py).

This class enables to deal with scattering in complex waveguides, two-dimensional (e.g. plates) or three-dimensional (arbitrarily shaped cross-section), inhomogeneous in the transverse directions, anisotropic. Complex-valued problems can be handled including the effects of non-propagating modes (evanescent, inhomogeneous), viscoelastic loss (complex material properties) or perfectly matched layers (PML) to simulate buried waveguides.

Transparent BCs are Waveguide objects, which must have been solved prior to the scattering problem solution, yielding the following object attributes: omega, eigenvalues, eigenvectors, eigenforces and traveling direction (see waveguide.py module for details). Transparent BCs are localized by their degrees of freedom in the global vector U. This means that the local degree of freedom i of a given eigenvector/eigenforce stored in a Waveguide object is located at the global degree of freedom dofs[i] of the FE model.

The user must supply the following inputs: - K, M, C, the global FE matrices (stiffness, mass and viscous damping) - F and F_spectrum, the global FE vector of internal excitation (i.e. sources inside the FE model), if any,

and its spectrum

• tbcs, a list of pairs (name, dofs) which characterize the transparent BCs, where name is a string specifying the attribute name of a given transparent BC (this attribute will be a Waveguide object) and dofs is a numpy array of the global degrees of freedom for this transparent BC.

• the ingoing mode coefficients, specified by the attribute coefficient in each transparent BC

Important convention: if the waveguide axis of a transparent BC is oriented outward (i.e. outside the FE box), dofs are positive, but if oriented inward (i.e. inside the FE box), a negative sign has to be assigned to dofs by the user.

The solution to the scattering problem yields: - the displacement U of the FE model for each angular frequency omega - the outgoing modal coefficients of every transparent BC and for each omega - the energy balance post-processed for each angular frequency

which enables to check the error due to the modal truncature introduced in the transparent BCs.

See Attributes below for more details.

Example:

# This simple example involves only one transparent boundary condition (e.g. waveguide scattering by a free edge)
# The tbc (the "inlet") is supposed to be at the left-hand side of the FE box so that its outward normal is negative

from waveguicsx.waveguide import Waveguide
from waveguicsx.scattering import Scattering

# Input parameters
omega = 2*np.sqrt(3)*np.linspace(1.48, 1.60, num=100) #normalized angular frequency range
nev = 30 #tbc number of eigenvalues requested at each frequency

# Scattering initialization
ws = Scattering(MPI.COMM_WORLD, M, K, 0*M, [('waveguide0', -tbc_dofs)]) #M and K are the mass and stiffness matrices of the FE box
#reminder: tbc_dofs are the global degrees of freedom, set negative by convention when the normal is negative (here, we suppose n=-ey)

#Solve waveguide problem associated with the tbc
ws.waveguide0 = Waveguide(MPI.COMM_WORLD, Ms, K0, K1, K2) #Ms, K0, K1 and K2 are SAFE matrices associated with the tbc (here, named 'waveguide0')
ws.waveguide0.set_parameters(omega=omega)
ws.waveguide0.solve(nev)
ws.waveguide0.compute_traveling_direction()
ws.waveguide0.compute_poynting_normalization()

# Solving scattering problem
mode = ws.waveguide0.track_mode(frequency_index, mode_index, threshold=0.98, plot=True) #track a mode, specified by its index at a given frequency, over the whole frequency range
ws.set_ingoing_mode('waveguide0', mode) #set mode as a single ingoing mode, coeff is 1 (here, power is also 1 thanks to poynting normalization)
ws.set_parameters()
ws.solve()

# Plot reflected power coefficients vs. angular frequency
ws.waveguide0.compute_complex_power()
ws.waveguide0.plot(y=('complex_power', lambda x:np.abs(np.real(x))), direction=-1)

Attributes

commmpi4py.MPI.Intracomm

MPI communicator (parallel processing)

M, K, Cpetsc4py.PETSc.Mat

FE matrices

Fpetsc4py.PETSc.Vec

FE force vector (internal excitation)

F_spectrumnumpy.ndarray

the modulation spectrum of F (size must be the same as omega)

tbclist of pairs (name, dofs)

name is a string corresponding to the desired attribute name of a tbc (Waveguide object), dofs is a numpy array of the degrees of freedom of the tbc (positive if outward negative if inward)

omeganumpy.ndarray

the angular frequency range, specified by the user in tbcs (Waveguide objects)

ksp: KSP object (PETSc object)

solver parameters

displacementlist of PETSc vectors

for each angular frequency, the FE displacement vector U (solution of the scattering problem)

energy_balancelist of 1d numpy array

for each angular frequency, energy_balance[i] gives the following three-component array, [Pin-i*omega/2*U^H*F, Ptot, -i*omega*U^H*D*U], where P=-i*omega/2*U^H*T: - the term Pin-i*omega/2*U^H*F represents the input power, supplied by ingoing modes and internal

forces (this term should have a negative real part)

• the term Ptot is the complex power flow of the sum of ingoing and ougoing modes

• the term -i*omega*U^H*D*U is related to the kinetic, potential and dissipated energies in the volume (the dissipated energy, defined as positive, is equal to the real part of this term divided by -2*omega).

• a perfect energy balance is achived if Ptot = -i*omega*U^H*D*U

name.coefficientlist of numpy.ndarray, where name is a Waveguide object associated with a given transparent BC

this transparent BC attribute stores modal coefficients at each frequency: - the coefficients of ingoing modes are considered as inputs, specified by the user (excitation in the scattering problem) - the coefficients of outgoing modes are considered as initially unknown, solution of the scattering problem Any non-zero outgoing amplitudes specified in name.coefficient prior to solving the scattering problem will hence be discarded and replaced with the scattering solution. If the attribute coefficient is empty prior to scattering (no specified ingoing modes), zero ingoing amplitudes will be considered by default.

Methods

__init__(comm:’_MPI.Comm’, M:PETSc.Mat, K:PETSc.Mat, C:PETSc.Mat, tbcs:list):

Constructor, initialization of scattering problem

set_ingoing_mode(tbc_name, mode, spectrum=None):

For a given tbc, specified by the string tbc_name, set the coefficient of a single mode to 1. mode is a list of indices identifying the mode position at each frequency.

set_internal_excitation(F, F_spectrum=None):

Set the internal excitation vector F and its spectrum F_spectrum

set_parameters(solver=’iterative’):

Set default parameters of KSP solver (stored in attribute ksp)

solve():

Solve the scattering problem repeatedly for the angular frequency range, solutions are stored as attributes (names: displacement, energy_balance)

plot_energy_balance():

Plot the three terms of energy_balance (complex modulus) at each frequency index for checking modal tbc truncature

set_ingoing_mode(tbc_name, mode, spectrum=None)

For a given tbc, specified by the string tbc_name, set the coefficient of a single mode to 1. The coefficients of all other modes are set to 0. mode is a list of indices identifying the mode position at each frequency (see also method track_mode() of class Waveguide). Please ensure that this mode is an ingoing mode (non-zero outgoing amplitudes will be discarded). When specified, spectrum is a vector of length omega (numpy.array) used to modulate the coefficient in terms of frequency (default: 1 for all frequencies).

set_internal_excitation(F, F_spectrum=None)

Set the internal excitation vector F (petsc4py.PETSc.Vec). When specified, F_spectrum is a vector of length omega (numpy.array) used to modulate F in terms of frequency (default: 1 for all frequencies). F and F_spectrum are stored as attributes.

set_parameters(solver='iterative')

Set default parameters of KSP solver (stored into the attribute ksp.) The preselected methods are CGS (iterative method) and MUMPS (direct method). CGS is used by default, set solver=’direct’ to use MUMPS instead. After calling this method, various PETSc parameters can be set by changing the attribute ksp manually.

solve()

Solve the scattering problem, i.e. the linear system repeatedly for the angular frequency range. The solutions are stored in the attributes displacement and energy_balance. The FE problem D*U=F is transformed into: Bout^T*(D*Bout-Tout)*Uout = Bout^T*(Tin-D*Bin)*Uin where: - D is the dynamic stiffness matrix (D=K-omega**2*M-1j*omega*C) - Bin and Bout (resp. Tin and Tout) are bases containing ingoing and outgoing modal

displacements (resp. forces) so that: U=Bin*Uin+Bout*Uout, F=Tin*Uin+Tout*Uout

• Uin contains internal forces and known ingoing amplitudes (zero by default, if not specified)

• Uout is the solution containing outgoing modal amplitudes and internal dofs

plot_energy_balance(ax=None, color='k', linewidth=1, linestyle='-', **kwargs)

Plot energy balance at each frequency index for checking modal tbc truncature. The energy balance is defined as the difference between tbc net flux and volume power divided by the input power (values close to zero indicate that enough modes have been included in the transparent BCs).

Parameters

ax: matplotlib axis

the matplotlib axis on which to plot data (created if None)

color: str, linewidth: int, linestyle: str, **kwargs are passed to ax.plot

Returns

ll: the matplotlib list of lines