PyMPDATA_examples.utils.nondivergent_vector_field_2d 

 1import numpy as np
 2
 3from PyMPDATA import VectorField
 4from PyMPDATA.boundary_conditions import Periodic
 5
 6
 7def nondivergent_vector_field_2d(grid, size, dt, stream_function: callable, halo):
 8    dx = size[0] / grid[0]
 9    dz = size[1] / grid[1]
10    dxX = 1 / grid[0]
11    dzZ = 1 / grid[1]
12
13    xX, zZ = x_vec_coord(grid)
14    rho_velocity_x = (
15        -(stream_function(xX, zZ + dzZ / 2) - stream_function(xX, zZ - dzZ / 2)) / dz
16    )
17
18    xX, zZ = z_vec_coord(grid)
19    rho_velocity_z = (
20        stream_function(xX + dxX / 2, zZ) - stream_function(xX - dxX / 2, zZ)
21    ) / dx
22
23    GC = [rho_velocity_x * dt / dx, rho_velocity_z * dt / dz]
24
25    # CFL condition
26    for val in GC:
27        np.testing.assert_array_less(np.abs(val), 1)
28
29    result = VectorField(GC, halo=halo, boundary_conditions=(Periodic(), Periodic()))
30
31    # nondivergence (of velocity field, hence dt)
32    assert np.amax(abs(result.div((dt, dt)).get())) < 5e-9
33
34    return result
35
36
37def x_vec_coord(grid):
38    nx = grid[0] + 1
39    nz = grid[1]
40    xX = np.repeat(np.linspace(0, grid[0], nx).reshape((nx, 1)), nz, axis=1) / grid[0]
41    assert np.amin(xX) == 0
42    assert np.amax(xX) == 1
43    assert xX.shape == (nx, nz)
44    zZ = (
45        np.repeat(np.linspace(1 / 2, grid[1] - 1 / 2, nz).reshape((1, nz)), nx, axis=0)
46        / grid[1]
47    )
48    assert np.amin(zZ) >= 0
49    assert np.amax(zZ) <= 1
50    assert zZ.shape == (nx, nz)
51    return xX, zZ
52
53
54def z_vec_coord(grid):
55    nx = grid[0]
56    nz = grid[1] + 1
57    xX = (
58        np.repeat(np.linspace(1 / 2, grid[0] - 1 / 2, nx).reshape((nx, 1)), nz, axis=1)
59        / grid[0]
60    )
61    assert np.amin(xX) >= 0
62    assert np.amax(xX) <= 1
63    assert xX.shape == (nx, nz)
64    zZ = np.repeat(np.linspace(0, grid[1], nz).reshape((1, nz)), nx, axis=0) / grid[1]
65    assert np.amin(zZ) == 0
66    assert np.amax(zZ) == 1
67    assert zZ.shape == (nx, nz)
68    return xX, zZ
def nondivergent_vector_field_2d(grid, size, dt, stream_function: <built-in function callable>, halo): 
 8def nondivergent_vector_field_2d(grid, size, dt, stream_function: callable, halo):
 9    dx = size[0] / grid[0]
10    dz = size[1] / grid[1]
11    dxX = 1 / grid[0]
12    dzZ = 1 / grid[1]
13
14    xX, zZ = x_vec_coord(grid)
15    rho_velocity_x = (
16        -(stream_function(xX, zZ + dzZ / 2) - stream_function(xX, zZ - dzZ / 2)) / dz
17    )
18
19    xX, zZ = z_vec_coord(grid)
20    rho_velocity_z = (
21        stream_function(xX + dxX / 2, zZ) - stream_function(xX - dxX / 2, zZ)
22    ) / dx
23
24    GC = [rho_velocity_x * dt / dx, rho_velocity_z * dt / dz]
25
26    # CFL condition
27    for val in GC:
28        np.testing.assert_array_less(np.abs(val), 1)
29
30    result = VectorField(GC, halo=halo, boundary_conditions=(Periodic(), Periodic()))
31
32    # nondivergence (of velocity field, hence dt)
33    assert np.amax(abs(result.div((dt, dt)).get())) < 5e-9
34
35    return result
def x_vec_coord(grid): 
38def x_vec_coord(grid):
39    nx = grid[0] + 1
40    nz = grid[1]
41    xX = np.repeat(np.linspace(0, grid[0], nx).reshape((nx, 1)), nz, axis=1) / grid[0]
42    assert np.amin(xX) == 0
43    assert np.amax(xX) == 1
44    assert xX.shape == (nx, nz)
45    zZ = (
46        np.repeat(np.linspace(1 / 2, grid[1] - 1 / 2, nz).reshape((1, nz)), nx, axis=0)
47        / grid[1]
48    )
49    assert np.amin(zZ) >= 0
50    assert np.amax(zZ) <= 1
51    assert zZ.shape == (nx, nz)
52    return xX, zZ
def z_vec_coord(grid): 
55def z_vec_coord(grid):
56    nx = grid[0]
57    nz = grid[1] + 1
58    xX = (
59        np.repeat(np.linspace(1 / 2, grid[0] - 1 / 2, nx).reshape((nx, 1)), nz, axis=1)
60        / grid[0]
61    )
62    assert np.amin(xX) >= 0
63    assert np.amax(xX) <= 1
64    assert xX.shape == (nx, nz)
65    zZ = np.repeat(np.linspace(0, grid[1], nz).reshape((1, nz)), nx, axis=0) / grid[1]
66    assert np.amin(zZ) == 0
67    assert np.amax(zZ) == 1
68    assert zZ.shape == (nx, nz)
69    return xX, zZ