lattice.py 11.4 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
from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException, NoSolutionException, InfiniteSolutionsException, UninitializedAttributeException
MischaD's avatar
MischaD committed
6
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, fill_rhs, approximate_ratio, timer, get_random_vector
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
        self._svd_tolerance = svd_tolerance
MischaD's avatar
MischaD committed
24

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

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

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

MischaD's avatar
MischaD committed
39
40
    def __str__(self):
        string = "D{} - Order: {} - Shells: {}".format(
MischaD's avatar
MischaD committed
41
            self._dimension, self._order, str(self._shell_list)
MischaD's avatar
MischaD committed
42
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
        )
        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
74

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

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

    @property
    def velocities(self):
        self._velocities

MischaD's avatar
MischaD committed
87
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
    @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
114
    @timer
MischaD's avatar
MischaD committed
115
116
117
118
119
120
    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:
MischaD's avatar
Bugfix    
MischaD committed
121
122
            roots = equation.real_roots()
            if not roots:
MischaD's avatar
MischaD committed
123
124
125
126
127
                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
MischaD's avatar
Bugfix    
MischaD committed
128
                if equation.eval(mesh_point) < 0:
MischaD's avatar
MischaD committed
129
130
131
                    continue
                # positive weights
                if i == 0:
MischaD's avatar
MischaD committed
132
                    cur_interval = sp.Union(cur_interval, sp.Interval(-oo, roots[i+1]))
MischaD's avatar
MischaD committed
133
                elif i == (len(roots) - 2):
MischaD's avatar
MischaD committed
134
135
136
137
                    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)
MischaD's avatar
MischaD committed
138
        self._interval = sp.Complement(interval, sp.FiniteSet(0))
MischaD's avatar
MischaD committed
139
140
        return interval

MischaD's avatar
MischaD committed
141
    @timer
MischaD's avatar
MischaD committed
142
143
144
    def calculate_weights(self):
        logger.debug("Entering Calculate Weights")
        logger.info(str(self))
MischaD's avatar
MischaD committed
145
146
147
        import random
        random.seed(self._seed)
        print("My Seed: {}".format(self._seed))
MischaD's avatar
MischaD committed
148
149
150
151
152
153
154
155
156
157
        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
158
        self._tensor_dimensions = [len(x) for x in self._possible_tensors]
MischaD's avatar
MischaD committed
159

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

MischaD's avatar
MischaD committed
165
        group = get_group(self._dimension)
MischaD's avatar
MischaD committed
166
        for shell in range(len(self._shell_list)):
MischaD's avatar
MischaD committed
167
            velocities = self.velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
MischaD's avatar
MischaD committed
168
            if len(velocities) == 0:
MischaD's avatar
MischaD committed
169
170
171
172
173
                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:
MischaD's avatar
MischaD committed
174
175
176
177
178
179
180
                # 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]
            grand_total_list.extend(velocities)
MischaD's avatar
MischaD committed
181
182
            velocities_amount += len(velocities)

MischaD's avatar
MischaD committed
183
184
        self._group  = group
        self._grand_total_list = grand_total_list
MischaD's avatar
MischaD committed
185
186
187
188
189
190
191
        self._velocities = velocities_amount
        self._shells = len(total_subshells) + 1
        logger.info("Velocities: {} Shells: {}".format(self._velocities, self._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
192
            #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
193

MischaD's avatar
MischaD committed
194
195
196
197
198
199
200
201
202
203
204
205
206
        # 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):
MischaD's avatar
MischaD committed
207
208
209
                vector = get_random_vector(self._dimension)
                #vector = 2 * np.random.rand(self._dimension) - 1
                #vector = vector / np.linalg.norm(vector)
MischaD's avatar
MischaD committed
210
211
212
213
214
215

                rows = []
                for shell_velocities in grand_total_list:
                    shell_sum = lattice_sum(vector, shell_velocities, rank)
                    rows.append(shell_sum)
                lhs.append(rows)
MischaD's avatar
MischaD committed
216
        self._lhs = np.array(lhs)
MischaD's avatar
MischaD committed
217
218
        A = np.array(lhs)
        b = fill_rhs(self._order, self._tensor_dimensions)
MischaD's avatar
MischaD committed
219
        self._rhs = np.copy(b)
MischaD's avatar
MischaD committed
220
221
222

        # --------------------------------------------------------------------------------------------------

MischaD's avatar
MischaD committed
223
224
225
226
        U, s, V = np.linalg.svd(A)
        self._U = np.copy(U)
        self._s = np.copy(s)
        self._V = np.copy(V)
MischaD's avatar
MischaD committed
227
228
229
230
231
232
233
234
235
236
237
        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
238
                raise NoSolutionException("Add shells or lower the order")  # change shells ?
MischaD's avatar
MischaD committed
239
240
            rows = rank
            s.resize(rows)
MischaD's avatar
MischaD committed
241
            rhs.resize((rows, rhs.shape[1]), refcheck=False)
MischaD's avatar
MischaD committed
242
243
244
245
246
247
248
249
250
251
252
253
254
255
            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
256
            raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
MischaD's avatar
MischaD committed
257
258
259
260
                                         "https://github.com/BDuenweg/Lattice-Boltzmann-weights")

        # ----------------------------------------------- Unique solution
        solution = np.dot(np.transpose(V), reduced_rhs)
MischaD's avatar
MischaD committed
261
        self._solution = solution
MischaD's avatar
MischaD committed
262
263
264
265
266
267
268
269
270
271
272
        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

MischaD's avatar
Bugfix    
MischaD committed
273
        self._coefficients = np.array(coefficients)
MischaD's avatar
MischaD committed
274
275
276
        # -------------------------------------------------- Prettify for sympy
        weight_polynomials = []
        c_s_sq = sp.Symbol("c_s_sq")
MischaD's avatar
Bugfix    
MischaD committed
277
278
279
280
281
        coeffs = np.zeros((solution.shape[0] + 1, solution.shape[1] + 1))
        coeffs[0, :] = coefficients[-1::-1]
        coeffs[1:, :-1] = solution[:,-1::-1]

        apporx_coeffs = [[approximate_ratio(x) for x in pol_coffs] for pol_coffs in coeffs]
MischaD's avatar
MischaD committed
282

MischaD's avatar
Bugfix    
MischaD committed
283
        weight_polynomials = [sp.Poly([approximate_ratio(x) for x in pol_coffs], c_s_sq) for pol_coffs in coeffs]
MischaD's avatar
MischaD committed
284
285
286
287
        # -------------------------------------------------- Calculate Roots

        self._weight_polynomials = weight_polynomials
        self.calculate_valid_interval()
MischaD's avatar
MischaD committed
288
289
        logger.info(self._interval.boundary)
        if self._interval == sp.EmptySet():
MischaD's avatar
Bugfix    
MischaD committed
290
            return []
MischaD's avatar
MischaD committed
291
292
        weights = [sp.N(x.eval(self._interval.boundary.inf)) for x in weight_polynomials]
        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
MischaD's avatar
MischaD committed
293
        return weights
MischaD's avatar
MischaD committed
294
295
296
297
298
299
300
301
302
303
304
305