Commit 64db569b authored by MischaD's avatar MischaD
Browse files

Finished first implementation

parent 483f9812
......@@ -17,7 +17,10 @@ def main():
def lbmweights(dimension, order, shells, seed):
shell_list = [int(x) for x in re.findall(u"\d+", shells)]
l = Lattice(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
l.calculate_weights()
weights = l.calculate_weights()
import sympy as sp
print(weights)
# TODO sp.N(weights[0]) etc
if __name__ == "__main__":
main()
......
......@@ -3,7 +3,7 @@ import sympy as sp
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
from .utils.utils import analyse_tensor_dimension, get_group, get_subshells, lattice_sum, fill_rhs, approximate_ratio, timer, get_random_vector
class Lattice:
......@@ -21,7 +21,6 @@ class Lattice:
self._shell_list = shell_list
self._seed = seed
self._svd_tolerance = svd_tolerance
np.random.seed(self._seed)
# LES
self._tensor_space_dimension = 0
......@@ -144,6 +143,9 @@ class Lattice:
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"
......@@ -161,7 +163,7 @@ class Lattice:
grand_total_list = []
total_subshells = []
group = get_group(self._dimension) # TODO test output
group = get_group(self._dimension)
for shell in range(len(self._shell_list)):
velocities = self.velocities_for_shell(dimension=self._dimension, shell=self._shell_list[shell])
if len(velocities) == 0:
......@@ -170,20 +172,26 @@ class Lattice:
subshells = get_subshells(velocities, group)
total_subshells.append(len(subshells))
if len(subshells) > 1:
# TODO l.119
raise ReducibleShellException("Shell {} is reducible")
grand_total_list.extend(subshells)
# 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)
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))
logger.info("Non zero shells: ")
for i, shell in enumerate(grand_total_list):
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: ...
......@@ -197,26 +205,27 @@ class Lattice:
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 = 2 * np.random.rand(self._dimension) - 1
vector = vector / np.linalg.norm(vector)
vector = get_random_vector(self._dimension)
#vector = 2 * np.random.rand(self._dimension) - 1
#vector = vector / np.linalg.norm(vector)
rows = []
for shell_velocities in grand_total_list:
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)
logger.info(b)
# --------------------------------------------------------------------------------------------------
logger.warning("A Overwritten for test purposes")
"""
A = np.array([[2., 4., 8.],
[0.50244727, 1.99021093, 8.03915626],
[0.47686341, 2.09254637, 7.6298145]])
"""
U, s, V = np.linalg.svd(A)
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]
......@@ -231,7 +240,7 @@ class Lattice:
raise NoSolutionException("Add shells or lower the order") # change shells ?
rows = rank
s.resize(rows)
rhs.resize((rows, rhs.shape[1]))
rhs.resize((rows, rhs.shape[1]), refcheck=False)
self._singular_values = len(s)
reduced_rhs = np.zeros((rows, rhs.shape[1]))
......@@ -251,6 +260,7 @@ class Lattice:
# ----------------------------------------------- Unique solution
solution = np.dot(np.transpose(V), reduced_rhs)
self._solution = solution
solution_columns = self._order // 2
assert solution_columns == solution.shape[1]
......@@ -262,6 +272,7 @@ class Lattice:
tmp -= solution[i, j] * len(grand_total_list[i])
coefficients[j+1] = tmp
self._coefficients = coefficients
# -------------------------------------------------- Prettify for sympy
weight_polynomials = []
c_s_sq = sp.Symbol("c_s_sq")
......@@ -274,41 +285,16 @@ class Lattice:
for j in range(solution_columns):
equation += approximate_ratio(solution[i, j]) * c_s_sq ** (j + 1)
weight_polynomials.append(equation)
for i in weight_polynomials:
print(i)
# -------------------------------------------------- Calculate Roots
"""
#interval = sp.Interval(-oo, oo)
#for equation in weight_polynomials:
# _interval = sp.solveset(equation > 0, c_s_sq, sp.Reals)
# interval = interval.intersect(_interval)
TOL = 1e-10
tmp = []
for i in range(len(coefficients)):
tmp.append(coefficients[len(coefficients) - 1 - i])
tmp_array = np.array(tmp)
roots = [root for root in np.roots(tmp_array) if abs(root.imag) < TOL and abs(root.real) > TOL]
rows = solution.shape[0]
cols = solution.shape[1]
for j in range(rows):
tmp = []
for i in range(cols):
tmp.append(solution[j, cols - 1 - i])
tmp_array = np.array(tmp)
roots.append([root for root in np.roots(tmp_array) if abs(root.imag) < TOL and abs(root.real) > TOL])
roots.sort()
if not len(roots):
return []
"""
self._weight_polynomials = weight_polynomials
self.calculate_valid_interval()
logger.info(self._interval)
logger.info(self._interval.boundary)
if self._interval == sp.EmptySet():
return None
weights = [x.subs(c_s_sq, self._interval.boundary.inf) for x in weight_polynomials]
return weights
......
......@@ -2,6 +2,7 @@ from .mylog import logger
import itertools
import numpy as np
import time
import random
import math
from fractions import Fraction
......@@ -131,6 +132,24 @@ def get_subshells(shell, group):
return subshells
def get_random_vector(dimension):
tmp = 2
while tmp > 1:
tmp = 0
vector = []
for dim in range(dimension):
val = 2 * random.random() - 1
vector.append(val)
tmp += val ** 2
for dim in range(dimension):
vector[dim] *= 1/ math.sqrt(tmp)
return vector
def double_factorial(n):
assert n >= 0
if n == 0 or n == 1:
......@@ -162,7 +181,7 @@ def fill_rhs(max_rank, shells):
def approximate_ratio(x):
if abs(x) < 1e-12:
if abs(x) < 1e-8:
return 0
denominator_limit = 1000000 if abs(x) >= 0.1 else int(1. / abs(x)) * 1000000
return Fraction(x).limit_denominator(denominator_limit)
......
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