lattice.py 12.5 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:
MischaD's avatar
MischaD committed
12

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

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

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

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

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

49
50
51
52
53
        # Final values
        self._discrete_velocities = []
        self._weights = []
        self._c_s_sq = None

MischaD's avatar
MischaD committed
54
55
    def __str__(self):
        string = "D{} - Order: {} - Shells: {}".format(
MischaD's avatar
MischaD committed
56
            self._dimension, self._order, str(self._shell_list)
MischaD's avatar
MischaD committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
        )
        if self._seed:
            string += " - Seed: {}".format(self._seed)
        return string

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

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

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

MischaD's avatar
MischaD committed
90
91
92
93
    @property
    def tensor_space_dimension(self):
        return np.array([len(x) for x in self._possible_tensors]).sum()

MischaD's avatar
MischaD committed
94
95
    @property
    def shells(self):
96
        return self._shells
MischaD's avatar
MischaD committed
97
98
99

    @property
    def velocities(self):
100
        return self._discrete_velocities
MischaD's avatar
MischaD committed
101

MischaD's avatar
MischaD committed
102
    def calculate_velocity_vectors(self):
MischaD's avatar
MischaD committed
103
        velocities_amount = 1
MischaD's avatar
MischaD committed
104
        velocity_vectors = []
MischaD's avatar
MischaD committed
105
        total_subshells = []
MischaD's avatar
MischaD committed
106

MischaD's avatar
MischaD committed
107
        group = get_group(self._dimension)
MischaD's avatar
MischaD committed
108
        for shell in range(len(self._shell_list)):
109
            velocities = velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
MischaD's avatar
MischaD committed
110
            if len(velocities) == 0:
MischaD's avatar
MischaD committed
111
112
                raise ImpossibleVelocityShellException("No velocity sets found for velocity shell {}.".format(shell))

113
            # If e.g. dimension = 2, then shell 25 in ambiguous. Both shell (3,4) and (5,0) are added
MischaD's avatar
MischaD committed
114
115
116
            subshells = get_subshells(velocities, group)
            total_subshells.append(len(subshells))
            if len(subshells) > 1:
MischaD's avatar
MischaD committed
117
118
119
120
121
122
                # 3 4, 0 5 can only be given together
                for i_subs, subshell in enumerate(subshells):
                    type = tuple(np.sort(abs(subshell[0])))
                velocities = subshells
            else:
                velocities = [velocities]
MischaD's avatar
MischaD committed
123
            velocity_vectors.extend(velocities)
MischaD's avatar
MischaD committed
124
125
126
            velocities_amount += len(velocities)

        logger.info("Velocities: {} Shells: {}".format(self._velocities, self._shells))
127
128
129
130
        self._velocities = velocities_amount
        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
131
132
        self._all_velocity_vectors = velocity_vectors
        return velocity_vectors
MischaD's avatar
MischaD committed
133

MischaD's avatar
MischaD committed
134
    def fill_lhs(self):
MischaD's avatar
MischaD committed
135
136
137
138
        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
139
                vector = get_random_vector(self._dimension)
MischaD's avatar
MischaD committed
140
141
                # vector = 2 * np.random.rand(self._dimension) - 1
                # vector = vector / np.linalg.norm(vector)
MischaD's avatar
MischaD committed
142
143

                rows = []
MischaD's avatar
MischaD committed
144
                for shell_velocities in self._all_velocity_vectors:
MischaD's avatar
MischaD committed
145
146
147
                    shell_sum = lattice_sum(vector, shell_velocities, rank)
                    rows.append(shell_sum)
                lhs.append(rows)
MischaD's avatar
MischaD committed
148
        self._lhs = np.array(lhs)
MischaD's avatar
MischaD committed
149
        return np.array(lhs)
MischaD's avatar
MischaD committed
150

MischaD's avatar
MischaD committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
    def fill_rhs(self):
        rhs = []
        cols = self._order // 2
        for k, rank in enumerate(np.arange(2, self._order + 2, 2)):
            for j in range(self._tensor_dimensions[k]):
                rows = []
                for i in range(cols):
                    c_s_power = 2 * i + 2
                    rows.append(1 if c_s_power == rank else 0)
                rhs.append(rows)
        self._rhs = np.array(rhs)
        return np.array(rhs)

    def svd(self):
        U, s, V = np.linalg.svd(self._lhs)
        rows = self._lhs.shape[0]
        cols = self._lhs.shape[1]
MischaD's avatar
MischaD committed
168
169
170
171

        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
172
        rhs = np.dot(np.transpose(U), self._rhs)
MischaD's avatar
MischaD committed
173
174
175
176
        if rank < rows:
            overlap = rows - rank
            if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
                # NO Solution
MischaD's avatar
MischaD committed
177
                raise NoSolutionException("Add shells or lower the order")  # change shells ?
MischaD's avatar
MischaD committed
178
179
            rows = rank
            s.resize(rows)
MischaD's avatar
MischaD committed
180
            rhs.resize((rows, rhs.shape[1]), refcheck=False)
MischaD's avatar
MischaD committed
181
182
183
184
185
186
187
188
189
190
191
192
193
            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
194
            raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
MischaD's avatar
MischaD committed
195
196
197
                                             "https://github.com/BDuenweg/Lattice-Boltzmann-weights")
        self._solution = np.dot(np.transpose(V), reduced_rhs)
        return self._solution
MischaD's avatar
MischaD committed
198

MischaD's avatar
MischaD committed
199
    def calculate_coefficients(self):
MischaD's avatar
MischaD committed
200
201
202
203
204
        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
205
206
207
            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
208

MischaD's avatar
MischaD committed
209
        coeffs = np.zeros((self._solution.shape[0] + 1, self._solution.shape[1] + 1))
MischaD's avatar
Bugfix    
MischaD committed
210
        coeffs[0, :] = coefficients[-1::-1]
MischaD's avatar
MischaD committed
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
        coeffs[1:, :-1] = self._solution[:, -1::-1]
        self._coefficients = coeffs
        return self._coefficients

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

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

241
    def calculate_polynomials(self):
MischaD's avatar
MischaD committed
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
        c_s_sq = sp.Symbol("c_s_sq")
        if self._order % 2 == 1:
            self._order -= 1
            logger.warning("Only odd order valid. Decrease by one")

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

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

        self.svd()

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

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

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

MischaD's avatar
MischaD committed
266
        self._interval = self.calculate_valid_interval()
MischaD's avatar
MischaD committed
267
268
        logger.info(self._interval.boundary)
        if self._interval == sp.EmptySet():
MischaD's avatar
Bugfix    
MischaD committed
269
            return []
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
        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)

    def calculate_weights(self, c_s_sq=None, boundary="inf"):
        """
        Calculate the weights for this lattice for a given speed of sound value c_s_sq.
        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:
        """
        if not self._weight_polynomials:
            self.calculate_polynomials()
        if not self._interval:  # Empty Interval
            return []
        if not c_s_sq:
            if boundary == "inf":
                c_s_sq = self._interval.inf
            else:
                c_s_sq = self._interval.sup
        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
307
        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
308
309
310
311
        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
312
        return weights
MischaD's avatar
MischaD committed
313

314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
    def velocity_set(self):
        """
        Returns the values that are important for the Lattice Boltzmann Model, i.e.
        a list of the velocity values their corresponding weights and the value for the speed of sound.
        :return: (c_s_sq, weights, velocities): c_s_sq
        """

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

        return self._c_s_sq, weights, velocities
MischaD's avatar
MischaD committed
336
337
338
339
340
341
342
343
344
345
346