Dynamic Mode Decomposition

Description

This processing computes a Dynamic Mode Decomposition (DMD) on 2D/3D field.

The DMD allows to characterize the evolution of a complex non-linear system by the evolution of a linear system of reduced dimension.

For the vector quantity of interest v(xn,t)Rn (for example the velocity field of a mesh), where t is the temporal variable and xn the spatial variable. The DMD provides a decomposition of this quantity in modes, each having a complex angular frequency ωi:

v(xn,t)i=1mciΦi(xn)ejωit

or in discrete form:

vki=1mλikciΦi(xn),with λik=ejωiΔt

with Φi are the DMD modes.

Using the DMD, the frequencies ωi are obtained with a least-square method.

Compared to a FFT, the DMD has thus the following advantages:

  • Less spectral leaking

  • Works even with very short signals (one or two periods of a signal can be enough)

  • Takes advantage of the large amount of spatial information

To illustrate the tool, the following reference case is chosen: a wave on a square 2D mesh of size 100×100 from (x,y)=(0,0) to (x,y)=(1,1). The unsteady field

p(x,y,t)=Asin(nxπx)sin(nyπy)cos(ωt+Φ(x))

is shown in the left figure below. The amplitude is shown in the middle figure and the phase on the right figure.

The following parameters are used:
  • A=100

  • nx=ny=1

  • f=215Hz

  • Φ(x)=4/3πx

DMD2D_pressure DMD2D_module DMD2D_phase

A database of N=100 snapshots is constructed using a time step Δt=2×104s, leading to a Nyquist frequency fe=1/2Δt=2500Hz and a resolution Δf=1/NΔ=50Hz.

Using the DMD algorithm, the spectrum and amplifications obtained are shown in the figures below:

DMD2D_Pmodulus DMD2D_Ampl

The exact value of 215 Hz expected from the frequency peak is present on the amplitude spectrum. The amplification is very close to zero, which is expected since the wave is neither amplified nor attenuated.

Construction

import antares
myt = antares.Treatment('dmd')

Parameters

  • base: Base

    The input base that contains one zone which contains one instant. The instant contains at least two 2-D arrays that represent the cell volumes and a scalar value.

  • time_step: float

    The time step between two instants Δt.

  • noiselevel: float, default= 1e-9

    The noise level should be increased when ill-conditioned and non-physical modes appear.

  • type: str, default= ‘mod/phi’

    The decomposition type of the DMD:

    • mod/phi for modulus/phase decomposition or

    • im/re for imaginery/real part decomposition.

  • variables: list(str)

    The variable names. The variable representing the volume of each cell must be the first variable in the list of variables. The other variables are the DMD variables.

  • memory_mode: bool, default= False

    If True, the initial base is deleted on the fly to limit memory usage.

  • list_freq: list(float), default= []

    Frequency peaks to detect in Hertz. If no specific frequency is given, the DMD will be achieved on all modes (costly).

  • delta_freq: float, default= 0

    The frequential resolution of the base, ie Δf=1/T, with T the total time of the database (T=number of instants×Δt)

Preconditions

The input base must be an unstructured mesh.

The variable representing the volume of each cell must be the first variable in the list of variables. The other variables are the DMD variables.

Moreover, the frequential resolution of the base must be Δf=1/T in Hertz, with T the total time of the database (T=number of instants×Δt).

Postconditions

The treatment returns one output base containing:

  • a zone named ‘spectrum’ with an instant named ‘spectrum’ with three variables ‘frequency’, ‘amplification’, and ‘p_modulus’ (magnitude)

  • a zone named ‘modes’. Each DMD variable modes will be stored in different instants. Each instant contains two variables depending on the type of DMD (e.g. ‘p_mod’ and ‘p_phi’)

The output zone ‘modes’ contains some keys ‘peak_modes_<variable name>’ in its Zone.attrs. These peak modes are the name of instants in the zone.

Each instant in the output zone ‘modes’ contains two keys ‘frequency’ and ‘amplification’ in its Instant.attrs.

References

P. J. Schmid: Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech., vol. 656, no. July 2010, pp. 5–28, 2010.

Example

import antares
treatment = antares.Treatment('dmd')
treatment['time_step'] = 0.1
treatment['variables'] = ['volumes', 'pressure', 'vx']
treatment['noiselevel'] = 1e-4
treatment['base'] = base
treatment['list_freq'] = [300.0, 600.0]
treatment['delta_freq'] = 50.0
treatment['type'] = 'im/re'
result = treatment.execute()

Main functions

class antares.treatment.TreatmentDmd.TreatmentDmd
execute()

Execute the treatment.

Returns:

a base that contains one zone with one instant. the instant has two variables.

Return type:

Base

Example

"""
This example illustrates the Dynamic Mode Decomposition.
It is a stand alone example : the Base used for the DMD
is not read in a file, it is created in the example
"""
import os
import numpy as np
from antares import Base, Instant, Treatment, Writer, Zone


def myfunc(x, y, t, A, w, Phi):
    return A * np.sin(np.pi*x)*np.sin(np.pi*y)*np.cos(w*t+Phi*x)


if not os.path.isdir('OUTPUT'):
    os.makedirs('OUTPUT')

# ----------------------
# Create the input base
# ----------------------

# Initial parameters
time_step = 2e-4
nb_snapshot = 50
T = nb_snapshot*time_step
A1 = 100.0
A2 = 10.0
f1 = 215.0
f2 = 120.0
w1 = f1*2.0*np.pi
w2 = f2*2.0*np.pi
Phi = 4.0/3.0*np.pi

# Create a 2D mesh
n = 50
mgrid = np.mgrid[0:n, 0:n]
xt = (mgrid[0]/(n-1.))
yt = (mgrid[1]/(n-1.))

# Initialise the base
base = Base()
base['zone_0'] = Zone()
base[0].shared['x'] = xt
base[0].shared['y'] = yt
base.unstructure()

# Define the cell volume
ncell = base[0].shared['x'].shape[0]
volume = np.ones((ncell))
base[0].shared['cell_volume'] = volume

# Extract coordinates and define time
x = base[0].shared['x']
y = base[0].shared['y']
time = np.linspace(0.0, T, nb_snapshot)

# Generate temporal database
for snapshot_idx in range(nb_snapshot):
    instant_name = 'snapshot_%s' % snapshot_idx
    base[0][instant_name] = Instant()
    base[0][instant_name]['pressure'] = myfunc(x, y, time[snapshot_idx], A1, w1, Phi)
    base[0][instant_name]['vx'] = myfunc(x, y, time[snapshot_idx], A2, w2, Phi)

# --------------
# DMD treatment
# --------------
treatment = Treatment('Dmd')
treatment['time_step'] = time_step
treatment['variables'] = ['cell_volume', 'pressure', 'vx']
treatment['noiselevel'] = 1e-9
treatment['base'] = base
treatment['list_freq'] = [f1, f2]
treatment['delta_freq'] = 50.0
treatment['type'] = 'im/re'
result = treatment.execute()

# -------------------
# Writing the result
# -------------------

# Write the spectrum (frequency, magnification and magnitude)
writer = Writer('column')
writer['filename'] = os.path.join('OUTPUT', 'ex_dmd_spectrum.dat')
writer['base'] = result[('spectrum', )]
writer.dump()

# Pick and write modes: here we retreive the peak modes of pressure for each variables
peak_modes = result['modes'].attrs['peak_modes_pressure']
newbase = result[('modes', ), peak_modes]
newbase[0].shared.connectivity = base[0][0].connectivity
newbase[0].shared['x'] = base[0][0]['x']
newbase[0].shared['y'] = base[0][0]['y']
writer = Writer('bin_tp')
writer['filename'] = os.path.join('OUTPUT', 'ex_dmd_pressure_peak_mode_<instant>.plt')
writer['base'] = newbase
writer.dump()

# Pick and write modes: here we retreive the peak modes of velocity for each variables
peak_modes = result['modes'].attrs['peak_modes_vx']
newbase = result[('modes', ), peak_modes]
newbase[0].shared.connectivity = base[0][0].connectivity
newbase[0].shared['x'] = base[0][0]['x']
newbase[0].shared['y'] = base[0][0]['y']
writer = Writer('bin_tp')
writer['filename'] = os.path.join('OUTPUT', 'ex_dmd_vx_peak_mode_<instant>.plt')
writer['base'] = newbase
writer.dump()

# Write the effective frequencies of identified modes
fic = open(os.path.join('OUTPUT', 'ex_dmd_list_modes.dat'), 'w')
for key in result[('modes', ), peak_modes][0].keys():
    fic.write(key)
    fic.write('\t')
    fic.write(str(result[('modes', ), peak_modes][0][key].attrs['frequency']))
    fic.write('\n')