Commit 91a0cfad authored by Michael Kuron's avatar Michael Kuron
Browse files

Merge branch 'IncompressibleCumulants' into 'master'

Fix minor issues and remove depricated usage of cumulant LB method

See merge request pycodegen/lbmpy!70
parents 847389c4 63b38380
Pipeline #31368 passed with stage
in 21 minutes and 8 seconds
......@@ -140,7 +140,7 @@
" "
],
"text/plain": [
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc00786c3d0>"
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f1be00e30d0>"
]
},
"execution_count": 2,
......@@ -326,7 +326,7 @@
" "
],
"text/plain": [
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc00774f640>"
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f1b8f8e0400>"
]
},
"execution_count": 5,
......@@ -404,7 +404,7 @@
" "
],
"text/plain": [
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc0075d7b50>"
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f1b8763a4f0>"
]
},
"execution_count": 6,
......@@ -489,7 +489,7 @@
"\n",
"x, y, z = MOMENT_SYMBOLS\n",
"\n",
"moments = mrt_orthogonal_modes_literature(get_stencil(\"D2Q9\"), is_weighted=True, is_cumulant=False)\n",
"moments = mrt_orthogonal_modes_literature(get_stencil(\"D2Q9\"), is_weighted=True)\n",
"moments"
]
},
......@@ -565,7 +565,7 @@
" "
],
"text/plain": [
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc04815aa30>"
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f1b874f5e20>"
]
},
"execution_count": 9,
......@@ -649,7 +649,7 @@
" "
],
"text/plain": [
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc0075b2e50>"
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f1b873e9cd0>"
]
},
"execution_count": 10,
......@@ -759,7 +759,7 @@
" "
],
"text/plain": [
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc007825f70>"
"<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f1b8766a520>"
]
},
"execution_count": 11,
......
......@@ -53,7 +53,6 @@ General:
- ``maxwellian_moments=True``: way to compute equilibrium moments/cumulants, if False the standard
discretized LBM equilibrium is used, otherwise the equilibrium moments are computed from the continuous Maxwellian.
This makes only a difference if sparse stencils are used e.g. D2Q9 and D3Q27 are not affected, D319 and DQ15 are
- ``cumulant=False``: use cumulants instead of moments (deprecated: use method=cumulant directly)
- ``initial_velocity=None``: initial velocity in domain, can either be a tuple (x,y,z) velocity to set a constant
velocity everywhere, or a numpy array with the same size of the domain, with a last coordinate of shape dim to set
velocities on cell level
......@@ -87,8 +86,6 @@ Cumulant methods:
- ``galilean_correction=False``: Special correction for D3Q27 cumulant LBMs. For Details see
:mod:`lbmpy.methods.centeredcumulant.galilean_correction`
- ``pre_simplification=True``: Simplifications applied during the derivaton of the collision rule for cumulant LBMs
For details see :mod:`lbmpy.methods.momentbased.moment_transforms`.
LES methods:
......@@ -112,7 +109,9 @@ Simplifications / Transformations:
load/store streams and thus speeds up the kernel on most architectures
- ``builtin_periodicity=(False,False,False)``: instead of handling periodicity by copying ghost layers, the periodicity
is built into the kernel. This parameters specifies if the domain is periodic in (x,y,z) direction. Even if the
periodicity is built into the kernel, the fields have one ghost layer to be consistent with other functions.
periodicity is built into the kernel, the fields have one ghost layer to be consistent with other functions.
- ``pre_simplification=True``: Simplifications applied during the derivaton of the collision rule for cumulant LBMs
For details see :mod:`lbmpy.methods.momentbased.moment_transforms`.
Field size information:
......@@ -190,6 +189,7 @@ For example, to modify the AST one can run::
func = create_lb_function(ast=ast, ...)
"""
from collections import OrderedDict
from copy import copy
import sympy as sp
......@@ -199,8 +199,9 @@ import lbmpy.methods.centeredcumulant.force_model as cumulant_force_model
from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor
from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule
from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc)
from lbmpy.methods.abstractlbmethod import RelaxationInfo
from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
from lbmpy.methods.momentbased.moment_transforms import PdfsToCentralMomentsByShiftMatrix
from lbmpy.methods.momentbased.moment_transforms import FastCentralMomentTransform
from lbmpy.methods.centeredcumulant.cumulant_transform import CentralMomentsToCumulantsByGeneratingFunc
from lbmpy.methods.creationfunctions import (
create_with_monomial_cumulants, create_with_polynomial_cumulants, create_with_default_polynomial_cumulants)
......@@ -416,7 +417,6 @@ def create_lb_method(**params):
'equilibrium_order': params['equilibrium_order'],
'force_model': force_model,
'maxwellian_moments': params['maxwellian_moments'],
'cumulant': params['cumulant'],
'c_s_sq': params['c_s_sq'],
}
......@@ -443,10 +443,6 @@ def create_lb_method(**params):
if all(get_order(m) < 2 for m in moments):
if params['entropic']:
return relaxation_rates[0]
elif params['cumulant']:
result = relaxation_rates[next_relaxation_rate[0]]
next_relaxation_rate[0] += 1
return result
else:
return 0
res = relaxation_rates[next_relaxation_rate[0]]
......@@ -494,11 +490,27 @@ def create_lb_method_from_existing(method, modification_function):
modification_function: function receiving (moment, equilibrium_value, relaxation_rate) as arguments,
i.e. one row of the relaxation table, returning a modified version
"""
relaxation_table = (modification_function(m, eq, rr)
for m, eq, rr in zip(method.moments, method.moment_equilibrium_values, method.relaxation_rates))
compressible = method.conserved_quantity_computation.compressible
cumulant = isinstance(method, CenteredCumulantBasedLbMethod)
return create_generic_mrt(method.stencil, relaxation_table, compressible, method.force_model, cumulant)
if isinstance(method, CenteredCumulantBasedLbMethod):
rr_dict = OrderedDict()
relaxation_table = (modification_function(m, eq, rr)
for m, eq, rr in
zip(method.cumulants, method.cumulant_equilibrium_values, method.relaxation_rates))
for cumulant, eq_value, rr in relaxation_table:
cumulant = sp.sympify(cumulant)
rr_dict[cumulant] = RelaxationInfo(eq_value, rr)
return CenteredCumulantBasedLbMethod(method.stencil, rr_dict, method.conserved_quantity_computation,
method.force_model,
galilean_correction=method.galilean_correction,
central_moment_transform_class=method.central_moment_transform_class,
cumulant_transform_class=method.cumulant_transform_class)
else:
relaxation_table = (modification_function(m, eq, rr)
for m, eq, rr in
zip(method.moments, method.moment_equilibrium_values, method.relaxation_rates))
return create_generic_mrt(method.stencil, relaxation_table, compressible, method.force_model)
# ----------------------------------------------------------------------------------------------------------------------
......@@ -564,11 +576,10 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'force_model': 'none',
'force': (0, 0, 0),
'maxwellian_moments': True,
'cumulant': False, # Depricated usage. Cumulant is now an own method
'initial_velocity': None,
'galilean_correction': False, # only available for D3Q27 cumulant methods
'central_moment_transform_class': PdfsToCentralMomentsByShiftMatrix,
'central_moment_transform_class': FastCentralMomentTransform,
'cumulant_transform_class': CentralMomentsToCumulantsByGeneratingFunc,
'entropic': False,
......
......@@ -159,21 +159,45 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
self.set_first_moment_relaxation_rate(self._force_model.override_momentum_relaxation_rate)
self.force_model_rr_override = True
@property
def central_moment_transform_class(self):
self._central_moment_transform_class
@property
def cumulants(self):
return tuple(self._cumulant_to_relaxation_info_dict.keys())
@property
def cumulant_equilibrium_values(self):
return tuple([e.equilibrium_value for e in self._cumulant_to_relaxation_info_dict.values()])
@property
def cumulant_transform_class(self):
self._cumulant_transform_class
@property
def first_order_equilibrium_moment_symbols(self, ):
return self._conserved_quantity_computation.first_order_moment_symbols
@property
def force_model(self):
return self._force_model
@property
def galilean_correction(self):
return self._galilean_correction
@property
def relaxation_info_dict(self):
return self._cumulant_to_relaxation_info_dict
@property
def zeroth_order_equilibrium_moment_symbol(self, ):
return self._conserved_quantity_computation.zeroth_order_moment_symbol
def relaxation_rates(self):
return tuple([e.relaxation_rate for e in self._cumulant_to_relaxation_info_dict.values()])
@property
def first_order_equilibrium_moment_symbols(self, ):
return self._conserved_quantity_computation.first_order_moment_symbols
def zeroth_order_equilibrium_moment_symbol(self, ):
return self._conserved_quantity_computation.zeroth_order_moment_symbol
def set_zeroth_moment_relaxation_rate(self, relaxation_rate):
e = sp.Rational(1, 1)
......@@ -200,18 +224,6 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
def set_force_model(self, force_model):
self._force_model = force_model
@property
def cumulants(self):
return sorted(self._cumulant_to_relaxation_info_dict.keys(), key=moment_sort_key)
@property
def cumulant_equilibrium_values(self):
return tuple([e.equilibrium_value for e in self._cumulant_to_relaxation_info_dict.values()])
@property
def relaxation_rates(self):
return tuple([e.relaxation_rate for e in self._cumulant_to_relaxation_info_dict.values()])
def _repr_html_(self):
table = """
<table style="border:none; width: 100%">
......
This diff is collapsed.
......@@ -142,7 +142,7 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
central_moments = self.forward_matrix * f_vec
main_assignments = [Assignment(sq_sym(moment_symbol_base, e), eq)
for e, eq in zip(self.moment_exponents, central_moments)]
symbol_gen = SymbolGen(subexpression_base, dtype=float)
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpression_symbol_generator=symbol_gen)
if simplification:
......@@ -157,7 +157,7 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
moment_vec = sp.Matrix(moments)
pdfs_from_moments = self.backward_matrix * moment_vec
main_assignments = [Assignment(f, eq) for f, eq in zip(pdf_symbols, pdfs_from_moments)]
symbol_gen = SymbolGen(subexpression_base, dtype=float)
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpression_symbol_generator=symbol_gen)
if simplification:
......@@ -228,7 +228,7 @@ class FastCentralMomentTransform(AbstractMomentTransform):
collect_partial_sums(e)
subexpressions = [Assignment(lhs, rhs) for lhs, rhs in subexpressions_dict.items()]
symbol_gen = SymbolGen(subexpression_base, dtype=float)
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification:
......@@ -244,7 +244,7 @@ class FastCentralMomentTransform(AbstractMomentTransform):
pdf_symbols, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT, simplification=False)
raw_equations = raw_equations.new_without_subexpressions()
symbol_gen = SymbolGen(subexpression_base, dtype=float)
symbol_gen = SymbolGen(subexpression_base)
ac = self._split_backward_equations(raw_equations, symbol_gen)
if simplification:
......@@ -441,7 +441,7 @@ class PdfsToRawMomentsTransform(AbstractMomentTransform):
collect_partial_sums(e)
subexpressions += [Assignment(lhs, rhs) for lhs, rhs in partial_sums_dict.items()]
symbol_gen = SymbolGen(subexpression_base, dtype=float)
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
ac.add_simplification_hint('cq_symbols_to_moments', self.get_cq_to_moment_symbols_dict(moment_symbol_base))
......@@ -457,7 +457,7 @@ class PdfsToRawMomentsTransform(AbstractMomentTransform):
post_collision_moments = [sq_sym(moment_symbol_base, e) for e in self.moment_exponents]
rm_to_f_vec = self.inv_moment_matrix * sp.Matrix(post_collision_moments)
main_assignments = [Assignment(f, eq) for f, eq in zip(pdf_symbols, rm_to_f_vec)]
symbol_gen = SymbolGen(subexpression_base, dtype=float)
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpression_symbol_generator=symbol_gen)
ac.add_simplification_hint('stencil', self.stencil)
......
......@@ -152,17 +152,7 @@ def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy,
iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
moment_matrix = mrt_method.moment_matrix
rel = mrt_method.relaxation_rates
eq = mrt_method.moment_equilibrium_values
eq = np.array(eq)
g_vals = [lb_velocity_field.center(i) for i, _ in enumerate(stencil)]
m0 = np.dot(moment_matrix.tolist(), g_vals)
m = m0 - eq
m = m * rel
non_equilibrium = np.dot(moment_matrix.inv().tolist(), m)
non_equilibrium = lb_velocity_field.center_vector - mrt_method.get_equilibrium_terms()
stress_tensor = [0] * 6
# Calculate Stress Tensor MRT
......
......@@ -3,7 +3,7 @@ import sympy as sp
from lbmpy.innerloopsplit import create_lbm_split_groups
from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
from lbmpy.methods.centeredcumulant.simplification import insert_aliases, insert_zeros
from lbmpy.methods.centeredcumulant.simplification import insert_aliases, insert_zeros, insert_constants
from lbmpy.methods.momentbased.momentbasedsimplifications import (
factor_density_after_factoring_relaxation_times, factor_relaxation_rates,
replace_common_quadratic_and_constant_term, replace_density_and_velocity, replace_second_order_velocity_products)
......@@ -36,5 +36,6 @@ def create_simplification_strategy(lb_method, split_inner_loop=False):
elif isinstance(lb_method, CenteredCumulantBasedLbMethod):
s.add(insert_zeros)
s.add(insert_aliases)
s.add(insert_constants)
s.add(lambda ac: ac.new_without_unused_subexpressions())
return s
import pytest
import sympy as sp
import numpy as np
from lbmpy.phasefield.analytical import (
analytic_interface_profile, chemical_potentials_from_free_energy, cosh_integral,
......@@ -7,6 +9,9 @@ from lbmpy.phasefield.analytical import (
symmetric_symbolic_surface_tension)
from pystencils.fd import evaluate_diffs, expand_diff_full
from lbmpy.phasefield.experiments2D import liquid_lens_setup
from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
from lbmpy.phasefield.post_processing import analytic_neumann_angles
def test_analytic_interface_solution():
"""Ensures that the tanh is an analytical solution for the prescribed free energy / chemical potential
......@@ -67,3 +72,27 @@ def test_pressure_tensor():
for f1_i, f2_i in zip(force_chem_pot, force_pressure_tensor):
assert sp.expand(f1_i - f2_i) == 0
def test_neumann_angle():
pytest.importorskip('skimage')
kappa3 = 0.03
alpha = 1
sc = liquid_lens_setup(domain_size=(150, 60), optimization={'target': 'cpu'},
kappas=(0.01, 0.02, kappa3),
cahn_hilliard_relaxation_rates=[np.nan, 1, 3 / 2],
cahn_hilliard_gammas=[1, 1, 1 / 3],
alpha=alpha)
sc.run(10000)
angles = liquid_lens_neumann_angles(sc.concentration[:, :, :])
assert sum(angles) == 360
analytic_angles = analytic_neumann_angles([0.01, 0.02, kappa3])
for ref, simulated in zip(analytic_angles, angles):
assert np.abs(ref - simulated) < 8
# to show the phasefield use:
# plt.phase_plot_for_step(sc)
import pytest
import numpy as np
from lbmpy.creationfunctions import create_lb_method, create_lb_method_from_existing
from lbmpy.methods import create_srt
from lbmpy.stencils import get_stencil
from lbmpy.methods.creationfunctions import create_with_default_polynomial_cumulants
......@@ -10,7 +11,7 @@ from lbmpy.scenarios import create_lid_driven_cavity
def test_weights(stencil_name):
stencil = get_stencil(stencil_name)
cumulant_method = create_with_default_polynomial_cumulants(stencil, [1])
moment_method = create_srt(stencil, 1, cumulant=False, compressible=True, maxwellian_moments=True)
moment_method = create_srt(stencil, 1, compressible=True, maxwellian_moments=True)
assert cumulant_method.weights == moment_method.weights
def test_cumulant_ldc():
......@@ -25,3 +26,25 @@ def test_cumulant_ldc():
sc_cumulant_3D.run(1000)
assert np.isfinite(np.max(sc_cumulant.velocity[:, :]))
assert np.isfinite(np.max(sc_cumulant_3D.velocity[:, :, :]))
def test_create_cumulant_method_from_existing():
method = create_lb_method(stencil='D2Q9', method='cumulant', relaxation_rate=1.5)
old_relaxation_info_dict = method.relaxation_info_dict
def modification_func(cumulant, eq, rate):
if rate == 0:
return cumulant, eq, 1.0
return cumulant, eq, rate
new_method = create_lb_method_from_existing(method, modification_func)
new_relaxation_info_dict = new_method.relaxation_info_dict
for i, (o, n) in enumerate(zip(old_relaxation_info_dict.items(), new_relaxation_info_dict.items())):
assert o[0] == n[0]
assert o[1].equilibrium_value == n[1].equilibrium_value
if o[1].relaxation_rate == 0:
assert n[1].relaxation_rate == 1.0
else:
assert o[1].relaxation_rate == n[1].relaxation_rate
......@@ -6,7 +6,7 @@ def test_gpu_block_size_limiting():
pytest.importorskip("pycuda")
too_large = 2048*2048
opt = {'target': 'gpu', 'gpu_indexing_params': {'block_size': (too_large, too_large, too_large)}}
ast = create_lb_ast(stencil='D3Q19', cumulant=True, relaxation_rate=1.8, optimization=opt,
ast = create_lb_ast(method='cumulant', stencil='D3Q19', relaxation_rate=1.8, optimization=opt,
compressible=True, force_model='schiller')
limited_block_size = ast.indexing.call_parameters((1024, 1024, 1024))
kernel = ast.compile()
......
......@@ -61,22 +61,22 @@ def test_relaxation_rate_setter():
def test_mrt_orthogonal():
m_ref = {}
moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), True, False)
moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), True)
m = create_lb_method(stencil=get_stencil("D2Q9"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_weighted_orthogonal
m_ref[("D2Q9", True)] = m
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q15"), True, False)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q15"), True)
m = create_lb_method(stencil=get_stencil("D3Q15"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_weighted_orthogonal
m_ref[("D3Q15", True)] = m
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True)
m = create_lb_method(stencil=get_stencil("D3Q19"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_weighted_orthogonal
m_ref[("D3Q19", True)] = m
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q27"), False, False)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q27"), False)
m = create_lb_method(stencil=get_stencil("D3Q27"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_orthogonal
m_ref[("D3Q27", False)] = m
......
......@@ -216,7 +216,7 @@ def test_ldc_mrt(action='Testing', plot="off"):
from lbmpy.stencils import get_stencil
if action == 'Testing' or action == 'Regenerate':
print("%s LidDrivenCavity MRT, compressible 0" % action)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True)
compare_lid_driven_cavity(domain_size=(16, 17, 18), lid_velocity=0.005, stencil='D3Q19',
method='MRT', nested_moments=moments, compressible=False, maxwellian_moments=False,
relaxation_rates=[1.3, 1.4, 1.5, 1.25, 1.36, 1.12], action=action, plot=plot)
......
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