PyMPDATA_examples.Jarecka_et_al_2015.formulae 

 1# adapted from:
 2# https://github.com/igfuw/shallow-water-elliptic-drop/blob/master/analytical/analytic_equations.py
 3# Code used in the paper of Jarecka, Jaruga, Smolarkiewicz -
 4# "A Spreading Drop of Shallow Water" (JCP 289, doi:10.1016/j.jcp.2015.02.003).
 5
 6import numba
 7import numpy as np
 8from scipy.integrate import odeint
 9
10
11def amplitude(x, y, lx, ly):
12    A = 1 / lx / ly
13    h = A * (1 - (x / lx) ** 2 - (y / ly) ** 2)
14    return np.where(h > 0, h, 0)
15
16
17@numba.njit()
18def deriv(y, _):
19    """
20    return derivatives of [lambda_x, dlambda_x/dt, lambda_y, dlambda_y/dt
21    four first-order ODEs based on  eq. 7  (Jarecka, Jaruga, Smolarkiewicz)
22    """
23    return np.array((y[1], 2.0 / y[0] ** 2 / y[2], y[3], 2.0 / y[0] / y[2] ** 2))
24
25
26def d2_el_lamb_lamb_t_evol(times, lamb_x0, lamb_y0):
27    """
28    solving coupled nonlinear second-order ODEs - eq. 7  (Jarecka, Jaruga, Smolarkiewicz)
29    returning array with first dim denoting time, second dim:
30    [lambda_x, dot{lambda_x}, lambda_y, dot{lambda_y}
31    """
32    assert times[0] == 0
33    yinit = np.array([lamb_x0, 0.0, lamb_y0, 0.0])  # initial values (dot_lamb = 0.)
34    result, info = odeint(deriv, yinit, times, full_output=True)
35    assert info["message"] == "Integration successful."
36    return result
def amplitude(x, y, lx, ly): 
12def amplitude(x, y, lx, ly):
13    A = 1 / lx / ly
14    h = A * (1 - (x / lx) ** 2 - (y / ly) ** 2)
15    return np.where(h > 0, h, 0)
@numba.njit()
def deriv(y, _): 
18@numba.njit()
19def deriv(y, _):
20    """
21    return derivatives of [lambda_x, dlambda_x/dt, lambda_y, dlambda_y/dt
22    four first-order ODEs based on  eq. 7  (Jarecka, Jaruga, Smolarkiewicz)
23    """
24    return np.array((y[1], 2.0 / y[0] ** 2 / y[2], y[3], 2.0 / y[0] / y[2] ** 2))

return derivatives of [lambda_x, dlambda_x/dt, lambda_y, dlambda_y/dt four first-order ODEs based on eq. 7 (Jarecka, Jaruga, Smolarkiewicz)

def d2_el_lamb_lamb_t_evol(times, lamb_x0, lamb_y0): 
27def d2_el_lamb_lamb_t_evol(times, lamb_x0, lamb_y0):
28    """
29    solving coupled nonlinear second-order ODEs - eq. 7  (Jarecka, Jaruga, Smolarkiewicz)
30    returning array with first dim denoting time, second dim:
31    [lambda_x, dot{lambda_x}, lambda_y, dot{lambda_y}
32    """
33    assert times[0] == 0
34    yinit = np.array([lamb_x0, 0.0, lamb_y0, 0.0])  # initial values (dot_lamb = 0.)
35    result, info = odeint(deriv, yinit, times, full_output=True)
36    assert info["message"] == "Integration successful."
37    return result

solving coupled nonlinear second-order ODEs - eq. 7 (Jarecka, Jaruga, Smolarkiewicz) returning array with first dim denoting time, second dim: [lambda_x, dot{lambda_x}, lambda_y, dot{lambda_y}