lattice.py 15.1 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
7
8
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, approximate_ratio, timer, get_random_vector, velocities_for_shell
from .shell import Shell
MischaD's avatar
MischaD committed
9
10


MischaD's avatar
MischaD committed
11
class Lattice:
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
    BY_NAME = {"D2Q9": {"dimension": 2, "order": 4, "shell_list": [1, 2, 4]},
               "D2V17": {"dimension": 2, "order": 6, "shell_list": [1, 2, 4, 8, 9]},
               "D2V37": {"dimension": 2, "order": 8, "shell_list": [1, 2, 4, 5, 8, 9, 10, 16]},
               "D3Q19": {"dimension": 3, "order": 4, "shell_list": [1, 2, 4]},
               "D3Q15": {"dimension": 3, "order": 4, "shell_list": [1, 3, 4]},
               "D3Q41-ZOT": {"dimension": 3, "order": 6, "shell_list": [1, 2, 3, 9, 16, 27], "boundary": "sup",
                             "unwanted_subshells": ["221", "511"]},
               }

    def __init__(self, dimension=2, order=4, shell_list=[1, 2, 4], seed=None, boundary=None, unwanted_subshells=[], svd_tolerance=1e-8):
        """Create an instance of the Lattice class without calculating the weights.
        Default values are chosen to lead to the standard D2Q9 lattice.
        Not all input values are valid. :func:`Lattice.init_by_order`

        Args:
            name: Initialize Lattice with the correct values for some known velocities. Possible names are
                D2Q9, D2V17, D2V37, D3Q17,
            dimension (int): Dimension of the resulting lattice. Default 2.
            order (int): Order of tensor isotropy. Only odd orders are
            shell_list (list): List of the Squared lengths of the velocities.
            unwanted_subshells (list) : List of strings that correspond to the type of a subshell that isn't supposed to
                to be part of lattice. Some shells are ambiguous and in the default case all possiblities are taken,
                but some lattice, like the D3Q41-ZOT need certain subshells to be remove. E.g. 2 Dimensions and
                shell 25 is added. Bother type = "50" and type = "43" could be used. Add one of them to the list of
                unwanted subshells to remove them.
            seed (float): Seed for the random number generator.
            svd_tolerance (float): tolerance to set the singular values to zero after svd.
MischaD's avatar
MischaD committed
39

MischaD's avatar
MischaD committed
40
        """
MischaD's avatar
MischaD committed
41
        self._dimension = dimension
MischaD's avatar
MischaD committed
42
43
44
        self._order = order
        self._shell_list = shell_list
        self._seed = seed
MischaD's avatar
MischaD committed
45
        self._svd_tolerance = svd_tolerance
MischaD's avatar
MischaD committed
46
        random.seed(self._seed)
MischaD's avatar
MischaD committed
47

MischaD's avatar
MischaD committed
48
        # Preperation
MischaD's avatar
MischaD committed
49
        self._tensor_space_dimension = 0
MischaD's avatar
MischaD committed
50
        self._all_velocity_vectors = []
MischaD's avatar
MischaD committed
51
        self._possible_tensors = []
MischaD's avatar
MischaD committed
52
        self._tensor_dimensions = []
MischaD's avatar
MischaD committed
53
        self._shells = 0
54
        self._unwanted_subshells = unwanted_subshells
MischaD's avatar
MischaD committed
55

MischaD's avatar
MischaD committed
56
57
58
59
        # LES
        self._lhs = []
        self._rhs = []

MischaD's avatar
MischaD committed
60
61
        # SVD
        self._singular_values = []
MischaD's avatar
MischaD committed
62
        self._solution = []
MischaD's avatar
MischaD committed
63

64
        # Polynomials
MischaD's avatar
MischaD committed
65
        self._coefficients = []
MischaD's avatar
MischaD committed
66
67
        self._weight_polynomials = []
        self._interval = None
68
        self._boundary = boundary if boundary == "sup" else "inf"
MischaD's avatar
MischaD committed
69

70
71
72
73
74
        # Final values
        self._discrete_velocities = []
        self._weights = []
        self._c_s_sq = None

MischaD's avatar
MischaD committed
75
76
    def __str__(self):
        string = "D{} - Order: {} - Shells: {}".format(
MischaD's avatar
MischaD committed
77
            self._dimension, self._order, str(self._shell_list)
MischaD's avatar
MischaD committed
78
79
80
81
82
83
        )
        if self._seed:
            string += " - Seed: {}".format(self._seed)
        return string

    @classmethod
84
85
86
87
88
89
90
91
92
93
    def from_name(cls, name, seed=None):
        """Initiate the Lattice with the correct values corresponding to a know Lattice in literature.
        A complete list can be seen in Lattice.BY_NAME.keys()"""
        kwargs = cls.BY_NAME.get(name)
        if not kwargs:
            raise ValueError("No Lattice with this name is known or implemented")
        return cls(**kwargs, seed=seed)

    @classmethod
    def from_order(cls, dimension, order, seed=None):
MischaD's avatar
MischaD committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
        """
        :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
119

MischaD's avatar
MischaD committed
120
121
122
123
    @property
    def tensor_space_dimension(self):
        return np.array([len(x) for x in self._possible_tensors]).sum()

MischaD's avatar
MischaD committed
124
125
    @property
    def shells(self):
126
        return self._shells
MischaD's avatar
MischaD committed
127
128
129

    @property
    def velocities(self):
130
131
        if not self._discrete_velocities and self._shells:
            self._discrete_velocities = [shell.type for shell in self._shells]
132
        return self._discrete_velocities
MischaD's avatar
MischaD committed
133

134
135
136
137
138
139
    @property
    def weights(self):
        if self._weights:
            return self._weights
        return self.calculate_weights()

MischaD's avatar
MischaD committed
140
141
    def calculate_velocity_vectors(self):
        velocity_vectors = []
MischaD's avatar
MischaD committed
142

MischaD's avatar
MischaD committed
143
        group = get_group(self._dimension)
MischaD's avatar
MischaD committed
144
        for shell in range(len(self._shell_list)):
145
            velocities = velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
MischaD's avatar
MischaD committed
146
            if len(velocities) == 0:
MischaD's avatar
MischaD committed
147
148
                raise ImpossibleVelocityShellException("No velocity sets found for velocity shell {}.".format(shell))

149
            # If e.g. dimension = 2, then shell 25 in ambiguous. Both shell (3,4) and (5,0) are added
MischaD's avatar
MischaD committed
150
151
            subshells = get_subshells(velocities, group)
            if len(subshells) > 1:
MischaD's avatar
MischaD committed
152
                # 3 4, 0 5 can only be given together
153
                velocities = []
MischaD's avatar
MischaD committed
154
                for i_subs, subshell in enumerate(subshells):
155
156
157
158
                    # get type for subshell. Example shell: all vectors of type [0, -1, 0] --> 100
                    type = ''.join([str(int(x)) for x in sorted(abs(subshell[0]), reverse=True)])
                    if type not in self._unwanted_subshells:
                        velocities.append(subshell)
MischaD's avatar
MischaD committed
159
160
            else:
                velocities = [velocities]
MischaD's avatar
MischaD committed
161
            velocity_vectors.extend(velocities)
MischaD's avatar
MischaD committed
162

163
164
165
        self._shells = [Shell([np.array([0]*self._dimension)])]
        for final_velocities in velocity_vectors:
            self._shells.append(Shell(final_velocities))
MischaD's avatar
MischaD committed
166
167
        self._all_velocity_vectors = velocity_vectors
        return velocity_vectors
MischaD's avatar
MischaD committed
168

MischaD's avatar
MischaD committed
169
    def fill_lhs(self):
MischaD's avatar
MischaD committed
170
171
172
173
        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
174
                vector = get_random_vector(self._dimension)
MischaD's avatar
MischaD committed
175
176
                # vector = 2 * np.random.rand(self._dimension) - 1
                # vector = vector / np.linalg.norm(vector)
MischaD's avatar
MischaD committed
177
178

                rows = []
MischaD's avatar
MischaD committed
179
                for shell_velocities in self._all_velocity_vectors:
MischaD's avatar
MischaD committed
180
181
182
                    shell_sum = lattice_sum(vector, shell_velocities, rank)
                    rows.append(shell_sum)
                lhs.append(rows)
MischaD's avatar
MischaD committed
183
        self._lhs = np.array(lhs)
MischaD's avatar
MischaD committed
184
        return np.array(lhs)
MischaD's avatar
MischaD committed
185

MischaD's avatar
MischaD committed
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    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
203
204
205
206

        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
207
        rhs = np.dot(np.transpose(U), self._rhs)
MischaD's avatar
MischaD committed
208
209
210
211
        if rank < rows:
            overlap = rows - rank
            if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
                # NO Solution
MischaD's avatar
MischaD committed
212
                raise NoSolutionException("Add shells or lower the order")  # change shells ?
MischaD's avatar
MischaD committed
213
214
            rows = rank
            s.resize(rows)
MischaD's avatar
MischaD committed
215
            rhs.resize((rows, rhs.shape[1]), refcheck=False)
MischaD's avatar
MischaD committed
216
217
218
219
220
221
222
223
224
225
226
227
228
            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
229
            raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
MischaD's avatar
MischaD committed
230
231
232
                                             "https://github.com/BDuenweg/Lattice-Boltzmann-weights")
        self._solution = np.dot(np.transpose(V), reduced_rhs)
        return self._solution
MischaD's avatar
MischaD committed
233

MischaD's avatar
MischaD committed
234
    def calculate_coefficients(self):
MischaD's avatar
MischaD committed
235
236
237
238
239
        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
240
241
242
            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
243

MischaD's avatar
MischaD committed
244
        coeffs = np.zeros((self._solution.shape[0] + 1, self._solution.shape[1] + 1))
MischaD's avatar
Bugfix    
MischaD committed
245
        coeffs[0, :] = coefficients[-1::-1]
MischaD's avatar
MischaD committed
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
        coeffs[1:, :-1] = self._solution[:, -1::-1]
        self._coefficients = coeffs
        return self._coefficients

    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

276
    def calculate_polynomials(self):
MischaD's avatar
MischaD committed
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
        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
294

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

MischaD's avatar
MischaD committed
297
        self._weight_polynomials = [sp.Poly([approximate_ratio(x) for x in pol_coffs], c_s_sq) for pol_coffs in self._coefficients]
298
299
        if not self._weight_polynomials:
            raise NoSolutionException("Something went wrong")
MischaD's avatar
MischaD committed
300

MischaD's avatar
MischaD committed
301
        self._interval = self.calculate_valid_interval()
MischaD's avatar
MischaD committed
302
        if self._interval == sp.EmptySet():
MischaD's avatar
Bugfix    
MischaD committed
303
            return []
304
305
306
307
308
309
310
311
312
313
314
315
316
317
        return self._weight_polynomials

    def calculate_reduced_weights(self, boundary="inf"):
        """
        Calculate the weights of this lattice for a reduced model.
        If the smallest or biggest value of the valid interval is taken, they reduce
        the amount of shells necessary by one.
        Overwrites the weight values in self.shell

        :param boundary: "inf" --> smallest possible value is taken, otherwise the largest
        :return: list of weights in the order of their shells.
        """
        return self.calculate_weights(boundary=boundary)

318
    def calculate_weights(self, c_s_sq=None, boundary=None):
319
        """
320
321
322
323
324
325
        Calculate the weights for this lattice for a given speed of sound value c_s_sq, or for
        a given boundary.
        It is first checked, whether a c_s_sq value is given and then boundary.
        If both are not giving, self._boundary is taken which is defined in the
        init function.

326
327
328
329
330
331
332
        If no value is given, the infimum is used, i.e. the smallest possible value.

        :param c_s_sq: Float value for the speed of sound, has to be inside of the valid interval
        :param boundary: Get reduced model by using the infimum or supremum of the interval. Fallback value if
            no c_s_sq is given
        :return:
        """
333

334
335
336
337
338
        if not self._weight_polynomials:
            self.calculate_polynomials()
        if not self._interval:  # Empty Interval
            return []
        if not c_s_sq:
339
340
            if not boundary:
                boundary = self._boundary
341
342
343
344
            if boundary == "inf":
                c_s_sq = self._interval.inf
            else:
                c_s_sq = self._interval.sup
345

346
347
348
349
        if self._interval.contains(c_s_sq) is sp.EmptySet:  # Invalid c_s_sq input
            return []

        weights = [sp.N(x.eval(c_s_sq)) for x in self._weight_polynomials]
MischaD's avatar
MischaD committed
350
        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
351
352
353
354
        for i, shell in enumerate(self.shells):
            shell.set_weight(c_s_sq, weights[i])
        self._c_s_sq = c_s_sq
        self._weights = weights
MischaD's avatar
MischaD committed
355
        return weights
MischaD's avatar
MischaD committed
356

357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
    def velocity_set(self):
        """
        Returns the values that are important for the Lattice Boltzmann Model, i.e.
        a list of the velocity values their corresponding weights and the value for the speed of sound.
        :return: (c_s_sq, weights, velocities): c_s_sq
        """

        weights = []
        for shell in self.shells:
            if shell.weight == 0:
                continue
            weights += [shell.weight] * shell.size
        velocities = np.zeros((self._dimension, len(weights)))
        i = 0
        for shell in self.shells:
            if shell.weight == 0:
                continue
            for velocity in shell.velocities:
                velocities[:,i] = velocity
                i += 1

        return self._c_s_sq, weights, velocities
MischaD's avatar
MischaD committed
379
380
381
382
383
384
385
386
387
388
389