Physical models

Summary

Package: dynamodels.physical (re-exported as romda.models.physical). All use IVPIntegrator (scipy.integrate.solve_ivp) except KS, which uses a discrete ETDRK4 map.

Class	Dim	Key parameters	Integrator
VdP	2	beta, zeta, kappa, law, omega	IVP
Lorenz63	3	rho, sigma, beta	IVP
Lorenz96	Nx	F, Nx	IVP
KS	Nx	nu, L, Nx	Discrete (ETDRK4)
Rijke	2Nm+Nc	beta, tau, C1, C2, kappa	IVP
Annular	4	omega, nu, c2beta, kappa, epsilon	IVP

---

dynamodels.physical.van_der_pol.VdP(**model_dict)

Bases: Model

Van der Pol oscillator — low-order model of a longitudinal thermoacoustic mode.

The acoustic pressure mode \(\eta\) evolves as

\[ \ddot{\eta} + \omega^2 \eta = \dot{\eta} \left( \beta - \zeta - \kappa \, g(\eta) \right), \]

with a cubic (\(g = \eta^2\), law='cubic') or arctangent-saturated (\(g = \eta^2 / (1 + \kappa \eta^2 / \beta)\), law='tan') heat-release law. The estimable parameters are the linear growth rate \(\beta\), the damping \(\zeta\) and the nonlinear saturation \(\kappa\).

References

Nóvoa & Magri (2022). Real-time thermoacoustic data assimilation. J. Fluid Mech., 948, A35. DOI: 10.1017/jfm.2022.653.

Source code in dynamodels/physical/van_der_pol.py

def __init__(self, **model_dict):

    psi0 = model_dict.pop('psi0', np.array([0.1, 0.1]))
    dt = model_dict.pop('dt', 1e-4)

    super().__init__(psi0=psi0, dt=dt, integrator_class=IVPIntegrator, **model_dict)

    #  Add fixed input_parameters
    self.alpha_labels = dict(beta='$\\beta$', zeta='$\\zeta$', kappa='$\\kappa$')
    self.alpha_lims = dict(zeta=(5, 120), kappa=(0.1, 20), beta=(5, 120))

dynamodels.physical.lorenz63.Lorenz63(**model_dict)

Bases: Model

Lorenz (1963) system — chaotic benchmark with three state variables.

\[ \dot{x} = \sigma (y - x), \qquad \dot{y} = x (\rho - z) - y, \qquad \dot{z} = x y - \beta z. \]

With the classical parameters (\(\sigma = 10\), \(\rho = 28\), \(\beta = 8/3\)) the system is chaotic with leading Lyapunov exponent \(\lambda_1 \approx 0.906\).

References

Lorenz (1963). Deterministic nonperiodic flow. J. Atmos. Sci., 20, 130–141.

Source code in dynamodels/physical/lorenz63.py

def __init__(self, **model_dict):


    psi0 = model_dict.pop('psi0', np.array([1.0, 1.0, 1.0]))
    dt = model_dict.pop('dt', 0.02)

    self.observe_dims = model_dict.pop('observe_dims', [0, 1, 2]) # Default to observing all dimensions if not specified
    self.Nq = len(self.observe_dims)

    super().__init__(psi0=psi0, dt=dt, integrator_class=IVPIntegrator, **model_dict)

    # measured 1/lambda1 at this rho (set once: params change by re-instantiation)
    self.t_lyap = self.t_lyap_from_table(self.rho, _LAM1_MEASURED, Lorenz63.t_lyap)

    self.alpha_labels = dict(rho='$\\rho$', sigma='$\\sigma$', beta='$\\beta$')

time_derivative(t, psi, sigma, rho, beta) staticmethod

Calculates the time derivative of the Lorenz 63 system. Note: This derivative must handle the augmented state vector (psi). The augmented parameters are stored after the core state (x, y, z).

Source code in dynamodels/physical/lorenz63.py

@staticmethod
def time_derivative(t, psi, sigma, rho, beta):
    """
    Calculates the time derivative of the Lorenz 63 system.
    Note: This derivative must handle the augmented state vector (psi).
    The augmented parameters are stored after the core state (x, y, z).
    """

    x1, x2, x3 = psi[:3]

    # The parameter values used for the current derivative calculation
    # These come from the 'params' dict passed by the integrator
    dx1 = sigma * (x2 - x1)
    dx2 = x1 * (rho - x3) - x2
    dx3 = x1 * x2 - beta * x3

    return (dx1, dx2, dx3) + (0,) * (len(psi) - 3)

visualize_attractor(psi_cases=None, **kwargs)

Visualizes the Lorenz attractor for given state trajectories.

Parameters:

Name	Type	Description	Default
psi_cases	list	State trajectories to plot, each of shape (Nt, 3) or (Nt, 3, Ne). Defaults to the model's own history.	None
**kwargs		Plotting options forwarded to the helper (color, figsize, ...).	{}

Source code in dynamodels/physical/lorenz63.py

def visualize_attractor(self, psi_cases=None, **kwargs):
    """
    Visualizes the Lorenz attractor for given state trajectories.
    Parameters
    ----------
    psi_cases : list, optional
        State trajectories to plot, each of shape ``(Nt, 3)`` or ``(Nt, 3, Ne)``.
        Defaults to the model's own history.
    **kwargs
        Plotting options forwarded to the helper (``color``, ``figsize``, ...).
    """
    if psi_cases is None:
        psi_cases = [self.hist[:, :3, :]]

    plot_attractor(psi_cases, **kwargs)

The Lorenz63 "butterfly" attractor at the chaotic point (ρ=28, σ=10, β=8/3), used as a twin-experiment test case for ensemble DA.

Ergodic exploration of the attractor over time.

Bifurcations of the long-term state as ρ varies.

dynamodels.physical.lorenz96.Lorenz96(**model_dict)

Bases: Model

Lorenz (1996) system — chaotic model of \(N_x\) variables on a periodic lattice.

\[ \dot{x}_i = \left(x_{i+1} - x_{i-2}\right) x_{i-1} - x_i + F, \qquad i = 0, \dots, N_x - 1, \]

with indices taken cyclically modulo \(N_x\) (i.e. \(x_{-1} = x_{N_x - 1}\), \(x_{-2} = x_{N_x - 2}\) and \(x_{N_x} = x_0\)). Each variable is coupled quadratically to its two upstream neighbours, damped linearly (\(-x_i\)) and driven by a constant forcing \(F\). With the classical parameters (\(N_x = 40\), \(F = 8\)) the system is chaotic.

References

Lorenz (1996). Predictability: a problem partly solved. Proc. Seminar on Predictability, Vol. 1, ECMWF, Reading, UK, 1-18.

Source code in dynamodels/physical/lorenz96.py

def __init__(self, **model_dict):

    self.Nx = model_dict.pop('Nx', 40)
    psi0 = model_dict.pop('psi0', np.array([1.6] + [1.0] * (self.Nx - 1)))
    dt = model_dict.pop('dt', 0.01)

    self.observed_idx = model_dict.pop('observed_idx', [0, self.Nx//2, self.Nx-1]) # Default to observing three dimensions if not specified
    self.Nq = len(self.observed_idx)

    super().__init__(psi0=psi0, dt=dt, integrator_class=IVPIntegrator, **model_dict)

    # measured 1/lambda1 at this F (Nx=10 table only; set once — params change
    # by re-instantiation)
    if self.Nx == 10:
        self.t_lyap = self.t_lyap_from_table(self.F, _LAM1_MEASURED_NX10, Lorenz96.t_lyap)

    self.alpha_labels = dict(F='$F$')

time_derivative(t, psi, Nx, F) staticmethod

Calculates the time derivative of the Lorenz 96 system (see class docstring).

Source code in dynamodels/physical/lorenz96.py

@staticmethod
def time_derivative(t, psi, Nx, F):
    """
    Calculates the time derivative of the Lorenz 96 system (see class docstring).
    """

    x = psi[:Nx]

    dx = (np.roll(x, -1) - np.roll(x, 2)) * np.roll(x, 1) - x + F

    if len(psi) > Nx:
         dx = np.concatenate((dx, np.zeros(len(psi) - Nx)))

    return dx

dynamodels.physical.kuramoto_sivashinsky.KS(**model_dict)

Bases: Model

Kuramoto-Sivashinsky equation.

\[ u_t + u_{xx} + \nu\, u_{xxxx} + u\,u_x = 0, \qquad x \in (0, L], \]

with periodic boundary conditions \(u(t, 0) = u(t, L)\), \(u_x(t, 0) = u_x(t, L)\). Solved with the ETDRK4 scheme in Fourier space, where the Fourier transform pair is

\[ \hat{u}(k) = \mathcal{F}[u(x)] = \frac{1}{L}\int_0^L u(x)\,e^{-ikx}\,\mathrm{d}x, \qquad u(x) = \mathcal{F}^{-1}[\hat{u}(k)] = \sum_k \hat{u}(k)\,e^{ikx}. \]

On the \(N_x\)-point grid the wavenumbers are \(\alpha_j = 2\pi j / L\), so the (diagonal) linear operator is \(\alpha_j^2 - \nu\,\alpha_j^4\) and the nonlinear term \(u\,u_x\) is computed in physical space and transformed back to Fourier space at every stage.

Parametrization (\(\nu\), \(L\)). The pair is independent: whichever of the two is given fixes that side of the operator.

• nu only (the default): the domain follows the standard nondimensionalization \(L = 2\pi/\sqrt{\nu}\), and the equation is then integrated in its \(\nu = 1\) form on that domain -- so self.nu is 1 afterwards and the stored \((N_x, \nu, L)\) always describes the operator that was actually integrated (this is what makes ntsa.respawn, which rebuilds from fixed_params, bit-faithful).
• L only: \(\nu = 1\) on the given domain.
• both: both are honoured as given, i.e. the genuine two-parameter system \(u_t + u_{xx} + \nu u_{xxxx} + u u_x = 0\) on \((0, L]\).

The one- and two-parameter forms are related by \(v(x', t') = \sqrt{\nu}\,u(x, t)\) with \(x = \sqrt{\nu}\,x'\) and \(t = \nu\,t'\): KS(Nx, L=Lx, nu=visc, dt=dt) and KS(Nx, L=Lx/sqrt(visc), dt=dt/visc) describe the same physical system.

Initialize the KS model.

Sets up the spatial grid, wavenumbers, sensor locations, and initial state.

Parameters:

Name	Type	Description	Default
**model_dict		Model parameters; supported keys are: Nx : int, number of spatial grid points (must be even). nu : float, viscosity parameter. L : float, domain length, domain is (0, L]. nu and L are independent -- see the class docstring for the resolution rules when only one of them is given. dt : float, time step size. initial_amplitude : float, amplitude of the initial condition. Nq : int, number of sensors. sensor_placement_method : str, 'grid' or 'random'. seed : int, random seed for sensor placement. psi0 : ndarray, initial state in Fourier space (optional).	{}

Source code in dynamodels/physical/kuramoto_sivashinsky.py

def __init__(self, **model_dict):
    """Initialize the KS model.

    Sets up the spatial grid, wavenumbers, sensor locations, and initial state.

    Parameters
    ----------
    **model_dict
        Model parameters; supported keys are:

        - ``Nx`` : int, number of spatial grid points (must be even).
        - ``nu`` : float, viscosity parameter.
        - ``L`` : float, domain length, domain is (0, L].
          ``nu`` and ``L`` are independent -- see the class docstring for the
          resolution rules when only one of them is given.
        - ``dt`` : float, time step size.
        - ``initial_amplitude`` : float, amplitude of the initial condition.
        - ``Nq`` : int, number of sensors.
        - ``sensor_placement_method`` : str, ``'grid'`` or ``'random'``.
        - ``seed`` : int, random seed for sensor placement.
        - ``psi0`` : ndarray, initial state in Fourier space (optional).
    """


    # 'nu' has a non-sentinel class default, so an EXPLICIT nu is what marks the
    # two-parameter form; 'L' uses its non-positive class default as the sentinel.
    nu_given = model_dict.get('nu') is not None

    for key in list(model_dict.keys()):
        if key in vars(KS):
            setattr(self, key, model_dict.pop(key))


    if self.Nx % 2 != 0:
        raise ValueError("Nx must be even.")

    L_given = self.L is not None and self.L > 0

    if not L_given and not nu_given and self.nu is None:
        raise ValueError("Either L or nu must be specified.")
    elif not L_given:
        # nu alone: standard nondimensionalization. The domain absorbs nu and the
        # equation is integrated in its nu = 1 form, so (Nx, nu, L) stays a faithful
        # description of the operator (respawn / filename keying).
        self.L = 2 * np.pi / np.sqrt(self.nu)
        self.nu = 1.
    elif not nu_given:
        self.nu = 1.
    # else: both given -- honour both (general two-parameter form).

    assert self.L is not None and self.L > 0, "L must be positive."
    assert self.nu is not None and self.nu > 0, "nu must be positive."

    # Fourier wavenumbers alpha_j = 2 pi j / L on the domain (0, L]
    self.k = 2 * np.pi * np.fft.rfftfreq(self.Nx, d=self.L / self.Nx)

    dt_requested = model_dict.pop('dt', 0.25)
    self.dt = dt_requested


    self.ETDRK4_f_terms = None  # This simply trigers the setter method.


    #  Select sensors ___________________________ #
    if self.sensor_placement_method not in ['grid', 'random']:
        raise NotImplementedError(f"sensor_placement_method '{self.sensor_placement_method}' not recognized.")

    if self.sensor_placement_method == 'grid':
        # Place sensors evenly spaced across the domain
        self.sensor_locations = np.linspace(0, self.Nx-1, self.Nq, endpoint=True, dtype=int)
    elif self.sensor_placement_method == 'random':
        # Place sensors at random locations in the domain
        self.sensor_locations = self.rng.integers(0, self.Nx-1, self.Nq)


    #   Init Model  #
    psi0 = model_dict.pop('psi0', None)
    if psi0 is None:
        # Initialize state in physical space and transform to spectral space
        u0 = self.initial_amplitude * self.rng.standard_normal(self.Nx)
        u0 -= np.mean(u0)  # Zero-mean initial condition
        u_hat = KS.physical_to_fourier(u0)[:, None]     # Transform to Fourier space
        psi0 = np.array(u_hat)


    super().__init__(psi0=psi0, dt=self.dt, integrator_class=DiscreteIntegrator, **model_dict)

    # Model's dt setter rounds to precision_t decimals, which silently perturbs
    # timesteps with more significant digits (e.g. dt = 0.1 * 71 / 16). Keep the
    # exact requested value for stepping (precision_t still governs time stamps)
    # and rebuild the ETDRK4 coefficients with it.
    self._dt = float(dt_requested)
    self.ETDRK4_f_terms = None

get_observables(Nt=1, loc=None, **kwargs)

Get the observable state in physical space at specified sensor locations.

Parameters:

Name	Type	Description	Default
Nt	int	Number of time steps to retrieve. Default is 1.	1
loc	array - like or str	Sensor locations to retrieve observables from. If 'all', returns observables at all spatial points. If None, returns observables at the predefined sensor locations.	None

Source code in dynamodels/physical/kuramoto_sivashinsky.py

def get_observables(self, Nt=1, loc=None, **kwargs):
    """
    Get the observable state in physical space at specified sensor locations.
    Parameters
    ----------
    Nt : int
        Number of time steps to retrieve. Default is 1.
    loc : array-like or str, optional
        Sensor locations to retrieve observables from. If 'all', returns observables at all spatial points.
        If None, returns observables at the predefined sensor locations.
    """
    if loc is None:
        loc = self.sensor_locations
    elif loc.lower() == 'all':
        loc = np.arange(self.Nx)

    if Nt == 1:
        return KS.fourier_to_physical(self.hist[-1, :self.Nk])[loc]
    else:
        return KS.fourier_to_physical(self.hist[-Nt:, :self.Nk])[:, loc]

ETDRK4_step(u_hat, nonlinear_operator, E, E2, Q, f1, f2, f3) staticmethod

Standard Kassam-Trefethen ETDRK4 step:

a_n = exp(L h / 2) u_n + Q N(u_n) b_n = exp(L h / 2) u_n + Q N(a_n) c_n = exp(L h / 2) a_n + Q (2 N(b_n) - N(u_n))

u_{n+1} = exp(L h) u_n + f1 N(u_n) + 2 f2 (N(a_n) + N(b_n)) + f3 N(c_n)

where h is the timestep, L the (diagonal) linear operator, N the nonlinear operator, and Q, f1, f2, f3 the contour-integrated phi-function coefficients (see ETDRK4_f_terms). The linear part is integrated exactly; the nonlinear terms with fourth-order accuracy.

Source code in dynamodels/physical/kuramoto_sivashinsky.py

@staticmethod
def ETDRK4_step(u_hat, nonlinear_operator, E, E2, Q, f1, f2, f3):
    """
    Standard Kassam-Trefethen ETDRK4 step:

    a_n = exp(L h / 2) u_n + Q N(u_n)
    b_n = exp(L h / 2) u_n + Q N(a_n)
    c_n = exp(L h / 2) a_n + Q (2 N(b_n) - N(u_n))

    u_{n+1} = exp(L h) u_n + f1 N(u_n) + 2 f2 (N(a_n) + N(b_n)) + f3 N(c_n)

    where h is the timestep, L the (diagonal) linear operator, N the
    nonlinear operator, and Q, f1, f2, f3 the contour-integrated
    phi-function coefficients (see ETDRK4_f_terms). The linear part is
    integrated exactly; the nonlinear terms with fourth-order accuracy.
    """

    N1 = nonlinear_operator(u_hat)
    a = E2 * u_hat + Q * N1
    N2 = nonlinear_operator(a)
    b = E2 * u_hat + Q * N2
    N3 = nonlinear_operator(b)
    c = E2 * a + Q * (2 * N3 - N1)
    N4 = nonlinear_operator(c)

    return E * u_hat + f1 * N1 + 2 * f2 * (N2 + N3) + f3 * N4

time_step(Nt=10, averaged=False, alpha=None)

Integrator for the KS model that supports ensembles and averaged ensemble propagation. Matches interface conventions of other models.

Parameters:

Name	Type	Description	Default
Nt	int	Number of time steps to integrate.	10
averaged	bool	If True, integrates the mean state and broadcasts ensemble deviations.	False
alpha	optional	Additional model parameters.	None

Returns:

Name	Type	Description
psi	ndarray	Forecasted state array of shape (Nt, Nphi, m).
t	ndarray	Time vector corresponding to each forecasted state.

Source code in dynamodels/physical/kuramoto_sivashinsky.py

def time_step(self, Nt=10, averaged=False, alpha=None):
    """
    Integrator for the KS model that supports ensembles and averaged ensemble propagation.
    Matches interface conventions of other models.

    Parameters
    ----------
    Nt : int
        Number of time steps to integrate.
    averaged : bool, optional
        If True, integrates the mean state and broadcasts ensemble deviations.
    alpha : optional
        Additional model parameters.

    Returns
    -------
    psi : np.ndarray
        Forecasted state array of shape (Nt, Nphi, m).
    t : np.ndarray
        Time vector corresponding to each forecasted state.
    """

    u0_hat = self.current_state

    if u0_hat.ndim == 1:  # reshape for non-ensemble
        u0_hat = u0_hat[:, None]

    t = np.round(self.current_time + np.arange(Nt + 1) * self.dt, self.precision_t)


    if averaged and self.ensemble:
        u0_hat_mean = np.mean(u0_hat, axis=1, keepdims=True)
        psi_deviation = u0_hat - u0_hat_mean

        psi_mean_arr = [u0_hat_mean[:, 0]]
        for _ in range(Nt):
            psi_mean_arr.append(KS.ETDRK4_step(psi_mean_arr[-1][:, None], **self.ETDRK4_f_terms)[:, 0])
        psi_mean_arr = np.stack(psi_mean_arr, axis=0)  # (Nt+1, N_x)

        # Broadcast deviations
        psi = np.array([psi_mean_arr[ii][:, None] + psi_deviation for ii in range(psi_mean_arr.shape[0])])  # (Nt+1, N_x, m)

    else:
        # Single member integration
        psi = [u0_hat]
        for _ in range(Nt):
            psi.append(KS.ETDRK4_step(psi[-1], **self.ETDRK4_f_terms))

        psi = np.stack(psi, axis=0)


    return psi, t

get_energy(Nt=0, u=None)

Compute the L2 energy of the solution: E = (1/L) * integral(u^2)dx

Returns:

float L2 energy

Source code in dynamodels/physical/kuramoto_sivashinsky.py

def get_energy(self, Nt=0, u=None):
    """
    Compute the L2 energy of the solution: E = (1/L) * integral(u^2)dx

    Returns:
    --------
    float
        L2 energy
    """

    if u is None:
        u = self.get_observable_hist(Nt=Nt, loc="all")


    if u.ndim == 2:
        u = u[np.newaxis, :]

    assert u.shape[1] == self.Nx

    return np.mean(u**2, axis=1)

get_enstrophy(Nt=0, u_hat=None)

Compute the enstrophy (integral of (u_x)^2).

Returns:

float Enstrophy

Source code in dynamodels/physical/kuramoto_sivashinsky.py

def get_enstrophy(self, Nt=0, u_hat=None):
    """
    Compute the enstrophy (integral of (u_x)^2).

    Returns:
    --------
    float
        Enstrophy
    """

    if u_hat is None:
        if Nt != 1:
            u_hat = self.hist[-Nt:]
        else:
            u_hat = self.current_state[np.newaxis, :]
    else:
        if u_hat.ndim == 2:
            u_hat = u_hat[np.newaxis, :]

    assert u_hat.shape[1] == self.k.shape[0]

    u_x_hat = self.first_derivative_x(u_hat)
    u_x = self.fourier_to_physical(u_x_hat)

    return np.mean(u_x**2, axis=1)

visualize_spatiotemporal_hist(y_hist=None, t=None, nrows=None, averaged=False, **kwargs)

Visualize the spatiotemporal evolution of the KS model in the physical space.

Source code in dynamodels/physical/kuramoto_sivashinsky.py

def visualize_spatiotemporal_hist(self, y_hist=None, t=None, nrows=None, averaged=False, **kwargs):
    """
    Visualize the spatiotemporal evolution of the KS model in the physical space.
    """

    if y_hist is None:
        y_hist = self.get_observable_hist(loc="all")

    if t is None:
        t = self.hist_t

    if not averaged:
        if nrows is None:
            nrows = min(10, y_hist.shape[-1])

        fig = plt.figure(figsize=(10, 1.5 * nrows))
        axs = fig.subplots(nrows=nrows, sharey=True, sharex=True)
        if nrows == 1:
            axs = [axs]

        lim = np.max(abs(y_hist))

        for mi, ax in enumerate(axs):
            im = ax.imshow(y_hist[:, :, mi].T,
                        aspect='auto', origin='lower',
                        cmap='RdBu_r', vmin=-lim, vmax=lim,
                        extent=[t[0], t[-1], self.x[0], self.x[-1]])  # TRANSPOSE


        axs[0].set(title=rf"KS spatiotemporal evolution. $L={self.L/np.pi:.2f}\pi, \nu={self.nu}$")
        axs[-1].set(xlabel="$t$")

        fig.colorbar(im, ax=axs, orientation='vertical', shrink=1/nrows)  #type: ignore
    else:
        # Averaged ensemble visualization
        y_mean_hist = np.mean(y_hist, axis=-1)

        fig, axs = plt.subplots(nrows=2, figsize=(10, 6), sharex=True)

        # Mean evolution
        lim_mean = np.max(abs(y_mean_hist))
        im0 = axs[0].imshow(y_mean_hist.T,
                            aspect='auto', origin='lower',
                            cmap='RdBu_r', vmin=-lim_mean, vmax=lim_mean,
                            extent=[t[0], t[-1], self.x[0], self.x[-1]])
        axs[0].set(title=rf"KS averaged spatiotemporal evolution (mean and std). $L={self.L/np.pi:.2f}\pi, \nu={self.nu}$") #type: ignore
        fig.colorbar(im0, ax=axs[0], orientation='vertical')

        # Deviation covariance evolution

        var_ensemble = np.var(y_hist, axis=-1, ddof=1).T            # (Nt, Nx)
        var_ensemble = np.sqrt(var_ensemble)                     # Standard deviation

        lim_dev = np.max(abs(var_ensemble))
        im1 = axs[1].imshow(var_ensemble,  # Plot covariance of deviations
                            aspect='auto', origin='lower',
                            cmap='magma', vmin=0, vmax=lim_dev,
                            extent=[t[0], t[-1], self.x[0], self.x[-1]])

        fig.colorbar(im1, ax=axs[1], orientation='vertical')
    # add the ticks and labels

    # Set spatial ticks as multiples of L
    assert self.L is not None, "L must be defined to set spatial ticks."
    ticks = (np.arange(4) + 1)* self.L/4
    tick_labels = [r"$L/4$", r"$L/2$", r"$3L/4$",r"$L$"]
    for ax in axs:
        ax.set(ylabel="$x$", yticks=ticks, yticklabels=tick_labels)

dynamodels.physical.rijke.Rijke(**model_dict)

Bases: Model

Rijke tube — longitudinal thermoacoustic low-order model.

The acoustic velocity and pressure perturbations are expanded on \(N_m\) Galerkin modes with wavenumbers \(k_j = j\pi/L\),

\[ u'(x, t) = \sum_{j=1}^{N_m} \eta_j(t) \cos(k_j x), \qquad p'(x, t) = -\sum_{j=1}^{N_m} \mu_j(t) \sin(k_j x), \]

giving the modal ODEs

\[ \dot{\eta}_j = \frac{k_j}{\bar\rho}\, \mu_j, \qquad \dot{\mu}_j = -k_j\, \bar\rho\, \bar{c}^2\, \eta_j - \frac{\bar{c}}{L}\, \zeta_j\, \mu_j + \dot{q}_j, \qquad \zeta_j = C_1\, j^2 + C_2\, \sqrt{j}, \]

where \(\bar\rho\), \(\bar{c}\) (and \(\bar{u}\), \(\bar{p}\), \(\bar\gamma\) below) are fixed mean-flow properties, weight-averaged across the temperature jump at the flame location \(x_f\), and \(\zeta_j\) is the modal damping. The heat release is projected onto the modes as

\[ \dot{q}_j = -\frac{2 (\bar\gamma - 1)}{L} \sin(k_j x_f)\, \dot{q}'(x_f, t), \]

with a gain–delay law relating \(\dot{q}'\) to the (time-delayed) acoustic velocity at the flame, \(u_f(t) \equiv u'(x_f, t - \tau)\): a square-root law (law='sqrt')

\[ \dot{q}'(x_f, t) = \bar{p}\, \bar{u}\, \beta \left[ \sqrt{\left| \tfrac{1}{3} + u_f(t) / \bar{u} \right|} - \sqrt{\tfrac{1}{3}} \right], \]

or a saturating arctangent law (law='tan')

\[ \dot{q}'(x_f, t) = \beta \sqrt{\beta / \kappa}\, \arctan\!\left( \sqrt{\beta / \kappa}\, u_f(t) \right). \]

The delay \(\tau\) is realized by advecting \(u'(x_f, t)\) along an auxiliary field discretized with \(N_c\) Chebyshev collocation points, and interpolating it at the point corresponding to the elapsed delay to obtain \(u_f(t)\).

The estimable parameters are \(\beta\), \(\tau\), the damping coefficients \(C_1\), \(C_2\), and \(\kappa\) (only active for law='tan'). The observables are the pressure at Nq microphone locations.

References

Nóvoa & Magri (2022). Real-time thermoacoustic data assimilation. J. Fluid Mech., 948, A35. DOI: 10.1017/jfm.2022.653.

Source code in dynamodels/physical/rijke.py

def __init__(self, **model_dict):

    if 'psi0' not in model_dict.keys():
        if 'Nm' in model_dict.keys():
            Nm = model_dict['Nm']
        else:
            Nm = self.Nm
        if 'Nc' in model_dict.keys():
            Nc = model_dict['Nc']
        else:
            Nc = self.Nc
        model_dict['psi0'] = .05 * np.hstack([np.ones(2 * Nm), np.zeros(Nc)])

    dt = model_dict.pop('dt', 1E-4)

    self.tau_adv = self.tau


    self.alpha_labels = dict(beta='$\\beta$', tau='$\\tau$', C1='$C_1$', C2='$C_2$', kappa='$\\kappa$')
    self.alpha_lims =  dict(beta=(0.01, 5), tau=(1E-6, self.tau_adv), C1=(0., 1.), C2=(0., 1.), kappa=(1E3, 1E8))


    # Chebyshev modes
    self.Dc, self.gc = Cheb(self.Nc, getg=True)

    # Microphone locations
    self.x_mic = np.linspace(self.xf, self.L, self.Nq + 1)[:-1]

    # Define modes frequency of each mode and sin cos etc
    jj = np.arange(1, self.Nm + 1)
    self.jpiL = jj * np.pi / self.L
    self.sinomjxf = np.sin(self.jpiL * self.xf)
    self.cosomjxf = np.cos(self.jpiL * self.xf)

    # Mean Flow Properties
    def weight_avg(y1, y2):
        return self.xf / self.L * y1 + (1. - self.xf / self.L) * y2

    self.meanFlow = dict(u=weight_avg(10, 11.1643), p=101300.,
                         gamma=1.4, T=weight_avg(300, 446.5282), R=287.1)
    self.meanFlow['rho'] = self.meanFlow['p'] / (self.meanFlow['R'] * self.meanFlow['T'])
    self.meanFlow['c'] = np.sqrt(self.meanFlow['gamma'] * self.meanFlow['R'] * self.meanFlow['T'])

    super().__init__(dt=dt, integrator_class=IVPIntegrator, **model_dict)

    # measured 1/lambda1 at this beta (set once: params change by re-instantiation)
    self.t_lyap = self.t_lyap_from_table(self.beta, _LAM1_MEASURED, Rijke.t_lyap)

time_derivative(t, psi, C1, C2, beta, kappa, tau, cosomjxf, Dc, gc, jpiL, L, law, meanFlow, Nc, Nm, tau_adv, sinomjxf) staticmethod

Time derivative of the Rijke tube governing equations (see class docstring).

Parameters:

Name	Type	Description	Default
t	float	Current time.	required
psi	ndarray	Augmented state vector; the first 2 * Nm + Nc entries are [eta (Nm,), mu (Nm,), v (Nc,)].	required
C1	float	Modal damping coefficients, \(\zeta_j = C_1 j^2 + C_2 \sqrt{j}\).	required
C2	float	Modal damping coefficients, \(\zeta_j = C_1 j^2 + C_2 \sqrt{j}\).	required
beta	float	Heat-release intensity.	required
kappa	float	Saturation parameter used by the 'tan' heat-release law.	required
tau	float	Time delay of the flame response.	required
cosomjxf	ndarray	Precomputed \(\cos(k_j x_f)\), \(\sin(k_j x_f)\) for each mode \(j\), used respectively to evaluate \(u'(x_f, t)\) and to project the heat release onto the \(\mu\) modes.	required
sinomjxf	ndarray	Precomputed \(\cos(k_j x_f)\), \(\sin(k_j x_f)\) for each mode \(j\), used respectively to evaluate \(u'(x_f, t)\) and to project the heat release onto the \(\mu\) modes.	required
Dc	ndarray	Chebyshev differentiation matrix and collocation points used to advect the delay line.	required
gc	ndarray	Chebyshev differentiation matrix and collocation points used to advect the delay line.	required
jpiL	ndarray	Modal wavenumbers \(k_j = j\pi/L\).	required
L	float	Tube length.	required
law	str	Heat-release law, 'sqrt' or 'tan'.	required
meanFlow	dict	Mean-flow properties at the flame (\(\bar\rho\), \(\bar u\), \(\bar p\), \(\bar c\), \(\bar\gamma\), \(\bar T\)).	required
Nc	int	Number of Chebyshev modes discretizing the delay line.	required
Nm	int	Number of Galerkin modes.	required
tau_adv	float	Reference advection time spanned by the delay line.	required

Returns:

Type	Description
ndarray	Concatenated time derivative of the augmented state vector.

Source code in dynamodels/physical/rijke.py

@staticmethod
def time_derivative(t, psi,
                    C1, C2, beta, kappa, tau,
                    cosomjxf, Dc, gc, jpiL, L, law, meanFlow, Nc, Nm, tau_adv, sinomjxf):
    r"""Time derivative of the Rijke tube governing equations (see class docstring).

    Parameters
    ----------
    t : float
        Current time.
    psi : ndarray
        Augmented state vector; the first ``2 * Nm + Nc`` entries are
        ``[eta (Nm,), mu (Nm,), v (Nc,)]``.
    C1, C2 : float
        Modal damping coefficients, $\zeta_j = C_1 j^2 + C_2 \sqrt{j}$.
    beta : float
        Heat-release intensity.
    kappa : float
        Saturation parameter used by the ``'tan'`` heat-release law.
    tau : float
        Time delay of the flame response.
    cosomjxf, sinomjxf : ndarray
        Precomputed $\cos(k_j x_f)$, $\sin(k_j x_f)$ for each mode $j$, used
        respectively to evaluate $u'(x_f, t)$ and to project the heat release
        onto the $\mu$ modes.
    Dc, gc : ndarray
        Chebyshev differentiation matrix and collocation points used to
        advect the delay line.
    jpiL : ndarray
        Modal wavenumbers $k_j = j\pi/L$.
    L : float
        Tube length.
    law : str
        Heat-release law, ``'sqrt'`` or ``'tan'``.
    meanFlow : dict
        Mean-flow properties at the flame ($\bar\rho$, $\bar u$, $\bar p$,
        $\bar c$, $\bar\gamma$, $\bar T$).
    Nc : int
        Number of Chebyshev modes discretizing the delay line.
    Nm : int
        Number of Galerkin modes.
    tau_adv : float
        Reference advection time spanned by the delay line.

    Returns
    -------
    ndarray
        Concatenated time derivative of the augmented state vector.
    """
    eta, mu, v = psi[:Nm], psi[Nm: 2 * Nm], psi[2 * Nm: 2 * Nm + Nc]

    # Advection equation boundary conditions
    v2 = np.hstack((np.dot(eta, cosomjxf), v))

    # Evaluate u(t_interp-tau) i.e. velocity at the flame at t_interp - tau
    x_tau = tau / tau_adv
    if x_tau < 1:
        f = splrep(gc, v2)
        u_tau = splev(x_tau, f)
    elif x_tau == 1:  # if no tau estimation, bypass interpolation to speed up code
        u_tau = v2[-1]
    else:
        raise Exception(f"tau = {tau} can't_interp be larger than tau_adv = {tau_adv}")

    # Compute damping and heat release law
    zeta = C1 * (jpiL * L / np.pi) ** 2 + C2 * (jpiL * L / np.pi) ** .5

    MF = meanFlow.copy()  # Physical properties
    if law == 'sqrt':
        q_dot = MF['p'] * MF['u'] * beta * (
                np.sqrt(abs(1. / 3 + u_tau / MF['u'])) - np.sqrt(1. / 3))  # [W/m2]=[m/s3]
    elif law == 'tan':
        q_dot = beta * np.sqrt(beta / kappa) * np.arctan(np.sqrt(beta / kappa) * u_tau)  # [m / s3]
    else:
        raise ValueError(f'Law "{law}" not defined')
    q_dot *= -2. * (MF['gamma'] - 1.) / L * sinomjxf  # [Pa/s]

    # governing equations
    deta_dt = jpiL / MF['rho'] * mu
    dmu_dt = - jpiL * MF['gamma'] * MF['p'] * eta - MF['c'] / L * zeta * mu + q_dot
    dv_dt = - 2. / tau_adv * np.dot(Dc, v2)

    return np.concatenate((deta_dt, dmu_dt, dv_dt[1:], np.zeros(len(psi) - (2 * Nm + Nc))))

visualize_spatiotemporal_hist(y_hist=None, t=None, nrows=None, averaged=False, reference_y=1.0, reference_t=1.0, **kwargs)

Visualize the spatiotemporal evolution of the Rijke tube model in the physical space.

Source code in dynamodels/physical/rijke.py

def visualize_spatiotemporal_hist(self, y_hist=None, t=None, nrows=None, averaged=False,
                                  reference_y=1.0, reference_t: float = 1.0, **kwargs):
    """
    Visualize the spatiotemporal evolution of the Rijke tube model in the physical space.
    """

    if y_hist is None:
        Nt = int(self.t_CR // self.dt)
        y_hist = self.get_observable_hist(loc="all", Nt=Nt)

    if t is None:
        t = self.hist_t

    (t,), t_lbl = normalized_time(reference_t, t)
    assert t is not None
    if reference_y != 1.0:
        y_hist = y_hist / reference_y


    # Set spatial ticks as multiples of L
    ticks = np.arange(5)* self.L/4
    # tick_labels = [r"$L/4$", r"$L/2$", r"$3L/4$",r"$L$"]

    if not averaged:
        if nrows is None:
            nrows = min(10, y_hist.shape[-1])

        fig = plt.figure(figsize=(10, 1.5 * nrows))
        axs = fig.subplots(nrows=nrows, sharey=True, sharex=True)
        if nrows == 1:
            axs = [axs]

        lim = np.max(abs(y_hist))
        for mi, ax in enumerate(axs):
            im = ax.imshow(y_hist[:, :, mi].T,
                        aspect='auto', origin='lower',
                        cmap='RdBu_r',
                        vmin=-lim, vmax=lim,
                        extent=[t[0], t[-1], 0, self.L])
            ax.set(ylabel="$x$")
            ax.set_yticks(ticks)
            # ax.set_yticklabels(tick_labels)
        fig.colorbar(im, ax=axs, orientation='vertical', shrink=1/nrows) #type: ignore


        axs[0].set(title=rf"Rijke spatiotemporal evolution $x_f={self.xf}$")
        axs[-1].set(xlabel=t_lbl)

    else:
        # Averaged ensemble visualization
        y_mean_hist = np.mean(y_hist, axis=-1)

        fig, axs = plt.subplots(nrows=2, figsize=(10, 6), sharex=True)

        # Mean evolution
        lim_mean = np.max(abs(y_mean_hist))
        im0 = axs[0].imshow(y_mean_hist.T,
                            aspect='auto', origin='lower',
                            cmap='RdBu_r', vmin=-lim_mean, vmax=lim_mean,
                            extent=[t[0], t[-1], 0, self.L])

        axs[0].set(title=rf"Rijke spatiotemporal evolution (mean and std) $x_f={self.xf}$")
        fig.colorbar(im0, ax=axs[0], orientation='vertical')
        # Deviation covariance evolution in percentage

        var_ensemble = np.var(y_hist, axis=-1, ddof=1).T            # (Nt, Nx)
        var_ensemble = np.sqrt(var_ensemble)                     # Standard deviation
        lim_dev = np.max(var_ensemble)
        im1 = axs[1].imshow(var_ensemble,  # Plot covariance of deviations
                            aspect='auto', origin='lower',
                            cmap='magma', vmin=0, vmax=lim_dev,
                            extent=[t[0], t[-1], 0, self.L])

        axs[1].set(xlabel=t_lbl)
        fig.colorbar(im1, ax=axs[1], orientation='vertical')


        for ax in axs:
            ax.set(yticks=ticks, ylabel="$x$")

Pressure field of the Rijke-tube low-order model.

dynamodels.physical.annular.Annular(**model_dict)

Bases: Model

Annular combustor — two coupled oscillators for the first azimuthal modes.

The acoustic pressure in the annulus obeys the wave equation with heat-release source and resistive/reactive asymmetries,

\[ \frac{\partial^2 p}{\partial t^2} + \zeta \frac{\partial p}{\partial t} - \left[ 1 + \epsilon \cos\!\big(2(\theta - \Theta_\epsilon)\big) \right] \frac{c^2}{r^2} \frac{\partial^2 p}{\partial \theta^2} = (\gamma - 1)\, \frac{\partial \dot{q}}{\partial t}, \qquad (\gamma - 1)\, \dot{q} = \beta \left[ 1 + c_2 \cos\!\big(2(\theta - \Theta_\beta)\big) \right] p - \kappa p^3 . \]

Decomposing the pressure field onto the first azimuthal mode pair (\(n = 1\)),

\[ p(\theta, t) = \eta_a(t) \cos(n\theta) + \eta_b(t) \sin(n\theta), \]

yields four coupled first-order ODEs for \((\eta_a, \dot{\eta}_a, \eta_b, \dot{\eta}_b)\):

\[ \ddot{\eta}_a = -\omega^2 \left[ \eta_a \big(1 + \tfrac{\epsilon}{2} \cos 2\Theta_\epsilon\big) + \eta_b \tfrac{\epsilon}{2} \sin 2\Theta_\epsilon \right] + \dot{\eta}_a \left[ 2\nu + \tfrac{c_2\beta}{2} \cos 2\Theta_\beta - \tfrac{3\kappa}{4} (3\eta_a^2 + \eta_b^2) \right] + \dot{\eta}_b \left[ \tfrac{c_2\beta}{2} \sin 2\Theta_\beta - \tfrac{3\kappa}{2} \eta_a \eta_b \right], \]

\[ \ddot{\eta}_b = -\omega^2 \left[ \eta_b \big(1 - \tfrac{\epsilon}{2} \cos 2\Theta_\epsilon\big) + \eta_a \tfrac{\epsilon}{2} \sin 2\Theta_\epsilon \right] + \dot{\eta}_b \left[ 2\nu - \tfrac{c_2\beta}{2} \cos 2\Theta_\beta - \tfrac{3\kappa}{4} (3\eta_b^2 + \eta_a^2) \right] + \dot{\eta}_a \left[ \tfrac{c_2\beta}{2} \sin 2\Theta_\beta - \tfrac{3\kappa}{2} \eta_a \eta_b \right]. \]

The estimable parameters are the growth rate \(\nu\), the resistive-asymmetry intensity \(c_2\beta\), the saturation \(\kappa\), the reactive-asymmetry amplitude \(\epsilon\) and phase \(\Theta_\epsilon\), the frequency \(\omega\) and the direction of maximum r.m.s. pressure \(\Theta_\beta\).

Example dynamical regimes:

• purely spinning mode: \((\nu, c_2\beta) = (30, 5)\);
• purely standing mode: \((\nu, c_2\beta) = (0, 50)\);
• mixed mode: \((\nu, c_2\beta) = (20, 18)\).

References

Nóvoa, Noiray, Dawson & Magri (2024). A real-time digital twin of azimuthal thermoacoustic instabilities. J. Fluid Mech., 1001, A49. DOI: 10.1017/jfm.2024.1052.

Source code in dynamodels/physical/annular.py

def __init__(self, **model_dict):

    dt = model_dict.pop('dt', 1. / 51200)
    psi0 = model_dict.pop('psi0', None)
    if psi0 is None:
        C0, X0, th0, ph0 = 10, 0, 0.63, 0  # %initial values
        # Conversion of the initial conditions from the quaternion formalism to the AB formalism
        Ai = C0 * np.sqrt(np.cos(th0) ** 2 * np.cos(X0) ** 2 + np.sin(th0) ** 2 * np.sin(X0) ** 2)
        Bi = C0 * np.sqrt(np.sin(th0) ** 2 * np.cos(X0) ** 2 + np.cos(th0) ** 2 * np.sin(X0) ** 2)
        phai = ph0 + np.arctan2(np.sin(th0) * np.sin(X0), np.cos(th0) * np.cos(X0))
        phbi = ph0 - np.arctan2(np.cos(th0) * np.sin(X0), np.sin(th0) * np.cos(X0))

        # %initial conditions for the fast oscillator equations
        psi0 = [Ai * np.cos(phai),
                -self.omega * Ai * np.sin(phai),
                Bi * np.cos(phbi),
                -self.omega * Bi * np.sin(phbi)]

        psi0 = np.array(psi0)  # initialise \eta_a, \dot{\eta_a}, \eta_b, \dot{\eta_b}

    super().__init__(psi0=psi0, dt=dt, integrator_class=IVPIntegrator, **model_dict)

    self.alpha_labels = dict(omega='$\\omega$', nu='$\\nu$', c2beta='$c_2\\beta $', kappa='$\\kappa$',
                             epsilon='$\\epsilon$', theta_b='$\\Theta_\\beta$', theta_e='$\\Theta_\\epsilon$')

    self.alpha_lims =  dict(omega=(1000 * 2 * np.pi, 1300 * 2 * np.pi),
                            nu=(-60., 100.), c2beta=(0., 100.),
                            theta_b=(0, 2 * np.pi), theta_e=(0, 2 * np.pi))

get_observables(Nt=1, loc=None, measure_modes=False, **kwargs)

pressure measurements at theta = [0º, 60º, 120º, 240º]

Source code in dynamodels/physical/annular.py

def get_observables(self, Nt=1, loc=None, measure_modes=False, **kwargs):
    """
    pressure measurements at theta = [0º, 60º, 120º, 240º]
    """
    if loc is None:
        loc = self.theta_mic

    if measure_modes:
        return self.hist[-Nt:, [0, 2], :]
    else:
        eta1, eta2 = self.hist[-Nt:, 0, :], self.hist[-Nt:, 2, :]
        if max(loc) > 2 * np.pi:
            raise ValueError('Theta must be in radians')

        p_mics = np.array([eta1 * np.cos(th) + eta2 * np.sin(th) for th in np.array(loc)])
        p_mics = p_mics.transpose(1, 0, 2)
        if Nt == 1:
            return p_mics.squeeze(axis=0)
        else:
            return p_mics

time_derivative(t, psi, nu, kappa, c2beta, theta_b, omega, epsilon, theta_e) staticmethod

Time derivative of the two coupled azimuthal oscillators (see class docstring).

Source code in dynamodels/physical/annular.py

@staticmethod
def time_derivative(t, psi, nu, kappa, c2beta, theta_b, omega, epsilon, theta_e):
    """Time derivative of the two coupled azimuthal oscillators (see class docstring)."""
    y_a, z_a, y_b, z_b = psi[:4]  # y = η, and z = dη/dt

    def k1(y1, y2, sign):
        return (2 * nu - 3. / 4 * kappa * (3 * y1 ** 2 + y2 ** 2) +
                sign * c2beta / 2. * np.cos(2. * theta_b))

    k2 = c2beta / 2. * np.sin(2. * theta_b) - 3. / 2 * kappa * y_a * y_b

    def k3(y1, y2, sign):
        return omega ** 2 * (y1 * (1 + sign * epsilon / 2. * np.cos(2. * theta_e)) +
                             y2 * epsilon / 2. * np.sin(2. * theta_e))

    dz_a = z_a * k1(y_a, y_b, sign=1) + z_b * k2 - k3(y_a, y_b, sign=1)
    dz_b = z_b * k1(y_b, y_a, sign=-1) + z_a * k2 - k3(y_b, y_a, sign=-1)

    return (z_a, dz_a, z_b, dz_b) + (0,) * (len(psi) - 4)

Azimuthal thermoacoustic simulation of the annular combustor model (ν=20, c₂β=10).