PySDM_examples.Szumowski_et_al_1998.fields

View Source

 1import numpy as np
 2
 3from PySDM.impl.arakawa_c import z_scalar_coord
 4
 5
 6def z_vec_coord(grid):
 7    nx = grid[0]
 8    nz = grid[1] + 1
 9    xX = (
10        np.repeat(np.linspace(1 / 2, grid[0] - 1 / 2, nx).reshape((nx, 1)), nz, axis=1)
11        / grid[0]
12    )
13    assert np.amin(xX) >= 0
14    assert np.amax(xX) <= 1
15    assert xX.shape == (nx, nz)
16    zZ = np.repeat(np.linspace(0, grid[1], nz).reshape((1, nz)), nx, axis=0) / grid[1]
17    assert np.amin(zZ) == 0
18    assert np.amax(zZ) == 1
19    assert zZ.shape == (nx, nz)
20    return xX, zZ
21
22
23def x_vec_coord(grid):
24    nx = grid[0] + 1
25    nz = grid[1]
26    xX = np.repeat(np.linspace(0, grid[0], nx).reshape((nx, 1)), nz, axis=1) / grid[0]
27    assert np.amin(xX) == 0
28    assert np.amax(xX) == 1
29    assert xX.shape == (nx, nz)
30    zZ = np.repeat(z_scalar_coord(grid).reshape((1, nz)), nx, axis=0) / grid[1]
31    assert np.amin(zZ) >= 0
32    assert np.amax(zZ) <= 1
33    assert zZ.shape == (nx, nz)
34    return xX, zZ
35
36
37def nondivergent_vector_field_2d(
38    grid: tuple, size: tuple, dt: float, stream_function: callable, t
39):
40    dx = size[0] / grid[0]
41    dz = size[1] / grid[1]
42    dxX = 1 / grid[0]
43    dzZ = 1 / grid[1]
44
45    xX, zZ = x_vec_coord(grid)
46    rho_velocity_x = (
47        -(stream_function(xX, zZ + dzZ / 2, t) - stream_function(xX, zZ - dzZ / 2, t))
48        / dz
49    )
50
51    xX, zZ = z_vec_coord(grid)
52    rho_velocity_z = (
53        stream_function(xX + dxX / 2, zZ, t) - stream_function(xX - dxX / 2, zZ, t)
54    ) / dx
55
56    rho_times_courant = [rho_velocity_x * dt / dx, rho_velocity_z * dt / dz]
57    return rho_times_courant

def z_vec_coord(grid): View Source

 7def z_vec_coord(grid):
 8    nx = grid[0]
 9    nz = grid[1] + 1
10    xX = (
11        np.repeat(np.linspace(1 / 2, grid[0] - 1 / 2, nx).reshape((nx, 1)), nz, axis=1)
12        / grid[0]
13    )
14    assert np.amin(xX) >= 0
15    assert np.amax(xX) <= 1
16    assert xX.shape == (nx, nz)
17    zZ = np.repeat(np.linspace(0, grid[1], nz).reshape((1, nz)), nx, axis=0) / grid[1]
18    assert np.amin(zZ) == 0
19    assert np.amax(zZ) == 1
20    assert zZ.shape == (nx, nz)
21    return xX, zZ

def x_vec_coord(grid): View Source

24def x_vec_coord(grid):
25    nx = grid[0] + 1
26    nz = grid[1]
27    xX = np.repeat(np.linspace(0, grid[0], nx).reshape((nx, 1)), nz, axis=1) / grid[0]
28    assert np.amin(xX) == 0
29    assert np.amax(xX) == 1
30    assert xX.shape == (nx, nz)
31    zZ = np.repeat(z_scalar_coord(grid).reshape((1, nz)), nx, axis=0) / grid[1]
32    assert np.amin(zZ) >= 0
33    assert np.amax(zZ) <= 1
34    assert zZ.shape == (nx, nz)
35    return xX, zZ

def nondivergent_vector_field_2d( grid: tuple, size: tuple, dt: float, stream_function: <built-in function callable>, t): View Source

38def nondivergent_vector_field_2d(
39    grid: tuple, size: tuple, dt: float, stream_function: callable, t
40):
41    dx = size[0] / grid[0]
42    dz = size[1] / grid[1]
43    dxX = 1 / grid[0]
44    dzZ = 1 / grid[1]
45
46    xX, zZ = x_vec_coord(grid)
47    rho_velocity_x = (
48        -(stream_function(xX, zZ + dzZ / 2, t) - stream_function(xX, zZ - dzZ / 2, t))
49        / dz
50    )
51
52    xX, zZ = z_vec_coord(grid)
53    rho_velocity_z = (
54        stream_function(xX + dxX / 2, zZ, t) - stream_function(xX - dxX / 2, zZ, t)
55    ) / dx
56
57    rho_times_courant = [rho_velocity_x * dt / dx, rho_velocity_z * dt / dz]
58    return rho_times_courant