PyMPDATA_examples.Smolarkiewicz_1984.settings

View Source

 1import numpy as np
 2from pystrict import strict
 3
 4
 5@strict
 6class Settings:
 7    def __init__(self, n: int, dt: float):
 8        self.grid = (n, n, n)
 9        self.dt = dt
10        self.L = 100
11        self.dx = self.L / n
12        self.dy = self.dx
13        self.dz = self.dx
14        self.h = 4
15        self.r = 15
16        d = 25 / np.sqrt(3)
17        self.x0 = 50 - d
18        self.y0 = 50 + d
19        self.z0 = 50 + d
20
21        self.omega = 0.1
22        self.xc = 50
23        self.yc = 50
24        self.zc = 50
25
26    @property
27    def advector(self):
28        """constant angular velocity rotational field"""
29
30        data = [None, None, None]
31        for index, letter in enumerate(("x", "y", "z")):
32            i, j, k = np.indices((g + (gi == index) for gi, g in enumerate(self.grid)))
33            if letter == "x":
34                data[index] = (
35                    -((j + 0.5) * self.dy - self.yc) + ((k + 0.5) * self.dz - self.zc)
36                ) / self.dx
37            elif letter == "y":
38                data[index] = (
39                    +((i + 0.5) * self.dx - self.xc) - ((k + 0.5) * self.dz - self.zc)
40                ) / self.dy
41            elif letter == "z":
42                data[index] = (
43                    -((i + 0.5) * self.dx - self.xc) + ((j + 0.5) * self.dy - self.yc)
44                ) / self.dz
45            data[index] *= self.omega / np.sqrt(3) * self.dt
46        return data
47
48    @property
49    def advectee(self):
50        i, j, k = np.indices(self.grid)
51        dist = (
52            ((i + 0.5) * self.dx - self.x0) ** 2
53            + ((j + 0.5) * self.dy - self.y0) ** 2
54            + ((k + 0.5) * self.dz - self.z0) ** 2
55        )
56        return np.where(dist - pow(self.r, 2) <= 0, self.h, 0)

@strict

class Settings: View Source

 6@strict
 7class Settings:
 8    def __init__(self, n: int, dt: float):
 9        self.grid = (n, n, n)
10        self.dt = dt
11        self.L = 100
12        self.dx = self.L / n
13        self.dy = self.dx
14        self.dz = self.dx
15        self.h = 4
16        self.r = 15
17        d = 25 / np.sqrt(3)
18        self.x0 = 50 - d
19        self.y0 = 50 + d
20        self.z0 = 50 + d
21
22        self.omega = 0.1
23        self.xc = 50
24        self.yc = 50
25        self.zc = 50
26
27    @property
28    def advector(self):
29        """constant angular velocity rotational field"""
30
31        data = [None, None, None]
32        for index, letter in enumerate(("x", "y", "z")):
33            i, j, k = np.indices((g + (gi == index) for gi, g in enumerate(self.grid)))
34            if letter == "x":
35                data[index] = (
36                    -((j + 0.5) * self.dy - self.yc) + ((k + 0.5) * self.dz - self.zc)
37                ) / self.dx
38            elif letter == "y":
39                data[index] = (
40                    +((i + 0.5) * self.dx - self.xc) - ((k + 0.5) * self.dz - self.zc)
41                ) / self.dy
42            elif letter == "z":
43                data[index] = (
44                    -((i + 0.5) * self.dx - self.xc) + ((j + 0.5) * self.dy - self.yc)
45                ) / self.dz
46            data[index] *= self.omega / np.sqrt(3) * self.dt
47        return data
48
49    @property
50    def advectee(self):
51        i, j, k = np.indices(self.grid)
52        dist = (
53            ((i + 0.5) * self.dx - self.x0) ** 2
54            + ((j + 0.5) * self.dy - self.y0) ** 2
55            + ((k + 0.5) * self.dz - self.z0) ** 2
56        )
57        return np.where(dist - pow(self.r, 2) <= 0, self.h, 0)

Settings(n: int, dt: float) View Source

 8    def __init__(self, n: int, dt: float):
 9        self.grid = (n, n, n)
10        self.dt = dt
11        self.L = 100
12        self.dx = self.L / n
13        self.dy = self.dx
14        self.dz = self.dx
15        self.h = 4
16        self.r = 15
17        d = 25 / np.sqrt(3)
18        self.x0 = 50 - d
19        self.y0 = 50 + d
20        self.z0 = 50 + d
21
22        self.omega = 0.1
23        self.xc = 50
24        self.yc = 50
25        self.zc = 50

grid

omega

advector View Source

27    @property
28    def advector(self):
29        """constant angular velocity rotational field"""
30
31        data = [None, None, None]
32        for index, letter in enumerate(("x", "y", "z")):
33            i, j, k = np.indices((g + (gi == index) for gi, g in enumerate(self.grid)))
34            if letter == "x":
35                data[index] = (
36                    -((j + 0.5) * self.dy - self.yc) + ((k + 0.5) * self.dz - self.zc)
37                ) / self.dx
38            elif letter == "y":
39                data[index] = (
40                    +((i + 0.5) * self.dx - self.xc) - ((k + 0.5) * self.dz - self.zc)
41                ) / self.dy
42            elif letter == "z":
43                data[index] = (
44                    -((i + 0.5) * self.dx - self.xc) + ((j + 0.5) * self.dy - self.yc)
45                ) / self.dz
46            data[index] *= self.omega / np.sqrt(3) * self.dt
47        return data

constant angular velocity rotational field

advectee View Source

49    @property
50    def advectee(self):
51        i, j, k = np.indices(self.grid)
52        dist = (
53            ((i + 0.5) * self.dx - self.x0) ** 2
54            + ((j + 0.5) * self.dy - self.y0) ** 2
55            + ((k + 0.5) * self.dz - self.z0) ** 2
56        )
57        return np.where(dist - pow(self.r, 2) <= 0, self.h, 0)