lattice.py 10.8 KB
Newer Older
MischaD's avatar
MischaD committed
1
import numpy as np
MischaD's avatar
MischaD committed
2
import sympy as sp
MischaD's avatar
MischaD committed
3
import random
MischaD's avatar
MischaD committed
4
from sympy import oo
Mischa Dombrowski's avatar
Mischa Dombrowski committed
5
from .utils.mylog import logger
MischaD's avatar
MischaD committed
6
from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException, NoSolutionException, InfiniteSolutionsException, UninitializedAttributeException
MischaD's avatar
MischaD committed
7
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, approximate_ratio, timer, get_random_vector
MischaD's avatar
MischaD committed
8
9


MischaD's avatar
MischaD committed
10
class Lattice:
MischaD's avatar
MischaD committed
11

MischaD's avatar
MischaD committed
12
    def __init__(self, dimension=2, order=4, shell_list=[1, 2, 4], seed=None, svd_tolerance=1e-8):
MischaD's avatar
MischaD committed
13
14
15
16
17
        """
        :param dimension:
        :param order:
        :param shell_list:
        :param seed:
MischaD's avatar
MischaD committed
18
        :param tolerance: tolerance to set the singular values to zero after svd
MischaD's avatar
MischaD committed
19
        """
MischaD's avatar
MischaD committed
20
        self._dimension = dimension
MischaD's avatar
MischaD committed
21
22
23
        self._order = order
        self._shell_list = shell_list
        self._seed = seed
MischaD's avatar
MischaD committed
24
        self._svd_tolerance = svd_tolerance
MischaD's avatar
MischaD committed
25
        random.seed(self._seed)
MischaD's avatar
MischaD committed
26

MischaD's avatar
MischaD committed
27
        # Preperation
MischaD's avatar
MischaD committed
28
        self._tensor_space_dimension = 0
MischaD's avatar
MischaD committed
29
        self._all_velocity_vectors = []
MischaD's avatar
MischaD committed
30
        self._possible_tensors = []
MischaD's avatar
MischaD committed
31
        self._tensor_dimensions = []
MischaD's avatar
MischaD committed
32
33
        self._velocities = 0
        self._shells = 0
MischaD's avatar
MischaD committed
34

MischaD's avatar
MischaD committed
35
36
37
38
        # LES
        self._lhs = []
        self._rhs = []

MischaD's avatar
MischaD committed
39
40
        # SVD
        self._singular_values = []
MischaD's avatar
MischaD committed
41
        self._solution = []
MischaD's avatar
MischaD committed
42

MischaD's avatar
MischaD committed
43
44
        # Output
        self._coefficients = []
MischaD's avatar
MischaD committed
45
46
47
        self._weight_polynomials = []
        self._interval = None

MischaD's avatar
MischaD committed
48
49
    def __str__(self):
        string = "D{} - Order: {} - Shells: {}".format(
MischaD's avatar
MischaD committed
50
            self._dimension, self._order, str(self._shell_list)
MischaD's avatar
MischaD committed
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
        )
        if self._seed:
            string += " - Seed: {}".format(self._seed)
        return string

    @classmethod
    def init_by_order(cls, dimension, order, seed=None):
        """
        :param dimension:
        :param order:
        :return:
        """
        shell_list = []
        if dimension == 2:
            if order == 4:
                shell_list = [1, 2, 4]
            elif order == 6:
                shell_list = [1, 2, 4, 8, 9]
            elif order == 8:
                shell_list = [1, 2, 4, 5, 8, 9, 10, 16]

        elif dimension == 3:
            if order == 4:
                shell_list = [1, 2, 4]

        if not shell_list:
            raise OrderNotImplementedException(
                "Cannot initialize by {} for dimension {}. Lower the order or try to come up with the shells yourself.".format(
                    order, dimension
                )
            )
        return cls(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
MischaD's avatar
MischaD committed
83

MischaD's avatar
MischaD committed
84
85
86
87
    @property
    def tensor_space_dimension(self):
        return np.array([len(x) for x in self._possible_tensors]).sum()

MischaD's avatar
MischaD committed
88
89
90
91
92
93
94
95
    @property
    def shells(self):
        self._shells

    @property
    def velocities(self):
        self._velocities

MischaD's avatar
MischaD committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
    @staticmethod
    def velocities_for_shell(dimension, shell):
        """
        :param dimension: Lattice dimension
        :param shell: Squared velocity of the shell
        :return: list of velocities with squared length equal to shell
        """
        velocities = []
        linear_lattice_size = 2 * shell + 1
        full_lattice_size = linear_lattice_size ** dimension
        for site in range(full_lattice_size):
            tmp = site
            cur_vel_squared = 0
            cur_vec = []

            for dim in range(dimension):
                coordinate = tmp % linear_lattice_size
                tmp = (tmp - coordinate) // linear_lattice_size
                shifted_coordinate = coordinate - shell
                cur_vec.append(shifted_coordinate)
                cur_vel_squared += shifted_coordinate ** 2

            if cur_vel_squared == shell:
                velocities.append(np.array(cur_vec, dtype=int))

        return velocities

MischaD's avatar
MischaD committed
123
    def calculate_velocity_vectors(self):
MischaD's avatar
MischaD committed
124
        velocities_amount = 1
MischaD's avatar
MischaD committed
125
        velocity_vectors = []
MischaD's avatar
MischaD committed
126
        total_subshells = []
MischaD's avatar
MischaD committed
127

MischaD's avatar
MischaD committed
128
        group = get_group(self._dimension)
MischaD's avatar
MischaD committed
129
        for shell in range(len(self._shell_list)):
MischaD's avatar
MischaD committed
130
            velocities = self.velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
MischaD's avatar
MischaD committed
131
            if len(velocities) == 0:
MischaD's avatar
MischaD committed
132
133
134
135
136
                raise ImpossibleVelocityShellException("No velocity sets found for velocity shell {}.".format(shell))

            subshells = get_subshells(velocities, group)
            total_subshells.append(len(subshells))
            if len(subshells) > 1:
MischaD's avatar
MischaD committed
137
138
139
140
141
142
                # 3 4, 0 5 can only be given together
                for i_subs, subshell in enumerate(subshells):
                    type = tuple(np.sort(abs(subshell[0])))
                velocities = subshells
            else:
                velocities = [velocities]
MischaD's avatar
MischaD committed
143
            velocity_vectors.extend(velocities)
MischaD's avatar
MischaD committed
144
145
146
147
            velocities_amount += len(velocities)

        self._velocities = velocities_amount
        self._shells = len(total_subshells) + 1
MischaD's avatar
MischaD committed
148

MischaD's avatar
MischaD committed
149
150
        logger.info("Velocities: {} Shells: {}".format(self._velocities, self._shells))

MischaD's avatar
MischaD committed
151
        for i, shell in enumerate(velocity_vectors):
MischaD's avatar
MischaD committed
152
153
            number_of_velocities = len(shell)
            type = tuple(np.sort(abs(shell[0])))
MischaD's avatar
MischaD committed
154
155
156
            # logger.info("Shell number {} with c_i^2 = {} and {} velocities of type {}".format(i+1, (shell[0]**2).sum(), number_of_velocities, type))
        self._all_velocity_vectors = velocity_vectors
        return velocity_vectors
MischaD's avatar
MischaD committed
157

MischaD's avatar
MischaD committed
158
    def fill_lhs(self):
MischaD's avatar
MischaD committed
159
160
161
162
        lhs = []
        for i, rank in enumerate(np.arange(2, self._order + 2, 2)):
            tensor_space_dimension = self._tensor_dimensions[i]
            for j in range(tensor_space_dimension):
MischaD's avatar
MischaD committed
163
                vector = get_random_vector(self._dimension)
MischaD's avatar
MischaD committed
164
165
                # vector = 2 * np.random.rand(self._dimension) - 1
                # vector = vector / np.linalg.norm(vector)
MischaD's avatar
MischaD committed
166
167

                rows = []
MischaD's avatar
MischaD committed
168
                for shell_velocities in self._all_velocity_vectors:
MischaD's avatar
MischaD committed
169
170
171
                    shell_sum = lattice_sum(vector, shell_velocities, rank)
                    rows.append(shell_sum)
                lhs.append(rows)
MischaD's avatar
MischaD committed
172
        self._lhs = np.array(lhs)
MischaD's avatar
MischaD committed
173
        return np.array(lhs)
MischaD's avatar
MischaD committed
174

MischaD's avatar
MischaD committed
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
    def fill_rhs(self):
        rhs = []
        cols = self._order // 2
        for k, rank in enumerate(np.arange(2, self._order + 2, 2)):
            for j in range(self._tensor_dimensions[k]):
                rows = []
                for i in range(cols):
                    c_s_power = 2 * i + 2
                    rows.append(1 if c_s_power == rank else 0)
                rhs.append(rows)
        self._rhs = np.array(rhs)
        return np.array(rhs)

    def svd(self):
        U, s, V = np.linalg.svd(self._lhs)
        rows = self._lhs.shape[0]
        cols = self._lhs.shape[1]
MischaD's avatar
MischaD committed
192
193
194
195

        self._singular_values = len(s)
        s = np.array([sv if sv > self._svd_tolerance else 0 for sv in s])
        rank = np.count_nonzero(s)
MischaD's avatar
MischaD committed
196
        rhs = np.dot(np.transpose(U), self._rhs)
MischaD's avatar
MischaD committed
197
198
199
200
        if rank < rows:
            overlap = rows - rank
            if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
                # NO Solution
MischaD's avatar
MischaD committed
201
                raise NoSolutionException("Add shells or lower the order")  # change shells ?
MischaD's avatar
MischaD committed
202
203
            rows = rank
            s.resize(rows)
MischaD's avatar
MischaD committed
204
            rhs.resize((rows, rhs.shape[1]), refcheck=False)
MischaD's avatar
MischaD committed
205
206
207
208
209
210
211
212
213
214
215
216
217
            self._singular_values = len(s)

        reduced_rhs = np.zeros((rows, rhs.shape[1]))
        for i, SingularValue in enumerate(s):
            for j in range(rhs.shape[1]):
                reduced_rhs[i, j] = rhs[i, j] / SingularValue

        reduced_rhs = np.zeros((rows, rhs.shape[1]))
        for i, singular_value in enumerate(s):
            reduced_rhs[i, :] = rhs[i, :] / singular_value

        if cols - rows > 0:
            # TODO Better
MischaD's avatar
MischaD committed
218
            raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
MischaD's avatar
MischaD committed
219
220
221
                                             "https://github.com/BDuenweg/Lattice-Boltzmann-weights")
        self._solution = np.dot(np.transpose(V), reduced_rhs)
        return self._solution
MischaD's avatar
MischaD committed
222

MischaD's avatar
MischaD committed
223
    def calculate_coefficients(self):
MischaD's avatar
MischaD committed
224
225
226
227
228
        solution_columns = self._order // 2
        coefficients = np.zeros((1 + solution_columns))
        coefficients[0] = 1
        for j in range(solution_columns):
            tmp = 0
MischaD's avatar
MischaD committed
229
230
231
            for i in range(len(self._all_velocity_vectors)):
                tmp -= self._solution[i, j] * len(self._all_velocity_vectors[i])
            coefficients[j + 1] = tmp
MischaD's avatar
MischaD committed
232

MischaD's avatar
MischaD committed
233
        coeffs = np.zeros((self._solution.shape[0] + 1, self._solution.shape[1] + 1))
MischaD's avatar
Bugfix    
MischaD committed
234
        coeffs[0, :] = coefficients[-1::-1]
MischaD's avatar
MischaD committed
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
        coeffs[1:, :-1] = self._solution[:, -1::-1]
        self._coefficients = coeffs
        return self._coefficients

    @timer
    def calculate_valid_interval(self):
        if not self._weight_polynomials:
            raise UninitializedAttributeException("Weight Polynomials not give. Calculate them first!")

        interval = sp.Interval(-oo, oo)
        for equation in self._weight_polynomials:
            roots = equation.real_roots()
            if not roots:
                continue
            roots = [roots[0] - 1] + roots + [roots[-1] + 1]
            cur_interval = sp.EmptySet()
            for i in range(len(roots) - 1):
                mesh_point = (roots[i + 1] + roots[i])/2
                if equation.eval(mesh_point) < 0:
                    continue
                # positive weights
                if i == 0:
                    cur_interval = sp.Union(cur_interval, sp.Interval(-oo, roots[i+1]))
                elif i == (len(roots) - 2):
                    cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], oo))
                else:
                    cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], roots[i+1]))
            interval = interval.intersect(cur_interval)
        self._interval = sp.Complement(interval, sp.FiniteSet(0))
        return self._interval

    @timer
    def calculate_weights(self):
        c_s_sq = sp.Symbol("c_s_sq")
        if self._order % 2 == 1:
            self._order -= 1
            logger.warning("Only odd order valid. Decrease by one")

        for k in np.arange(2, self._order + 2, 2):
            possible_tensors = analyse_tensor_dimension(k)
            self._possible_tensors.append(possible_tensors)
        self._tensor_dimensions = [len(x) for x in self._possible_tensors]

        self.calculate_velocity_vectors()
        self.fill_lhs()
        self.fill_rhs()

        self.svd()

        self.calculate_coefficients()
MischaD's avatar
Bugfix    
MischaD committed
285

MischaD's avatar
MischaD committed
286
        apporx_coeffs = [[approximate_ratio(x) for x in pol_coffs] for pol_coffs in self._coefficients]
MischaD's avatar
MischaD committed
287

MischaD's avatar
MischaD committed
288
        self._weight_polynomials = [sp.Poly([approximate_ratio(x) for x in pol_coffs], c_s_sq) for pol_coffs in self._coefficients]
MischaD's avatar
MischaD committed
289

MischaD's avatar
MischaD committed
290
        self._interval = self.calculate_valid_interval()
MischaD's avatar
MischaD committed
291
292
        logger.info(self._interval.boundary)
        if self._interval == sp.EmptySet():
MischaD's avatar
Bugfix    
MischaD committed
293
            return []
MischaD's avatar
MischaD committed
294
        weights = [sp.N(x.eval(self._interval.boundary.inf)) for x in self._weight_polynomials]
MischaD's avatar
MischaD committed
295
        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
MischaD's avatar
MischaD committed
296
        return weights
MischaD's avatar
MischaD committed
297
298
299
300
301
302
303
304
305
306
307
308