dapper.mods.Lorenz96.extras

Extra functionality (not necessary for the EnKF or the particle filter).

 1"""Extra functionality (not necessary for the EnKF or the particle filter)."""
 2
 3import numpy as np
 4
 5import dapper.tools.liveplotting as LP
 6from dapper.mods.integration import integrate_TLM
 7
 8
 9def d2x_dtdx(x):
10    """Tangent linear model (TLM). I.e. the Jacobian of dxdt(x)."""
11    M = len(x)
12    F = np.zeros((M, M))
13
14    def md(i):
15        return np.mod(i, M)  # modulo
16
17    for i in range(M):
18        F[i, i] = -1.0
19        F[i, i - 2] = -x[i - 1]
20        F[i, md(i + 1)] = +x[i - 1]
21        F[i, i - 1] = x[md(i + 1)] - x[i - 2]
22
23    return F
24
25
26# dstep_dx = FD_Jac(step)
27def dstep_dx(x, t, dt):
28    """Compute resolvent (propagator) of the TLM. I.e. the Jacobian of `step(x)`."""
29    # For L96, method='analytic' >> 'approx'
30    return integrate_TLM(d2x_dtdx(x), dt, method="analytic")
31
32
33# Add some non-default liveplotters
34def LPs(jj=None):
35    return [
36        (1, LP.spatial1d(jj)),
37        (1, LP.correlations),
38        (0, LP.spectral_errors),
39        (0, LP.phase_particles(True, jj)),
40        (0, LP.sliding_marginals(jj)),
41    ]
def d2x_dtdx(x):
10def d2x_dtdx(x):
11    """Tangent linear model (TLM). I.e. the Jacobian of dxdt(x)."""
12    M = len(x)
13    F = np.zeros((M, M))
14
15    def md(i):
16        return np.mod(i, M)  # modulo
17
18    for i in range(M):
19        F[i, i] = -1.0
20        F[i, i - 2] = -x[i - 1]
21        F[i, md(i + 1)] = +x[i - 1]
22        F[i, i - 1] = x[md(i + 1)] - x[i - 2]
23
24    return F

Tangent linear model (TLM). I.e. the Jacobian of dxdt(x).

def dstep_dx(x, t, dt):
28def dstep_dx(x, t, dt):
29    """Compute resolvent (propagator) of the TLM. I.e. the Jacobian of `step(x)`."""
30    # For L96, method='analytic' >> 'approx'
31    return integrate_TLM(d2x_dtdx(x), dt, method="analytic")

Compute resolvent (propagator) of the TLM. I.e. the Jacobian of step(x).

def LPs(jj=None):
35def LPs(jj=None):
36    return [
37        (1, LP.spatial1d(jj)),
38        (1, LP.correlations),
39        (0, LP.spectral_errors),
40        (0, LP.phase_particles(True, jj)),
41        (0, LP.sliding_marginals(jj)),
42    ]