Page 222 - Process Modelling and Simulation With Finite Element Methods

Simulation and Nonlinear Dynamics 209

This piece of code deserves several comments:

There is no direct way of calculating eigenpairs using f emeig for transient
models. Instead, we use assemble to create the three sparse matrices of equation
(5.14) needed for the generalized eigenpair solution. f emeig does accept these
three matrices as inputs, and does not further reference the fem structure that
created them. This recipe also works for parametric continuation, which has a p-
list of parameters and an array of solutions for each parameter. Transient models
have a fem.sol with a t-list of times and an associated array of the same length of
solution vectors.

f erne ig creates a solution structure (see fem.sol) with subfields solhmbda (list
of eigenvalues) and so1.u (array of eigenvectors of same length as sol.lambda).
eigenvectors have the same structure as solution vectors, but only satisfy the
linearized, homogeneous boundary conditions and do not satisfy the pde itself.
In particular, there is no requirement that the load vector (L) or the constraint
load vector (M) be satisfied. Indeed, the code above does not use L or M in
computing eigenpairs.

femeig sometimes has difficulty finding large decay rates. Even though I
requested twenty ‘SM’ eigenpairs, after t=0.04, it can only find 19, and after
t=0.42, only 18. femeig uses the sparse eigenanalysis routines of ARPACK,
which is essentially iterative, to compute eigenvalues and eigenvectors. This
package has difficulty in finding and distinguishing zero eigenvalues (associated
with singular systems). Since [I31 and [15] show theoretically and numerically
that the linear stability theory has a neutral mode at zero wavenumber and at a
finite cut-off wavenumber of the longwave unstable wave packet, &, the linear
system is nearly singular and will have difficulty resolving these neutral (or
numerically near-neutral) modes.
Figure 5.19 shows the eigenvalue with least real part at each instant in time, as
computed from vf-eig.m. Due to the computational intensity of computing the
eigenvalues of these large sparse matrices, it is recommended to execute this m-
file script as a background job in MATLAB:
matlab -nojvm <vf-eig.m >err 2>err &
on the UNIX command shell.

From Figure 5.19, it can be seen that for a range of times shortly after t=O, up
until t=O. 15, the largest growth rate is roughly constant, but with substantial
scatter. Our animations showed steady growth of the fingered instability during
this interval. That there is substantial scatter is not surprising to me, as
Zimmerman and Homsy [15] computed average growth rates versus
wavenumber for 1282 Fourier modes in a similar model but with anisotropic
dispersion. They found growth of power in each mode was sporadic from time