Skip to content
GitLab
Commit 957923e3 authored by Markus Holzer's avatar Markus Holzer
Browse files

Fixed some problemns

parent 718a1f35
Pipeline #33618 failed with stages
in 13 minutes and 49 seconds
Expand all Hide whitespace changes
Inline Side-by-side
lbmpy/methods/centeredcumulant/centered_cumulants.py
View file @ 957923e3
......@@ -42,6 +42,29 @@ def get_default_polynomial_cumulants_for_stencil(stencil):
[x**2 * y**2]
]
elif have_same_entries(stencil, get_stencil("D3Q15")):
# D3Q15 central moments by Premnath et al. https://arxiv.org/pdf/1202.6081.pdf.
return [
[sp.sympify(1)], # density is conserved
[x, y, z], # momentum might be affected by forcing
[x * y,
x * z,
y * z,
x ** 2 - y ** 2,
x ** 2 - z ** 2], # shear
[x ** 2 + y ** 2 + z ** 2], # bulk
[x * (x**2 + y**2 + z**2),
y * (x**2 + y**2 + z**2),
z * (x**2 + y**2 + z**2)],
[x * y * z],
[x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2]
]
elif have_same_entries(stencil, get_stencil("D3Q19")):
# D3Q19 cumulants are obtained by pruning the D3Q27 cumulant set as
# described by Coreixas, 2019.
......
lbmpy/moment_transforms/centralmomenttransforms.py
View file @ 957923e3
......@@ -7,7 +7,7 @@ from pystencils.simp.assignment_collection import SymbolGen
from pystencils.sympyextensions import subs_additive, fast_subs
from lbmpy.moments import moment_matrix, monomial_to_polynomial_transformation_matrix,\
set_up_shift_matrix, contained_moments, moments_up_to_order
set_up_shift_matrix, contained_moments, moments_up_to_order, non_aliased_polynomial_moments
from lbmpy.moments import statistical_quantity_symbol as sq_sym
from .abstractmomenttransform import (
......@@ -28,6 +28,10 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
**kwargs):
assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
if moment_polynomials:
# Remove aliases
moment_polynomials = non_aliased_polynomial_moments(moment_polynomials, stencil)
super(PdfsToCentralMomentsByMatrix, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
moment_polynomials=moment_polynomials, **kwargs)
......@@ -101,6 +105,10 @@ class FastCentralMomentTransform(AbstractMomentTransform):
**kwargs):
assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
if moment_polynomials:
# Remove aliases
moment_polynomials = non_aliased_polynomial_moments(moment_polynomials, stencil)
super(FastCentralMomentTransform, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations,
......@@ -318,6 +326,10 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractMomentTransform):
**kwargs):
assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
if moment_polynomials:
# Remove aliases
moment_polynomials = non_aliased_polynomial_moments(moment_polynomials, stencil)
super(PdfsToCentralMomentsByShiftMatrix, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations,
......
lbmpy_tests/centeredcumulant/test_central_moment_transform.py
View file @ 957923e3
import sympy as sp
from pystencils.stencil import have_same_entries
from lbmpy.moments import get_default_moment_set_for_stencil, extract_monomials, statistical_quantity_symbol
from lbmpy.methods.centeredcumulant.centered_cumulants import get_default_polynomial_cumulants_for_stencil
import pytest
from lbmpy.stencils import get_stencil
......@@ -9,20 +11,38 @@ from lbmpy.moment_transforms import (
)
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
def test_forward_transform(stencil):
@pytest.mark.parametrize('type', ['monomial', 'polynomial'])
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
def test_forward_transform(type, stencil):
stencil = get_stencil(stencil)
dim = len(stencil[0])
q = len(stencil)
moment_exponents = list(extract_monomials(get_default_moment_set_for_stencil(stencil)))
moment_exponents = sorted(moment_exponents, key=sum)
if type == 'monomial':
moment_polynomials = get_default_moment_set_for_stencil(stencil)
elif type == 'polynomial':
moment_polynomials = [item for sublist in get_default_polynomial_cumulants_for_stencil(stencil)
for item in sublist]
pdfs = sp.symbols(f"f_:{q}")
rho = sp.Symbol('rho')
u = sp.symbols(f"u_:{dim}")
matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_exponents, rho, u)
fast_transform = FastCentralMomentTransform(stencil, moment_exponents, rho, u)
shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_exponents, rho, u)
matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_polynomials, rho, u)
fast_transform = FastCentralMomentTransform(stencil, moment_polynomials, rho, u)
shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_polynomials, rho, u)
assert matrix_transform.moment_exponents == fast_transform.moment_exponents
assert shift_transform.moment_exponents == fast_transform.moment_exponents
if type == 'monomial' and not have_same_entries(stencil, get_stencil('D3Q15')):
assert matrix_transform.mono_to_poly_matrix == sp.eye(q)
assert fast_transform.mono_to_poly_matrix == sp.eye(q)
assert shift_transform.mono_to_poly_matrix == sp.eye(q)
else:
assert not matrix_transform.mono_to_poly_matrix == sp.eye(q)
assert not fast_transform.mono_to_poly_matrix == sp.eye(q)
assert not shift_transform.mono_to_poly_matrix == sp.eye(q)
moment_exponents = matrix_transform.moment_exponents
f_to_k_matrix = matrix_transform.forward_transform(pdfs)
f_to_k_matrix = f_to_k_matrix.new_without_subexpressions().main_assignments_dict
......@@ -41,21 +61,36 @@ def test_forward_transform(stencil):
assert (rhs_matrix - rhs_fast) == 0, f"Mismatch between matrix and fast transform at {moment_symbol}."
assert (rhs_matrix - rhs_shift) == 0, f"Mismatch between matrix and shift-matrix transform at {moment_symbol}."
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
def test_backward_transform(stencil):
@pytest.mark.parametrize('type', ['monomial', 'polynomial'])
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
def test_backward_transform(type, stencil):
stencil = get_stencil(stencil)
dim = len(stencil[0])
q = len(stencil)
moment_exponents = list(extract_monomials(get_default_moment_set_for_stencil(stencil)))
moment_exponents = sorted(moment_exponents, key=sum)
if type == 'monomial':
moment_polynomials = get_default_moment_set_for_stencil(stencil)
elif type == 'polynomial':
moment_polynomials = [item for sublist in get_default_polynomial_cumulants_for_stencil(stencil)
for item in sublist]
pdfs = sp.symbols(f"f_:{q}")
rho = sp.Symbol('rho')
u = sp.symbols(f"u_:{dim}")
matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_exponents, rho, u)
fast_transform = FastCentralMomentTransform(stencil, moment_exponents, rho, u)
shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_exponents, rho, u)
matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_polynomials, rho, u)
fast_transform = FastCentralMomentTransform(stencil, moment_polynomials, rho, u)
shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_polynomials, rho, u)
assert matrix_transform.moment_exponents == fast_transform.moment_exponents
assert shift_transform.moment_exponents == fast_transform.moment_exponents
if type == 'monomial' and not have_same_entries(stencil, get_stencil('D3Q15')):
assert matrix_transform.mono_to_poly_matrix == sp.eye(q)
assert fast_transform.mono_to_poly_matrix == sp.eye(q)
assert shift_transform.mono_to_poly_matrix == sp.eye(q)
else:
assert not matrix_transform.mono_to_poly_matrix == sp.eye(q)
assert not fast_transform.mono_to_poly_matrix == sp.eye(q)
assert not shift_transform.mono_to_poly_matrix == sp.eye(q)
k_to_f_matrix = matrix_transform.backward_transform(pdfs)
k_to_f_matrix = k_to_f_matrix.new_without_subexpressions().main_assignments_dict
......
lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn_2d.ipynb
View file @ 957923e3
%% Cell type:code id: tags:
``` python
from collections import OrderedDict
import math
from pystencils.session import *
from lbmpy.session import *
from pystencils.boundaries.boundaryhandling import BoundaryHandling
from lbmpy.moments import MOMENT_SYMBOLS
from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti
from lbmpy.phasefield_allen_cahn.kernel_equations import *
from lbmpy.phasefield_allen_cahn.force_model import CentralMomentMultiphaseForceModel
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
```
%% Cell type:code id: tags:
``` python
try:
import pycuda
except ImportError:
pycuda = None
gpu = False
target = 'cpu'
print('No pycuda installed')
if pycuda:
gpu = True
target = 'gpu'
```
%% Cell type:code id: tags:
``` python
stencil_phase = get_stencil("D2Q9")
stencil_hydro = get_stencil("D2Q9")
fd_stencil = get_stencil("D2Q9")
assert(len(stencil_phase[0]) == len(stencil_hydro[0]))
dimensions = len(stencil_phase[0])
```
%% Cell type:code id: tags:
``` python
# domain
L0 = 256
domain_size = (L0, 4 * L0)
# time step
timesteps = 1000
# reference time
reference_time = 8000
# density of the heavier fluid
rho_H = 1.0
# calculate the parameters for the RTI
parameters = calculate_parameters_rti(reference_length=L0,
reference_time=reference_time,
density_heavy=rho_H,
capillary_number=0.44,
reynolds_number=3000,
atwood_number=0.998,
peclet_number=1000,
density_ratio=1000,
viscosity_ratio=100)
# get the parameters
rho_L = parameters.get("density_light")
mu_H = parameters.get("dynamic_viscosity_heavy")
mu_L = parameters.get("dynamic_viscosity_light")
tau_H = parameters.get("relaxation_time_heavy")
tau_L = parameters.get("relaxation_time_light")
sigma = parameters.get("surface_tension")
M = parameters.get("mobility")
gravitational_acceleration = parameters.get("gravitational_acceleration")
drho3 = (rho_H - rho_L)/3
# interface thickness
W = 5
# coeffcient related to surface tension
beta = 12.0 * (sigma/W)
# coeffcient related to surface tension
kappa = 1.5 * sigma*W
# relaxation rate allen cahn (h)
w_c = 1.0/(0.5 + (3.0 * M))
```
%% Cell type:code id: tags:
``` python
dh = ps.create_data_handling((domain_size), periodicity = (True, False), parallel=False, default_target=target)
g = dh.add_array("g", values_per_cell=len(stencil_hydro))
dh.fill("g", 0.0, ghost_layers=True)
h = dh.add_array("h",values_per_cell=len(stencil_phase))
dh.fill("h", 0.0, ghost_layers=True)
g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
dh.fill("g_tmp", 0.0, ghost_layers=True)
h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
dh.fill("h_tmp", 0.0, ghost_layers=True)
u = dh.add_array("u", values_per_cell=dh.dim)
dh.fill("u", 0.0, ghost_layers=True)
C = dh.add_array("C")
dh.fill("C", 0.0, ghost_layers=True)
C_tmp = dh.add_array("C_tmp")
dh.fill("C_tmp", 0.0, ghost_layers=True)
```
%% Cell type:code id: tags:
``` python
tau = 0.5 + tau_L + (C.center) * (tau_H - tau_L)
s8 = 1/(tau)
rho = rho_L + (C.center) * (rho_H - rho_L)
body_force = [0, 0, 0]
body_force[1] = gravitational_acceleration * rho
```
%% Cell type:code id: tags:
``` python
method_phase = create_lb_method(stencil=stencil_phase, method="central_moment", compressible=True,
relaxation_rates=[0, w_c, 1, 1, 1], equilibrium_order=4)
relaxation_rate=1, equilibrium_order=4)
method_phase.set_first_moment_relaxation_rate(w_c)
method_hydro = create_lb_method(stencil=stencil_phase, method="central_moment", compressible=False,
relaxation_rates=[s8, 1, 1], equilibrium_order=4)
relaxation_rate=s8, equilibrium_order=4)
```
%% Cell type:code id: tags:
``` python
# initialize the domain
def Initialize_distributions():
Nx = domain_size[0]
Ny = domain_size[1]
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
x = np.zeros_like(block.midpoint_arrays[0])
x[:, :] = block.midpoint_arrays[0]
y = np.zeros_like(block.midpoint_arrays[1])
y[:, :] = block.midpoint_arrays[1]
y -= 2 * L0
tmp = 0.1 * Nx * np.cos((2 * math.pi * x) / Nx)
init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (W / 2))
block["C"][:, :] = init_values
if gpu:
dh.all_to_gpu()
dh.run_kernel(h_init)
dh.run_kernel(g_init)
```
%% Cell type:code id: tags:
``` python
h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W, fd_stencil=fd_stencil)
g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
```
%% Cell type:code id: tags:
``` python
force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W, fd_stencil=stencil_phase)]
force_model_h = CentralMomentMultiphaseForceModel(force=force_h)
force_g = hydrodynamic_force(g, C, method_hydro, tau, rho_H, rho_L, kappa, beta, body_force)
force_model_g = CentralMomentMultiphaseForceModel(force=force_g, rho=rho)
```
%% Cell type:code id: tags:
``` python
allen_cahn_lb = get_collision_assignments_phase(lb_method=method_phase,
velocity_input=u,
output={'density': C_tmp},
force_model=force_model_h,
symbolic_fields={"symbolic_field": h,
"symbolic_temporary_field": h_tmp},
kernel_type='stream_pull_collide')
allen_cahn_lb = sympy_cse(allen_cahn_lb)
ast_allen_cahn_lb = ps.create_kernel(allen_cahn_lb, target=dh.default_target, cpu_openmp=True)
kernel_allen_cahn_lb = ast_allen_cahn_lb.compile()
```
%% Cell type:code id: tags:
``` python
hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
density=rho,
velocity_input=u,
force_model=force_model_g,
sub_iterations=2,
symbolic_fields={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
kernel_type='collide_stream_push')
hydro_lb_update_rule = sympy_cse(hydro_lb_update_rule)
ast_hydro_lb = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_hydro_lb = ast_hydro_lb.compile()
```
%% Cell type:code id: tags:
``` python
# periodic Boundarys for g, h and C
periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,
streaming_pattern='push')
periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,
streaming_pattern='pull')
# No slip boundary for the phasefield lbm
bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
target=dh.default_target, name='boundary_handling_h',
streaming_pattern='pull')
# No slip boundary for the velocityfield lbm
bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
target=dh.default_target, name='boundary_handling_g',
streaming_pattern='push')
contact_angle = BoundaryHandling(dh, C.name, fd_stencil, target=dh.default_target)
contact = ContactAngle(45, W)
wall = NoSlip()
if dimensions == 2:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
bh_hydro.set_boundary(wall, make_slice[:, 0])
bh_hydro.set_boundary(wall, make_slice[:, -1])
contact_angle.set_boundary(contact, make_slice[:, 0])
contact_angle.set_boundary(contact, make_slice[:, -1])
else:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
bh_hydro.set_boundary(wall, make_slice[:, 0, :])
bh_hydro.set_boundary(wall, make_slice[:, -1, :])
contact_angle.set_boundary(contact, make_slice[:, 0, :])
contact_angle.set_boundary(contact, make_slice[:, -1, :])
bh_allen_cahn.prepare()
bh_hydro.prepare()
```
%% Cell type:code id: tags:
``` python
# definition of the timestep for the immiscible fluids model
def Improved_PhaseField_h_g():
periodic_BC_h()
bh_allen_cahn()
# run the phase-field LB
dh.run_kernel(kernel_allen_cahn_lb)
dh.swap("C", "C_tmp")
contact_angle()
# periodic BC of the phase-field
periodic_BC_C()
dh.run_kernel(kernel_hydro_lb)
periodic_BC_g()
bh_hydro()
# field swaps
dh.swap("h", "h_tmp")
dh.swap("g", "g_tmp")
```
%% Cell type:code id: tags:
``` python
Initialize_distributions()
for i in range(0, timesteps):
Improved_PhaseField_h_g()
```
%% Cell type:code id: tags:
``` python
if gpu:
dh.to_cpu("C")
Ny = domain_size[1]
pos_liquid_front = (np.argmax(dh.cpu_arrays["C"][L0//2, :] > 0.5) - Ny//2)/L0
pos_bubble_front = (np.argmax(dh.cpu_arrays["C"][0, :] > 0.5) - Ny//2)/L0
```
%% Cell type:code id: tags:
``` python
assert(np.isclose(pos_liquid_front, -1e-1, atol=0.01))
assert(np.isclose(pos_bubble_front, 9e-2, atol=0.01))
```
......
lbmpy_tests/test_central_moment.py
View file @ 957923e3
......@@ -4,7 +4,8 @@ import sympy as sp
from lbmpy.creationfunctions import create_lb_method
from lbmpy.forcemodels import Luo
from lbmpy.maxwellian_equilibrium import get_weights
from lbmpy.moments import get_default_moment_set_for_stencil, moment_matrix, set_up_shift_matrix
from lbmpy.moments import moment_matrix, set_up_shift_matrix
from lbmpy.methods.centeredcumulant.centered_cumulants import get_default_polynomial_cumulants_for_stencil
from lbmpy.scenarios import create_lid_driven_cavity
from lbmpy.stencils import get_stencil
......@@ -13,11 +14,11 @@ from lbmpy.moment_transforms import PdfsToCentralMomentsByShiftMatrix
def test_central_moment_ldc():
sc_central_moment = create_lid_driven_cavity((20, 20), method='central_moment',
relaxation_rates=[1.8, 1, 1], equilibrium_order=4,
relaxation_rate=1.8, equilibrium_order=4,
compressible=True, force=(-1e-10, 0))
sc_central_mometn_3D = create_lid_driven_cavity((20, 20, 3), method='central_moment',
relaxation_rates=[1.8, 1, 1, 1, 1], equilibrium_order=4,
relaxation_rate=1.8, equilibrium_order=4,
compressible=True, force=(-1e-10, 0, 0))
sc_central_moment.run(1000)
......@@ -30,7 +31,7 @@ def test_central_moment_class():
stencil = get_stencil("D2Q9")
method = create_lb_method(stencil=stencil, method='central_moment',
relaxation_rates=[1.2, 1.3, 1.4], equilibrium_order=4, compressible=True)
relaxation_rates=[1.2, 1.3, 1.4, 1.5], equilibrium_order=4, compressible=True)
srt = create_lb_method(stencil=stencil, method='srt',
relaxation_rate=1.2, equilibrium_order=4, compressible=True)
......@@ -41,9 +42,10 @@ def test_central_moment_class():
force_model = Luo(force=sp.symbols(f"F_:{2}"))
eq = (rho, 0, 0, rho * cs_sq, rho * cs_sq, 0, 0, 0, rho * cs_sq ** 2)
eq = (rho, 0, 0, 0, 0, 2 * rho * cs_sq, 0, 0, rho * cs_sq ** 2)
default_moments = get_default_moment_set_for_stencil(stencil)
default_moments = get_default_polynomial_cumulants_for_stencil(stencil)
default_moments = [item for sublist in default_moments for item in sublist]
assert method.central_moment_transform_class == PdfsToCentralMomentsByShiftMatrix
assert method.conserved_quantity_computation.zeroth_order_moment_symbol == rho
......@@ -90,11 +92,11 @@ def test_central_moment_class():
method = create_lb_method(stencil="D2Q9", method='central_moment',
relaxation_rates=[1.7, 1.8, 1.2, 1.3, 1.4], equilibrium_order=4, compressible=True)
assert method.relaxation_matrix[0, 0] == 1.7
assert method.relaxation_matrix[1, 1] == 1.8
assert method.relaxation_matrix[2, 2] == 1.8
assert method.relaxation_matrix[0, 0] == 0
assert method.relaxation_matrix[1, 1] == 0
assert method.relaxation_matrix[2, 2] == 0
method = create_lb_method(stencil="D2Q9", method='central_moment',
relaxation_rates=[1.3] * 9, equilibrium_order=4, compressible=True)
assert sum(method.relaxation_rates) == 1.3 * 9
\ No newline at end of file
np.testing.assert_almost_equal(sum(method.relaxation_rates), 1.3 * 6)
\ No newline at end of file
Supports Markdown
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