Commit 7928603a authored by mischa's avatar mischa
Browse files

Assignment Collection in 2D and 3D

parent e2728aff
......@@ -35,11 +35,10 @@ def lbmweights(dimension, order, shells, seed):
lbmweights --dimension=2 --order=4 --shells="1,2,4"
"""
lattice = Lattice.from_name("D2Q9")
lattice = Lattice.from_name("D2V17")
lattice.weights
lattice.assignment_collection
if __name__ == "__main__":
main()
......@@ -60,6 +60,7 @@ class Lattice:
seed=None,
boundary=None,
unwanted_subshells=[],
eq_order=None,
svd_tolerance=1e-8,
):
"""Create an instance of the Lattice class without calculating the weights.
......@@ -77,6 +78,8 @@ class Lattice:
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.
eq_order (int): Explicitly specify the order of accuracy for the equilibrium distribution function.
Only relevant if AssignmentCollections are used. For 2D resonable defaults are used.
seed (float): Seed for the random number generator.
svd_tolerance (float): tolerance to set the singular values to zero after svd.
......@@ -85,6 +88,7 @@ class Lattice:
self._order = order
self._shell_list = shell_list
self._seed = seed
self._eq_order = eq_order
self._svd_tolerance = svd_tolerance
random.seed(self._seed)
......@@ -113,6 +117,7 @@ class Lattice:
# Final values
self._discrete_velocities = []
self._weights = []
self._q = None
self._c_s_sq = None
def __str__(self):
......@@ -164,12 +169,17 @@ class Lattice:
return self._shells
@property
def velocities(self):
"""Get a list of all types of velocities. Velocities with zero weight are included
For D2O9 this equals [00, 10, 11, 20]"""
if not self._discrete_velocities and self._shells:
self._discrete_velocities = [shell.type for shell in self._shells]
return self._discrete_velocities
def q(self):
if self._q:
return self._q
self._q = len(self.velocity_set()[1])
return self._q
@property
def c_s_sq(self):
if self._c_s_sq:
return self._c_s_sq
return self.velocity_set()[0]
@property
def weights(self):
......@@ -178,25 +188,43 @@ class Lattice:
return self.calculate_weights()
@property
def c_s_sq(self):
if self._c_s_sq:
return self._c_s_sq
return self.velocity_set()[0]
def velocities(self):
"""Get a list of all types of velocities. Velocities with zero weight are included
For D2O9 this equals [00, 10, 11, 20]"""
if not self._discrete_velocities and self._shells:
self._discrete_velocities = [shell.type for shell in self._shells]
return self._discrete_velocities
@property
@ps_installed
def assignment_collection(self):
print("okay")
return
def eq_order(self):
"""Calculate the eq order based on the amount of velocities given and dimension.
According to Philippi et. al, 17 and 37 velocities are the minimum amount respectively
to use for third and fourth order accurate models. """
if self._eq_order:
return self._eq_order
if self._dimension == 2:
if self.q >= 37:
self._eq_order = 4
elif self.q >= 17:
self._eq_order = 3
else:
self._eq_order = 2
if self._dimension == 3:
# Should be done explicitly for three dimensions using the constructor
self._eq_order = 2
return self._eq_order
@property
@ps_installed
def assignment_collection(self, ps):
if self._dimension == 3:
return self._assignment_collection3d()
Q = len(self.weights)
weights = self.weights
c_s_sq = float(self.c_s_sq)
c_ix = self.velocities[0]
c_iy = self.velocities[1]
src, dst = ps.fields("src, dst({Q}): double[2D]".format(Q=Q))
c_s_sq, weights, velocities = self.velocity_set()
Q = self.q
c_ix = velocities[0]
c_iy = velocities[1]
src, dst = ps.fields("src({Q}), dst({Q}): double[2D]".format(Q=Q))
density, u_x, u_y, omega = sp.symbols("density, u_x, u_y, omega")
density_formula = sum([src[0, 0](i) for i in range(Q)])
......@@ -208,35 +236,77 @@ class Lattice:
ps.Assignment(u_y, vel_y_formula), ]
for i in range(Q):
feq_formula = weights[i] * density * (
feq_formula = weights[i] * density[0, 0] * (
1
+ (u_x * c_ix[i] + u_y * c_iy[i]) / c_s_sq
+ (u_x * c_ix[i] + u_y * c_iy[i]) ** 2 / (2 * c_s_sq ** 2)
- (u_x ** 2 + u_y ** 2) / (2 * c_s_sq)
)
if len(c_ix) >= 17:
feq_formula += weights[i] * density * (
-2 * (1 / (2 * c_s_sq)) ** 2 * (
c_ix[i] * u_x ** 3 + c_ix[i] * u_x * u_y ** 2 + c_iy[i] * u_y * u_x ** 2 + c_iy[
i] * u_y ** 3)
+ 4 / 3 * (1 / (2 * c_s_sq)) ** 3 * (u_x * c_ix[i] + u_y * c_iy[i]) ** 3
if self.eq_order >= 3:
feq_formula += weights[i] * density[0, 0] * (
-1 / (2 * c_s_sq**2) * (u_x * c_ix[i] + u_y * c_iy[i]) * (u_x ** 2 + u_y ** 2)
+ 4 / 3 * (1 / (2 * c_s_sq)) ** 3 * (u_x * c_ix[i] + u_y * c_iy[i]) ** 3
)
if len(c_ix) >= 37:
feq_formula += weights[i] * density * \
if self.eq_order >= 4:
feq_formula += weights[i] * density[0, 0] * \
(1 / (6 * c_s_sq ** 2) * (
(u_x * c_ix[i] + u_y * c_iy[i]) ** 4 / (4 * c_s_sq ** 2) - (
3 * (u_x ** 2 + u_y ** 2) *
(u_x * c_ix[i] + u_y * c_iy[i]) ** 2) / (2 * c_s_sq) + 3 * (
u_x ** 4 + u_y ** 4)))
for i in range(Q):
symbolic_description.append(
ps.Assignment(dst[-1 * c_iy[i], c_ix[i]](i), omega * feq_formula + (1 - omega) * src[0, 0](i))
)
print("okakychanp")
return ps.AssignmentCollection(symbolic_description)
def _assignment_collection3d(self):
pass
@property
@ps_installed
def _assignment_collection3d(self, ps):
c_s_sq, weights, velocities = self.velocity_set()
Q = self.q
c_ix = velocities[0]
c_iy = velocities[1]
c_iz = velocities[2]
src, dst = ps.fields("src({Q}), dst({Q}): double[2D]".format(Q=Q))
density, u_x, u_y, u_z, omega = sp.symbols("density, u_x, u_y, u_z, omega")
density_formula = sum([src[0, 0](i) for i in range(Q)])
vel_x_formula = (sum([src[0, 0](i) * c_ix[i] for i in range(Q)])) / density
vel_y_formula = (sum([src[0, 0](i) * c_iy[i] for i in range(Q)])) / density
vel_z_formula = (sum([src[0, 0](i) * c_iz[i] for i in range(Q)])) / density
symbolic_description = [ps.Assignment(density, density_formula),
ps.Assignment(u_x, vel_x_formula),
ps.Assignment(u_y, vel_y_formula),
ps.Assignment(u_z, vel_z_formula),]
for i in range(Q):
feq_formula = weights[i] * density[0, 0] * (
1
+ (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) / c_s_sq
+ (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 2 / (2 * c_s_sq ** 2)
- (u_x ** 2 + u_y ** 2 + u_z ** 2) / (2 * c_s_sq)
)
if self.eq_order >= 3:
feq_formula += weights[i] * density[0, 0] * (
-1 / (2 * c_s_sq ** 2) * (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) * (u_x ** 2 + u_y ** 2 + u_z ** 2)
+ 4 / 3 * (1 / (2 * c_s_sq)) ** 3 * (u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 3
)
if self.eq_order >= 4:
feq_formula += weights[i] * density[0, 0] * \
(1 / (6 * c_s_sq ** 2) * (
(u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 4 / (4 * c_s_sq ** 2) - (
3 * (u_x ** 2 + u_y ** 2 + u_z ** 2) *
(u_x * c_ix[i] + u_y * c_iy[i] + u_z * c_iz[i]) ** 2) / (2 * c_s_sq) + 3 * (
u_x ** 4 + u_y ** 4 + u_z ** 4)))
for i in range(Q):
symbolic_description.append(
ps.Assignment(dst[-1 * c_iy[i], c_ix[i]](i), omega * feq_formula + (1 - omega) * src[0, 0](i))
)
return ps.AssignmentCollection(symbolic_description)
def shell_from_type(self, type):
for shell in self._shells:
......
......@@ -63,6 +63,7 @@ def ps_installed(func, *args, **kwargs):
ps = __import__("pystencils")
except ModuleNotFoundError:
raise ModuleNotFoundError("Pystencils not found. Please install it to get AssignmentCollections")
kwargs["ps"] = ps
ret = func(*args, **kwargs)
return ret
......
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