Lattice Symmetries

.codeclimate.yml

@@ -0,0 +1,36 @@
+
+engines:
+  # ... CONFIG CONTENT ...
+  pep8:
+    enabled: true
+    # ... CONFIG CONTENT ...
+    checks:
+        E731:  # do not assign a lambda expression, use a def
+            enabled: false
+        E402:  # module level import not at top of file
+            enabled: false
+        E501: # maximum line length
+            enabled: false
+
+checks:
+    argument-count:
+        enabled: false
+    complex-logic:
+        enabled: false
+    file-lines:
+        enabled: false
+    method-complexity:
+        enabled: false
+    method-count:
+        enabled: false
+    method-lines:
+        enabled: false
+    nested-control-flow:
+        enabled: false
+    return-statements:
+        enabled: false
+    similar-code:
+        enabled: false
+    identical-code:
+        enabled: true
+

.pep8speaks.yml

@@ -0,0 +1,10 @@
+
+scanner:
+    diff_only: True  # If False, the entire file touched by the Pull Request is scanned for errors. If True, only the diff is scanned.
+    linter: pycodestyle  # Other option is flake8
+
+pycodestyle:  # Same as scanner.linter value. Other option is flake8
+    max-line-length: 121  # Default is 79 in PEP 8
+    ignore:  # Errors and warnings to ignore
+        - E731  # do not assign a lambda expression, use a def
+        - E402  # module level import not at top of file

.travis.yml

@@ -35,7 +35,7 @@ install:
 # Command to run tests, e.g. python setup.py test
 script:
   # run unit tests
-  - py.test
+  - py.test --runslow
   # run integration tests
   - lettuce --no-cuda convergence
   - lettuce --no-cuda benchmark

lettuce/__init__.py

@@ -10,6 +10,7 @@
 __version__ = get_versions()['version']
 del get_versions
 
+
 from lettuce.util import *
 from lettuce.unit import *
 from lettuce.lattices import *
@@ -25,5 +26,7 @@
 from lettuce.simulation import *
 from lettuce.force import *
 from lettuce.observables import *
+from lettuce.symmetry import *
+from lettuce.neural import *
 
 from lettuce.flows import *

lettuce/boundary.py

@@ -117,9 +117,9 @@ def __call__(self, f):
         u = self.lattice.u(f)
         u_w = u[[slice(None)] + self.index] + 0.5 * (u[[slice(None)] + self.index] - u[[slice(None)] + self.neighbor])
         f[[np.array(self.lattice.stencil.opposite)[self.velocities]] + self.index] = (
-                - f[[self.velocities] + self.index] + self.w * self.lattice.rho(f)[[slice(None)] + self.index] *
-                (2 + torch.einsum(self.dims, self.lattice.e[self.velocities], u_w) ** 2 / self.lattice.cs ** 4
-                 - (torch.norm(u_w, dim=0) / self.lattice.cs) ** 2)
+            - f[[self.velocities] + self.index] + self.w * self.lattice.rho(f)[[slice(None)] + self.index] *
+            (2 + torch.einsum(self.dims, self.lattice.e[self.velocities], u_w) ** 2 / self.lattice.cs ** 4
+             - (torch.norm(u_w, dim=0) / self.lattice.cs) ** 2)
         )
         return f
 

lettuce/collision.py

@@ -3,17 +3,26 @@
 """
 
 import torch
+import numpy as np
 
 from lettuce.equilibrium import QuadraticEquilibrium
-from lettuce.util import LettuceException
+from lettuce.moments import DEFAULT_TRANSFORM
+from lettuce.util import LettuceCollisionNotDefined
+from lettuce.stencils import D2Q9, D3Q27
 
 __all__ = [
-    "BGKCollision", "KBCCollision2D", "KBCCollision3D", "MRTCollision", "RegularizedCollision",
-    "SmagorinskyCollision", "TRTCollision", "BGKInitialization"
+    "Collision",
+    "BGKCollision", "KBCCollision2D", "KBCCollision3D", "MRTCollision",
+    "RegularizedCollision", "SmagorinskyCollision", "TRTCollision", "BGKInitialization",
 ]
 
 
-class BGKCollision:
+class Collision:
+    def __call__(self, f):
+        return NotImplemented
+
+
+class BGKCollision(Collision):
     def __init__(self, lattice, tau, force=None):
         self.force = force
         self.lattice = lattice
@@ -28,17 +37,27 @@ def __call__(self, f):
         return f - 1.0 / self.tau * (f - feq) + Si
 
 
-class MRTCollision:
+class MRTCollision(Collision):
     """Multiple relaxation time collision operator
 
     This is an MRT operator in the most general sense of the word.
     The transform does not have to be linear and can, e.g., be any moment or cumulant transform.
     """
 
-    def __init__(self, lattice, transform, relaxation_parameters):
+    def __init__(self, lattice, relaxation_parameters, transform=None):
         self.lattice = lattice
-        self.transform = transform
-        self.relaxation_parameters = lattice.convert_to_tensor(relaxation_parameters)
+        if transform is None:
+            try:
+                self.transform = DEFAULT_TRANSFORM[lattice.stencil](lattice)
+            except KeyError:
+                raise LettuceCollisionNotDefined("No entry for stencil {lattice.stencil} in moments.DEFAULT_TRANSFORM")
+        else:
+            self.transform = transform
+        if isinstance(relaxation_parameters, float):
+            tau = relaxation_parameters
+            self.relaxation_parameters = lattice.convert_to_tensor(tau * np.ones(lattice.stencil.Q()))
+        else:
+            self.relaxation_parameters = lattice.convert_to_tensor(relaxation_parameters)
 
     def __call__(self, f):
         m = self.transform.transform(f)
@@ -48,7 +67,7 @@ def __call__(self, f):
         return f
 
 
-class TRTCollision:
+class TRTCollision(Collision):
     """Two relaxation time collision model - standard implementation (cf. Krüger 2017)
         """
 
@@ -62,14 +81,14 @@ def __call__(self, f):
         u = self.lattice.u(f)
         feq = self.lattice.equilibrium(rho, u)
         f_diff_neq = ((f + f[self.lattice.stencil.opposite]) - (feq + feq[self.lattice.stencil.opposite])) / (
-                2.0 * self.tau_plus)
+            2.0 * self.tau_plus)
         f_diff_neq += ((f - f[self.lattice.stencil.opposite]) - (feq - feq[self.lattice.stencil.opposite])) / (
-                2.0 * self.tau_minus)
+            2.0 * self.tau_minus)
         f = f - f_diff_neq
         return f
 
 
-class RegularizedCollision:
+class RegularizedCollision(Collision):
     """Regularized LBM according to Jonas Latt and Bastien Chopard (2006)"""
 
     def __init__(self, lattice, tau):
@@ -100,12 +119,13 @@ def __call__(self, f):
         return f
 
 
-class KBCCollision2D:
+class KBCCollision2D(Collision):
     """Entropic multi-relaxation time model according to Karlin et al. in two dimensions"""
 
     def __init__(self, lattice, tau):
         self.lattice = lattice
-        assert lattice.Q == 9, LettuceException("KBC2D only realized for D2Q9")
+        if not lattice.stencil == D2Q9:
+            raise LettuceCollisionNotDefined("This implementation only works for the D2Q9 stencil.")
         self.tau = tau
         self.beta = 1. / (2 * tau)
 
@@ -178,12 +198,13 @@ def __call__(self, f):
         return f
 
 
-class KBCCollision3D:
+class KBCCollision3D(Collision):
     """Entropic multi-relaxation time-relaxation time model according to Karlin et al. in three dimensions"""
 
     def __init__(self, lattice, tau):
         self.lattice = lattice
-        assert lattice.Q == 27, LettuceException("KBC only realized for D3Q27")
+        if not lattice.stencil == D3Q27:
+            raise LettuceCollisionNotDefined("This implementation only works for the D3Q27 stencil.")
         self.tau = tau
         self.beta = 1. / (2 * tau)
 
@@ -269,7 +290,7 @@ def __call__(self, f):
         return f
 
 
-class SmagorinskyCollision:
+class SmagorinskyCollision(Collision):
     """Smagorinsky large eddy simulation (LES) collision model with BGK operator."""
 
     def __init__(self, lattice, tau, smagorinsky_constant=0.17, force=None):
@@ -324,3 +345,4 @@ def __call__(self, f):
         mnew[self.momentum_indices] = rho * self.u
         f = self.moment_transformation.inverse_transform(mnew)
         return f
+

lettuce/lattices.py

@@ -8,6 +8,7 @@
 Its stencil is still accessible trough Lattice.stencil.
 """
 
+import copy
 import warnings
 import numpy as np
 import torch
@@ -102,15 +103,16 @@ def mv(self, m, v):
     def einsum(self, equation, fields, **kwargs):
         """Einstein summation on local fields."""
         input, output = equation.split("->")
+        out = copy.copy(output)
         inputs = input.split(",")
+        xyz = "xyz"[:self.D]
         for i, inp in enumerate(inputs):
-            if len(inp) == len(fields[i].shape):
+            if "..." in inp:
+                raise LettuceException("... not allowed in lattice.einsum")
+            elif len(inp) == len(fields[i].shape):
                 pass
-            elif len(inp) == len(fields[i].shape) - self.D:
-                inputs[i] += "..."
-                if not output.endswith("..."):
-                    output += "..."
             else:
-                raise LettuceException("Bad dimension.")
-        equation = ",".join(inputs) + "->" + output
+                inputs[i] = f"...{inp}{xyz}"
+                out = f"...{output}{xyz}"
+        equation = ",".join(inputs) + "->" + out
         return torch.einsum(equation, fields, **kwargs)

lettuce/moments.py

@@ -11,7 +11,7 @@
 
 __all__ = [
     "moment_tensor", "get_default_moment_transform", "Moments", "Transform", "D1Q3Transform",
-    "D2Q9Lallemand", "D2Q9Dellar", "D3Q27Hermite"
+    "D2Q9Lallemand", "D2Q9Dellar", "D3Q27Hermite", "DEFAULT_TRANSFORM"
 ]
 
 _ALL_STENCILS = get_subclasses(Stencil, module=lettuce)
@@ -24,15 +24,6 @@ def moment_tensor(e, multiindex):
         return np.prod(np.power(e, multiindex[..., None, :]), axis=-1)
 
 
-def get_default_moment_transform(lattice):
-    if lattice.stencil == D1Q3:
-        return D1Q3Transform(lattice)
-    if lattice.stencil == D2Q9:
-        return D2Q9Lallemand(lattice)
-    else:
-        raise LettuceException(f"No default moment transform for lattice {lattice}.")
-
-
 class Moments:
     def __init__(self, lattice):
         self.rho = moment_tensor(lattice.e, lattice.convert_to_tensor(np.zeros(lattice.D)))
@@ -54,10 +45,10 @@ def __getitem__(self, moment_names):
         return [self.names.index(name) for name in moment_names]
 
     def transform(self, f):
-        return f
+        return self.lattice.einsum("ij,j->i", (self.matrix, f))
 
     def inverse_transform(self, m):
-        return m
+        return self.lattice.einsum("ij,j->i", (self.inverse, m))
 
     def equilibrium(self, m):
         """A very inefficient and basic implementation of the equilibrium moments.
@@ -83,20 +74,14 @@ class D1Q3Transform(Transform):
         [0, 1 / 2, 1 / 2],
         [0, -1 / 2, 1 / 2]
     ])
-    names = ["rho", "j", "e"]
+    names = ["rho", "j_x", "e_xx"]
     supported_stencils = [D1Q3]
 
     def __init__(self, lattice):
         super(D1Q3Transform, self).__init__(lattice, self.names)
         self.matrix = self.lattice.convert_to_tensor(self.matrix)
         self.inverse = self.lattice.convert_to_tensor(self.inverse)
 
-    def transform(self, f):
-        return self.lattice.mv(self.matrix, f)
-
-    def inverse_transform(self, m):
-        return self.lattice.mv(self.inverse, m)
-
     # def equilibrium(self, m):
     #    # TODO
     #    raise NotImplementedError
@@ -125,7 +110,7 @@ class D2Q9Dellar(Transform):
          [1 / 36, -1 / 12, -1 / 12, 1 / 54, 1 / 36, 1 / 54, 1 / 36, -1 / 24, -1 / 24],
          [1 / 36, 1 / 12, -1 / 12, 1 / 54, -1 / 36, 1 / 54, 1 / 36, 1 / 24, -1 / 24]]
     )
-    names = ['rho', 'jx', 'jy', 'Pi_xx', 'Pi_xy', 'PI_yy', 'N', 'Jx', 'Jy']
+    names = ['rho', 'jx', 'jy', 'Pi_xx', 'Pi_xy', 'PI_yy', 'N_xxyy', 'J_xxy', 'J_xyy']
     supported_stencils = [D2Q9]
 
     def __init__(self, lattice):
@@ -135,14 +120,7 @@ def __init__(self, lattice):
         self.matrix = self.lattice.convert_to_tensor(self.matrix)
         self.inverse = self.lattice.convert_to_tensor(self.inverse)
 
-    def transform(self, f):
-        return self.lattice.mv(self.matrix, f)
-
-    def inverse_transform(self, m):
-        return self.lattice.mv(self.inverse, m)
-
     def equilibrium(self, m):
-        warnings.warn("I am not 100% sure if this equilibrium is correct.", ExperimentalWarning)
         meq = torch.zeros_like(m)
         rho = m[0]
         jx = m[1]
@@ -192,12 +170,6 @@ def __init__(self, lattice):
         self.matrix = self.lattice.convert_to_tensor(self.matrix)
         self.inverse = self.lattice.convert_to_tensor(self.inverse)
 
-    def transform(self, f):
-        return self.lattice.mv(self.matrix, f)
-
-    def inverse_transform(self, m):
-        return self.lattice.mv(self.inverse, m)
-
     def equilibrium(self, m):
         """From Lallemand and Luo"""
         warnings.warn("I am not 100% sure if this equilibrium is correct.", ExperimentalWarning)
@@ -414,12 +386,6 @@ def __init__(self, lattice):
         self.matrix = self.lattice.convert_to_tensor(self.matrix)
         self.inverse = self.lattice.convert_to_tensor(self.inverse)
 
-    def transform(self, f):
-        return self.lattice.mv(self.matrix, f)
-
-    def inverse_transform(self, m):
-        return self.lattice.mv(self.inverse, m)
-
     def equilibrium(self, m):
         meq = torch.zeros_like(m)
         rho = m[0]
@@ -454,3 +420,18 @@ def equilibrium(self, m):
         meq[25] = jx * jy * jy * jz * jz / rho ** 4
         meq[26] = jx * jy * jx * jz * jy * jz / rho ** 5
         return meq
+
+
+DEFAULT_TRANSFORM = {
+    D1Q3: D1Q3Transform,
+    D2Q9: D2Q9Dellar,
+    D3Q27: D3Q27Hermite
+}
+
+
+def get_default_moment_transform(lattice):
+    try:
+        transform_class = DEFAULT_TRANSFORM[lattice.stencil]
+    except KeyError:
+        raise LettuceException(f"No default moment transform for lattice {lattice}.")
+    return transform_class(lattice)