lattice.py 16.9 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
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
from .exceptions import (
    OrderNotImplementedException,
    ImpossibleVelocityShellException,
    ReducibleShellException,
    NoSolutionException,
    InfiniteSolutionsException,
    UninitializedAttributeException,
)
from .utils.utils import (
    analyse_tensor_dimension,
    get_group,
    get_subshells,
    lattice_sum,
    approximate_ratio,
    timer,
    get_random_vector,
    velocities_for_shell,
)
24
from .shell import Shell
MischaD's avatar
MischaD committed
25
26


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

    def __init__(
        self,
        dimension=2,
        order=4,
        shell_list=[1, 2, 4],
        seed=None,
        boundary=None,
        unwanted_subshells=[],
        svd_tolerance=1e-8,
    ):
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
        """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
81

MischaD's avatar
MischaD committed
82
        """
MischaD's avatar
MischaD committed
83
        self._dimension = dimension
MischaD's avatar
MischaD committed
84
85
86
        self._order = order
        self._shell_list = shell_list
        self._seed = seed
MischaD's avatar
MischaD committed
87
        self._svd_tolerance = svd_tolerance
MischaD's avatar
MischaD committed
88
        random.seed(self._seed)
MischaD's avatar
MischaD committed
89

MischaD's avatar
MischaD committed
90
        # Preperation
MischaD's avatar
MischaD committed
91
        self._tensor_space_dimension = 0
MischaD's avatar
MischaD committed
92
        self._all_velocity_vectors = []
MischaD's avatar
MischaD committed
93
        self._possible_tensors = []
MischaD's avatar
MischaD committed
94
        self._tensor_dimensions = []
MischaD's avatar
MischaD committed
95
        self._shells = 0
96
        self._unwanted_subshells = unwanted_subshells
MischaD's avatar
MischaD committed
97

MischaD's avatar
MischaD committed
98
99
100
101
        # LES
        self._lhs = []
        self._rhs = []

MischaD's avatar
MischaD committed
102
103
        # SVD
        self._singular_values = []
MischaD's avatar
MischaD committed
104
        self._solution = []
MischaD's avatar
MischaD committed
105

106
        # Polynomials
MischaD's avatar
MischaD committed
107
        self._coefficients = []
MischaD's avatar
MischaD committed
108
109
        self._weight_polynomials = []
        self._interval = None
110
        self._boundary = boundary if boundary == "sup" else "inf"
MischaD's avatar
MischaD committed
111

112
113
114
115
116
        # Final values
        self._discrete_velocities = []
        self._weights = []
        self._c_s_sq = None

MischaD's avatar
MischaD committed
117
118
    def __str__(self):
        string = "D{} - Order: {} - Shells: {}".format(
MischaD's avatar
MischaD committed
119
            self._dimension, self._order, str(self._shell_list)
MischaD's avatar
MischaD committed
120
        )
MischaD's avatar
MischaD committed
121
122
123
        if self._boundary:
            string += " - Boundary: {}".format(self._boundary)
        if self._unwanted_subshells:
MischaD's avatar
MischaD committed
124
125
126
            string += " - Unwanted Subshells: {}".format(
                "".join([s + ", " for s in self._unwanted_subshells])
            )
MischaD's avatar
MischaD committed
127
128
129
        return string

    @classmethod
MischaD's avatar
MischaD committed
130
    def from_name(cls, name):
131
132
133
134
135
        """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")
MischaD's avatar
MischaD committed
136
        return cls(**kwargs)
137
138

    @classmethod
MischaD's avatar
MischaD committed
139
    def from_order(cls, dimension, order):
MischaD's avatar
MischaD committed
140
141
142
143
144
        """
        :param dimension:
        :param order:
        :return:
        """
MischaD's avatar
MischaD committed
145
146
147
148
149
150
151
152
        if order % 2 != 0:
            order -= 1

        for key, value in cls.BY_NAME.items():
            if value.get("dimension") == dimension and value.get("order") == order:
                return cls(**cls.BY_NAME.get(key))

        raise ValueError(
MischaD's avatar
MischaD committed
153
            "No Lattice known with order {} and dimension {}.".format(order, dimension)
MischaD's avatar
MischaD committed
154
        )
MischaD's avatar
MischaD committed
155
        return cls(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
MischaD's avatar
MischaD committed
156

MischaD's avatar
MischaD committed
157
158
159
160
    @property
    def tensor_space_dimension(self):
        return np.array([len(x) for x in self._possible_tensors]).sum()

MischaD's avatar
MischaD committed
161
162
    @property
    def shells(self):
163
        return self._shells
MischaD's avatar
MischaD committed
164
165
166

    @property
    def velocities(self):
MischaD's avatar
MischaD committed
167
168
        """Get a list of all types of velocities. Velocities with zero weight are included
        For D2O9 this equals [00, 10, 11, 20]"""
169
170
        if not self._discrete_velocities and self._shells:
            self._discrete_velocities = [shell.type for shell in self._shells]
171
        return self._discrete_velocities
MischaD's avatar
MischaD committed
172

173
174
175
176
177
178
    @property
    def weights(self):
        if self._weights:
            return self._weights
        return self.calculate_weights()

MischaD's avatar
MischaD committed
179
180
181
182
183
184
    @property
    def c_s_sq(self):
        if self._c_s_sq:
            return self._c_s_sq
        return self.velocity_set()[0]

MischaD's avatar
MischaD committed
185
186
187
188
189
    def shell_from_type(self, type):
        for shell in self._shells:
            if shell.type == type:
                return shell

MischaD's avatar
MischaD committed
190
    def _calculate_velocity_vectors(self):
MischaD's avatar
MischaD committed
191
        velocity_vectors = []
MischaD's avatar
MischaD committed
192

MischaD's avatar
MischaD committed
193
        group = get_group(self._dimension)
MischaD's avatar
MischaD committed
194
        for shell in range(len(self._shell_list)):
MischaD's avatar
MischaD committed
195
196
197
            velocities = velocities_for_shell(
                dimension=self._dimension, shell=self._shell_list[shell]
            )
MischaD's avatar
MischaD committed
198
            if len(velocities) == 0:
MischaD's avatar
MischaD committed
199
200
201
                raise ImpossibleVelocityShellException(
                    "No velocity sets found for velocity shell {}.".format(shell)
                )
MischaD's avatar
MischaD committed
202

203
            # If e.g. dimension = 2, then shell 25 in ambiguous. Both shell (3,4) and (5,0) are added
MischaD's avatar
MischaD committed
204
205
            subshells = get_subshells(velocities, group)
            if len(subshells) > 1:
MischaD's avatar
MischaD committed
206
                # 3 4, 0 5 can only be given together
207
                velocities = []
MischaD's avatar
MischaD committed
208
                for i_subs, subshell in enumerate(subshells):
209
                    # get type for subshell. Example shell: all vectors of type [0, -1, 0] --> 100
MischaD's avatar
MischaD committed
210
211
212
                    type = "".join(
                        [str(int(x)) for x in sorted(abs(subshell[0]), reverse=True)]
                    )
213
214
                    if type not in self._unwanted_subshells:
                        velocities.append(subshell)
MischaD's avatar
MischaD committed
215
216
            else:
                velocities = [velocities]
MischaD's avatar
MischaD committed
217
            velocity_vectors.extend(velocities)
MischaD's avatar
MischaD committed
218

MischaD's avatar
MischaD committed
219
        self._shells = [Shell([np.array([0] * self._dimension)])]
220
221
        for final_velocities in velocity_vectors:
            self._shells.append(Shell(final_velocities))
MischaD's avatar
MischaD committed
222
223
        self._all_velocity_vectors = velocity_vectors
        return velocity_vectors
MischaD's avatar
MischaD committed
224

MischaD's avatar
MischaD committed
225
    def _fill_lhs(self):
MischaD's avatar
MischaD committed
226
227
228
229
        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
230
                vector = get_random_vector(self._dimension)
MischaD's avatar
MischaD committed
231
232
                # vector = 2 * np.random.rand(self._dimension) - 1
                # vector = vector / np.linalg.norm(vector)
MischaD's avatar
MischaD committed
233
234

                rows = []
MischaD's avatar
MischaD committed
235
                for shell_velocities in self._all_velocity_vectors:
MischaD's avatar
MischaD committed
236
237
238
                    shell_sum = lattice_sum(vector, shell_velocities, rank)
                    rows.append(shell_sum)
                lhs.append(rows)
MischaD's avatar
MischaD committed
239
        self._lhs = np.array(lhs)
MischaD's avatar
MischaD committed
240
        return np.array(lhs)
MischaD's avatar
MischaD committed
241

MischaD's avatar
MischaD committed
242
    def _fill_rhs(self):
MischaD's avatar
MischaD committed
243
244
245
246
247
248
249
250
251
252
253
254
        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)

MischaD's avatar
MischaD committed
255
    def _svd(self):
MischaD's avatar
MischaD committed
256
257
258
        U, s, V = np.linalg.svd(self._lhs)
        rows = self._lhs.shape[0]
        cols = self._lhs.shape[1]
MischaD's avatar
MischaD committed
259
260
261
262

        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
263
        rhs = np.dot(np.transpose(U), self._rhs)
MischaD's avatar
MischaD committed
264
265
266
267
        if rank < rows:
            overlap = rows - rank
            if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
                # NO Solution
MischaD's avatar
MischaD committed
268
269
270
                raise NoSolutionException(
                    "Add shells or lower the order"
                )  # change shells ?
MischaD's avatar
MischaD committed
271
272
            rows = rank
            s.resize(rows)
MischaD's avatar
MischaD committed
273
            rhs.resize((rows, rhs.shape[1]), refcheck=False)
MischaD's avatar
MischaD committed
274
275
276
277
278
279
280
281
282
283
284
285
            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:
MischaD's avatar
MischaD committed
286
287
288
289
            raise InfiniteSolutionsException(
                "Reduce order, remove shells, or use continue.py in "
                "https://github.com/BDuenweg/Lattice-Boltzmann-weights"
            )
MischaD's avatar
MischaD committed
290
291
        self._solution = np.dot(np.transpose(V), reduced_rhs)
        return self._solution
MischaD's avatar
MischaD committed
292

MischaD's avatar
MischaD committed
293
    def _calculate_coefficients(self):
MischaD's avatar
MischaD committed
294
295
296
297
        """Use _all_velocity_vectors and _solution to calculate the coefficients of the
        weigth polynomials. Writes them into _coefficients
        """

MischaD's avatar
MischaD committed
298
299
300
301
302
        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
303
304
305
            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
306

MischaD's avatar
MischaD committed
307
        coeffs = np.zeros((self._solution.shape[0] + 1, self._solution.shape[1] + 1))
MischaD's avatar
Bugfix    
MischaD committed
308
        coeffs[0, :] = coefficients[-1::-1]
MischaD's avatar
MischaD committed
309
310
311
312
        coeffs[1:, :-1] = self._solution[:, -1::-1]
        self._coefficients = coeffs
        return self._coefficients

MischaD's avatar
MischaD committed
313
    def _calculate_valid_interval(self):
MischaD's avatar
MischaD committed
314
        if not self._weight_polynomials:
MischaD's avatar
MischaD committed
315
316
317
            raise UninitializedAttributeException(
                "Weight Polynomials not give. Calculate them first!"
            )
MischaD's avatar
MischaD committed
318
319
320
321
322
323
324
325
326

        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):
MischaD's avatar
MischaD committed
327
                mesh_point = (roots[i + 1] + roots[i]) / 2
MischaD's avatar
MischaD committed
328
329
330
331
                if equation.eval(mesh_point) < 0:
                    continue
                # positive weights
                if i == 0:
MischaD's avatar
MischaD committed
332
333
334
                    cur_interval = sp.Union(
                        cur_interval, sp.Interval(-oo, roots[i + 1])
                    )
MischaD's avatar
MischaD committed
335
336
337
                elif i == (len(roots) - 2):
                    cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], oo))
                else:
MischaD's avatar
MischaD committed
338
339
340
                    cur_interval = sp.Union(
                        cur_interval, sp.Interval(roots[i], roots[i + 1])
                    )
MischaD's avatar
MischaD committed
341
342
343
344
            interval = interval.intersect(cur_interval)
        self._interval = sp.Complement(interval, sp.FiniteSet(0))
        return self._interval

MischaD's avatar
MischaD committed
345
346
347
348
349
350
351
    def solution_expected(self):
        """Evaluate if the order, dimension and amount of shells entered are expected to retrieve a solution.
        This is in no way a guarantee on the existence of a solution or whether or not the program
        will finish successful, but it can be used as a quick hint on the existence of a solution.

        Changes the state of self._solution and self._velocity vectors for future calculations."""

MischaD's avatar
MischaD committed
352
353
354
355
356
357
358
359
360
        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]

MischaD's avatar
MischaD committed
361
362
363
        self._calculate_velocity_vectors()
        self._fill_lhs()
        self._fill_rhs()
MischaD's avatar
MischaD committed
364

MischaD's avatar
MischaD committed
365
        self._svd()
MischaD's avatar
MischaD committed
366

MischaD's avatar
MischaD committed
367
        return self._solution is not None
MischaD's avatar
Bugfix    
MischaD committed
368

MischaD's avatar
MischaD committed
369
370
    def calculate_polynomials(self, approximate=True):
        """Main algorithm as proposed by D. Spiller and B Duenweg.
MischaD's avatar
MischaD committed
371

MischaD's avatar
MischaD committed
372
373
374
        Divided into the part where the existence of a solution and calculated together with all the velocity
        vectors. If one is expected the coefficients are calculated, the interval evaluated and the weight
        polynomials returned.
MischaD's avatar
MischaD committed
375

MischaD's avatar
MischaD committed
376
377
        :param approximate: Use python Fraction to approximate the calculated polynomial coefficients as rational
         number. This should give better results, since they should have a rational representation.
378

MischaD's avatar
MischaD committed
379
        :return: List of Sympy Polynomials."""
380

MischaD's avatar
MischaD committed
381
382
383
384
385
        self.solution_expected()
        self._calculate_coefficients()

        c_s_sq = sp.Symbol("c_s_sq")
        if approximate:
MischaD's avatar
MischaD committed
386
387
388
389
            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
390
        else:
MischaD's avatar
MischaD committed
391
392
393
            self._weight_polynomials = [
                sp.Poly(pol_coffs, c_s_sq) for pol_coffs in self._coefficients
            ]
MischaD's avatar
MischaD committed
394
395
396

        self._interval = self._calculate_valid_interval()
        return [] if self._interval == sp.EmptySet else self._weight_polynomials
397

398
    def calculate_weights(self, c_s_sq=None, boundary=None):
399
        """
400
401
402
403
404
405
        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.

406
407
408
409
410
        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
MischaD's avatar
MischaD committed
411
        :return: weights With zero weight(s) included
412
        """
413

414
415
416
417
418
        if not self._weight_polynomials:
            self.calculate_polynomials()
        if not self._interval:  # Empty Interval
            return []
        if not c_s_sq:
419
420
            if not boundary:
                boundary = self._boundary
421
422
423
424
            if boundary == "inf":
                c_s_sq = self._interval.inf
            else:
                c_s_sq = self._interval.sup
425

426
427
428
        if self._interval.contains(c_s_sq) is sp.EmptySet:  # Invalid c_s_sq input
            return []

MischaD's avatar
MischaD committed
429
        weights = [x.eval(c_s_sq) for x in self._weight_polynomials]
MischaD's avatar
MischaD committed
430
        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
431
432
433
434
        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
435
        return weights
MischaD's avatar
MischaD committed
436

437
438
    def velocity_set(self):
        """
MischaD's avatar
MischaD committed
439
440
        Removes unwanted shells, i.e. weight == = and
        returns the values that are important for the Lattice Boltzmann Model, i.e.
441
        a list of the velocity values their corresponding weights and the value for the speed of sound.
MischaD's avatar
MischaD committed
442
443

        :return: (c_s_sq, weights, velocities)
444
        """
MischaD's avatar
MischaD committed
445
446
        if not self.weights:
            self.calculate_weights()
447
448
449
450
451
452

        weights = []
        for shell in self.shells:
            if shell.weight == 0:
                continue
            weights += [shell.weight] * shell.size
MischaD's avatar
MischaD committed
453
        velocities = np.zeros((self._dimension, len(weights)), dtype=np.int32)
454
455
456
457
458
        i = 0
        for shell in self.shells:
            if shell.weight == 0:
                continue
            for velocity in shell.velocities:
MischaD's avatar
MischaD committed
459
                velocities[:, i] = velocity
460
461
462
                i += 1

        return self._c_s_sq, weights, velocities