Skip to content
GitLab
Commit 89051627 authored by MischaD's avatar MischaD
Browse files

Refactoring

parent 826f9ed0
Hide whitespace changes
Inline Side-by-side
lbmweights/lattice.py
View file @ 89051627
import numpy as np
import sympy as sp
import random
from sympy import oo
from lbmweights.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, fill_rhs, approximate_ratio, timer, get_random_vector
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, approximate_ratio, timer, get_random_vector
class Lattice:
......@@ -21,18 +22,26 @@ class Lattice:
self._shell_list = shell_list
self._seed = seed
self._svd_tolerance = svd_tolerance
random.seed(self._seed)
# LES
# Preperation
self._tensor_space_dimension = 0
self._all_velocity_vectors = []
self._possible_tensors = []
self._tensor_dimensions = []
self._velocities = 0
self._shells = 0
# LES
self._lhs = []
self._rhs = []
# SVD
self._amount_singular_values = 0
self._singular_values = []
self._solution = []
# Output
self._coefficients = []
self._weight_polynomials = []
self._interval = None
......@@ -111,55 +120,9 @@ class Lattice:
return velocities
@timer
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:
roots = equation.real_roots()
if not roots:
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.eval(mesh_point) < 0:
continue
# positive weights
if i == 0:
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]))
interval = interval.intersect(cur_interval)
self._interval = sp.Complement(interval, sp.FiniteSet(0))
return interval
@timer
def calculate_weights(self):
logger.debug("Entering Calculate Weights")
logger.info(str(self))
import random
random.seed(self._seed)
print("My Seed: {}".format(self._seed))
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)
self._tensor_dimensions = [len(x) for x in self._possible_tensors]
# FROM HERE ...
def calculate_velocity_vectors(self):
velocities_amount = 1
grand_total_list = []
velocity_vectors = []
total_subshells = []
group = get_group(self._dimension)
......@@ -177,60 +140,60 @@ class Lattice:
velocities = subshells
else:
velocities = [velocities]
grand_total_list.extend(velocities)
velocity_vectors.extend(velocities)
velocities_amount += len(velocities)
self._group = group
self._grand_total_list = grand_total_list
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):
for i, shell in enumerate(velocity_vectors):
number_of_velocities = len(shell)
type = tuple(np.sort(abs(shell[0])))
#logger.info("Shell number {} with c_i^2 = {} and {} velocities of type {}".format(i+1, (shell[0]**2).sum(), number_of_velocities, type))
# 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':}
# logger.info("Shell number {} with c_i^2 = {} and {} velocities of type {}".format(i+1, (shell[0]**2).sum(), number_of_velocities, type))
self._all_velocity_vectors = velocity_vectors
return velocity_vectors
def fill_lhs(self):
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 = get_random_vector(self._dimension)
#vector = 2 * np.random.rand(self._dimension) - 1
#vector = vector / np.linalg.norm(vector)
# vector = 2 * np.random.rand(self._dimension) - 1
# vector = vector / np.linalg.norm(vector)
rows = []
for shell_velocities in grand_total_list:
for shell_velocities in self._all_velocity_vectors:
shell_sum = lattice_sum(vector, shell_velocities, rank)
rows.append(shell_sum)
lhs.append(rows)
self._lhs = np.array(lhs)
A = np.array(lhs)
b = fill_rhs(self._order, self._tensor_dimensions)
self._rhs = np.copy(b)
return np.array(lhs)
# --------------------------------------------------------------------------------------------------
U, s, V = np.linalg.svd(A)
self._U = np.copy(U)
self._s = np.copy(s)
self._V = np.copy(V)
rows = A.shape[0]
cols = A.shape[1]
def fill_rhs(self):
rhs = []
cols = self._order // 2
for k, rank in enumerate(np.arange(2, self._order + 2, 2)):
for j in range(self._tensor_dimensions[k]):
rows = []
for i in range(cols):
c_s_power = 2 * i + 2
rows.append(1 if c_s_power == rank else 0)
rhs.append(rows)
self._rhs = np.array(rhs)
return np.array(rhs)
def svd(self):
U, s, V = np.linalg.svd(self._lhs)
rows = self._lhs.shape[0]
cols = self._lhs.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)
rhs = np.dot(np.transpose(U), self._rhs)
if rank < rows:
overlap = rows - rank
if np.linalg.norm(rhs[-overlap:]) > self._svd_tolerance:
......@@ -246,7 +209,6 @@ class Lattice:
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
......@@ -254,41 +216,82 @@ class Lattice:
if cols - rows > 0:
# TODO Better
raise InfiniteSolutionsException("Reduce order, remove shells, or use continue.py in "
"https://github.com/BDuenweg/Lattice-Boltzmann-weights")
"https://github.com/BDuenweg/Lattice-Boltzmann-weights")
self._solution = np.dot(np.transpose(V), reduced_rhs)
return self._solution
# ----------------------------------------------- Unique solution
solution = np.dot(np.transpose(V), reduced_rhs)
self._solution = solution
def calculate_coefficients(self):
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
for i in range(len(self._all_velocity_vectors)):
tmp -= self._solution[i, j] * len(self._all_velocity_vectors[i])
coefficients[j + 1] = tmp
self._coefficients = np.array(coefficients)
# -------------------------------------------------- Prettify for sympy
weight_polynomials = []
c_s_sq = sp.Symbol("c_s_sq")
coeffs = np.zeros((solution.shape[0] + 1, solution.shape[1] + 1))
coeffs = np.zeros((self._solution.shape[0] + 1, self._solution.shape[1] + 1))
coeffs[0, :] = coefficients[-1::-1]
coeffs[1:, :-1] = solution[:,-1::-1]
coeffs[1:, :-1] = self._solution[:, -1::-1]
self._coefficients = coeffs
return self._coefficients
@timer
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:
roots = equation.real_roots()
if not roots:
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.eval(mesh_point) < 0:
continue
# positive weights
if i == 0:
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]))
interval = interval.intersect(cur_interval)
self._interval = sp.Complement(interval, sp.FiniteSet(0))
return self._interval
@timer
def calculate_weights(self):
c_s_sq = sp.Symbol("c_s_sq")
if self._order % 2 == 1:
self._order -= 1
logger.warning("Only odd order valid. Decrease by one")
for k in np.arange(2, self._order + 2, 2):
possible_tensors = analyse_tensor_dimension(k)
self._possible_tensors.append(possible_tensors)
self._tensor_dimensions = [len(x) for x in self._possible_tensors]
self.calculate_velocity_vectors()
self.fill_lhs()
self.fill_rhs()
self.svd()
self.calculate_coefficients()
apporx_coeffs = [[approximate_ratio(x) for x in pol_coffs] for pol_coffs in coeffs]
apporx_coeffs = [[approximate_ratio(x) for x in pol_coffs] for pol_coffs in self._coefficients]
weight_polynomials = [sp.Poly([approximate_ratio(x) for x in pol_coffs], c_s_sq) for pol_coffs in coeffs]
# -------------------------------------------------- Calculate Roots
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 = weight_polynomials
self.calculate_valid_interval()
self._interval = self.calculate_valid_interval()
logger.info(self._interval.boundary)
if self._interval == sp.EmptySet():
return []
weights = [sp.N(x.eval(self._interval.boundary.inf)) for x in weight_polynomials]
weights = [sp.N(x.eval(self._interval.boundary.inf)) for x in self._weight_polynomials]
weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
return weights
......
lbmweights/utils/utils.py
View file @ 89051627
......@@ -167,17 +167,6 @@ def lattice_sum(vector, velocities, rank):
return ret / double_factorial(rank - 1)
def fill_rhs(max_rank, shells):
rhs = []
cols = max_rank // 2
for k, rank in enumerate(np.arange(2, max_rank + 2, 2)):
for j in range(shells[k]):
rows = []
for i in range(cols):
c_s_power = 2 * i + 2
rows.append(1 if c_s_power == rank else 0)
rhs.append(rows)
return np.array(rhs)
def approximate_ratio(x):
......
Supports Markdown
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