Skip to content

GitLab

Commit 073ef30c authored by MischaD's avatar MischaD
Browse files

Black Codeformatting

parent 3135fec2
Hide whitespace changes
Inline Side-by-side
cli.py
View file @ 073ef30c
......@@ -9,6 +9,7 @@ from lbmweights.utils.mylog import logger
def main():
pass
@timer
def benchmark(x=10):
for i in range(x):
......@@ -17,18 +18,27 @@ def benchmark(x=10):
@main.command()
@click.option('--dimension', type=int, default=2, help='Dimension of the lattice, default 2')
@click.option('--order', type=int, default=4, help='Order of isotropy, default is 4')
@click.option('--shells', type=str, default="1,2,4", help='String of velocity shell size. Squared length of the shell')
@click.option('--seed', type=int, help='Random number to be used as seed')
@click.option(
"--dimension", type=int, default=2, help="Dimension of the lattice, default 2"
)
@click.option("--order", type=int, default=4, help="Order of isotropy, default is 4")
@click.option(
"--shells",
type=str,
default="1,2,4",
help="String of velocity shell size. Squared length of the shell",
)
@click.option("--seed", type=int, help="Random number to be used as seed")
def lbmweights(dimension, order, shells, seed):
shell_list = [int(x) for x in re.findall(u"\d+", shells)]
lattice = Lattice(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
lattice = Lattice(
dimension=dimension, order=order, shell_list=shell_list, seed=seed
)
#benchmark(10)
# benchmark(10)
print(lattice.velocity_set())
if __name__ == "__main__":
main()
lbmweights/exceptions.py
View file @ 073ef30c
class OrderNotImplementedException(Exception):
"""Order in classmethod of Lattice not implemented
"""
pass
class ImpossibleVelocityShellException(Exception):
"""No velocity sets could be constructed. Only positive integer possible
"""
pass
class ReducibleShellException(Exception):
"""Shell is reducible although it should be irreducible"""
pass
class NoSolutionException(Exception):
"""Equation system cannot be solved"""
pass
class InfiniteSolutionsException(Exception):
"""Infinite Solutions. May be solveable with Continue.py in https://github.com/BDuenweg/Lattice-Boltzmann-weights"""
pass
class UninitializedAttributeException(Exception):
"""Calling a function that is not supposed to be called yet"""
pass
\ No newline at end of file
pass
lbmweights/lattice.py
View file @ 073ef30c
......@@ -3,22 +3,64 @@ import sympy as sp
import random
from sympy import oo
from .utils.mylog import logger
from .exceptions import OrderNotImplementedException, ImpossibleVelocityShellException, ReducibleShellException, NoSolutionException, InfiniteSolutionsException, UninitializedAttributeException
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, approximate_ratio, timer, get_random_vector, velocities_for_shell
from .exceptions import (
OrderNotImplementedException,
ImpossibleVelocityShellException,
ReducibleShellException,
NoSolutionException,
InfiniteSolutionsException,
UninitializedAttributeException,
)
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
class Lattice:
BY_NAME = {"D2Q9": {"dimension": 2, "order": 4, "shell_list": [1, 2, 4], "seed": 44},
"D2V17": {"dimension": 2, "order": 6, "shell_list": [1, 2, 4, 8, 9], "seed": 44},
"D2V37": {"dimension": 2, "order": 8, "shell_list": [1, 2, 4, 5, 8, 9, 10, 16], "seed": 44},
"D3Q19": {"dimension": 3, "order": 4, "shell_list": [1, 2, 4], "seed": 44},
"D3Q15": {"dimension": 3, "order": 4, "shell_list": [1, 3, 4], "seed": 44},
"D3Q41-ZOT": {"dimension": 3, "order": 6, "shell_list": [1, 2, 3, 9, 16, 27], "boundary": "sup",
"unwanted_subshells": ["221", "511"], "seed": 44},
}
def __init__(self, dimension=2, order=4, shell_list=[1, 2, 4], seed=None, boundary=None, unwanted_subshells=[], svd_tolerance=1e-8):
BY_NAME = {
"D2Q9": {"dimension": 2, "order": 4, "shell_list": [1, 2, 4], "seed": 44},
"D2V17": {
"dimension": 2,
"order": 6,
"shell_list": [1, 2, 4, 8, 9],
"seed": 44,
},
"D2V37": {
"dimension": 2,
"order": 8,
"shell_list": [1, 2, 4, 5, 8, 9, 10, 16],
"seed": 44,
},
"D3Q19": {"dimension": 3, "order": 4, "shell_list": [1, 2, 4], "seed": 44},
"D3Q15": {"dimension": 3, "order": 4, "shell_list": [1, 3, 4], "seed": 44},
"D3Q41-ZOT": {
"dimension": 3,
"order": 6,
"shell_list": [1, 2, 3, 9, 16, 27],
"boundary": "sup",
"unwanted_subshells": ["221", "511"],
"seed": 44,
},
}
def __init__(
self,
dimension=2,
order=4,
shell_list=[1, 2, 4],
seed=None,
boundary=None,
unwanted_subshells=[],
svd_tolerance=1e-8,
):
"""Create an instance of the Lattice class without calculating the weights.
Default values are chosen to lead to the standard D2Q9 lattice.
Not all input values are valid. :func:`Lattice.init_by_order`
......@@ -79,7 +121,9 @@ class Lattice:
if self._boundary:
string += " - Boundary: {}".format(self._boundary)
if self._unwanted_subshells:
string += " - Unwanted Subshells: {}".format(''.join([s + ', ' for s in self._unwanted_subshells]))
string += " - Unwanted Subshells: {}".format(
"".join([s + ", " for s in self._unwanted_subshells])
)
return string
@classmethod
......@@ -106,9 +150,7 @@ class Lattice:
return cls(**cls.BY_NAME.get(key))
raise ValueError(
"No Lattice known with order {} and dimension {}.".format(
order, dimension
)
"No Lattice known with order {} and dimension {}.".format(order, dimension)
)
return cls(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
......@@ -144,9 +186,13 @@ class Lattice:
group = get_group(self._dimension)
for shell in range(len(self._shell_list)):
velocities = velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
velocities = velocities_for_shell(
dimension=self._dimension, shell=self._shell_list[shell]
)
if len(velocities) == 0:
raise ImpossibleVelocityShellException("No velocity sets found for velocity shell {}.".format(shell))
raise ImpossibleVelocityShellException(
"No velocity sets found for velocity shell {}.".format(shell)
)
# If e.g. dimension = 2, then shell 25 in ambiguous. Both shell (3,4) and (5,0) are added
subshells = get_subshells(velocities, group)
......@@ -155,14 +201,16 @@ class Lattice:
velocities = []
for i_subs, subshell in enumerate(subshells):
# get type for subshell. Example shell: all vectors of type [0, -1, 0] --> 100
type = ''.join([str(int(x)) for x in sorted(abs(subshell[0]), reverse=True)])
type = "".join(
[str(int(x)) for x in sorted(abs(subshell[0]), reverse=True)]
)
if type not in self._unwanted_subshells:
velocities.append(subshell)
else:
velocities = [velocities]
velocity_vectors.extend(velocities)
self._shells = [Shell([np.array([0]*self._dimension)])]
self._shells = [Shell([np.array([0] * self._dimension)])]
for final_velocities in velocity_vectors:
self._shells.append(Shell(final_velocities))
self._all_velocity_vectors = velocity_vectors
......@@ -211,7 +259,9 @@ class Lattice:
overlap = rows - rank
if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
# NO Solution
raise NoSolutionException("Add shells or lower the order") # change shells ?
raise NoSolutionException(
"Add shells or lower the order"
) # change shells ?
rows = rank
s.resize(rows)
rhs.resize((rows, rhs.shape[1]), refcheck=False)
......@@ -227,8 +277,10 @@ class Lattice:
reduced_rhs[i, :] = rhs[i, :] / singular_value
if cols - rows > 0:
raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
"https://github.com/BDuenweg/Lattice-Boltzmann-weights")
raise InfiniteSolutionsException(
"Reduce order, remove shells, or use continue.py in "
"https://github.com/BDuenweg/Lattice-Boltzmann-weights"
)
self._solution = np.dot(np.transpose(V), reduced_rhs)
return self._solution
......@@ -254,7 +306,9 @@ class Lattice:
def _calculate_valid_interval(self):
if not self._weight_polynomials:
raise UninitializedAttributeException("Weight Polynomials not give. Calculate them first!")
raise UninitializedAttributeException(
"Weight Polynomials not give. Calculate them first!"
)
interval = sp.Interval(-oo, oo)
for equation in self._weight_polynomials:
......@@ -264,16 +318,20 @@ class Lattice:
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
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]))
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]))
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
......@@ -319,10 +377,14 @@ class Lattice:
c_s_sq = sp.Symbol("c_s_sq")
if approximate:
self._weight_polynomials = [sp.Poly([approximate_ratio(x) for x in pol_coffs], c_s_sq) for pol_coffs in self._coefficients]
self._weight_polynomials = [
sp.Poly([approximate_ratio(x) for x in pol_coffs], c_s_sq)
for pol_coffs in self._coefficients
]
else:
self._weight_polynomials = [sp.Poly(pol_coffs, c_s_sq) for pol_coffs in
self._coefficients]
self._weight_polynomials = [
sp.Poly(pol_coffs, c_s_sq) for pol_coffs in self._coefficients
]
self._interval = self._calculate_valid_interval()
return [] if self._interval == sp.EmptySet else self._weight_polynomials
......@@ -388,18 +450,7 @@ class Lattice:
if shell.weight == 0:
continue
for velocity in shell.velocities:
velocities[:,i] = velocity
velocities[:, i] = velocity
i += 1
return self._c_s_sq, weights, velocities
lbmweights/shell.py
View file @ 073ef30c
......@@ -6,13 +6,17 @@ class Shell:
def __init__(self, velocities):
self._weight = None
self._c_s_sq = None
self.shell = abs(velocities[0]**2).sum()
self.type = ''.join([str(int(x)) for x in sorted(abs(velocities[0]), reverse=True)])
self.shell = abs(velocities[0] ** 2).sum()
self.type = "".join(
[str(int(x)) for x in sorted(abs(velocities[0]), reverse=True)]
)
self.size = len(velocities)
self.velocities = velocities
def __repr__(self):
repr_string = "Shell: {} - Type: {} - Amount of velocities: {}".format(self.shell, self.type, self.size)
repr_string = "Shell: {} - Type: {} - Amount of velocities: {}".format(
self.shell, self.type, self.size
)
if self._weight is not None:
repr_string += " - Weight: {}".format(self.weight)
return repr_string
......
lbmweights/utils/mylog.py
View file @ 073ef30c
......@@ -19,5 +19,3 @@ ch.setFormatter(formatter)
logger = logging.getLogger("LBMweights")
logger.setLevel(logging.DEBUG)
log_levels = dict(debug=logging.DEBUG, info=logging.INFO, warning=logging.WARNING)
lbmweights/utils/utils.py
View file @ 073ef30c
......@@ -93,11 +93,13 @@ def get_group(dimension):
unit_vectors = np.arange(1, dimension + 1)
permutations = itertools.permutations(unit_vectors, dimension)
sign_combinations = np.zeros((2**dimension, dimension))
sign_combinations = np.zeros((2 ** dimension, dimension))
format_ = "0{}b".format(dimension)
for i in range(2**dimension):
for i in range(2 ** dimension):
binary = format(i, format_)
sign_combinations[i] = np.array([2*int(binary[j]) - 1 for j in range(dimension)])
sign_combinations[i] = np.array(
[2 * int(binary[j]) - 1 for j in range(dimension)]
)
group = []
for permutation in permutations:
......@@ -105,7 +107,7 @@ def get_group(dimension):
for sign_combination in sign_combinations:
group.append(np.multiply(sign_combination, M))
transformation_amount = np.math.factorial(dimension) * 2 ** dimension
assert len(group) == transformation_amount # Todo Exception
assert len(group) == transformation_amount
return group
......@@ -134,6 +136,8 @@ def velocities_for_shell(dimension, shell):
velocities.append(np.array(cur_vec, dtype=int))
return velocities
def contains(arr, list_of_arr):
return any((arr == x).all() for x in list_of_arr)
......@@ -167,7 +171,7 @@ def get_random_vector(dimension):
vector.append(val)
tmp += val ** 2
for dim in range(dimension):
vector[dim] *= 1/ math.sqrt(tmp)
vector[dim] *= 1 / math.sqrt(tmp)
return vector
......@@ -175,7 +179,7 @@ def double_factorial(n):
assert n >= 0
if n == 0 or n == 1:
return 1
return n * double_factorial(n-2)
return n * double_factorial(n - 2)
def lattice_sum(vector, velocities, rank):
......@@ -188,14 +192,8 @@ def lattice_sum(vector, velocities, rank):
return ret / double_factorial(rank - 1)
def approximate_ratio(x):
if abs(x) < 1e-8:
return 0
denominator_limit = 1000000 if abs(x) >= 0.1 else int(1. / abs(x)) * 1000000
denominator_limit = 1000000 if abs(x) >= 0.1 else int(1.0 / abs(x)) * 1000000
return Fraction(x).limit_denominator(denominator_limit)
test/test_lattice.py
View file @ 073ef30c
......@@ -17,21 +17,24 @@ class SolutionExpected(unittest.TestCase):
class TestQ9(unittest.TestCase):
def setUp(self):
self.seed = 20
self.lattice = Lattice(dimension=2, order=4, shell_list=[1, 2, 4], seed=self.seed)
self.lattice = Lattice(
dimension=2, order=4, shell_list=[1, 2, 4], seed=self.seed
)
def testOutput(self):
weights = self.lattice.calculate_weights()
self.assertEqual(len(weights), 4)
interval = self.lattice._interval
self.assertAlmostEqual(sp.N(interval.inf), 1/3, places=6)
self.assertAlmostEqual(sp.N(interval.sup), 2/3, places=6)
self.assertAlmostEqual(sp.N(interval.inf), 1 / 3, places=6)
self.assertAlmostEqual(sp.N(interval.sup), 2 / 3, places=6)
class TestV17(unittest.TestCase):
def setUp(self):
self.seed = 44
self.lattice = Lattice(dimension=2, order=6, shell_list=[1,2,4,8,9], seed=self.seed)
self.lattice = Lattice(
dimension=2, order=6, shell_list=[1, 2, 4, 8, 9], seed=self.seed
)
def testOutput(self):
weights = self.lattice.calculate_weights()
......@@ -44,20 +47,24 @@ class TestV17(unittest.TestCase):
class Test3D(unittest.TestCase):
def setUp(self):
self.seed = 20
self.lattice = Lattice(dimension=3, order=6, shell_list=[1,2,3,4,12,16], seed=self.seed)
self.lattice = Lattice(
dimension=3, order=6, shell_list=[1, 2, 3, 4, 12, 16], seed=self.seed
)
def testOutput(self):
weights = self.lattice.calculate_weights()
self.assertEqual(len(weights), 7)
interval = self.lattice._interval
self.assertAlmostEqual(sp.N(interval.inf), 0.3510760, places=6)
self.assertAlmostEqual(sp.N(interval.sup), 4/9, places=6)
self.assertAlmostEqual(sp.N(interval.sup), 4 / 9, places=6)
class TestV37(unittest.TestCase):
def setUp(self):
self.seed = 20
self.lattice = Lattice(dimension=2, order=8, shell_list=[1,2,4,5,8,9,10,16], seed=self.seed)
self.lattice = Lattice(
dimension=2, order=8, shell_list=[1, 2, 4, 5, 8, 9, 10, 16], seed=self.seed
)
def testOutput(self):
weights = self.lattice.calculate_weights()
......@@ -71,7 +78,12 @@ class TestReducibleShell(unittest.TestCase):
def setUp(self):
self.seed = 20
# 25 can be shell 05 or shell 34
self.lattice = Lattice(dimension=2, order=10, shell_list=[1,2,4,5,8,9,10,13,16,25], seed=44)
self.lattice = Lattice(
dimension=2,
order=10,
shell_list=[1, 2, 4, 5, 8, 9, 10, 13, 16, 25],
seed=44,
)
def testOutput(self):
weights = self.lattice.calculate_weights()
......@@ -84,7 +96,9 @@ class TestReducibleShell(unittest.TestCase):
class TestSupremum(unittest.TestCase):
def setUp(self):
self.seed = 20
self.lattice = Lattice(dimension=2, order=4, shell_list=[1, 2, 4], seed=self.seed, boundary="sup")
self.lattice = Lattice(
dimension=2, order=4, shell_list=[1, 2, 4], seed=self.seed, boundary="sup"
)
def testOutput(self):
weights = self.lattice.calculate_weights()
......@@ -96,7 +110,13 @@ class TestSupremum(unittest.TestCase):
class TestUnwantedSubshells(unittest.TestCase):
def setUp(self):
self.seed = 20
self.lattice = Lattice(dimension=3, order=6, shell_list=[1,2,3,9,16,27], boundary="sup", unwanted_subshells=["221", "511"])
self.lattice = Lattice(
dimension=3,
order=6,
shell_list=[1, 2, 3, 9, 16, 27],
boundary="sup",
unwanted_subshells=["221", "511"],
)
def testOutput(self):
weights = self.lattice.calculate_weights()
......@@ -109,7 +129,9 @@ class TestUnwantedSubshells(unittest.TestCase):
class TestInitSchemes(unittest.TestCase):
def setUp(self):
self.seed = 20
self.lattice = Lattice(dimension=2, order=4, shell_list=[1, 2, 4], seed=self.seed)
self.lattice = Lattice(
dimension=2, order=4, shell_list=[1, 2, 4], seed=self.seed
)
self.lattice_from_order = Lattice.from_order(dimension=2, order=4)
self.lattice_from_name = Lattice.from_name("D2Q9")
......@@ -118,6 +140,6 @@ class TestInitSchemes(unittest.TestCase):
self.assertTrue(self.lattice.weights == self.lattice_from_order.weights)
if __name__ == '__main__':
if __name__ == "__main__":
logger.disabled = True
unittest.main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment