moments.py 26 KB
 Martin Bauer committed Apr 10, 2018 1 # -*- coding: utf-8 -*-  Martin Bauer committed Oct 29, 2018 2 r"""  Martin Bauer committed Feb 09, 2017 3 4 Module Overview ~~~~~~~~~~~~~~~  Martin Bauer committed Oct 04, 2016 5   Martin Bauer committed Feb 09, 2017 6 This module provides functions to  Martin Bauer committed Oct 04, 2016 7   Martin Bauer committed Feb 09, 2017 8 9 10 11 - generate moments according to certain patterns - compute moments of discrete probability distribution functions - create transformation rules into moment space - orthogonalize moment bases  Martin Bauer committed Oct 04, 2016 12 13   Martin Bauer committed Feb 09, 2017 14 15 Moment Representation ~~~~~~~~~~~~~~~~~~~~~  Martin Bauer committed Oct 04, 2016 16   Martin Bauer committed Feb 09, 2017 17 Moments can be represented in two ways:  Martin Bauer committed Oct 04, 2016 18   Martin Bauer committed Feb 09, 2017 19 20 21 - by an index :math:i,j: defined as :math:m_{ij} := \sum_{\mathbf{d} \in stencil} <\mathbf{d}, \mathbf{x}> f_i - or by a polynomial in the variables x,y and z. For example the polynomial :math:x^2 y^1 z^3 + x + 1 is describing the linear combination of moments: :math:m_{213} + m_{100} + m_{000}  Martin Bauer committed Oct 04, 2016 22   Martin Bauer committed Feb 09, 2017 23 24 The polynomial description is more powerful, since one polynomial can express a linear combination of single moments. All moment polynomials have to use MOMENT_SYMBOLS (which is a module variable) as degrees of freedom.  Martin Bauer committed Oct 04, 2016 25   Martin Bauer committed Feb 09, 2017 26 Example ::  Martin Bauer committed Oct 04, 2016 27   Martin Bauer committed Apr 10, 2018 28 29  >>> from lbmpy.moments import MOMENT_SYMBOLS >>> x, y, z = MOMENT_SYMBOLS  Martin Bauer committed Apr 10, 2018 30  >>> second_order_moment = x*y + y*z  Martin Bauer committed Oct 04, 2016 31   Martin Bauer committed Feb 09, 2017 32 33 34 35  Functions ~~~~~~~~~  Martin Bauer committed Feb 09, 2017 36 37 """ import itertools  Martin Bauer committed Sep 26, 2017 38 import math  Martin Bauer committed Apr 10, 2018 39 from collections import Counter, defaultdict  Martin Bauer committed Jul 11, 2019 40 41 42 from copy import copy from typing import Iterable, List, Optional, Sequence, Tuple, TypeVar  Martin Bauer committed Sep 26, 2017 43 import sympy as sp  Martin Bauer committed Oct 04, 2016 44   Martin Bauer committed Sep 26, 2017 45 from pystencils.cache import memorycache  Martin Bauer committed Apr 10, 2018 46 from pystencils.sympyextensions import remove_higher_order_terms  Martin Bauer committed Feb 09, 2017 47   Martin Bauer committed Apr 10, 2018 48 49 MOMENT_SYMBOLS = sp.symbols('x y z') T = TypeVar('T')  Martin Bauer committed Oct 04, 2016 50 51   Martin Bauer committed Feb 09, 2017 52 53 54 # ------------------------------ Discrete (Exponent Tuples) ------------------------------------------------------------  Martin Bauer committed Apr 10, 2018 55 def moment_multiplicity(exponent_tuple):  Martin Bauer committed Feb 09, 2017 56 57 58 59  """ Returns number of permutations of the given moment tuple Example:  Martin Bauer committed Apr 10, 2018 60  >>> moment_multiplicity((2,0,0))  Martin Bauer committed Feb 09, 2017 61  3  Martin Bauer committed Apr 10, 2018 62  >>> list(moment_permutations((2,0,0)))  Martin Bauer committed Feb 09, 2017 63 64  [(0, 0, 2), (0, 2, 0), (2, 0, 0)] """  Martin Bauer committed Apr 10, 2018 65 66  c = Counter(exponent_tuple) result = math.factorial(len(exponent_tuple))  Martin Bauer committed Oct 04, 2016 67 68 69 70 71  for d in c.values(): result //= math.factorial(d) return result  Martin Bauer committed Apr 10, 2018 72 def pick_representative_moments(moments):  Martin Bauer committed Feb 09, 2017 73  """Picks the representative i.e. of each permutation group only one is kept"""  Martin Bauer committed Apr 10, 2018 74  to_remove = []  Martin Bauer committed Feb 09, 2017 75  for m in moments:  Martin Bauer committed Apr 10, 2018 76 77 78  permutations = list(moment_permutations(m)) to_remove += permutations[1:] return set(moments) - set(to_remove)  Martin Bauer committed Feb 09, 2017 79 80   Martin Bauer committed Apr 10, 2018 81 def moment_permutations(exponent_tuple):  Martin Bauer committed Feb 09, 2017 82  """Returns all (unique) permutations of the given tuple"""  Martin Bauer committed Apr 10, 2018 83  return __unique_permutations(exponent_tuple)  Martin Bauer committed Oct 04, 2016 84 85   Martin Bauer committed Apr 10, 2018 86 def moments_of_order(order, dim=3, include_permutations=True):  Martin Bauer committed Oct 04, 2016 87  """All tuples of length 'dim' which sum equals 'order'"""  Martin Bauer committed Apr 10, 2018 88  for item in __fixed_sum_tuples(dim, order, ordered=not include_permutations):  Michael Kuron committed Feb 13, 2017 89 90 91  assert(len(item) == dim) assert(sum(item) == order) yield item  Martin Bauer committed Oct 04, 2016 92 93   Martin Bauer committed Apr 10, 2018 94 def moments_up_to_order(order, dim=3, include_permutations=True):  Martin Bauer committed Oct 04, 2016 95  """All tuples of length 'dim' which sum is smaller than 'order' """  Martin Bauer committed Apr 10, 2018 96 97  single_moment_iterators = [moments_of_order(o, dim, include_permutations) for o in range(order + 1)] return tuple(itertools.chain(*single_moment_iterators))  Martin Bauer committed Oct 04, 2016 98 99   Martin Bauer committed Apr 10, 2018 100 def moments_up_to_component_order(order, dim=3):  Martin Bauer committed Oct 04, 2016 101  """All tuples of length 'dim' where each entry is smaller or equal to 'order' """  Martin Bauer committed Feb 09, 2017 102  return tuple(itertools.product(*[range(order + 1)] * dim))  Martin Bauer committed Oct 04, 2016 103 104   Martin Bauer committed Apr 10, 2018 105 def extend_moments_with_permutations(exponent_tuples):  Martin Bauer committed Oct 04, 2016 106  """Returns all permutations of the given exponent tuples"""  Martin Bauer committed Apr 10, 2018 107 108 109 110  all_moments = [] for i in exponent_tuples: all_moments += list(moment_permutations(i)) return __unique(all_moments)  Martin Bauer committed Oct 04, 2016 111 112   Martin Bauer committed Feb 09, 2017 113 114 115 # ------------------------------ Representation Conversions ------------------------------------------------------------  Martin Bauer committed Apr 10, 2018 116 def exponent_to_polynomial_representation(exponent_tuple):  Martin Bauer committed Feb 09, 2017 117 118 119 120  """ Converts an exponent tuple to corresponding polynomial representation Example:  Martin Bauer committed Apr 10, 2018 121  >>> exponent_to_polynomial_representation( (2,1,3) )  Martin Bauer committed Feb 09, 2017 122 123 124  x**2*y*z**3 """ poly = 1  Martin Bauer committed Apr 10, 2018 125 126  for sym, tuple_entry in zip(MOMENT_SYMBOLS[:len(exponent_tuple)], exponent_tuple): poly *= sym ** tuple_entry  Martin Bauer committed Feb 09, 2017 127 128 129  return poly  Martin Bauer committed Apr 10, 2018 130 131 132 def exponents_to_polynomial_representations(sequence_of_exponent_tuples): """Applies :func:exponent_to_polynomial_representation to given sequence""" return tuple([exponent_to_polynomial_representation(t) for t in sequence_of_exponent_tuples])  Martin Bauer committed Feb 09, 2017 133 134   Martin Bauer committed Apr 10, 2018 135 def polynomial_to_exponent_representation(polynomial, dim=3):  Martin Bauer committed Feb 09, 2017 136 137 138 139  """ Converts a linear combination of moments in polynomial representation into exponent representation :returns list of tuples where the first element is the coefficient and the second element is the exponent tuple  Martin Bauer committed Oct 04, 2016 140   Martin Bauer committed Feb 09, 2017 141 142  Example: >>> x , y, z = MOMENT_SYMBOLS  Martin Bauer committed Apr 10, 2018 143  >>> set(polynomial_to_exponent_representation(1 + (42 * x**2 * y**2 * z) )) == {(42, (2, 2, 1)), (1, (0, 0, 0))}  Martin Bauer committed Feb 10, 2017 144  True  Martin Bauer committed Feb 09, 2017 145 146 147 148  """ assert dim <= 3 x, y, z = MOMENT_SYMBOLS polynomial = polynomial.expand()  Martin Bauer committed Apr 10, 2018 149  coeff_exp_tuple_representation = []  Martin Bauer committed Jan 31, 2018 150 151 152  summands = [polynomial] if polynomial.func != sp.Add else polynomial.args for expr in summands:  Martin Bauer committed Feb 09, 2017 153 154  if len(expr.atoms(sp.Symbol) - set(MOMENT_SYMBOLS)) > 0: raise ValueError("Invalid moment polynomial: " + str(expr))  Martin Bauer committed Jan 31, 2018 155  c, x_exp, y_exp, z_exp = sp.Wild('c'), sp.Wild('xexp'), sp.Wild('yexp'), sp.Wild('zc')  Martin Bauer committed Apr 10, 2018 156 157 158  match_res = expr.match(c * x**x_exp * y**y_exp * z**z_exp) assert match_res[x_exp].is_integer and match_res[y_exp].is_integer and match_res[z_exp].is_integer exp_tuple = (int(match_res[x_exp]), int(match_res[y_exp]), int(match_res[z_exp]),)  Martin Bauer committed Feb 09, 2017 159 160  if dim < 3: for i in range(dim, 3):  Martin Bauer committed Apr 10, 2018 161 162 163 164  assert exp_tuple[i] == 0, "Used symbols in polynomial are not representable in that dimension" exp_tuple = exp_tuple[:dim] coeff_exp_tuple_representation.append((match_res[c], exp_tuple)) return coeff_exp_tuple_representation  Martin Bauer committed Oct 04, 2016 165   Martin Bauer committed Feb 09, 2017 166   Martin Bauer committed Apr 10, 2018 167 def moment_sort_key(moment):  Martin Bauer committed Oct 21, 2017 168 169  """Sort key function for moments to sort them by (in decreasing priority) order, number of occuring symbols, length of string representation, string representation"""  Martin Bauer committed Apr 10, 2018 170 171  mom_str = str(moment) return get_order(moment), len(moment.atoms(sp.Symbol)), len(mom_str), mom_str  Martin Bauer committed Oct 21, 2017 172 173   Martin Bauer committed Apr 10, 2018 174 def sort_moments_into_groups_of_same_order(moments):  Martin Bauer committed Mar 05, 2017 175 176 177  """Returns a dictionary mapping the order (int) to a list of moments with that order.""" result = defaultdict(list) for i, moment in enumerate(moments):  Martin Bauer committed Apr 10, 2018 178  order = get_order(moment)  Martin Bauer committed Mar 05, 2017 179 180 181  result[order].append(moment) return result  Martin Bauer committed Oct 04, 2016 182 183 184 # -------------------- Common Function working with exponent tuples and polynomial moments -----------------------------  Martin Bauer committed Apr 10, 2018 185 def is_even(moment):  Martin Bauer committed Feb 09, 2017 186 187 188 189  """ A moment is considered even when under sign reversal nothing changes i.e. :math:m(-x,-y,-z) = m(x,y,z) For the exponent tuple representation that means that the exponent sum is even e.g.  Martin Bauer committed Feb 09, 2017 190  >>> x , y, z = MOMENT_SYMBOLS  Martin Bauer committed Apr 10, 2018 191  >>> is_even(x**2 * y**2)  Martin Bauer committed Feb 09, 2017 192  True  Martin Bauer committed Apr 10, 2018 193  >>> is_even(x**2 * y)  Martin Bauer committed Feb 09, 2017 194  False  Martin Bauer committed Apr 10, 2018 195  >>> is_even((1,0,0))  Martin Bauer committed Feb 09, 2017 196  False  Martin Bauer committed Apr 10, 2018 197  >>> is_even(1)  Martin Bauer committed Feb 09, 2017 198  True  Martin Bauer committed Oct 04, 2016 199 200 201 202  """ if type(moment) is tuple: return sum(moment) % 2 == 0 else:  Martin Bauer committed Feb 09, 2017 203  moment = sp.sympify(moment)  Martin Bauer committed Oct 04, 2016 204 205 206  opposite = moment for s in MOMENT_SYMBOLS: opposite = opposite.subs(s, -s)  Martin Bauer committed Feb 09, 2017 207  return sp.expand(moment - opposite) == 0  Martin Bauer committed Oct 04, 2016 208 209   Martin Bauer committed Apr 10, 2018 210 def get_moment_indices(moment_exponent_tuple):  Martin Bauer committed Jun 09, 2017 211 212 213  """Returns indices for a given exponent tuple: Example:  Martin Bauer committed Apr 10, 2018 214  >>> get_moment_indices((2,1,0))  Martin Bauer committed Jun 09, 2017 215  [0, 0, 1]  Martin Bauer committed Apr 10, 2018 216  >>> get_moment_indices((0,0,3))  Martin Bauer committed Jun 09, 2017 217 218 219  [2, 2, 2] """ result = []  Martin Bauer committed Apr 10, 2018 220  for i, element in enumerate(moment_exponent_tuple):  Martin Bauer committed Jun 09, 2017 221 222 223 224  result += [i] * element return result  Martin Bauer committed Apr 10, 2018 225 def get_exponent_tuple_from_indices(moment_indices, dim):  Martin Bauer committed Mar 22, 2018 226  result = [0] * dim  Martin Bauer committed Apr 10, 2018 227  for i in moment_indices:  Martin Bauer committed Mar 22, 2018 228 229 230 231  result[i] += 1 return tuple(result)  Martin Bauer committed Apr 10, 2018 232 def get_order(moment):  Martin Bauer committed Apr 10, 2018 233  """Computes polynomial order of given moment.  Martin Bauer committed Oct 10, 2016 234   Martin Bauer committed Feb 09, 2017 235  Examples:  Martin Bauer committed Feb 09, 2017 236  >>> x , y, z = MOMENT_SYMBOLS  Martin Bauer committed Apr 10, 2018 237  >>> get_order(x**2 * y + x)  Martin Bauer committed Feb 09, 2017 238  3  Martin Bauer committed Apr 10, 2018 239  >>> get_order(z**4 * x**2)  Martin Bauer committed Feb 09, 2017 240  6  Martin Bauer committed Apr 10, 2018 241  >>> get_order((2,1,0))  Martin Bauer committed Feb 09, 2017 242  3  Martin Bauer committed Feb 09, 2017 243  """  Martin Bauer committed Feb 09, 2017 244 245  if isinstance(moment, tuple): return sum(moment)  Martin Bauer committed Oct 10, 2016 246 247  if len(moment.atoms(sp.Symbol)) == 0: return 0  Martin Bauer committed Apr 10, 2018 248 249 250 251 252  leading_coefficient = sp.polys.polytools.LM(moment) symbols_in_leading_coefficient = leading_coefficient.atoms(sp.Symbol) return sum([sp.degree(leading_coefficient, gen=m) for m in symbols_in_leading_coefficient])  Martin Bauer committed Apr 10, 2018 253 def non_aliased_moment(moment_tuple: Sequence[int]) -> Tuple[int, ...]:  Martin Bauer committed Apr 10, 2018 254 255 256 257 258 259 260  """Takes a moment exponent tuple and returns the non-aliased version of it. For first neighborhood stencils, all moments with exponents 3, 5, 7... are equal to exponent 1, and exponents 4, 6, 8... are equal to exponent 2. This is because first neighborhood stencils only have values d ∈ {+1, 0, -1}. So for example d**5 is always the same as d**3 and d, and d**6 == d**4 == d**2 Example:  Martin Bauer committed Apr 10, 2018 261  >>> non_aliased_moment((5, 4, 2))  Martin Bauer committed Apr 10, 2018 262  (1, 2, 2)  Martin Bauer committed Apr 10, 2018 263  >>> non_aliased_moment((9, 1, 2))  Martin Bauer committed Apr 10, 2018 264 265 266 267 268 269 270 271 272 273  (1, 1, 2) """ moment = list(moment_tuple) result = [] for element in moment: if element > 2: result.append(2 - (element % 2)) else: result.append(element) return tuple(result)  Martin Bauer committed Oct 10, 2016 274 275   Martin Bauer committed Apr 10, 2018 276 def is_shear_moment(moment):  Martin Bauer committed Feb 09, 2017 277 278  """Shear moments in 3D are: x*y, x*z and y*z - in 2D its only x*y""" if type(moment) is tuple:  Martin Bauer committed Apr 10, 2018 279  moment = exponent_to_polynomial_representation(moment)  Martin Bauer committed Apr 10, 2018 280  return moment in is_shear_moment.shear_moments  Martin Bauer committed Apr 10, 2018 281 282   Martin Bauer committed Apr 10, 2018 283 is_shear_moment.shear_moments = set([c[0] * c[1] for c in itertools.combinations(MOMENT_SYMBOLS, 2)])  Martin Bauer committed Feb 09, 2017 284 285   Michael Kuron committed Feb 13, 2017 286 @memorycache(maxsize=512)  Martin Bauer committed Nov 13, 2019 287 def discrete_moment(func, moment, stencil, shift_velocity=None):  Martin Bauer committed Oct 29, 2018 288  r"""  Martin Bauer committed Feb 09, 2017 289 290 291 292 293 294 295  Computes discrete moment of given distribution function .. math :: \sum_{d \in stencil} p(d) f_i where :math:p(d) is the moment polynomial where :math:x, y, z have been replaced with the components of the stencil direction, and :math:f_i is the i'th entry in the passed function sequence  Martin Bauer committed Feb 09, 2017 296   Martin Bauer committed Sep 05, 2018 297 298 299 300  Args: func: list of distribution functions for each direction moment: can either be a exponent tuple, or a sympy polynomial expression stencil: sequence of directions  Martin Bauer committed Nov 13, 2019 301 302  shift_velocity: velocity vector u to compute central moments, the lattice velocity is replaced by (lattice_velocity - shift_velocity)  Martin Bauer committed Oct 04, 2016 303  """  Martin Bauer committed Apr 10, 2018 304  assert len(stencil) == len(func)  Martin Bauer committed Nov 13, 2019 305 306  if shift_velocity is None: shift_velocity = (0,) * len(stencil[0])  Martin Bauer committed Oct 04, 2016 307  res = 0  Martin Bauer committed Apr 10, 2018 308  for factor, e in zip(func, stencil):  Martin Bauer committed Oct 04, 2016 309  if type(moment) is tuple:  Martin Bauer committed Nov 17, 2019 310  for vel, shift, exponent in zip(e, shift_velocity, moment):  Martin Bauer committed Nov 13, 2019 311  factor *= (vel - shift) ** exponent  Martin Bauer committed Oct 04, 2016 312 313  res += factor else:  Martin Bauer committed Oct 05, 2016 314  weight = moment  Martin Bauer committed Nov 13, 2019 315 316  for variable, e_i, shift in zip(MOMENT_SYMBOLS, e, shift_velocity): weight = weight.subs(variable, e_i - shift)  Martin Bauer committed Oct 05, 2016 317  res += weight * factor  Martin Bauer committed Oct 04, 2016 318   Martin Bauer committed Feb 09, 2017 319  return res  Martin Bauer committed Oct 04, 2016 320 321   Martin Bauer committed Nov 13, 2019 322 def moment_matrix(moments, stencil, shift_velocity=None):  Martin Bauer committed Feb 09, 2017 323 324 325 326 327 328  """ Returns transformation matrix to moment space each row corresponds to a moment, each column to a direction of the stencil The entry i,j is the i'th moment polynomial evaluated at direction j """  Martin Bauer committed Nov 13, 2019 329 330  if shift_velocity is None: shift_velocity = (0,) * len(stencil[0])  Martin Bauer committed Oct 04, 2016 331 332 333 334  if type(moments[0]) is tuple: def generator(row, column): result = sp.Rational(1, 1)  Martin Bauer committed Nov 13, 2019 335 336  for exponent, stencil_entry, shift in zip(moments[row], stencil[column], shift_velocity): result *= (sp.sympify(stencil_entry) - shift) ** exponent  Martin Bauer committed Oct 04, 2016 337 338 339 340  return result else: def generator(row, column): evaluated = moments[row]  Martin Bauer committed Nov 13, 2019 341 342  for var, stencil_entry, shift in zip(MOMENT_SYMBOLS, stencil[column], shift_velocity): evaluated = evaluated.subs(var, stencil_entry - shift)  Martin Bauer committed Oct 04, 2016 343 344 345 346 347  return evaluated return sp.Matrix(len(moments), len(stencil), generator)  Martin Bauer committed Apr 10, 2018 348 def gram_schmidt(moments, stencil, weights=None):  Martin Bauer committed Oct 29, 2018 349  r"""  Martin Bauer committed Feb 09, 2017 350 351  Computes orthogonal set of moments using the method by Gram-Schmidt  Martin Bauer committed Oct 29, 2018 352 353 354 355 356 357 358 359  Args: moments: sequence of moments, either in tuple or polynomial form stencil: stencil as sequence of directions weights: optional weights, that define the scalar product which is used for normalization. Scalar product :math:< a,b > = \sum a_i b_i w_i with weights :math:w_i. Passing no weights sets all weights to 1. Returns: set of orthogonal moments in polynomial form  Martin Bauer committed Feb 09, 2017 360  """  Martin Bauer committed Oct 04, 2016 361 362 363 364 365 366 367  if weights is None: weights = sp.eye(len(stencil)) if type(weights) is list: assert len(weights) == len(stencil) weights = sp.diag(*weights) if type(moments[0]) is tuple:  Martin Bauer committed Apr 10, 2018 368  moments = list(exponents_to_polynomial_representations(moments))  Martin Bauer committed Oct 04, 2016 369  else:  Martin Bauer committed Feb 09, 2017 370  moments = list(copy(moments))  Martin Bauer committed Oct 04, 2016 371   Martin Bauer committed Apr 10, 2018 372 373 374 375 376  pdfs_to_moments = moment_matrix(moments, stencil).transpose() moment_matrix_columns = [pdfs_to_moments.col(i) for i in range(pdfs_to_moments.cols)] orthogonalized_vectors = [] for i in range(len(moment_matrix_columns)): current_element = moment_matrix_columns[i]  Martin Bauer committed Oct 04, 2016 377  for j in range(i):  Martin Bauer committed Apr 10, 2018 378 379 380  prev_element = orthogonalized_vectors[j] denominator = prev_element.dot(weights * prev_element) if denominator == 0:  Martin Bauer committed Oct 04, 2016 381  raise ValueError("Not an independent set of vectors given: "  Martin Bauer committed Oct 18, 2016 382  "vector %d is dependent on previous vectors" % (i,))  Martin Bauer committed Apr 10, 2018 383 384  overlap = current_element.dot(weights * prev_element) / denominator current_element -= overlap * prev_element  Martin Bauer committed Oct 04, 2016 385  moments[i] -= overlap * moments[j]  Martin Bauer committed Apr 10, 2018 386  orthogonalized_vectors.append(current_element)  Martin Bauer committed Feb 09, 2017 387 388  return moments  Martin Bauer committed Oct 04, 2016 389 390   Martin Bauer committed Apr 10, 2018 391 def get_default_moment_set_for_stencil(stencil):  Martin Bauer committed Feb 09, 2017 392 393 394  """ Returns a sequence of moments that are commonly used to construct a LBM equilibrium for the given stencil """  Martin Bauer committed Sep 19, 2018 395  from lbmpy.stencils import get_stencil  Martin Bauer committed May 05, 2019 396  from pystencils.stencil import have_same_entries  Martin Bauer committed Feb 09, 2017 397   Martin Bauer committed Apr 10, 2018 398  to_poly = exponents_to_polynomial_representations  Martin Bauer committed Feb 09, 2017 399   Martin Bauer committed May 05, 2019 400  if have_same_entries(stencil, get_stencil("D2Q9")):  Martin Bauer committed Apr 10, 2018 401  return sorted(to_poly(moments_up_to_component_order(2, dim=2)), key=moment_sort_key)  Martin Bauer committed Feb 09, 2017 402   Martin Bauer committed Apr 10, 2018 403  all27_moments = moments_up_to_component_order(2, dim=3)  Martin Bauer committed May 05, 2019 404  if have_same_entries(stencil, get_stencil("D3Q27")):  Martin Bauer committed Apr 10, 2018 405  return to_poly(all27_moments)  Martin Bauer committed May 05, 2019 406  if have_same_entries(stencil, get_stencil("D3Q19")):  Martin Bauer committed Apr 10, 2018 407 408 409  non_matched_moments = [(1, 2, 2), (1, 1, 2), (2, 2, 2), (1, 1, 1)] moments19 = set(all27_moments) - set(extend_moments_with_permutations(non_matched_moments)) return sorted(to_poly(moments19), key=moment_sort_key)  Martin Bauer committed May 05, 2019 410  if have_same_entries(stencil, get_stencil("D3Q15")):  Martin Bauer committed Feb 09, 2017 411  x, y, z = MOMENT_SYMBOLS  Martin Bauer committed Apr 10, 2018 412 413 414 415 416 417 418 419  non_matched_moments = [(1, 2, 0), (2, 2, 0), (1, 1, 2), (1, 2, 2), (2, 2, 2)] additional_moments = (6 * (x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2), 3 * (x * (y ** 2 + z ** 2)), 3 * (y * (x ** 2 + z ** 2)), 3 * (z * (x ** 2 + y ** 2)), ) to_remove = set(extend_moments_with_permutations(non_matched_moments)) return sorted(to_poly(set(all27_moments) - to_remove) + additional_moments, key=moment_sort_key)  Martin Bauer committed Feb 09, 2017 420 421 422 423  raise NotImplementedError("No default moment system available for this stencil - define matched moments yourself")  Martin Bauer committed Apr 10, 2018 424 425 426 def extract_monomials(sequence_of_polynomials, dim=3): """ Returns a set of exponent tuples of all monomials contained in the given set of polynomials  Martin Bauer committed Aug 09, 2019 427 428 429 430  Args: sequence_of_polynomials: sequence of polynomials in the MOMENT_SYMBOLS dim: length of returned exponent tuples  Martin Bauer committed Apr 10, 2018 431 432 433 434 435 436 437 438 439  >>> x, y, z = MOMENT_SYMBOLS >>> extract_monomials([x**2 + y**2 + y, y + y**2]) {(0, 2, 0), (0, 1, 0), (2, 0, 0)} >>> extract_monomials([x**2 + y**2 + y, y + y**2], dim=2) {(0, 1), (2, 0), (0, 2)} """ monomials = set() for polynomial in sequence_of_polynomials:  Martin Bauer committed Apr 10, 2018 440 441  for factor, exponent_tuple in polynomial_to_exponent_representation(polynomial): monomials.add(exponent_tuple[:dim])  Martin Bauer committed Apr 10, 2018 442 443 444 445 446 447  return monomials def monomial_to_polynomial_transformation_matrix(monomials, polynomials): """ Returns a transformation matrix from a monomial to a polynomial representation  Martin Bauer committed Aug 09, 2019 448 449 450 451  Args: monomials: sequence of exponent tuples polynomials: sequence of polynomials in the MOMENT_SYMBOLS  Martin Bauer committed Apr 10, 2018 452 453 454 455 456 457 458 459 460 461 462 463 464  >>> x, y, z = MOMENT_SYMBOLS >>> polys = [7 * x**2 + 3 * x + 2 * y **2, \ 9 * x**2 - 5 * x] >>> mons = list(extract_monomials(polys, dim=2)) >>> monomial_to_polynomial_transformation_matrix(mons, polys) Matrix([ [7, 3, 2], [9, -5, 0]]) """ dim = len(monomials[0]) result = sp.zeros(len(polynomials), len(monomials))  Martin Bauer committed Apr 10, 2018 465  for polynomial_idx, polynomial in enumerate(polynomials):  Martin Bauer committed Apr 10, 2018 466 467  for factor, exponent_tuple in polynomial_to_exponent_representation(polynomial): exponent_tuple = exponent_tuple[:dim]  Martin Bauer committed Apr 10, 2018 468  result[polynomial_idx, monomials.index(exponent_tuple)] = factor  Martin Bauer committed Apr 10, 2018 469 470 471  return result  Martin Bauer committed Feb 09, 2017 472 473 474 # ---------------------------------- Visualization ---------------------------------------------------------------------  Martin Bauer committed Apr 10, 2018 475 def moment_equality_table(stencil, discrete_equilibrium=None, continuous_equilibrium=None,  Martin Bauer committed Apr 10, 2018 476  max_order=4, truncate_order=3):  Martin Bauer committed Oct 04, 2016 477 478 479  """ Creates a table showing which moments of a discrete stencil/equilibrium coincide with the corresponding continuous moments  Martin Bauer committed Feb 09, 2017 480   Martin Bauer committed Jun 25, 2018 481 482 483 484 485 486 487 488 489 490  Args: stencil: list of stencil velocities discrete_equilibrium: list of sympy expr to compute discrete equilibrium for each direction, if left to default the standard discrete Maxwellian equilibrium is used continuous_equilibrium: continuous equilibrium, if left to default, the continuous Maxwellian is used max_order: compare moments up to this order (number of rows in table) truncate_order: moments are considered equal if they match up to this order Returns: Object to display in an Jupyter notebook  Martin Bauer committed Oct 04, 2016 491 492  """ import ipy_table  Martin Bauer committed Apr 10, 2018 493  from lbmpy.continuous_distribution_measures import continuous_moment  Martin Bauer committed Jun 09, 2017 494   Martin Bauer committed Oct 04, 2016 495  dim = len(stencil[0])  Martin Bauer committed Apr 27, 2018 496  u = sp.symbols("u_:{dim}".format(dim=dim))  Martin Bauer committed Apr 10, 2018 497  if discrete_equilibrium is None:  Martin Bauer committed Apr 10, 2018 498 499 500  from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium discrete_equilibrium = discrete_maxwellian_equilibrium(stencil, c_s_sq=sp.Rational(1, 3), compressible=True, u=u, order=truncate_order)  Martin Bauer committed Apr 10, 2018 501  if continuous_equilibrium is None:  Martin Bauer committed Apr 10, 2018 502 503  from lbmpy.maxwellian_equilibrium import continuous_maxwellian_equilibrium continuous_equilibrium = continuous_maxwellian_equilibrium(dim=dim, u=u, c_s_sq=sp.Rational(1, 3))  Martin Bauer committed Feb 09, 2017 504   Martin Bauer committed Oct 04, 2016 505  table = []  Martin Bauer committed Apr 10, 2018 506 507  matched_moments = 0 non_matched_moments = 0  Martin Bauer committed Oct 04, 2016 508   Martin Bauer committed Apr 10, 2018 509  moments_list = [list(moments_of_order(o, dim, include_permutations=False)) for o in range(max_order + 1)]  Martin Bauer committed Oct 04, 2016 510 511  colors = dict()  Martin Bauer committed Apr 10, 2018 512  nr_of_columns = max([len(v) for v in moments_list]) + 1  Martin Bauer committed Oct 04, 2016 513   Martin Bauer committed Apr 10, 2018 514 515 516  header_row = [' '] * nr_of_columns header_row[0] = 'order' table.append(header_row)  Martin Bauer committed Oct 04, 2016 517   Martin Bauer committed Apr 10, 2018 518 519  for order, moments in enumerate(moments_list): row = [' '] * nr_of_columns  Martin Bauer committed Oct 04, 2016 520  row[0] = '%d' % (order,)  Martin Bauer committed Apr 10, 2018 521  for moment, col_idx in zip(moments, range(1, len(row))):  Martin Bauer committed Apr 10, 2018 522 523 524  multiplicity = moment_multiplicity(moment) dm = discrete_moment(discrete_equilibrium, moment, stencil) cm = continuous_moment(continuous_equilibrium, moment, symbols=sp.symbols("v_0 v_1 v_2")[:dim])  Martin Bauer committed Oct 04, 2016 525  difference = sp.simplify(dm - cm)  Martin Bauer committed Apr 10, 2018 526 527  if truncate_order: difference = sp.simplify(remove_higher_order_terms(difference, symbols=u, order=truncate_order))  Martin Bauer committed Oct 04, 2016 528  if difference != 0:  Martin Bauer committed Apr 10, 2018 529  colors[(order + 1, col_idx)] = 'Orange'  Martin Bauer committed Apr 10, 2018 530  non_matched_moments += multiplicity  Martin Bauer committed Oct 04, 2016 531  else:  Martin Bauer committed Apr 10, 2018 532  colors[(order + 1, col_idx)] = 'lightGreen'  Martin Bauer committed Apr 10, 2018 533  matched_moments += multiplicity  Martin Bauer committed Oct 04, 2016 534   Martin Bauer committed Apr 10, 2018 535  row[col_idx] = '%s x %d' % (moment, moment_multiplicity(moment))  Martin Bauer committed Oct 04, 2016 536 537 538  table.append(row)  Martin Bauer committed Apr 10, 2018 539  table_display = ipy_table.make_table(table)  Martin Bauer committed Oct 04, 2016 540  ipy_table.set_row_style(0, color='#ddd')  Martin Bauer committed Apr 10, 2018 541 542  for cell_idx, color in colors.items(): ipy_table.set_cell_style(cell_idx[0], cell_idx[1], color=color)  Martin Bauer committed Oct 04, 2016 543 544  print("Matched moments %d - non matched moments %d - total %d" %  Martin Bauer committed Apr 10, 2018 545  (matched_moments, non_matched_moments, matched_moments + non_matched_moments))  Martin Bauer committed Oct 04, 2016 546   Martin Bauer committed Apr 10, 2018 547  return table_display  Martin Bauer committed Feb 09, 2017 548 549   Martin Bauer committed Apr 10, 2018 550 def moment_equality_table_by_stencil(name_to_stencil_dict, moments, truncate_order=3):  Martin Bauer committed Feb 09, 2017 551 552 553 554  """ Creates a table for display in IPython notebooks that shows which moments agree between continuous and discrete equilibrium, group by stencils  Martin Bauer committed Jun 25, 2018 555 556 557 558 559  Args: name_to_stencil_dict: dict from stencil name to stencil moments: sequence of moments to compare - assumes that permutations have similar properties so just one representative is shown labeled with its multiplicity truncate_order: compare up to this order  Martin Bauer committed Feb 09, 2017 560 561  """ import ipy_table  Martin Bauer committed Apr 10, 2018 562 563 564  from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium from lbmpy.maxwellian_equilibrium import continuous_maxwellian_equilibrium from lbmpy.continuous_distribution_measures import continuous_moment  Martin Bauer committed Feb 09, 2017 565   Martin Bauer committed Apr 10, 2018 566  stencil_names = []  Martin Bauer committed Feb 09, 2017 567  stencils = []  Martin Bauer committed Apr 10, 2018 568  for key, value in name_to_stencil_dict.items():  Martin Bauer committed Apr 10, 2018 569  stencil_names.append(key)  Martin Bauer committed Feb 09, 2017 570 571  stencils.append(value)  Martin Bauer committed Apr 10, 2018 572  moments = list(pick_representative_moments(moments))  Martin Bauer committed Feb 09, 2017 573 574  colors = {}  Martin Bauer committed Apr 10, 2018 575  for stencil_idx, stencil in enumerate(stencils):  Martin Bauer committed Feb 09, 2017 576  dim = len(stencil[0])  Martin Bauer committed Apr 27, 2018 577  u = sp.symbols("u_:{dim}".format(dim=dim))  Martin Bauer committed Apr 10, 2018 578 579 580  discrete_equilibrium = discrete_maxwellian_equilibrium(stencil, c_s_sq=sp.Rational(1, 3), compressible=True, u=u, order=truncate_order) continuous_equilibrium = continuous_maxwellian_equilibrium(dim=dim, u=u, c_s_sq=sp.Rational(1, 3))  Martin Bauer committed Feb 09, 2017 581   Martin Bauer committed Apr 10, 2018 582  for moment_idx, moment in enumerate(moments):  Martin Bauer committed Feb 09, 2017 583  moment = moment[:dim]  Martin Bauer committed Apr 10, 2018 584 585  dm = discrete_moment(discrete_equilibrium, moment, stencil) cm = continuous_moment(continuous_equilibrium, moment, symbols=sp.symbols("v_0 v_1 v_2")[:dim])  Martin Bauer committed Feb 09, 2017 586  difference = sp.simplify(dm - cm)  Martin Bauer committed Apr 10, 2018 587 588  if truncate_order: difference = sp.simplify(remove_higher_order_terms(difference, symbols=u, order=truncate_order))  Martin Bauer committed Apr 10, 2018 589  colors[(moment_idx + 1, stencil_idx + 2)] = 'Orange' if difference != 0 else 'lightGreen'  Martin Bauer committed Feb 09, 2017 590 591  table = []  Martin Bauer committed Apr 10, 2018 592 593  header_row = [' ', '#'] + stencil_names table.append(header_row)  Martin Bauer committed Feb 09, 2017 594  for moment in moments:  Martin Bauer committed Apr 10, 2018 595  row = [str(moment), str(moment_multiplicity(moment))] + [' '] * len(stencils)  Martin Bauer committed Feb 09, 2017 596 597  table.append(row)  Martin Bauer committed Apr 10, 2018 598  table_display = ipy_table.make_table(table)  Martin Bauer committed Feb 09, 2017 599  ipy_table.set_row_style(0, color='#ddd')  Martin Bauer committed Apr 10, 2018 600 601  for cell_idx, color in colors.items(): ipy_table.set_cell_style(cell_idx[0], cell_idx[1], color=color)  Martin Bauer committed Feb 09, 2017 602   Martin Bauer committed Apr 10, 2018 603  return table_display  Martin Bauer committed Feb 09, 2017 604 605   Martin Bauer committed Apr 10, 2018 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 # --------------------------------------- Internal Functions ----------------------------------------------------------- def __unique(seq: Sequence[T]) -> List[T]: """Removes duplicates from a sequence in an order preserving way. Example: >>> __unique((1, 2, 3, 1, 4, 6, 3)) [1, 2, 3, 4, 6] """ seen = {} result = [] for item in seq: if item in seen: continue seen[item] = 1 result.append(item) return result def __unique_permutations(elements: Sequence[T]) -> Iterable[T]: """Generates all unique permutations of the passed sequence. Example: >>> list(__unique_permutations([1, 1, 2])) [(1, 1, 2), (1, 2, 1), (2, 1, 1)] """ if len(elements) == 1: yield (elements[0],) else: unique_elements = set(elements) for first_element in unique_elements: remaining_elements = list(elements) remaining_elements.remove(first_element) for sub_permutation in __unique_permutations(remaining_elements): yield (first_element,) + sub_permutation def __fixed_sum_tuples(tuple_length: int, tuple_sum: int, allowed_values: Optional[Sequence[int]] = None, ordered: bool = False) -> Iterable[Tuple[int, ...]]: """Generates all possible tuples of positive integers with a fixed sum of all entries. Args: tuple_length: length of the returned tuples tuple_sum: summing over the entries of a yielded tuple should give this number allowed_values: optional sequence of positive integers that are considered as tuple entries zero has to be in the set of allowed values if None, all possible positive integers are allowed ordered: if True, only tuples are returned where the entries are descending, i.e. where the entries are ordered Examples: Generate all 2-tuples where the sum of entries is 3 >>> list(__fixed_sum_tuples(tuple_length=2, tuple_sum=3)) [(0, 3), (1, 2), (2, 1), (3, 0)] Same with ordered tuples only >>> list(__fixed_sum_tuples(tuple_length=2, tuple_sum=3, ordered=True)) [(2, 1), (3, 0)] Restricting the allowed values, note that zero has to be in the allowed values! >>> list(__fixed_sum_tuples(tuple_length=3, tuple_sum=4, allowed_values=(0, 1, 3))) [(0, 1, 3), (0, 3, 1), (1, 0, 3), (1, 3, 0), (3, 0, 1), (3, 1, 0)] """ if not allowed_values: allowed_values = set(range(0, tuple_sum + 1)) assert 0 in allowed_values and all(i >= 0 for i in allowed_values) def recursive_helper(current_list, position, total_sum): new_position = position + 1 if new_position < len(current_list): for i in allowed_values: current_list[position] = i new_sum = total_sum - i if new_sum < 0: continue for item in recursive_helper(current_list, new_position, new_sum): yield item else: if total_sum in allowed_values: current_list[-1] = total_sum if not ordered: yield tuple(current_list) if ordered and current_list == sorted(current_list, reverse=True): yield tuple(current_list)  Martin Bauer committed Apr 10, 2018 693  return recursive_helper([0] * tuple_length, 0, tuple_sum)