Commit 7189a62e authored by Martin Bauer's avatar Martin Bauer
Browse files

flake8 linter

- removed warnings
- added flake8 as CI target
parent 8e9c74e7
from lbmpy.boundaries.boundaryconditions import NoSlip, UBB, FixedDensity, NeumannByCopy
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
__all__ = ['NoSlip', 'UBB', 'FixedDensity', 'NeumannByCopy', 'LatticeBoltzmannBoundaryHandling']
......@@ -181,7 +181,8 @@ class FixedDensity(Boundary):
density_symbol = cqc.defined_symbols()['density']
density = self._density
equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density).new_without_subexpressions()
equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density)
equilibrium_input = equilibrium_input.new_without_subexpressions()
density_eq = equilibrium_input.main_assignments[0]
assert density_eq.lhs == density_symbol
transformed_density = density_eq.rhs
......
......@@ -544,7 +544,7 @@ def chapman_enskog_ansatz(equation, time_derivative_orders=(1, 3), spatial_deriv
max_expansion_order = max(max_expansion_order, stop_order)
equation = equation.subs(subs_dict)
equation = expand_using_linearity(equation, functions=expanded_pdf_symbols).expand().collect(eps)
result = {eps_order: equation.coeff(eps ** eps_order) for eps_order in range(1, 2*max_expansion_order)}
result = {eps_order: equation.coeff(eps ** eps_order) for eps_order in range(1, 2 * max_expansion_order)}
result[0] = equation.subs(eps, 0)
return result
......
import sympy as sp
from pystencils.derivative import full_diff_expand
from lbmpy.chapman_enskog.chapman_enskog import CeMoment, Diff, take_moments
from lbmpy.chapman_enskog.derivative import full_diff_expand
from lbmpy.chapman_enskog.chapman_enskog import expanded_symbol
from lbmpy.moments import moments_up_to_order, moments_of_order
from collections import OrderedDict, namedtuple
......@@ -35,7 +35,7 @@ def get_solvability_conditions(dim, order):
solvability_conditions = {}
for name in ["\\Pi", "\\Upsilon"]:
for moment_tuple in moments_up_to_order(1, dim=dim):
for k in range(order+1):
for k in range(order + 1):
solvability_conditions[CeMoment(name, superscript=k, moment_tuple=moment_tuple)] = 0
return solvability_conditions
......@@ -72,14 +72,15 @@ def determine_higher_order_moments(epsilon_hierarchy, relaxation_rates, moment_c
for eps_order in range(1, order):
eps_eq = epsilon_hierarchy[eps_order]
for order in range(order+1):
for order in range(order + 1):
eqs = sp.Matrix([full_expand(tm(eps_eq * m)) for m in poly_moments(order, dim)])
unknown_moments = [m for m in eqs.atoms(CeMoment)
if m.superscript == eps_order and sum(m.moment_tuple) == order]
if len(unknown_moments) == 0:
for eq in eqs:
time_diffs_in_expr = [d for d in eq.atoms(Diff) if
(d.target == 't' or d.target == sp.Symbol("t")) and d.superscript == eps_order]
t = sp.Symbol("t")
time_diffs_in_expr = [d for d in eq.atoms(Diff)
if (d.target == 't' or d.target == t) and d.superscript == eps_order]
if len(time_diffs_in_expr) == 0:
continue
assert len(time_diffs_in_expr) == 1, \
......
import sympy as sp
import functools
from lbmpy.chapman_enskog.chapman_enskog import LbMethodEqMoments, CeMoment, take_moments, insert_moments
from pystencils.derivative import Diff, DiffOperator, expand_using_linearity, normalize_diff_order
from pystencils.derivative import collect_derivatives, create_nested_diff
from pystencils.sympyextensions import normalize_product, multidimensional_sum, kronecker_delta
from lbmpy.chapman_enskog.derivative import Diff, DiffOperator, expand_using_linearity, normalize_diff_order
from lbmpy.chapman_enskog.chapman_enskog import LbMethodEqMoments, CeMoment, take_moments, insert_moments
from lbmpy.chapman_enskog.chapman_enskog import expanded_symbol, chapman_enskog_ansatz, remove_higher_order_u
......
from pystencils.derivative import *
import sympy as sp
from pystencils.derivative import Diff
def chapman_enskog_derivative_expansion(expr, label, eps=sp.Symbol("epsilon"), start_order=1, stop_order=4):
......
......@@ -168,8 +168,8 @@ from lbmpy.relaxationrates import relaxation_rate_from_magic_number
from lbmpy.stencils import get_stencil, stencils_have_same_entries
import lbmpy.forcemodels as forcemodels
from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.updatekernels import StreamPullTwoFieldsAccessor, PeriodicTwoFieldsAccessor, CollideOnlyInplaceAccessor, \
create_lbm_kernel, create_stream_pull_with_output_kernel
from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor, PeriodicTwoFieldsAccessor, CollideOnlyInplaceAccessor
from lbmpy.updatekernels import create_lbm_kernel, create_stream_pull_with_output_kernel
def create_lb_function(ast=None, optimization={}, **kwargs):
......
......@@ -23,7 +23,7 @@ def __get_indexed_symbols(passed_symbols, prefix, indices):
if passed_symbols is not None:
return passed_symbols
else:
format_string = "%s_" + "_".join(["%d"]*dim)
format_string = "%s_" + "_".join(["%d"] * dim)
tuple_indices = []
for i in indices:
tuple_indices.append(i if type(i) is tuple else (i,))
......@@ -65,7 +65,7 @@ def __cumulant_raw_moment_transform(index, dependent_var_dict, outer_function, d
subs_dict[result_symbol] = dependent_var_dict[idx]
return result_symbol
zeroth_moment = create_moment_symbol((0,)*dim)
zeroth_moment = create_moment_symbol((0,) * dim)
def outer_function_derivative(n):
x = zeroth_moment
......
......@@ -55,9 +55,10 @@ class StreamPullTwoFieldsAccessor(PdfFieldAccessor):
class PeriodicTwoFieldsAccessor(PdfFieldAccessor):
"""Access scheme that builds periodicity into the kernel.
Introduces a condition on every load, such that at the borders the periodic value is loaded. The periodicity is specified as a tuple of booleans, one for
each direction. The second parameter `ghost_layers` specifies the number of assumed ghost layers of the field.
For the periodic kernel itself no ghost layers are required, however other kernels might need them.
Introduces a condition on every load, such that at the borders the periodic value is loaded. The periodicity is
specified as a tuple of booleans, one for each direction. The second parameter `ghost_layers` specifies the number
of assumed ghost layers of the field. For the periodic kernel itself no ghost layers are required,
however other kernels might need them.
"""
def __init__(self, periodicity, ghost_layers=0):
self._periodicity = periodicity
......@@ -158,6 +159,3 @@ def visualize_pdf_field_accessor(pdf_field_accessor, figure=None):
ax_left.set_title("Read")
ax_right.set_title("Write")
......@@ -6,8 +6,8 @@ r"""
Get started:
-----------
This module offers different models to introduce a body force in the lattice Boltzmann scheme.
If you don't know which model to choose, use :class:`lbmpy.forcemodels.Guo`.
This module offers different models to introduce a body force in the lattice Boltzmann scheme.
If you don't know which model to choose, use :class:`lbmpy.forcemodels.Guo`.
For incompressible collision models the :class:`lbmpy.forcemodels.Buick` model can be better.
......@@ -26,10 +26,10 @@ to the acceleration :math:`\pmb{a}` for all force models.
.. math ::
\sum_q \pmb{c}_q F_q = \pmb{a}
\sum_q \pmb{c}_q F_q = \pmb{a}
The second order moment is different for the forcing models - if it is zero the model is suited for
The second order moment is different for the forcing models - if it is zero the model is suited for
incompressible flows. For weakly compressible collision operators a force model with a corrected second order moment
should be chosen.
......@@ -52,7 +52,7 @@ Models with nonzero second moment have:
F_q = \frac{w_q}{c_s^2} c_{qi} \; a_i + \frac{w_q}{c_s^4} (c_{qi} c_{qj} - c_s^2 \delta_{ij} ) u_j \, a_i
For all force models the computation of the macroscopic velocity has to be adapted (shifted) by adding a term
For all force models the computation of the macroscopic velocity has to be adapted (shifted) by adding a term
:math:`S_{macro}` that we call "macroscopic velocity shift"
.. math ::
......@@ -62,7 +62,7 @@ For all force models the computation of the macroscopic velocity has to be adapt
S_{macro} = \frac{\Delta t}{2} \sum_q F_q
Some models also shift the velocity entering the equilibrium distribution.
Some models also shift the velocity entering the equilibrium distribution.
Comparison
----------
......@@ -70,7 +70,8 @@ Comparison
Force models can be distinguished by 2 options:
**Option 1**:
:math:`C_F = 1` and equilibrium velocity is not shifted, or :math:`C_F=(1 - \frac{\omega}{2})` and equilibrum is shifted.
:math:`C_F = 1` and equilibrium velocity is not shifted, or :math:`C_F=(1 - \frac{\omega}{2})`
and equilibrum is shifted.
**Option 2**:
second velocity moment is zero or :math:`F_i u_j + F_j u_i`
......@@ -212,4 +213,3 @@ class EDM(object):
def default_velocity_shift(density, force):
return [f_i / (2 * density) for f_i in force]
......@@ -141,7 +141,7 @@ def get_pipe_velocity_field(domain_size, u_max, flow_direction=0, diameter=None)
radius = int(round(diameter / 2))
params = [np.arange(s) + 0.5 for s in domain_size]
grid = np.meshgrid(*params, indexing='ij') # type: Tuple[np.array]
grid: Tuple[np.array] = np.meshgrid(*params, indexing='ij')
dist = 0
for i in range(len(domain_size)):
......@@ -202,8 +202,8 @@ def add_black_and_white_image(boundary_handling, image_file, target_slice=None,
# resize necessary if aspect ratio should be constant
if zoomed_image.shape != target_size:
resized_image = np.zeros(target_size, dtype=np.bool)
mid = [(ts - s)//2 for ts, s in zip(target_size, zoomed_image.shape)]
resized_image[mid[0]:zoomed_image.shape[0]+mid[0], mid[1]:zoomed_image.shape[1]+mid[1]] = zoomed_image
mid = [(ts - s) // 2 for ts, s in zip(target_size, zoomed_image.shape)]
resized_image[mid[0]:zoomed_image.shape[0] + mid[0], mid[1]:zoomed_image.shape[1] + mid[1]] = zoomed_image
zoomed_image = resized_image
def callback(*coordinates):
......
......@@ -46,7 +46,7 @@ def create_lbm_split_groups(cr: LbmCollisionRule):
dim = len(stencil[0])
for direction, eq in zip(stencil, cr.main_assignments):
if direction == tuple([0]*dim):
if direction == tuple([0] * dim):
split_groups[0].append(eq.lhs)
continue
......
from lbmpy.methods.abstractlbmethod import AbstractLbMethod
from lbmpy.methods.momentbased import MomentBasedLbMethod, RelaxationInfo
from lbmpy.methods.creationfunctions import create_srt, create_trt, create_trt_with_magic_number, create_mrt_orthogonal, \
create_with_continuous_maxwellian_eq_moments, create_with_discrete_maxwellian_eq_moments, create_trt_kbc, create_mrt_raw, \
create_mrt3
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation, DensityVelocityComputation
from lbmpy.methods.momentbased import RelaxationInfo
from lbmpy.methods.creationfunctions import create_srt, create_trt, create_trt_with_magic_number, create_trt_kbc, \
create_mrt_orthogonal, create_mrt_raw, create_mrt3, \
create_with_continuous_maxwellian_eq_moments, create_with_discrete_maxwellian_eq_moments
__all__ = ['RelaxationInfo',
'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
'create_mrt_orthogonal', 'create_mrt_raw', 'create_mrt3',
'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments']
......@@ -58,4 +58,3 @@ class AbstractLbMethod(abc.ABC):
def get_collision_rule(self):
"""Returns an LbmCollisionRule i.e. an equation collection with a reference to the method.
This collision rule defines the collision operator."""
......@@ -2,7 +2,7 @@ import abc
import sympy as sp
from collections import OrderedDict
from pystencils.assignment_collection import AssignmentCollection
from pystencils.field import Field, Assignment
from pystencils import Field, Assignment
class AbstractConservedQuantityComputation(abc.ABC):
......@@ -63,8 +63,8 @@ class AbstractConservedQuantityComputation(abc.ABC):
be slightly different that the ones used as input for the equilibrium e.g. due to a force model.
:param pdfs: values for the pdf entries
:param output_quantity_names_to_symbols: dict mapping of conserved quantity names (See :func:`conserved_quantities`)
to symbols or field accesses where they should be written to
:param output_quantity_names_to_symbols: dict mapping of conserved quantity names
(See :func:`conserved_quantities`) to symbols or field accesses where they should be written to
"""
@abc.abstractmethod
......@@ -131,11 +131,12 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
def equilibrium_input_equations_from_pdfs(self, pdfs):
dim = len(self._stencil[0])
eq_coll = get_equations_for_zeroth_and_first_order_moment(self._stencil, pdfs, self._symbolOrder0,
self._symbolsOrder1[:dim])
self._symbolsOrder1[:dim])
if self._compressible:
eq_coll = divide_first_order_moments_by_rho(eq_coll, dim)
eq_coll = apply_force_model_shift('equilibrium_velocity_shift', dim, eq_coll, self._forceModel, self._compressible)
eq_coll = apply_force_model_shift('equilibrium_velocity_shift', dim, eq_coll,
self._forceModel, self._compressible)
return eq_coll
def equilibrium_input_equations_from_init_values(self, density=1, velocity=(0, 0, 0)):
......@@ -284,8 +285,8 @@ def divide_first_order_moments_by_rho(assignment_collection, dim):
"""
old_eqs = assignment_collection.main_assignments
rho = old_eqs[0].lhs
new_first_order_moment_eq = [Assignment(eq.lhs, eq.rhs / rho) for eq in old_eqs[1:dim+1]]
new_eqs = [old_eqs[0]] + new_first_order_moment_eq + old_eqs[dim+1:]
new_first_order_moment_eq = [Assignment(eq.lhs, eq.rhs / rho) for eq in old_eqs[1:dim + 1]]
new_eqs = [old_eqs[0]] + new_first_order_moment_eq + old_eqs[dim + 1:]
return assignment_collection.copy(new_eqs)
......
......@@ -310,26 +310,27 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
shear_tensor_trace = sum(shear_tensor_diagonal)
shear_tensor_trace_free_diagonal = [dim * d - shear_tensor_trace for d in shear_tensor_diagonal]
energy_transport_tensor = list(exponents_to_polynomial_representations([a for a in moments_of_order(3, dim, True)
if 3 not in a]))
poly_repr = exponents_to_polynomial_representations
energy_transport_tensor = list(poly_repr([a for a in moments_of_order(3, dim, True)
if 3 not in a]))
explicitly_defined = set(rho + velocity + shear_tensor_off_diagonal +
shear_tensor_diagonal + energy_transport_tensor)
rest = list(set(exponents_to_polynomial_representations(moments_up_to_component_order(2, dim))) - explicitly_defined)
rest = list(set(poly_repr(moments_up_to_component_order(2, dim))) - explicitly_defined)
assert len(rest) + len(explicitly_defined) == 3**dim
# naming according to paper Karlin 2015: Entropic multirelaxation lattice Boltzmann models for turbulent flows
# naming according to paper Karlin 2015: Entropic multi relaxation lattice Boltzmann models for turbulent flows
d = shear_tensor_off_diagonal + shear_tensor_trace_free_diagonal[:-1]
t = [shear_tensor_trace]
q = energy_transport_tensor
if method_name == 'KBC-N1':
decomposition = [d, t+q+rest]
decomposition = [d, t + q + rest]
elif method_name == 'KBC-N2':
decomposition = [d+t, q+rest]
decomposition = [d + t, q + rest]
elif method_name == 'KBC-N3':
decomposition = [d+q, t+rest]
decomposition = [d + q, t + rest]
elif method_name == 'KBC-N4':
decomposition = [d+t+q, rest]
decomposition = [d + t + q, rest]
else:
raise ValueError("Unknown model. Supported models KBC-Nx where x in (1,2,3,4)")
......
......@@ -134,7 +134,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
post_collision_subexpressions=False, include_force_terms=True):
def tuple_to_symbol(exp, prefix):
dim = len(exp)
format_string = prefix + "_" + "_".join(["%d"]*dim)
format_string = prefix + "_" + "_".join(["%d"] * dim)
return sp.Symbol(format_string % exp)
def substitute_conserved_quantities(expressions, cqe):
......
......@@ -60,7 +60,8 @@ def add_entropy_condition(collision_rule, omega_output_field=None):
const_part = decomposition.constant_exprs()
for update_eq in collision_rule.main_assignments:
index = collision_rule.method.post_collision_pdf_symbols.index(update_eq.lhs)
new_eq = Assignment(update_eq.lhs, const_part[index] + omega_s * ds_symbols[index] + omega_h * dh_symbols[index])
new_eq = Assignment(update_eq.lhs,
const_part[index] + omega_s * ds_symbols[index] + omega_h * dh_symbols[index])
new_update_equations.append(new_eq)
new_collision_rule = collision_rule.copy(new_update_equations, collision_rule.subexpressions + subexpressions)
new_collision_rule.simplification_hints['entropic'] = True
......@@ -118,26 +119,26 @@ def add_iterative_entropy_condition(collision_rule, free_omega=None, newton_iter
# 2) get equilibrium from method and define subexpressions for it
eq_terms = [eq.rhs for eq in collision_rule.method.get_equilibrium().main_assignments]
eq_symbols = sp.symbols("entropicFeq_:%d" % (len(eq_terms,)))
eq_symbols = sp.symbols("entropicFeq_:%d" % (len(eq_terms, )))
eq_subexpressions = [Assignment(a, b) for a, b in zip(eq_symbols, eq_terms)]
# 3) find coefficients of entropy derivatives
entropy_diff = sp.diff(discrete_approx_entropy(rr_polynomials, eq_symbols), free_omega)
coefficients_first_diff = [c.expand() for c in reversed(sp.poly(entropy_diff, free_omega).all_coeffs())]
sym_coeff_diff1 = sp.symbols("entropicDiffCoeff_:%d" % (len(coefficients_first_diff,)))
sym_coeff_diff1 = sp.symbols("entropicDiffCoeff_:%d" % (len(coefficients_first_diff, )))
coefficient_eqs = [Assignment(a, b) for a, b in zip(sym_coeff_diff1, coefficients_first_diff)]
sym_coeff_diff2 = [(i+1) * coeff for i, coeff in enumerate(sym_coeff_diff1[1:])]
sym_coeff_diff2 = [(i + 1) * coeff for i, coeff in enumerate(sym_coeff_diff1[1:])]
# 4) define Newtons method update iterations
newton_iteration_equations = []
intermediate_omegas = [sp.Symbol("omega_iter_%i" % (i,)) for i in range(newton_iterations + 1)]
intermediate_omegas[0] = initial_value
intermediate_omegas[-1] = free_omega
for omega_idx in range(len(intermediate_omegas)-1):
for omega_idx in range(len(intermediate_omegas) - 1):
rhs_omega = intermediate_omegas[omega_idx]
lhs_omega = intermediate_omegas[omega_idx+1]
diff1_poly = sum([coeff * rhs_omega**i for i, coeff in enumerate(sym_coeff_diff1)])
diff2_poly = sum([coeff * rhs_omega**i for i, coeff in enumerate(sym_coeff_diff2)])
lhs_omega = intermediate_omegas[omega_idx + 1]
diff1_poly = sum([coeff * rhs_omega ** i for i, coeff in enumerate(sym_coeff_diff1)])
diff2_poly = sum([coeff * rhs_omega ** i for i, coeff in enumerate(sym_coeff_diff2)])
newton_eq = Assignment(lhs_omega, rhs_omega - diff1_poly / diff2_poly)
newton_iteration_equations.append(newton_eq)
......@@ -179,7 +180,7 @@ def discrete_approx_entropy(func, reference):
.. math ::
S = - \sum_i f_i \left( \frac{f_i}{r_i} - 1 \right)
"""
return -sum([f_i * ((f_i / r_i)-1) for f_i, r_i in zip(func, reference)])
return -sum([f_i * ((f_i / r_i) - 1) for f_i, r_i in zip(func, reference)])
def _get_entropy_maximizing_omega(omega_s, f_eq, ds, dh):
......
......@@ -76,9 +76,9 @@ class MomentBasedLbMethod(AbstractLbMethod):
def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False):
relaxation_matrix = sp.eye(len(self.relaxation_rates))
return self._get_collision_rule_with_relaxation_matrix(relaxation_matrix,
conserved_quantity_equations=conserved_quantity_equations,
include_force_terms=include_force_terms)
return self._collision_rule_with_relaxation_matrix(relaxation_matrix,
conserved_quantity_equations=conserved_quantity_equations,
include_force_terms=include_force_terms)
def get_equilibrium_terms(self):
equilibrium = self.get_equilibrium()
......@@ -87,8 +87,8 @@ class MomentBasedLbMethod(AbstractLbMethod):
def get_collision_rule(self, conserved_quantity_equations=None):
d = sp.diag(*self.relaxation_rates)
relaxation_rate_sub_expressions, d = self._generate_relaxation_matrix(d)
ac = self._get_collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
True, conserved_quantity_equations)
ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
True, conserved_quantity_equations)
return ac
def set_zeroth_moment_relaxation_rate(self, relaxation_rate):
......@@ -170,8 +170,8 @@ class MomentBasedLbMethod(AbstractLbMethod):
weights.append(value)
return weights
def _get_collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
conserved_quantity_equations=None):
def _collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
conserved_quantity_equations=None):
f = sp.Matrix(self.pre_collision_pdf_symbols)
pdf_to_moment = self.moment_matrix
m_eq = sp.Matrix(self.moment_equilibrium_values)
......
......@@ -42,4 +42,3 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
rho = sp.Symbol("rho")
equilibrium[0] = rho - sp.expand(sum(equilibrium[1:]))
return create_from_equilibrium(stencil, tuple(equilibrium), relaxation_rate, compressible=True)
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