Commit 21ff9da8 authored by MischaD's avatar MischaD
Browse files

Added velocity set calculation - now ready for usage in lbm-bigger-stencil-neighbourhoods

parent 2839d405
......@@ -17,7 +17,12 @@ 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)
weights = l.calculate_weights()
weights_polynomials = l.calculate_polynomials()
#weights = l.calculate_weights(boundary="inf")
c_s_sq, weights, velocities = l.velocity_set()
print(weights)
print(velocities[0])
print(velocities[1])
if __name__ == "__main__":
......
......@@ -4,7 +4,8 @@ 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
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:
......@@ -40,11 +41,16 @@ class Lattice:
self._singular_values = []
self._solution = []
# Output
# Polynomials
self._coefficients = []
self._weight_polynomials = []
self._interval = None
# Final values
self._discrete_velocities = []
self._weights = []
self._c_s_sq = None
def __str__(self):
string = "D{} - Order: {} - Shells: {}".format(
self._dimension, self._order, str(self._shell_list)
......@@ -87,38 +93,11 @@ class Lattice:
@property
def shells(self):
self._shells
return self._shells
@property
def velocities(self):
self._velocities
@staticmethod
def velocities_for_shell(dimension, shell):
"""
:param dimension: Lattice dimension
:param shell: Squared velocity of the shell
:return: list of velocities with squared length equal to shell
"""
velocities = []
linear_lattice_size = 2 * shell + 1
full_lattice_size = linear_lattice_size ** dimension
for site in range(full_lattice_size):
tmp = site
cur_vel_squared = 0
cur_vec = []
for dim in range(dimension):
coordinate = tmp % linear_lattice_size
tmp = (tmp - coordinate) // linear_lattice_size
shifted_coordinate = coordinate - shell
cur_vec.append(shifted_coordinate)
cur_vel_squared += shifted_coordinate ** 2
if cur_vel_squared == shell:
velocities.append(np.array(cur_vec, dtype=int))
return velocities
return self._discrete_velocities
def calculate_velocity_vectors(self):
velocities_amount = 1
......@@ -127,10 +106,11 @@ class Lattice:
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])
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))
# If e.g. dimension = 2, then shell 25 in ambiguous. Both shell (3,4) and (5,0) are added
subshells = get_subshells(velocities, group)
total_subshells.append(len(subshells))
if len(subshells) > 1:
......@@ -143,15 +123,11 @@ class Lattice:
velocity_vectors.extend(velocities)
velocities_amount += len(velocities)
self._velocities = velocities_amount
self._shells = len(total_subshells) + 1
logger.info("Velocities: {} Shells: {}".format(self._velocities, self._shells))
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))
self._velocities = velocities_amount
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
return velocity_vectors
......@@ -236,7 +212,6 @@ class Lattice:
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!")
......@@ -263,8 +238,7 @@ class Lattice:
self._interval = sp.Complement(interval, sp.FiniteSet(0))
return self._interval
@timer
def calculate_weights(self):
def calculate_polynomials(self):
c_s_sq = sp.Symbol("c_s_sq")
if self._order % 2 == 1:
self._order -= 1
......@@ -286,15 +260,79 @@ class Lattice:
apporx_coeffs = [[approximate_ratio(x) for x in pol_coffs] 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]
if not self._weight_polynomials:
raise NoSolutionException("Something went wrong")
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 self._weight_polynomials]
return self._weight_polynomials
def calculate_reduced_weights(self, boundary="inf"):
"""
Calculate the weights of this lattice for a reduced model.
If the smallest or biggest value of the valid interval is taken, they reduce
the amount of shells necessary by one.
Overwrites the weight values in self.shell
:param boundary: "inf" --> smallest possible value is taken, otherwise the largest
:return: list of weights in the order of their shells.
"""
return self.calculate_weights(boundary=boundary)
def calculate_weights(self, c_s_sq=None, boundary="inf"):
"""
Calculate the weights for this lattice for a given speed of sound value c_s_sq.
If no value is given, the infimum is used, i.e. the smallest possible value.
:param c_s_sq: Float value for the speed of sound, has to be inside of the valid interval
:param boundary: Get reduced model by using the infimum or supremum of the interval. Fallback value if
no c_s_sq is given
:return:
"""
if not self._weight_polynomials:
self.calculate_polynomials()
if not self._interval: # Empty Interval
return []
if not c_s_sq:
if boundary == "inf":
c_s_sq = self._interval.inf
else:
c_s_sq = self._interval.sup
if self._interval.contains(c_s_sq) is sp.EmptySet: # Invalid c_s_sq input
return []
weights = [sp.N(x.eval(c_s_sq)) for x in self._weight_polynomials]
weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
for i, shell in enumerate(self.shells):
shell.set_weight(c_s_sq, weights[i])
self._c_s_sq = c_s_sq
self._weights = weights
return weights
def velocity_set(self):
"""
Returns the values that are important for the Lattice Boltzmann Model, i.e.
a list of the velocity values their corresponding weights and the value for the speed of sound.
:return: (c_s_sq, weights, velocities): c_s_sq
"""
weights = []
for shell in self.shells:
if shell.weight == 0:
continue
weights += [shell.weight] * shell.size
velocities = np.zeros((self._dimension, len(weights)))
i = 0
for shell in self.shells:
if shell.weight == 0:
continue
for velocity in shell.velocities:
velocities[:,i] = velocity
i += 1
return self._c_s_sq, weights, velocities
......
......@@ -109,6 +109,31 @@ def get_group(dimension):
return group
def velocities_for_shell(dimension, shell):
"""
:param dimension: Lattice dimension
:param shell: Squared velocity of the shell
:return: list of velocities with squared length equal to shell
"""
velocities = []
linear_lattice_size = 2 * shell + 1
full_lattice_size = linear_lattice_size ** dimension
for site in range(full_lattice_size):
tmp = site
cur_vel_squared = 0
cur_vec = []
for dim in range(dimension):
coordinate = tmp % linear_lattice_size
tmp = (tmp - coordinate) // linear_lattice_size
shifted_coordinate = coordinate - shell
cur_vec.append(shifted_coordinate)
cur_vel_squared += shifted_coordinate ** 2
if cur_vel_squared == 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)
......
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