lattice.py 15.4 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
        )
MischaD's avatar
MischaD committed
79
80
81
82
        if self._boundary:
            string += " - Boundary: {}".format(self._boundary)
        if self._unwanted_subshells:
            string += " - Unwanted Subshells: {}".format(''.join([s + ', ' for s in self._unwanted_subshells]))
MischaD's avatar
MischaD committed
83
84
85
        return string

    @classmethod
86
87
88
89
90
91
92
93
94
95
    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
96
97
98
99
100
        """
        :param dimension:
        :param order:
        :return:
        """
MischaD's avatar
MischaD committed
101
102
103
104
105
106
107
108
109
110
        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(
            "No Lattice known with order {} and dimension {}.".format(
                order, dimension
MischaD's avatar
MischaD committed
111
            )
MischaD's avatar
MischaD committed
112
        )
MischaD's avatar
MischaD committed
113
        return cls(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
MischaD's avatar
MischaD committed
114

MischaD's avatar
MischaD committed
115
116
117
118
    @property
    def tensor_space_dimension(self):
        return np.array([len(x) for x in self._possible_tensors]).sum()

MischaD's avatar
MischaD committed
119
120
    @property
    def shells(self):
121
        return self._shells
MischaD's avatar
MischaD committed
122
123
124

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

131
132
133
134
135
136
    @property
    def weights(self):
        if self._weights:
            return self._weights
        return self.calculate_weights()

MischaD's avatar
MischaD committed
137
138
139
140
141
    def shell_from_type(self, type):
        for shell in self._shells:
            if shell.type == type:
                return shell

MischaD's avatar
MischaD committed
142
    def _calculate_velocity_vectors(self):
MischaD's avatar
MischaD committed
143
        velocity_vectors = []
MischaD's avatar
MischaD committed
144

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

151
            # If e.g. dimension = 2, then shell 25 in ambiguous. Both shell (3,4) and (5,0) are added
MischaD's avatar
MischaD committed
152
153
            subshells = get_subshells(velocities, group)
            if len(subshells) > 1:
MischaD's avatar
MischaD committed
154
                # 3 4, 0 5 can only be given together
155
                velocities = []
MischaD's avatar
MischaD committed
156
                for i_subs, subshell in enumerate(subshells):
157
158
159
160
                    # 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
161
162
            else:
                velocities = [velocities]
MischaD's avatar
MischaD committed
163
            velocity_vectors.extend(velocities)
MischaD's avatar
MischaD committed
164

165
166
167
        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
168
169
        self._all_velocity_vectors = velocity_vectors
        return velocity_vectors
MischaD's avatar
MischaD committed
170

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

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

MischaD's avatar
MischaD committed
188
    def _fill_rhs(self):
MischaD's avatar
MischaD committed
189
190
191
192
193
194
195
196
197
198
199
200
        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
201
    def _svd(self):
MischaD's avatar
MischaD committed
202
203
204
        U, s, V = np.linalg.svd(self._lhs)
        rows = self._lhs.shape[0]
        cols = self._lhs.shape[1]
MischaD's avatar
MischaD committed
205
206
207
208

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

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

MischaD's avatar
MischaD committed
246
        coeffs = np.zeros((self._solution.shape[0] + 1, self._solution.shape[1] + 1))
MischaD's avatar
Bugfix    
MischaD committed
247
        coeffs[0, :] = coefficients[-1::-1]
MischaD's avatar
MischaD committed
248
249
250
251
        coeffs[1:, :-1] = self._solution[:, -1::-1]
        self._coefficients = coeffs
        return self._coefficients

MischaD's avatar
MischaD committed
252
    def _calculate_valid_interval(self):
MischaD's avatar
MischaD committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
        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

278
    def calculate_polynomials(self):
MischaD's avatar
MischaD committed
279
280
281
282
283
284
285
286
287
288
        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]

MischaD's avatar
MischaD committed
289
290
291
        self._calculate_velocity_vectors()
        self._fill_lhs()
        self._fill_rhs()
MischaD's avatar
MischaD committed
292

MischaD's avatar
MischaD committed
293
        self._svd()
MischaD's avatar
MischaD committed
294

MischaD's avatar
MischaD committed
295
        self._calculate_coefficients()
MischaD's avatar
Bugfix    
MischaD committed
296

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

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

MischaD's avatar
MischaD committed
303
        self._interval = self._calculate_valid_interval()
MischaD's avatar
MischaD committed
304
        if self._interval == sp.EmptySet():
MischaD's avatar
Bugfix    
MischaD committed
305
            return []
306
307
308
309
310
311
312
313
314
315
316
317
318
319
        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)

320
    def calculate_weights(self, c_s_sq=None, boundary=None):
321
        """
322
323
324
325
326
327
        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.

328
329
330
331
332
333
334
        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:
        """
335

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

348
349
350
351
        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
352
        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
353
354
355
356
        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
357
        return weights
MischaD's avatar
MischaD committed
358

359
360
361
362
363
364
    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
        """
MischaD's avatar
MischaD committed
365
366
        if not self.weights:
            self.calculate_weights()
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382

        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
383
384
385
386
387
388
389
390
391
392
393