forked from lettucecfd/lettuce
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtest_moments.py
100 lines (77 loc) · 3.96 KB
/
test_moments.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
import numpy as np
import pytest
from lettuce.moments import *
from lettuce.stencils import *
from lettuce.lattices import Lattice
def test_moments_density_array(stencil):
rho_tensor = moment_tensor(stencil.e, np.array([0]*stencil.D()))
assert rho_tensor == pytest.approx(np.ones((stencil.Q())))
def test_more_moments_density_array(stencil):
rho_tensor = moment_tensor(stencil.e, np.array([[0]*stencil.D()]))
assert rho_tensor == pytest.approx(np.ones((1,stencil.Q())))
def test_moments_density_tensor(lattice):
rho_tensor = moment_tensor(lattice.e, lattice.convert_to_tensor(([0]*lattice.D)))
assert rho_tensor.shape == (lattice.Q,)
assert rho_tensor.cpu().numpy() == pytest.approx(np.ones((lattice.Q)))
def test_more_moments_density_tensor(lattice):
rho_tensor = moment_tensor(lattice.e, lattice.convert_to_tensor(([[0]*lattice.D])))
assert rho_tensor.shape == (1,lattice.Q)
assert rho_tensor.cpu().numpy() == pytest.approx(np.ones((1,lattice.Q)))
@pytest.mark.parametrize("MomentSet", (D2Q9Dellar, D2Q9Lallemand))
def test_conserved_moments_d2q9(MomentSet):
multiindices = np.array([
[0, 0], [1,0], [0,1]
])
m = moment_tensor(D2Q9.e, multiindices)
assert m == pytest.approx(MomentSet.matrix[:3,:])
def test_inverse_transform(f_transform):
f, transform = f_transform
lattice = transform.lattice
retransformed = lattice.convert_to_numpy(transform.inverse_transform(transform.transform(f)))
original = lattice.convert_to_numpy(f)
assert retransformed == pytest.approx(original, abs=1e-5)
def test_getitem(dtype_device):
dtype, device = dtype_device
moments = D2Q9Lallemand(Lattice(D2Q9, device, dtype))
assert moments["jx", "jy"] == [1,2]
assert moments["rho"] == [0]
def test_moment_equilibrium_dellar(dtype_device):
dtype, device = dtype_device
lattice = Lattice(D2Q9, device, dtype)
moments = D2Q9Dellar(lattice)
np.random.seed(1)
f = lattice.convert_to_tensor(np.random.random([lattice.Q] + [3] * lattice.D))
meq1 = lattice.convert_to_numpy(moments.transform(lattice.equilibrium(lattice.rho(f), lattice.u(f))))
meq2 = lattice.convert_to_numpy(moments.equilibrium(moments.transform(f)))
assert meq1 == pytest.approx(meq2, abs=1e-5)
def test_moment_equilibrium_lallemand(dtype_device):
dtype, device = dtype_device
lattice = Lattice(D2Q9, device, dtype)
moments = D2Q9Lallemand(lattice)
np.random.seed(1)
f = lattice.convert_to_tensor(np.random.random([lattice.Q] + [3] * lattice.D))
meq1 = lattice.convert_to_numpy(moments.transform(lattice.equilibrium(lattice.rho(f), lattice.u(f))))
meq2 = lattice.convert_to_numpy(moments.equilibrium(moments.transform(f)))
same_moments = moments["rho", "jx", "jy", "qx", "qy"]
assert meq1[same_moments] == pytest.approx(meq2[same_moments], abs=1e-5)
def test_moment_equilibrium_D3Q27Hermite(dtype_device):
dtype, device = dtype_device
lattice = Lattice(D3Q27, device, dtype)
moments = D3Q27Hermite(lattice)
np.random.seed(1)
f = lattice.convert_to_tensor(np.random.random([lattice.Q] + [3] * lattice.D))
meq1 = lattice.convert_to_numpy(moments.transform(lattice.equilibrium(lattice.rho(f), lattice.u(f))))
meq2 = lattice.convert_to_numpy(moments.equilibrium(moments.transform(f)))
same_moments = moments['rho', 'jx', 'jy', 'jz', 'Pi_xx', 'Pi_xy', 'PI_xz', 'PI_yy', 'PI_yz', 'PI_zz']
assert meq1[same_moments] == pytest.approx(meq2[same_moments], abs=1e-5)
@pytest.mark.parametrize("MomentSet", (D2Q9Dellar, D2Q9Lallemand, D3Q27Hermite))
def test_orthogonality(dtype_device,MomentSet):
dtype, device = dtype_device
lattice = Lattice(MomentSet.supported_stencils[0], device, dtype)
moments = MomentSet(lattice)
M = Lattice.convert_to_numpy(moments.matrix)
if MomentSet == D2Q9Lallemand:
Md = np.round(M @ M.T, 4)
else:
Md = np.round(M @ np.diag(lattice.stencil.w) @ M.T, 4)
assert np.where(np.diag(np.ones(lattice.stencil.Q())),Md!=0.0,Md==0.0).all()