Commit b50d3e74 authored by MischaD's avatar MischaD
Browse files

Supremum Boundary + UnwantedSubshells

parent 40ef02fc
......@@ -15,14 +15,26 @@ def main():
@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)]
l = Lattice(dimension=dimension, order=order, shell_list=shell_list, seed=seed)
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])
# 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_polynomials = l.calculate_polynomials()
# weights = l.calculate_weights(boundary="inf")
#lattice1 = Lattice.from_name("D2V37")
#lattice2 = Lattice.from_order(dimension=2, order=6)
#lattice3 = Lattice(dimension=3, order=4, shell_list=[1,2,4])
lattice = Lattice.from_name("D3Q41-ZOT")
lattice.calculate_weights()
cs, weights, velocities = lattice.velocity_set()
print(lattice.weights)
print(lattice.velocities)
#print(velocities.shape[1])
#print(cs)
#print(weights)
#print(velocities)
#print(Lattice.from_name("D3Q41-ZOT").weights)
if __name__ == "__main__":
......
......@@ -9,14 +9,34 @@ from .shell import Shell
class Lattice:
BY_NAME = {"D2Q9": {"dimension": 2, "order": 4, "shell_list": [1, 2, 4]},
"D2V17": {"dimension": 2, "order": 6, "shell_list": [1, 2, 4, 8, 9]},
"D2V37": {"dimension": 2, "order": 8, "shell_list": [1, 2, 4, 5, 8, 9, 10, 16]},
"D3Q19": {"dimension": 3, "order": 4, "shell_list": [1, 2, 4]},
"D3Q15": {"dimension": 3, "order": 4, "shell_list": [1, 3, 4]},
"D3Q41-ZOT": {"dimension": 3, "order": 6, "shell_list": [1, 2, 3, 9, 16, 27], "boundary": "sup",
"unwanted_subshells": ["221", "511"]},
}
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`
Args:
name: Initialize Lattice with the correct values for some known velocities. Possible names are
D2Q9, D2V17, D2V37, D3Q17,
dimension (int): Dimension of the resulting lattice. Default 2.
order (int): Order of tensor isotropy. Only odd orders are
shell_list (list): List of the Squared lengths of the velocities.
unwanted_subshells (list) : List of strings that correspond to the type of a subshell that isn't supposed to
to be part of lattice. Some shells are ambiguous and in the default case all possiblities are taken,
but some lattice, like the D3Q41-ZOT need certain subshells to be remove. E.g. 2 Dimensions and
shell 25 is added. Bother type = "50" and type = "43" could be used. Add one of them to the list of
unwanted subshells to remove them.
seed (float): Seed for the random number generator.
svd_tolerance (float): tolerance to set the singular values to zero after svd.
def __init__(self, dimension=2, order=4, shell_list=[1, 2, 4], seed=None, svd_tolerance=1e-8):
"""
:param dimension:
:param order:
:param shell_list:
:param seed:
:param tolerance: tolerance to set the singular values to zero after svd
"""
self._dimension = dimension
self._order = order
......@@ -30,8 +50,8 @@ class Lattice:
self._all_velocity_vectors = []
self._possible_tensors = []
self._tensor_dimensions = []
self._velocities = 0
self._shells = 0
self._unwanted_subshells = unwanted_subshells
# LES
self._lhs = []
......@@ -45,6 +65,7 @@ class Lattice:
self._coefficients = []
self._weight_polynomials = []
self._interval = None
self._boundary = boundary if boundary == "sup" else "inf"
# Final values
self._discrete_velocities = []
......@@ -60,7 +81,16 @@ class Lattice:
return string
@classmethod
def init_by_order(cls, dimension, order, seed=None):
def from_name(cls, name, seed=None):
"""Initiate the Lattice with the correct values corresponding to a know Lattice in literature.
A complete list can be seen in Lattice.BY_NAME.keys()"""
kwargs = cls.BY_NAME.get(name)
if not kwargs:
raise ValueError("No Lattice with this name is known or implemented")
return cls(**kwargs, seed=seed)
@classmethod
def from_order(cls, dimension, order, seed=None):
"""
:param dimension:
:param order:
......@@ -97,12 +127,18 @@ class Lattice:
@property
def velocities(self):
if not self._discrete_velocities and self._shells:
self._discrete_velocities = [shell.type for shell in self._shells]
return self._discrete_velocities
@property
def weights(self):
if self._weights:
return self._weights
return self.calculate_weights()
def calculate_velocity_vectors(self):
velocities_amount = 1
velocity_vectors = []
total_subshells = []
group = get_group(self._dimension)
for shell in range(len(self._shell_list)):
......@@ -112,19 +148,18 @@ class Lattice:
# 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:
# 3 4, 0 5 can only be given together
velocities = []
for i_subs, subshell in enumerate(subshells):
type = tuple(np.sort(abs(subshell[0])))
velocities = 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)])
if type not in self._unwanted_subshells:
velocities.append(subshell)
else:
velocities = [velocities]
velocity_vectors.extend(velocities)
velocities_amount += len(velocities)
logger.info("Velocities: {} Shells: {}".format(self._velocities, self._shells))
self._velocities = velocities_amount
self._shells = [Shell([np.array([0]*self._dimension)])]
for final_velocities in velocity_vectors:
self._shells.append(Shell(final_velocities))
......@@ -264,7 +299,6 @@ class Lattice:
raise NoSolutionException("Something went wrong")
self._interval = self.calculate_valid_interval()
logger.info(self._interval.boundary)
if self._interval == sp.EmptySet():
return []
return self._weight_polynomials
......@@ -281,9 +315,14 @@ class Lattice:
"""
return self.calculate_weights(boundary=boundary)
def calculate_weights(self, c_s_sq=None, boundary="inf"):
def calculate_weights(self, c_s_sq=None, boundary=None):
"""
Calculate the weights for this lattice for a given speed of sound value c_s_sq.
Calculate the weights for this lattice for a given speed of sound value c_s_sq, or for
a given boundary.
It is first checked, whether a c_s_sq value is given and then boundary.
If both are not giving, self._boundary is taken which is defined in the
init function.
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
......@@ -291,15 +330,19 @@ class Lattice:
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 not boundary:
boundary = self._boundary
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 []
......
......@@ -7,7 +7,7 @@ class Shell:
self._weight = None
self._c_s_sq = None
self.shell = abs(velocities[0]**2).sum()
self.type = tuple(np.sort(abs(velocities[0])))
self.type = ''.join([str(int(x)) for x in sorted(abs(velocities[0]), reverse=True)])
self.size = len(velocities)
self.velocities = velocities
......
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