lattice.py 11.9 KB
Newer Older
MischaD's avatar
MischaD committed
1
import numpy as np
MischaD's avatar
MischaD committed
2
3
import sympy as sp
from sympy import oo
MischaD's avatar
MischaD committed
4
from lbmweights.utils.mylog import logger
MischaD's avatar
MischaD committed
5
6
from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException, NoSolutionException, InfiniteSolutionsException, UninitializedAttributeException
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, fill_rhs, approximate_ratio, timer
MischaD's avatar
MischaD committed
7
8


MischaD's avatar
MischaD committed
9
class Lattice:
MischaD's avatar
MischaD committed
10

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

MischaD's avatar
MischaD committed
26
        # LES
MischaD's avatar
MischaD committed
27
28
        self._tensor_space_dimension = 0
        self._possible_tensors = []
MischaD's avatar
MischaD committed
29
        self._tensor_dimensions = []
MischaD's avatar
MischaD committed
30
31
        self._velocities = 0
        self._shells = 0
MischaD's avatar
MischaD committed
32

MischaD's avatar
MischaD committed
33
34
35
36
        # SVD
        self._amount_singular_values = 0
        self._singular_values = []

MischaD's avatar
MischaD committed
37
38
39
        self._weight_polynomials = []
        self._interval = None

MischaD's avatar
MischaD committed
40
41
    def __str__(self):
        string = "D{} - Order: {} - Shells: {}".format(
MischaD's avatar
MischaD committed
42
            self._dimension, self._order, str(self._shell_list)
MischaD's avatar
MischaD committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
        )
        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
75

MischaD's avatar
MischaD committed
76
77
78
79
    @property
    def tensor_space_dimension(self):
        return np.array([len(x) for x in self._possible_tensors]).sum()

MischaD's avatar
MischaD committed
80
81
82
83
84
85
86
87
    @property
    def shells(self):
        self._shells

    @property
    def velocities(self):
        self._velocities

MischaD's avatar
MischaD committed
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
    @staticmethod
    def velocities_for_shell(dimension, shell):
        """
        :param dimension: Lattice dimension
        :param shell: Squared velocity of the shell
        :return: list of velocities with squared length equal to shell
        """
        velocities = []
        linear_lattice_size = 2 * shell + 1
        full_lattice_size = linear_lattice_size ** dimension
        for site in range(full_lattice_size):
            tmp = site
            cur_vel_squared = 0
            cur_vec = []

            for dim in range(dimension):
                coordinate = tmp % linear_lattice_size
                tmp = (tmp - coordinate) // linear_lattice_size
                shifted_coordinate = coordinate - shell
                cur_vec.append(shifted_coordinate)
                cur_vel_squared += shifted_coordinate ** 2

            if cur_vel_squared == shell:
                velocities.append(np.array(cur_vec, dtype=int))

        return velocities

MischaD's avatar
MischaD committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
    def calculate_valid_interval(self):
        if not self._weight_polynomials:
            raise UninitializedAttributeException("Weight Polynomials not give. Calculate them first!")

        interval = sp.Interval(-oo, oo)
        root = []
        c_s_sq = sp.Symbol("c_s_sq")
        for equation in self._weight_polynomials:
            roots = list(sp.roots(equation, c_s_sq, filter="R").keys())
            roots.sort()
            if len(roots) == 0:
                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.subs(c_s_sq, mesh_point) < 0:
                    continue
                # positive weights
                if i == 0:
                    cur_interval = sp.Union(cur_interval, sp.Interval(-oo, roots[i]))
                elif i == (len(roots) - 1):
                    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 = interval
        return interval

MischaD's avatar
MischaD committed
144
145
146
147
148
149
150
151
152
153
154
155
156
    def calculate_weights(self):
        logger.debug("Entering Calculate Weights")
        logger.info(str(self))
        if self._order % 2 == 1:
            self._order -= 1
            logger.warning("Odd order for isotropy entered. This package only handles symmetric lattices, therefore"
                           "the equations for odd orders are automatically fulfilled. "
                           "Reduction of order to {}".format(self._order))

        # Do tensor dimension analysis for every order up to the desired
        for k in np.arange(2, self._order + 2, 2):
            possible_tensors = analyse_tensor_dimension(k)
            self._possible_tensors.append(possible_tensors)
MischaD's avatar
MischaD committed
157
        self._tensor_dimensions = [len(x) for x in self._possible_tensors]
MischaD's avatar
MischaD committed
158

MischaD's avatar
MischaD committed
159
        # FROM HERE  ...
MischaD's avatar
MischaD committed
160
161
        velocities_amount = 1
        grand_total_list = []
MischaD's avatar
MischaD committed
162
        total_subshells = []
MischaD's avatar
MischaD committed
163

MischaD's avatar
MischaD committed
164
        group = get_group(self._dimension) # TODO test output
MischaD's avatar
MischaD committed
165
        for shell in range(len(self._shell_list)):
MischaD's avatar
MischaD committed
166
            velocities = self.velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
MischaD's avatar
MischaD committed
167
            if len(velocities) == 0:
MischaD's avatar
MischaD committed
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
                raise ImpossibleVelocityShellException("No velocity sets found for velocity shell {}.".format(shell))

            subshells = get_subshells(velocities, group)
            total_subshells.append(len(subshells))
            if len(subshells) > 1:
                # TODO l.119
                raise ReducibleShellException("Shell {} is reducible")
            grand_total_list.extend(subshells)
            velocities_amount += len(velocities)

        self._velocities = velocities_amount
        self._shells = len(total_subshells) + 1
        logger.info("Velocities: {} Shells: {}".format(self._velocities, self._shells))

        logger.info("Non zero shells: ")
        for i, shell in enumerate(grand_total_list):
            number_of_velocities = len(shell)
            type = tuple(np.sort(abs(shell[0])))
MischaD's avatar
MischaD committed
186
            #logger.info("Shell number {} with c_i^2 = {} and {} velocities of type {}".format(i+1, (shell[0]**2).sum(), number_of_velocities, type))
MischaD's avatar
MischaD committed
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
        # TO HERE: Only Calculate GrandTotalList

        # Here: ...
        # SpacialDimension = self._dimension
        # MaxTensorRank = self._order
        # ListOfTensordimensions = [len(x) for x in self._possible_tensors]
        # GrandTotalList = grand_total_list == list shells which are list of arrays with the actual values
        # Arguments = <class 'dict'>: {'d': 2, 'm': 4, 'c': [1, 2, 4], 's': 44, 'y': 'test': 'quiet': 'write_latex':}

        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):
                vector = 2 * np.random.rand(self._dimension) - 1
                vector = vector / np.linalg.norm(vector)

                rows = []
                for shell_velocities in grand_total_list:
                    shell_sum = lattice_sum(vector, shell_velocities, rank)
                    rows.append(shell_sum)
                lhs.append(rows)
        A = np.array(lhs)
        b = fill_rhs(self._order, self._tensor_dimensions)

        # --------------------------------------------------------------------------------------------------
        logger.warning("A Overwritten for test purposes")
MischaD's avatar
MischaD committed
213
        """
MischaD's avatar
MischaD committed
214
215
216
        A = np.array([[2., 4., 8.],
                      [0.50244727, 1.99021093, 8.03915626],
                      [0.47686341, 2.09254637,  7.6298145]])
MischaD's avatar
MischaD committed
217
                      """
MischaD's avatar
MischaD committed
218
219
220
221
222
223
224
225
226
227
228
229
230
        U, s, V = np.linalg.svd(A)

        rows = A.shape[0]
        cols = A.shape[1]

        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)
        rhs = np.dot(np.transpose(U), b)
        if rank < rows:
            overlap = rows - rank
            if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
                # NO Solution
MischaD's avatar
MischaD committed
231
                raise NoSolutionException("Add shells or lower the order")  # change shells ?
MischaD's avatar
MischaD committed
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
            rows = rank
            s.resize(rows)
            rhs.resize((rows, rhs.shape[1]))
            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

        x = np.copy(reduced_rhs)
        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
249
            raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
MischaD's avatar
MischaD committed
250
251
252
253
                                         "https://github.com/BDuenweg/Lattice-Boltzmann-weights")

        # ----------------------------------------------- Unique solution
        solution = np.dot(np.transpose(V), reduced_rhs)
MischaD's avatar
MischaD committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
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
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
        solution_columns = self._order // 2
        assert solution_columns == solution.shape[1]

        coefficients = np.zeros((1 + solution_columns))
        coefficients[0] = 1
        for j in range(solution_columns):
            tmp = 0
            for i in range(len(grand_total_list)):
                tmp -= solution[i, j] * len(grand_total_list[i])
            coefficients[j+1] = tmp

        # -------------------------------------------------- Prettify for sympy
        weight_polynomials = []
        c_s_sq = sp.Symbol("c_s_sq")
        equation = 0
        for j in range(len(coefficients)):
            equation += approximate_ratio(coefficients[j]) * c_s_sq ** j
        weight_polynomials.append(equation)
        for i in range(len(grand_total_list)):
            equation = 0
            for j in range(solution_columns):
                equation += approximate_ratio(solution[i, j]) * c_s_sq ** (j + 1)
            weight_polynomials.append(equation)
        for i in weight_polynomials:
            print(i)

        # -------------------------------------------------- Calculate Roots
        """
        #interval = sp.Interval(-oo, oo)
        #for equation in weight_polynomials:
        #    _interval = sp.solveset(equation > 0, c_s_sq, sp.Reals)
        #    interval = interval.intersect(_interval)

        TOL = 1e-10
        tmp = []
        for i in range(len(coefficients)):
            tmp.append(coefficients[len(coefficients) - 1 - i])
        tmp_array = np.array(tmp)
        roots = [root for root in np.roots(tmp_array) if abs(root.imag) < TOL and abs(root.real) > TOL]

        rows = solution.shape[0]
        cols = solution.shape[1]
        for j in range(rows):
            tmp = []
            for i in range(cols):
                tmp.append(solution[j, cols - 1 - i])
            tmp_array = np.array(tmp)
            roots.append([root for root in np.roots(tmp_array) if abs(root.imag) < TOL and abs(root.real) > TOL])
        roots.sort()
        if not len(roots):
            return []
        """

        self._weight_polynomials = weight_polynomials
        self.calculate_valid_interval()
        logger.info(self._interval)