lattice.py 38.3 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
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,
MischaD's avatar
MischaD committed
21
    ps_installed,
MischaD's avatar
MischaD committed
22
23
24
    get_random_vector,
    velocities_for_shell,
)
25
from .shell import Shell
MischaD's avatar
MischaD committed
26
27


MischaD's avatar
MischaD committed
28
class Lattice:
MischaD's avatar
MischaD committed
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
    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=[],
mischa's avatar
mischa committed
63
        eq_order=None,
MischaD's avatar
MischaD committed
64
65
        svd_tolerance=1e-8,
    ):
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.
mischa's avatar
mischa committed
81
82
            eq_order (int): Explicitly specify the order of accuracy for the equilibrium distribution function.
                Only relevant if AssignmentCollections are used. For 2D resonable defaults are used.
83
84
            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
85

MischaD's avatar
MischaD committed
86
        """
MischaD's avatar
MischaD committed
87
        self._dimension = dimension
MischaD's avatar
MischaD committed
88
89
90
        self._order = order
        self._shell_list = shell_list
        self._seed = seed
mischa's avatar
mischa committed
91
        self._eq_order = eq_order
MischaD's avatar
MischaD committed
92
        self._svd_tolerance = svd_tolerance
MischaD's avatar
MischaD committed
93
        random.seed(self._seed)
94
        self._debug = False
MischaD's avatar
MischaD committed
95

MischaD's avatar
MischaD committed
96
        # Preperation
MischaD's avatar
MischaD committed
97
        self._tensor_space_dimension = 0
MischaD's avatar
MischaD committed
98
        self._all_velocity_vectors = []
MischaD's avatar
MischaD committed
99
        self._possible_tensors = []
MischaD's avatar
MischaD committed
100
        self._tensor_dimensions = []
MischaD's avatar
MischaD committed
101
        self._shells = 0
102
        self._unwanted_subshells = unwanted_subshells
MischaD's avatar
MischaD committed
103

MischaD's avatar
MischaD committed
104
105
106
107
        # LES
        self._lhs = []
        self._rhs = []

MischaD's avatar
MischaD committed
108
109
        # SVD
        self._singular_values = []
MischaD's avatar
MischaD committed
110
        self._solution = []
MischaD's avatar
MischaD committed
111

112
        # Polynomials
MischaD's avatar
MischaD committed
113
        self._coefficients = []
MischaD's avatar
MischaD committed
114
115
        self._weight_polynomials = []
        self._interval = None
116
        self._boundary = boundary if boundary == "sup" else "inf"
MischaD's avatar
MischaD committed
117

118
119
        # Final values
        self._discrete_velocities = []
120
        self._reduced_weights = []
121
        self._weights = []
122
        self._velocities = tuple()
mischa's avatar
mischa committed
123
        self._q = None
124
125
        self._c_s_sq = None

MischaD's avatar
MischaD committed
126
127
    def __str__(self):
        string = "D{} - Order: {} - Shells: {}".format(
MischaD's avatar
MischaD committed
128
            self._dimension, self._order, str(self._shell_list)
MischaD's avatar
MischaD committed
129
        )
MischaD's avatar
MischaD committed
130
131
132
        if self._boundary:
            string += " - Boundary: {}".format(self._boundary)
        if self._unwanted_subshells:
MischaD's avatar
MischaD committed
133
134
135
            string += " - Unwanted Subshells: {}".format(
                "".join([s + ", " for s in self._unwanted_subshells])
            )
MischaD's avatar
MischaD committed
136
137
138
        return string

    @classmethod
MischaD's avatar
MischaD committed
139
    def from_name(cls, name):
140
141
142
143
144
        """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
145
        return cls(**kwargs)
146
147

    @classmethod
MischaD's avatar
MischaD committed
148
    def from_order(cls, dimension, order):
MischaD's avatar
MischaD committed
149
150
151
152
153
        """
        :param dimension:
        :param order:
        :return:
        """
MischaD's avatar
MischaD committed
154
155
156
157
158
159
160
161
        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
162
            "No Lattice known with order {} and dimension {}.".format(order, dimension)
MischaD's avatar
MischaD committed
163
        )
MischaD's avatar
MischaD committed
164
        return cls(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
MischaD's avatar
MischaD committed
165

MischaD's avatar
MischaD committed
166
167
168
169
    @property
    def tensor_space_dimension(self):
        return np.array([len(x) for x in self._possible_tensors]).sum()

MischaD's avatar
MischaD committed
170
171
    @property
    def shells(self):
172
        return self._shells
MischaD's avatar
MischaD committed
173

mischa's avatar
mischa committed
174
175
176
177
178
179
180
181
182
183
184
    @property
    def q(self):
        if self._q:
            return self._q
        self._q = len(self.velocity_set()[1])
        return self._q

    @property
    def c_s_sq(self):
        if self._c_s_sq:
            return self._c_s_sq
mischa's avatar
mischa committed
185
186
        self._c_s_sq = self.velocity_set()[0]
        return self.c_s_sq
mischa's avatar
mischa committed
187
188
189

    @property
    def weights(self):
190
191
192
        if not self._weights:
            self.velocity_set()
        return self._weights
mischa's avatar
mischa committed
193

MischaD's avatar
MischaD committed
194
195
    @property
    def velocities(self):
196
197
198
199
200
201
202
203
204
205
206
        """Return weights and velocities in walberla fashion"""
        if not self._velocities:
            self.velocity_set()
        return self._velocities

    @property
    def reduced_weights(self):
        if self._reduced_weights:
            return self._reduced_weights
        return self.calculate_weights()

MischaD's avatar
MischaD committed
207

208
    @property
mischa's avatar
mischa committed
209
210
211
212
213
214
    def eq_order(self):
        """Calculate the eq order based on the amount of velocities given and dimension.
        According to Philippi et. al, 17 and 37 velocities are the minimum amount respectively
        to use for third and fourth order accurate models. """
        if self._eq_order:
            return self._eq_order
215

mischa's avatar
mischa committed
216
217
218
219
220
221
222
223
224
225
226
        if self._dimension == 2:
            if self.q >= 37:
                self._eq_order = 4
            elif self.q >= 17:
                self._eq_order = 3
            else:
                self._eq_order = 2
        if self._dimension == 3:
            # Should be done explicitly for three dimensions using the constructor
            self._eq_order = 2
        return self._eq_order
MischaD's avatar
MischaD committed
227

mischa's avatar
mischa committed
228
229
230
231
232
    @property
    def necessary_ghost_layers(self):
        velocities = self.velocities
        return np.array(velocities).max()

MischaD's avatar
MischaD committed
233
234
    @property
    @ps_installed
mischa's avatar
mischa committed
235
    def assignment_collection(self, ps):
MischaD's avatar
MischaD committed
236
237
        if self._dimension == 3:
            return self._assignment_collection3d()
mischa's avatar
mischa committed
238
239
240
241
        symbolic_description = []
        c_s_sq = self.c_s_sq
        weights = self.weights
        velocities = self.velocities
mischa's avatar
mischa committed
242
        Q = self.q
mischa's avatar
mischa committed
243
244
245
246
        weights = [complex(x).real for x in weights]
        c_ix = [vel[0] for vel in velocities]
        c_iy = [vel[1] for vel in velocities]

mischa's avatar
mischa committed
247
        src, dst = ps.fields("src({Q}), dst({Q}): double[2D]".format(Q=Q))
MischaD's avatar
MischaD committed
248
249
        density, u_x, u_y, omega = sp.symbols("density, u_x, u_y, omega")

mischa's avatar
mischa committed
250
251
252
253
        for i in range(Q):
            symbolic_description.append(
                ps.Assignment(dst[0, 0](i), src[c_ix[i], c_iy[i]](i), )
            )
MischaD's avatar
MischaD committed
254

mischa's avatar
mischa committed
255
256
257
258
259
260
261
        density_formula = sum([dst[0, 0](i) for i in range(Q)])
        vel_x_formula = (sum([dst[0, 0](i) * c_ix[i] for i in range(Q)])) / density
        vel_y_formula = (sum([dst[0, 0](i) * c_iy[i] for i in range(Q)])) / density

        symbolic_description.append(ps.Assignment(density, density_formula))
        symbolic_description.append(ps.Assignment(u_x, vel_x_formula))
        symbolic_description.append(ps.Assignment(u_y, vel_y_formula))
MischaD's avatar
MischaD committed
262

263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
        if not self._debug:
            for i in range(Q):
                feq_formula = weights[i] * density * (
                        1
                        + (u_x * c_ix[i] + u_y * c_iy[i]) / c_s_sq
                        + (u_x * c_ix[i] + u_y * c_iy[i]) ** 2 / (2 * c_s_sq ** 2)
                        - (u_x ** 2 + u_y ** 2) / (2 * c_s_sq)
                )
                if self.eq_order >= 3:
                    feq_formula += weights[i] * density * (
                            -1 / (2 * c_s_sq ** 2) * (u_x * c_ix[i] + u_y * c_iy[i]) * (u_x ** 2 + u_y ** 2)
                            + 4 / 3 * (1 / (2 * c_s_sq)) ** 3 * (u_x * c_ix[i] + u_y * c_iy[i]) ** 3
                    )
                if self.eq_order >= 4:
                    feq_formula += weights[i] * density * \
                                   (1 / (6 * c_s_sq ** 2) * (
                                           (u_x * c_ix[i] + u_y * c_iy[i]) ** 4 / (4 * c_s_sq ** 2) - (
                                           3 * (u_x ** 2 + u_y ** 2) *
mischa's avatar
mischa committed
281
                                           (u_x * c_ix[i] + u_y * c_iy[i]) ** 2) / (2 * c_s_sq) + 3/4 * (
282
283
284
285
                                                   u_x ** 4 + u_y ** 4)))

                symbolic_description.append(
                    ps.Assignment(dst[0, 0](i), omega * feq_formula + (1 - omega) * dst[0, 0](i))
MischaD's avatar
MischaD committed
286
                )
Mischa Dombrowski's avatar
Cleanup    
Mischa Dombrowski committed
287

mischa's avatar
mischa committed
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
        if self._debug:
            raise AssertionError("debug in lbmweights not valid")
            print("debug")
            sym = symbolic_description
            sym.append(ps.Assignment(dst[0, 0](0),
                                     -5 * np.sqrt(193) * density * u_x ** 2 / 48 + 193 * density * u_x
                                          ** 2 / 720 - 5 * np.sqrt(
                                         193) * density * u_y ** 2 / 48 + 193 * density * u_y
                                          ** 2 / 720 + 23 * density / 324 + 193 * np.sqrt(193) * density
                                          / 8100))
            sym.append(ps.Assignment(dst[0, 0](1),
                                     density * u_x ** 2 * u_y / 2 - 1973 * density * u_x
                                          ** 2 / 4800 + 43 * np.sqrt(193) * density * u_x
                                          ** 2 / 960 - 3 * np.sqrt(
                                         193) * density * u_y ** 3 / 80 - density * u_y
                                          ** 3 / 8 - 773 * density * u_y ** 2 / 4800 + 11 * np.sqrt(
                                         193) * density * u_y ** 2 / 192 - 307 * density * u_y
                                          / 1200 - np.sqrt(193) * density * u_y / 240 - 91 * np.sqrt(
                                         193) * density / 18000 + 671 * density / 3600))
            sym.append(ps.Assignment(dst[0, 0](2),
                                     -3 * np.sqrt(193) * density * u_x ** 3 / 80 - density * u_x
                                          ** 3 / 8 - 773 * density * u_x ** 2 / 4800 + 11 * np.sqrt(
                                         193) * density * u_x ** 2 / 192 + density * u_x * u_y
                                          ** 2 / 2 - 307 * density * u_x / 1200 - np.sqrt(
                                         193) * density * u_x / 240 - 1973 * density * u_y
                                          ** 2 / 4800 + 43 * np.sqrt(193) * density * u_y
                                          ** 2 / 960 - 91 * np.sqrt(
                                         193) * density / 18000 + 671 * density / 3600))
            sym.append(ps.Assignment(dst[0, 0](3),
                                     density * u_x ** 3 / 8 + 3 * np.sqrt(193) * density * u_x
                                          ** 3 / 80 - 773 * density * u_x ** 2 / 4800 + 11 * np.sqrt(
                                         193) * density * u_x ** 2 / 192 - density * u_x * u_y
                                          ** 2 / 2 + np.sqrt(
                                         193) * density * u_x / 240 + 307 * density * u_x
                                          / 1200 - 1973 * density * u_y ** 2 / 4800 + 43 * np.sqrt(
                                         193) * density * u_y ** 2 / 960 - 91 * np.sqrt(193) * density
                                          / 18000 + 671 * density / 3600))
            sym.append(ps.Assignment(dst[0, 0](4),
                                     -density * u_x ** 2 * u_y / 2 - 1973 * density * u_x
                                          ** 2 / 4800 + 43 * np.sqrt(
                                         193) * density * u_x ** 2 / 960 + density * u_y
                                          ** 3 / 8 + 3 * np.sqrt(
                                         193) * density * u_y ** 3 / 80 - 773 * density * u_y
                                          ** 2 / 4800 + 11 * np.sqrt(193) * density * u_y
                                          ** 2 / 192 + np.sqrt(
                                         193) * density * u_y / 240 + 307 * density * u_y
                                          / 1200 - 91 * np.sqrt(193) * density / 18000 + 671 * density
                                          / 3600))
            sym.append(ps.Assignment(dst[0, 0](5), density * u_x ** 3 / 9 + np.sqrt(193) * density * u_x
                 ** 3 / 45 - density * u_x ** 2 * u_y / 4 - np.sqrt(
                193) * density * u_x ** 2 / 40 + 41 * density * u_x ** 2 / 200 - density *
                                     u_x * u_y ** 2 / 4 + np.sqrt(
                193) * density * u_x * u_y / 120 + 391 * density * u_x * u_y
                                          / 3000 - 91 * density * u_x / 1800 - np.sqrt(
                193) * density * u_x / 360 + density * u_y ** 3 / 9 + np.sqrt(193) * density *
                                     u_y ** 3 / 45 - np.sqrt(
                193) * density * u_y ** 2 / 40 + 41 * density * u_y ** 2 / 200 - 91 * density
                                          * u_y / 1800 - np.sqrt(
                193) * density * u_y / 360 + 17 * np.sqrt(193) * density / 27000 + 131 * density
                                          / 5400))
            sym.append(ps.Assignment(dst[0, 0](6),
                                     -np.sqrt(193) * density * u_x ** 3 / 45 - density * u_x
                                          ** 3 / 9 - density * u_x ** 2 * u_y / 4 - np.sqrt(
                                         193) * density * u_x ** 2 / 40 + 41 * density * u_x
                                          ** 2 / 200 + density * u_x * u_y ** 2 / 4 - 391 * density
                                          * u_x * u_y / 3000 - np.sqrt(
                                         193) * density * u_x * u_y / 120 + np.sqrt(193) * density *
                                     u_x / 360 + 91 * density * u_x / 1800 + density * u_y
                                          ** 3 / 9 + np.sqrt(
                                         193) * density * u_y ** 3 / 45 - np.sqrt(
                                         193) * density * u_y ** 2 / 40 + 41 * density * u_y
                                          ** 2 / 200 - 91 * density * u_y / 1800 - np.sqrt(
                                         193) * density * u_y / 360 + 17 * np.sqrt(193) * density
                                          / 27000 + 131 * density / 5400))
            sym.append(ps.Assignment(dst[0, 0](7), density * u_x ** 3 / 9 + np.sqrt(193) * density * u_x
                 ** 3 / 45 + density * u_x ** 2 * u_y / 4 - np.sqrt(
                193) * density * u_x ** 2 / 40 + 41 * density * u_x ** 2 / 200 - density *
                                     u_x * u_y ** 2 / 4 - 391 * density * u_x * u_y
                                          / 3000 - np.sqrt(
                193) * density * u_x * u_y / 120 - 91 * density * u_x / 1800 - np.sqrt(193) *
                                     density * u_x / 360 - np.sqrt(
                193) * density * u_y ** 3 / 45 - density * u_y ** 3 / 9 - np.sqrt(
                193) * density * u_y ** 2 / 40 + 41 * density * u_y ** 2 / 200 + np.sqrt(
                193) * density * u_y / 360 + 91 * density * u_y / 1800 + 17 * np.sqrt(193) *
                                     density / 27000 + 131 * density / 5400))
            sym.append(ps.Assignment(dst[0, 0](8),
                                     -np.sqrt(193) * density * u_x ** 3 / 45 - density * u_x
                                          ** 3 / 9 + density * u_x ** 2 * u_y / 4 - np.sqrt(
                                         193) * density * u_x ** 2 / 40 + 41 * density * u_x
                                          ** 2 / 200 + density * u_x * u_y ** 2 / 4 + np.sqrt(
                                         193) * density * u_x * u_y / 120 + 391 * density * u_x
                                          * u_y / 3000 + np.sqrt(
                                         193) * density * u_x / 360 + 91 * density * u_x
                                          / 1800 - np.sqrt(193) * density * u_y ** 3 / 45 - density *
                                     u_y ** 3 / 9 - np.sqrt(
                                         193) * density * u_y ** 2 / 40 + 41 * density * u_y
                                          ** 2 / 200 + np.sqrt(
                                         193) * density * u_y / 360 + 91 * density * u_y
                                          / 1800 + 17 * np.sqrt(193) * density / 27000 + 131 * density
                                          / 5400))
            sym.append(ps.Assignment(dst[0, 0](9),
                                     -np.sqrt(193) * density * u_x ** 3 / 360 - density * u_x
                                          ** 3 / 72 + density * u_x ** 2 / 4800 + np.sqrt(
                                         193) * density * u_x ** 2 / 960 - np.sqrt(
                                         193) * density * u_x * u_y / 480 + 359 * density * u_x
                                          * u_y / 12000 - 71 * density * u_x / 3600 + np.sqrt(
                                         193) * density * u_x / 720 - np.sqrt(193) * density * u_y
                                          ** 3 / 360 - density * u_y ** 3 / 72 + density * u_y
                                          ** 2 / 4800 + np.sqrt(
                                         193) * density * u_y ** 2 / 960 - 71 * density * u_y
                                          / 3600 + np.sqrt(193) * density * u_y / 720 - 49 * np.sqrt(
                                         193) * density / 54000 + 137 * density / 10800))
            sym.append(ps.Assignment(dst[0, 0](10), density * u_x ** 3 / 72 + np.sqrt(193) * density * u_x
                 ** 3 / 360 + density * u_x ** 2 / 4800 + np.sqrt(
                193) * density * u_x ** 2 / 960 - 359 * density * u_x * u_y / 12000 + np.sqrt(
                193) * density * u_x * u_y / 480 - np.sqrt(
                193) * density * u_x / 720 + 71 * density * u_x / 3600 - np.sqrt(
                193) * density * u_y ** 3 / 360 - density * u_y ** 3 / 72 + density * u_y
                                          ** 2 / 4800 + np.sqrt(
                193) * density * u_y ** 2 / 960 - 71 * density * u_y / 3600 + np.sqrt(193) *
                                     density * u_y / 720 - 49 * np.sqrt(
                193) * density / 54000 + 137 * density / 10800))
            sym.append(ps.Assignment(dst[0, 0](11),
                                     -np.sqrt(193) * density * u_x ** 3 / 360 - density * u_x
                                          ** 3 / 72 + density * u_x ** 2 / 4800 + np.sqrt(
                                         193) * density * u_x ** 2 / 960 - 359 * density * u_x *
                                     u_y / 12000 + np.sqrt(
                                         193) * density * u_x * u_y / 480 - 71 * density * u_x
                                          / 3600 + np.sqrt(193) * density * u_x / 720 + density *
                                     u_y ** 3 / 72 + np.sqrt(
                                         193) * density * u_y ** 3 / 360 + density * u_y
                                          ** 2 / 4800 + np.sqrt(193) * density * u_y ** 2 / 960 - np.sqrt(
                                         193) * density * u_y / 720 + 71 * density * u_y
                                          / 3600 - 49 * np.sqrt(193) * density / 54000 + 137 * density
                                          / 10800))
            sym.append(ps.Assignment(dst[0, 0](12), density * u_x ** 3 / 72 + np.sqrt(193) * density * u_x
                 ** 3 / 360 + density * u_x ** 2 / 4800 + np.sqrt(
                193) * density * u_x ** 2 / 960 - np.sqrt(193) * density * u_x * u_y
                                          / 480 + 359 * density * u_x * u_y / 12000 - np.sqrt(
                193) * density * u_x / 720 + 71 * density * u_x / 3600 + density * u_y
                                          ** 3 / 72 + np.sqrt(
                193) * density * u_y ** 3 / 360 + density * u_y ** 2 / 4800 + np.sqrt(193) *
                                     density * u_y ** 2 / 960 - np.sqrt(
                193) * density * u_y / 720 + 71 * density * u_y / 3600 - 49 * np.sqrt(193) *
                                     density / 54000 + 137 * density / 10800))
            sym.append(ps.Assignment(dst[0, 0](13),
                                     -np.sqrt(193) * density * u_x ** 2 / 2880 - density * u_x
                                          ** 2 / 14400 - density * u_y ** 3 / 72 + np.sqrt(
                                         193) * density * u_y ** 3 / 720 - np.sqrt(
                                         193) * density * u_y ** 2 / 576 + 133 * density * u_y
                                          ** 2 / 4800 - 77 * density * u_y / 3600 + np.sqrt(
                                         193) * density * u_y / 720 - 101 * np.sqrt(193) * density
                                          / 162000 + 289 * density / 32400))
            sym.append(ps.Assignment(dst[0, 0](14),
                                     -density * u_x ** 3 / 72 + np.sqrt(193) * density * u_x
                                          ** 3 / 720 - np.sqrt(
                                         193) * density * u_x ** 2 / 576 + 133 * density * u_x
                                          ** 2 / 4800 - 77 * density * u_x / 3600 + np.sqrt(
                                         193) * density * u_x / 720 - np.sqrt(193) * density * u_y
                                          ** 2 / 2880 - density * u_y ** 2 / 14400 - 101 * np.sqrt(
                                         193) * density / 162000 + 289 * density / 32400))
            sym.append(ps.Assignment(dst[0, 0](15),
                                     -np.sqrt(193) * density * u_x ** 3 / 720 + density * u_x
                                          ** 3 / 72 - np.sqrt(
                                         193) * density * u_x ** 2 / 576 + 133 * density * u_x
                                          ** 2 / 4800 - np.sqrt(
                                         193) * density * u_x / 720 + 77 * density * u_x
                                          / 3600 - np.sqrt(
                                         193) * density * u_y ** 2 / 2880 - density * u_y
                                          ** 2 / 14400 - 101 * np.sqrt(193) * density / 162000 + 289 * density
                                          / 32400))
            sym.append(ps.Assignment(dst[0, 0](16),
                                     -np.sqrt(193) * density * u_x ** 2 / 2880 - density * u_x
                                          ** 2 / 14400 - np.sqrt(
                                         193) * density * u_y ** 3 / 720 + density * u_y
                                          ** 3 / 72 - np.sqrt(
                                         193) * density * u_y ** 2 / 576 + 133 * density * u_y
                                          ** 2 / 4800 - np.sqrt(
                                         193) * density * u_y / 720 + 77 * density * u_y
                                          / 3600 - 101 * np.sqrt(193) * density / 162000 + 289 * density
                                          / 32400))

MischaD's avatar
MischaD committed
470
471
        return ps.AssignmentCollection(symbolic_description)

mischa's avatar
mischa committed
472
473
474
475
476
    @property
    @ps_installed
    def _assignment_collection3d(self, ps):
        c_s_sq, weights, velocities = self.velocity_set()
        Q = self.q
477
478
479
        c_ix = [vel[0] for vel in velocities]
        c_iy = [vel[1] for vel in velocities]
        c_iz = [vel[2] for vel in velocities]
mischa's avatar
mischa committed
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
        src, dst = ps.fields("src({Q}), dst({Q}): double[2D]".format(Q=Q))
        density, u_x, u_y, u_z, omega = sp.symbols("density, u_x, u_y, u_z, omega")

        density_formula = sum([src[0, 0](i) for i in range(Q)])
        vel_x_formula = (sum([src[0, 0](i) * c_ix[i] for i in range(Q)])) / density
        vel_y_formula = (sum([src[0, 0](i) * c_iy[i] for i in range(Q)])) / density
        vel_z_formula = (sum([src[0, 0](i) * c_iz[i] for i in range(Q)])) / density

        symbolic_description = [ps.Assignment(density, density_formula),
                                ps.Assignment(u_x, vel_x_formula),
                                ps.Assignment(u_y, vel_y_formula),
                                ps.Assignment(u_z, vel_z_formula),]

        for i in range(Q):
            feq_formula = weights[i] * density[0, 0] * (
                    1
                    + (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) / c_s_sq
                    + (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 2 / (2 * c_s_sq ** 2)
                    - (u_x ** 2 + u_y ** 2 + u_z ** 2) / (2 * c_s_sq)
            )
            if self.eq_order >= 3:
                feq_formula += weights[i] * density[0, 0] * (
                        -1 / (2 * c_s_sq ** 2) * (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) * (u_x ** 2 + u_y ** 2 + u_z ** 2)
                        + 4 / 3 * (1 / (2 * c_s_sq)) ** 3 * (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 3
                )
            if self.eq_order >= 4:
                feq_formula += weights[i] * density[0, 0] * \
                               (1 / (6 * c_s_sq ** 2) * (
                                       (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 4 / (4 * c_s_sq ** 2) - (
                                           3 * (u_x ** 2 + u_y ** 2 + u_z ** 2) *
                                           (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 2) / (2 * c_s_sq) + 3 * (
                                                   u_x ** 4 + u_y ** 4 + u_z ** 4)))

        for i in range(Q):
            symbolic_description.append(
                ps.Assignment(dst[-1 * c_iy[i], c_ix[i]](i), omega * feq_formula + (1 - omega) * src[0, 0](i))
            )
        return ps.AssignmentCollection(symbolic_description)
MischaD's avatar
MischaD committed
518

MischaD's avatar
MischaD committed
519
520
521
522
523
    def shell_from_type(self, type):
        for shell in self._shells:
            if shell.type == type:
                return shell

MischaD's avatar
MischaD committed
524
    def _calculate_velocity_vectors(self):
MischaD's avatar
MischaD committed
525
        velocity_vectors = []
MischaD's avatar
MischaD committed
526

MischaD's avatar
MischaD committed
527
        group = get_group(self._dimension)
MischaD's avatar
MischaD committed
528
        for shell in range(len(self._shell_list)):
MischaD's avatar
MischaD committed
529
530
531
            velocities = velocities_for_shell(
                dimension=self._dimension, shell=self._shell_list[shell]
            )
MischaD's avatar
MischaD committed
532
            if len(velocities) == 0:
MischaD's avatar
MischaD committed
533
534
535
                raise ImpossibleVelocityShellException(
                    "No velocity sets found for velocity shell {}.".format(shell)
                )
MischaD's avatar
MischaD committed
536

537
            # If e.g. dimension = 2, then shell 25 in ambiguous. Both shell (3,4) and (5,0) are added
MischaD's avatar
MischaD committed
538
539
            subshells = get_subshells(velocities, group)
            if len(subshells) > 1:
MischaD's avatar
MischaD committed
540
                # 3 4, 0 5 can only be given together
541
                velocities = []
MischaD's avatar
MischaD committed
542
                for i_subs, subshell in enumerate(subshells):
543
                    # get type for subshell. Example shell: all vectors of type [0, -1, 0] --> 100
MischaD's avatar
MischaD committed
544
545
546
                    type = "".join(
                        [str(int(x)) for x in sorted(abs(subshell[0]), reverse=True)]
                    )
547
548
                    if type not in self._unwanted_subshells:
                        velocities.append(subshell)
MischaD's avatar
MischaD committed
549
550
            else:
                velocities = [velocities]
MischaD's avatar
MischaD committed
551
            velocity_vectors.extend(velocities)
MischaD's avatar
MischaD committed
552

MischaD's avatar
MischaD committed
553
        self._shells = [Shell([np.array([0] * self._dimension)])]
554
555
        for final_velocities in velocity_vectors:
            self._shells.append(Shell(final_velocities))
MischaD's avatar
MischaD committed
556
557
        self._all_velocity_vectors = velocity_vectors
        return velocity_vectors
MischaD's avatar
MischaD committed
558

MischaD's avatar
MischaD committed
559
    def _fill_lhs(self):
MischaD's avatar
MischaD committed
560
561
562
563
        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
564
                vector = get_random_vector(self._dimension)
MischaD's avatar
MischaD committed
565
566
                # vector = 2 * np.random.rand(self._dimension) - 1
                # vector = vector / np.linalg.norm(vector)
MischaD's avatar
MischaD committed
567
568

                rows = []
MischaD's avatar
MischaD committed
569
                for shell_velocities in self._all_velocity_vectors:
MischaD's avatar
MischaD committed
570
571
572
                    shell_sum = lattice_sum(vector, shell_velocities, rank)
                    rows.append(shell_sum)
                lhs.append(rows)
MischaD's avatar
MischaD committed
573
        self._lhs = np.array(lhs)
MischaD's avatar
MischaD committed
574
        return np.array(lhs)
MischaD's avatar
MischaD committed
575

MischaD's avatar
MischaD committed
576
    def _fill_rhs(self):
MischaD's avatar
MischaD committed
577
578
579
580
581
582
583
584
585
586
587
588
        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
589
    def _svd(self):
MischaD's avatar
MischaD committed
590
591
592
        U, s, V = np.linalg.svd(self._lhs)
        rows = self._lhs.shape[0]
        cols = self._lhs.shape[1]
MischaD's avatar
MischaD committed
593
594
595
596

        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
597
        rhs = np.dot(np.transpose(U), self._rhs)
MischaD's avatar
MischaD committed
598
599
600
601
        if rank < rows:
            overlap = rows - rank
            if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
                # NO Solution
MischaD's avatar
MischaD committed
602
603
604
                raise NoSolutionException(
                    "Add shells or lower the order"
                )  # change shells ?
MischaD's avatar
MischaD committed
605
606
            rows = rank
            s.resize(rows)
MischaD's avatar
MischaD committed
607
            rhs.resize((rows, rhs.shape[1]), refcheck=False)
MischaD's avatar
MischaD committed
608
609
610
611
612
613
614
615
616
617
618
619
            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
620
621
622
623
            raise InfiniteSolutionsException(
                "Reduce order, remove shells, or use continue.py in "
                "https://github.com/BDuenweg/Lattice-Boltzmann-weights"
            )
MischaD's avatar
MischaD committed
624
625
        self._solution = np.dot(np.transpose(V), reduced_rhs)
        return self._solution
MischaD's avatar
MischaD committed
626

MischaD's avatar
MischaD committed
627
    def _calculate_coefficients(self):
MischaD's avatar
MischaD committed
628
629
630
631
        """Use _all_velocity_vectors and _solution to calculate the coefficients of the
        weigth polynomials. Writes them into _coefficients
        """

MischaD's avatar
MischaD committed
632
633
634
635
636
        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
637
638
639
            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
640

MischaD's avatar
MischaD committed
641
        coeffs = np.zeros((self._solution.shape[0] + 1, self._solution.shape[1] + 1))
MischaD's avatar
Bugfix    
MischaD committed
642
        coeffs[0, :] = coefficients[-1::-1]
MischaD's avatar
MischaD committed
643
644
645
646
        coeffs[1:, :-1] = self._solution[:, -1::-1]
        self._coefficients = coeffs
        return self._coefficients

MischaD's avatar
MischaD committed
647
    def _calculate_valid_interval(self):
MischaD's avatar
MischaD committed
648
        if not self._weight_polynomials:
MischaD's avatar
MischaD committed
649
650
651
            raise UninitializedAttributeException(
                "Weight Polynomials not give. Calculate them first!"
            )
MischaD's avatar
MischaD committed
652
653
654
655
656
657
658
659
660

        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
661
                mesh_point = (roots[i + 1] + roots[i]) / 2
MischaD's avatar
MischaD committed
662
663
664
665
                if equation.eval(mesh_point) < 0:
                    continue
                # positive weights
                if i == 0:
MischaD's avatar
MischaD committed
666
667
668
                    cur_interval = sp.Union(
                        cur_interval, sp.Interval(-oo, roots[i + 1])
                    )
MischaD's avatar
MischaD committed
669
670
671
                elif i == (len(roots) - 2):
                    cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], oo))
                else:
MischaD's avatar
MischaD committed
672
673
674
                    cur_interval = sp.Union(
                        cur_interval, sp.Interval(roots[i], roots[i + 1])
                    )
MischaD's avatar
MischaD committed
675
676
677
678
            interval = interval.intersect(cur_interval)
        self._interval = sp.Complement(interval, sp.FiniteSet(0))
        return self._interval

MischaD's avatar
MischaD committed
679
680
681
682
683
684
685
    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
686
687
688
689
690
691
692
693
694
        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
695
696
697
        self._calculate_velocity_vectors()
        self._fill_lhs()
        self._fill_rhs()
MischaD's avatar
MischaD committed
698

MischaD's avatar
MischaD committed
699
        self._svd()
MischaD's avatar
MischaD committed
700

MischaD's avatar
MischaD committed
701
        return self._solution is not None
MischaD's avatar
Bugfix    
MischaD committed
702

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

MischaD's avatar
MischaD committed
706
707
708
        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
709

MischaD's avatar
MischaD committed
710
711
        :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.
712

MischaD's avatar
MischaD committed
713
        :return: List of Sympy Polynomials."""
714

MischaD's avatar
MischaD committed
715
716
717
        self.solution_expected()
        self._calculate_coefficients()

mischa's avatar
mischa committed
718
        c_s_sq = sp.Symbol("c_s_sq", real=True)
MischaD's avatar
MischaD committed
719
        if approximate:
MischaD's avatar
MischaD committed
720
721
722
723
            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
724
        else:
MischaD's avatar
MischaD committed
725
726
727
            self._weight_polynomials = [
                sp.Poly(pol_coffs, c_s_sq) for pol_coffs in self._coefficients
            ]
MischaD's avatar
MischaD committed
728
729
730

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

732
    def calculate_weights(self, c_s_sq=None, boundary=None):
733
        """
734
735
736
737
738
739
        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.

740
741
742
743
744
        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
745
        :return: weights With zero weight(s) included
746
        """
747

748
749
750
751
752
        if not self._weight_polynomials:
            self.calculate_polynomials()
        if not self._interval:  # Empty Interval
            return []
        if not c_s_sq:
753
754
            if not boundary:
                boundary = self._boundary
755
756
757
758
            if boundary == "inf":
                c_s_sq = self._interval.inf
            else:
                c_s_sq = self._interval.sup
759

760
761
762
        if self._interval.contains(c_s_sq) is sp.EmptySet:  # Invalid c_s_sq input
            return []

MischaD's avatar
MischaD committed
763
        weights = [x.eval(c_s_sq) for x in self._weight_polynomials]
MischaD's avatar
MischaD committed
764
        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
765
766
        for i, shell in enumerate(self.shells):
            shell.set_weight(c_s_sq, weights[i])
mischa's avatar
mischa committed
767
768
769
770
771
        if isinstance(c_s_sq, sp.CRootOf):
            c_s_sq = complex(c_s_sq)
            if c_s_sq.imag != 0:
                raise ValueError("Imaginary root returned. This should not happen")
            c_s_sq = float(c_s_sq.real)
772
773
774
        if c_s_sq < 0:
            raise ValueError("Negative Value for c_s_sq")

775
        self._c_s_sq = c_s_sq
776
        self._reduced_weights = weights
MischaD's avatar
MischaD committed
777
        return weights
MischaD's avatar
MischaD committed
778

779
780
    def velocity_set(self):
        """
MischaD's avatar
MischaD committed
781
782
        Removes unwanted shells, i.e. weight == = and
        returns the values that are important for the Lattice Boltzmann Model, i.e.
783
        a list of the velocity values their corresponding weights and the value for the speed of sound.
MischaD's avatar
MischaD committed
784
785

        :return: (c_s_sq, weights, velocities)
786
        """
787
        if not self.reduced_weights:
MischaD's avatar
MischaD committed
788
            self.calculate_weights()
789
790
791
792
793
794

        weights = []
        for shell in self.shells:
            if shell.weight == 0:
                continue
            weights += [shell.weight] * shell.size
MischaD's avatar
MischaD committed
795
        velocities = np.zeros((self._dimension, len(weights)), dtype=np.int32)
796
797
798
799
800
        i = 0
        for shell in self.shells:
            if shell.weight == 0:
                continue
            for velocity in shell.velocities:
MischaD's avatar
MischaD committed
801
                velocities[:, i] = velocity
802
                i += 1
803
804
        self._weights = weights
        self._velocities = tuple(e for e in zip(velocities[0], velocities[1]))
805
806

        return self._c_s_sq, weights, velocities