Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • Sparse
  • WallLaw
  • improved_comm
  • master
  • mr_refactor_wfb
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
  • release/1.3.7
43 results

Target

Select target project
  • ravi.k.ayyala/lbmpy
  • brendan-waters/lbmpy
  • anirudh.jonnalagadda/lbmpy
  • jbadwaik/lbmpy
  • alexander.reinauer/lbmpy
  • itischler/lbmpy
  • he66coqe/lbmpy
  • ev81oxyl/lbmpy
  • Bindgen/lbmpy
  • da15siwa/lbmpy
  • holzer/lbmpy
  • RudolfWeeber/lbmpy
  • pycodegen/lbmpy
13 results
Select Git revision
  • GetterSetterAPI
  • HRR
  • HydroPressure
  • InplaceConfig
  • Outflow
  • PhaseField
  • Sparse
  • UBBVelocity
  • UpdateAPISparse
  • WallLaw
  • WetNodeBoundaries
  • csebug
  • feature/sparse
  • feature/try
  • improved_comm
  • install_requires
  • master
  • phaseField
  • relaxationrates
  • test_martin
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
57 results
Show changes
Showing
with 682 additions and 121 deletions
File moved
......@@ -10,11 +10,18 @@ from lbmpy.macroscopic_value_kernels import (
create_advanced_velocity_setter_collision_rule, pdf_initialization_assignments)
from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.stencils import LBStencil
from pystencils import create_data_handling, create_kernel, make_slice, Target, Backend
from pystencils import CreateKernelConfig
from pystencils import create_data_handling, create_kernel, make_slice, Target
from pystencils.slicing import SlicedGetter
from pystencils.timeloop import TimeLoop
from ._compat import IS_PYSTENCILS_2
if not IS_PYSTENCILS_2:
from pystencils import Backend
class LatticeBoltzmannStep:
def __init__(self, domain_size=None, lbm_kernel=None, periodicity=False,
......@@ -24,7 +31,9 @@ class LatticeBoltzmannStep:
velocity_input_array_name=None, time_step_order='stream_collide', flag_interface=None,
alignment_if_vectorized=64, fixed_loop_sizes=True,
timeloop_creation_function=TimeLoop,
lbm_config=None, lbm_optimisation=None, config=None, **method_parameters):
lbm_config=None, lbm_optimisation=None,
config: CreateKernelConfig | None = None,
**method_parameters):
if optimization is None:
optimization = {}
......@@ -36,7 +45,10 @@ class LatticeBoltzmannStep:
raise ValueError("When passing a data_handling, the domain_size parameter can not be specified")
if config is not None:
target = config.target
if IS_PYSTENCILS_2:
target = config.get_target()
else:
target = config.target
else:
target = optimization.get('target', Target.CPU)
......@@ -49,6 +61,9 @@ class LatticeBoltzmannStep:
default_target=target,
parallel=False)
if lbm_config:
method_parameters['stencil'] = lbm_config.stencil
if 'stencil' not in method_parameters:
method_parameters['stencil'] = LBStencil(Stencil.D2Q9) \
if data_handling.dim == 2 else LBStencil(Stencil.D3Q27)
......@@ -57,7 +72,11 @@ class LatticeBoltzmannStep:
lbm_config, lbm_optimisation, config)
# the parallel datahandling understands only numpy datatypes. Strings lead to an errors
field_dtype = config.data_type.default_factory().numpy_dtype
if IS_PYSTENCILS_2:
from pystencils import create_type
field_dtype = create_type(config.get_option("default_dtype")).numpy_dtype
else:
field_dtype = config.data_type.default_factory().numpy_dtype
if lbm_kernel:
q = lbm_kernel.method.stencil.Q
......@@ -72,12 +91,21 @@ class LatticeBoltzmannStep:
self.density_data_name = name + "_density" if density_data_name is None else density_data_name
self.density_data_index = density_data_index
self._gpu = target == Target.GPU
if IS_PYSTENCILS_2:
self._gpu = target.is_gpu()
else:
self._gpu = target == Target.GPU
layout = lbm_optimisation.field_layout
alignment = False
if config.backend == Backend.C and config.cpu_vectorize_info:
alignment = alignment_if_vectorized
if IS_PYSTENCILS_2:
if config.get_target().is_vector_cpu() and config.cpu.vectorize.enable:
alignment = alignment_if_vectorized
else:
if config.backend == Backend.C and config.cpu_vectorize_info:
alignment = alignment_if_vectorized
self._data_handling.add_array(self._pdf_arr_name, values_per_cell=q, gpu=self._gpu, layout=layout,
latex_name='src', dtype=field_dtype, alignment=alignment)
......@@ -147,10 +175,14 @@ class LatticeBoltzmannStep:
self._sync_tmp = data_handling.synchronization_function([self._tmp_arr_name], stencil_name, target,
stencil_restricted=True)
self._boundary_handling = LatticeBoltzmannBoundaryHandling(self.method, self._data_handling, self._pdf_arr_name,
name=name + "_boundary_handling",
flag_interface=flag_interface,
target=target, openmp=config.cpu_openmp)
self._boundary_handling = LatticeBoltzmannBoundaryHandling(
self.method, self._data_handling, self._pdf_arr_name,
name=name + "_boundary_handling",
flag_interface=flag_interface,
target=target,
openmp=config.cpu_openmp,
**({"default_dtype": field_dtype} if IS_PYSTENCILS_2 else dict())
)
self._lbm_config = lbm_config
self._lbm_optimisation = lbm_optimisation
......@@ -220,7 +252,7 @@ class LatticeBoltzmannStep:
@property
def config(self):
"""Configutation of pystencils parameters"""
return self.config
return self._config
def _get_slice(self, data_name, slice_obj, masked):
if slice_obj is None:
......@@ -436,12 +468,19 @@ class LatticeBoltzmannStep:
rho_field = rho_field.center if self.density_data_index is None else rho_field(self.density_data_index)
vel_field = self._data_handling.fields[self.velocity_data_name]
if IS_PYSTENCILS_2:
gen_config = CreateKernelConfig(target=Target.CPU)
gen_config.cpu.openmp.enable = self._config.cpu.openmp.get_option("enable")
gen_config.default_dtype = self._config.get_option("default_dtype")
else:
gen_config = CreateKernelConfig(target=Target.CPU, cpu_openmp=self._config.cpu_openmp)
getter_eqs = cqc.output_equations_from_pdfs(pdf_field.center_vector,
{'density': rho_field, 'velocity': vel_field})
getter_kernel = create_kernel(getter_eqs, target=Target.CPU, cpu_openmp=self._config.cpu_openmp).compile()
getter_kernel = create_kernel(getter_eqs, config=gen_config).compile()
setter_eqs = pdf_initialization_assignments(lb_method, rho_field,
vel_field.center_vector, pdf_field.center_vector)
setter_eqs = create_simplification_strategy(lb_method)(setter_eqs)
setter_kernel = create_kernel(setter_eqs, target=Target.CPU, cpu_openmp=self._config.cpu_openmp).compile()
setter_kernel = create_kernel(setter_eqs, config=gen_config).compile()
return getter_kernel, setter_kernel
from typing import Sequence, Any
from abc import ABC, abstractmethod
import numpy as np
import sympy as sp
from ._compat import IS_PYSTENCILS_2
if not IS_PYSTENCILS_2:
raise ImportError("`lbmpy.lookup_tables` is only available when running with pystencils 2.x")
from pystencils import Assignment
from pystencils.sympyextensions import TypedSymbol
from pystencils.types.quick import Arr
from pystencils.types import UserTypeSpec, create_type
class LookupTables(ABC):
@abstractmethod
def get_array_declarations(self) -> list[Assignment]:
pass
class NeighbourOffsetArrays(LookupTables):
@staticmethod
def neighbour_offset(dir_idx, stencil):
if isinstance(sp.sympify(dir_idx), sp.Integer):
return stencil[dir_idx]
else:
return tuple(
[
sp.IndexedBase(symbol, shape=(1,))[dir_idx]
for symbol in NeighbourOffsetArrays._offset_symbols(stencil)
]
)
@staticmethod
def _offset_symbols(stencil):
q = len(stencil)
dim = len(stencil[0])
return [
TypedSymbol(f"neighbour_offset_{d}", Arr(create_type("int32"), q))
for d in ["x", "y", "z"][:dim]
]
def __init__(self, stencil, offsets_dtype: UserTypeSpec = np.int32):
self._offsets_dtype = create_type(
offsets_dtype
) # TODO: Currently, this has no effect
self._stencil = stencil
self._dim = len(stencil[0])
def get_array_declarations(self) -> list[Assignment]:
array_symbols = NeighbourOffsetArrays._offset_symbols(self._stencil)
return [
Assignment(arrsymb, tuple((d[i] for d in self._stencil)))
for i, arrsymb in enumerate(array_symbols)
]
class MirroredStencilDirections(LookupTables):
@staticmethod
def mirror_stencil(direction, mirror_axis):
assert mirror_axis <= len(
direction
), f"only {len(direction)} axis available for mirage"
direction = list(direction)
direction[mirror_axis] = -direction[mirror_axis]
return tuple(direction)
@staticmethod
def _mirrored_symbol(mirror_axis, stencil):
axis = ["x", "y", "z"]
q = len(stencil)
return TypedSymbol(
f"{axis[mirror_axis]}_axis_mirrored_stencil_dir", Arr(create_type("int32"), q)
)
def __init__(self, stencil, mirror_axis, dtype=np.int32):
self._offsets_dtype = create_type(dtype) # TODO: Currently, this has no effect
self._mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(
mirror_axis, stencil
)
self._mirrored_directions = tuple(
stencil.index(
MirroredStencilDirections.mirror_stencil(direction, mirror_axis)
)
for direction in stencil
)
def get_array_declarations(self) -> list[Assignment]:
return [Assignment(self._mirrored_stencil_symbol, self._mirrored_directions)]
class LbmWeightInfo(LookupTables):
def __init__(self, lb_method, data_type="double"):
self._weights = lb_method.weights
self._weights_array = TypedSymbol("weights", Arr(create_type(data_type), len(self._weights)))
def weight_of_direction(self, dir_idx, lb_method=None):
if isinstance(sp.sympify(dir_idx), sp.Integer):
assert lb_method is not None
return lb_method.weights[dir_idx].evalf(17)
else:
return sp.IndexedBase(self._weights_array, shape=(1,))[dir_idx]
def get_array_declarations(self) -> list[Assignment]:
return [Assignment(self._weights_array, tuple(self._weights))]
class TranslationArraysNode(LookupTables):
def __init__(self, array_content: Sequence[tuple[TypedSymbol, Sequence[Any]]]):
self._decls = [
Assignment(symb, tuple(content)) for symb, content in array_content
]
def __str__(self):
return "Variable PDF Access Translation Arrays"
def __repr__(self):
return "Variable PDF Access Translation Arrays"
def get_array_declarations(self) -> list[Assignment]:
return self._decls
import functools
import sympy as sp
from copy import deepcopy
from lbmpy.simplificationfactory import create_simplification_strategy
from pystencils import create_kernel, CreateKernelConfig
from pystencils import create_kernel, CreateKernelConfig, Assignment
from pystencils.field import Field, get_layout_of_array
from pystencils.enums import Target
from pystencils import Target
from lbmpy.advanced_streaming.utility import get_accessor, Timestep
from lbmpy.relaxationrates import get_shear_relaxation_rate
from lbmpy.utils import second_order_moment_tensor
def pdf_initialization_assignments(lb_method, density, velocity, pdfs,
streaming_pattern='pull', previous_timestep=Timestep.BOTH,
set_pre_collision_pdfs=False):
"""Assignments to initialize the pdf field with equilibrium"""
def get_field_accesses(lb_method, pdfs, streaming_pattern, previous_timestep, pre_collision_pdfs):
if isinstance(pdfs, Field):
accessor = get_accessor(streaming_pattern, previous_timestep)
if set_pre_collision_pdfs:
if pre_collision_pdfs:
field_accesses = accessor.read(pdfs, lb_method.stencil)
else:
field_accesses = accessor.write(pdfs, lb_method.stencil)
elif streaming_pattern == 'pull' and not set_pre_collision_pdfs:
elif streaming_pattern == 'pull' and not pre_collision_pdfs:
field_accesses = pdfs
else:
raise ValueError("Invalid value of pdfs: A PDF field reference is required to derive "
+ f"initialization assignments for streaming pattern {streaming_pattern}.")
return field_accesses
def get_individual_or_common_fraction_field(psm_config):
if psm_config.individual_fraction_field is not None:
return psm_config.individual_fraction_field
else:
return psm_config.fraction_field
def pdf_initialization_assignments(lb_method, density, velocity, pdfs, psm_config=None,
streaming_pattern='pull', previous_timestep=Timestep.BOTH,
set_pre_collision_pdfs=False):
"""Assignments to initialize the pdf field with equilibrium"""
field_accesses = get_field_accesses(lb_method, pdfs, streaming_pattern, previous_timestep, set_pre_collision_pdfs)
if isinstance(density, Field):
density = density.center
if isinstance(velocity, Field):
velocity = velocity.center_vector
......@@ -35,23 +49,40 @@ def pdf_initialization_assignments(lb_method, density, velocity, pdfs,
setter_eqs = lb_method.get_equilibrium(conserved_quantity_equations=inp_eqs)
setter_eqs = setter_eqs.new_with_substitutions({sym: field_accesses[i]
for i, sym in enumerate(lb_method.post_collision_pdf_symbols)})
if lb_method.fraction_field is not None:
if psm_config is None:
raise ValueError("If setting up LBM with PSM, please specify a PSM config in the macroscopic setter")
else:
for equ in setter_eqs:
if equ.lhs in lb_method.first_order_equilibrium_moment_symbols:
pos = lb_method.first_order_equilibrium_moment_symbols.index(equ.lhs)
new_rhs = 0
if isinstance(equ.rhs, sp.core.Add):
for summand in equ.rhs.args:
if summand in velocity:
new_rhs += (1.0 - psm_config.fraction_field.center) * summand
else:
new_rhs += summand.subs(lb_method.fraction_field, psm_config.fraction_field.center)
else:
new_rhs += (1.0 - psm_config.fraction_field.center) * equ.rhs
fraction_field = get_individual_or_common_fraction_field(psm_config)
for p in range(psm_config.max_particles_per_cell):
new_rhs += psm_config.object_velocity_field(p * lb_method.dim + pos) * fraction_field.center(p)
setter_eqs.subexpressions.remove(equ)
setter_eqs.subexpressions.append(Assignment(equ.lhs, new_rhs))
return setter_eqs
def macroscopic_values_getter(lb_method, density, velocity, pdfs,
def macroscopic_values_getter(lb_method, density, velocity, pdfs, psm_config=None,
streaming_pattern='pull', previous_timestep=Timestep.BOTH,
use_pre_collision_pdfs=False):
if isinstance(pdfs, Field):
accessor = get_accessor(streaming_pattern, previous_timestep)
if use_pre_collision_pdfs:
field_accesses = accessor.read(pdfs, lb_method.stencil)
else:
field_accesses = accessor.write(pdfs, lb_method.stencil)
elif streaming_pattern == 'pull' and not use_pre_collision_pdfs:
field_accesses = pdfs
else:
raise ValueError("Invalid value of pdfs: A PDF field reference is required to derive "
+ f"getter assignments for streaming pattern {streaming_pattern}.")
field_accesses = get_field_accesses(lb_method, pdfs, streaming_pattern, previous_timestep, use_pre_collision_pdfs)
cqc = lb_method.conserved_quantity_computation
assert not (velocity is None and density is None)
output_spec = {}
......@@ -59,12 +90,54 @@ def macroscopic_values_getter(lb_method, density, velocity, pdfs,
output_spec['velocity'] = velocity
if density is not None:
output_spec['density'] = density
return cqc.output_equations_from_pdfs(field_accesses, output_spec)
getter_equ = cqc.output_equations_from_pdfs(field_accesses, output_spec)
if lb_method.fraction_field is not None:
if psm_config.fraction_field is None:
raise ValueError("If setting up LBM with PSM, please specify a PSM config in the macroscopic getter")
else:
if lb_method.force_model is not None:
for equ in getter_equ:
if equ.lhs in lb_method.force_model.symbolic_force_vector:
new_rhs = equ.rhs.subs(lb_method.fraction_field, psm_config.fraction_field.center)
getter_equ.subexpressions.remove(equ)
getter_equ.subexpressions.append(Assignment(equ.lhs, new_rhs))
for i, equ in enumerate(getter_equ.main_assignments[-lb_method.dim:]):
new_rhs = (1.0 - psm_config.fraction_field.center) * equ.rhs
fraction_field = get_individual_or_common_fraction_field(psm_config)
for p in range(psm_config.max_particles_per_cell):
new_rhs += psm_config.object_velocity_field(p * lb_method.dim + i) * fraction_field.center(p)
getter_equ.main_assignments.remove(equ)
getter_equ.main_assignments.append(Assignment(equ.lhs, new_rhs))
getter_equ.topological_sort()
return getter_equ
macroscopic_values_setter = pdf_initialization_assignments
def strain_rate_tensor_getter(lb_method, strain_rate_tensor, pdfs, streaming_pattern='pull',
previous_timestep=Timestep.BOTH, use_pre_collision_pdfs=False):
field_accesses = get_field_accesses(lb_method, pdfs, streaming_pattern, previous_timestep, use_pre_collision_pdfs)
if isinstance(strain_rate_tensor, Field):
strain_rate_tensor = strain_rate_tensor.center_vector
omega_s = get_shear_relaxation_rate(lb_method)
equilibrium = lb_method.equilibrium_distribution
rho = equilibrium.density if equilibrium.compressible else equilibrium.background_density
f_neq = sp.Matrix([field_accesses[i] for i in range(lb_method.stencil.Q)]) - lb_method.get_equilibrium_terms()
pi = second_order_moment_tensor(f_neq, lb_method.stencil)
strain_rate_tensor_equ = - 1.5 * (omega_s / rho) * pi
assignments = [Assignment(strain_rate_tensor[i * lb_method.stencil.D + j], strain_rate_tensor_equ[i, j])
for i in range(lb_method.stencil.D) for j in range(lb_method.stencil.D)]
return assignments
def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None,
ghost_layers=1, iteration_slice=None,
field_layout='numpy', target=Target.CPU,
......
......@@ -27,16 +27,16 @@ import warnings
import numpy as np
# Optional packages cpuinfo, pycuda and psutil for hardware queries
# Optional packages cpuinfo, cupy and psutil for hardware queries
try:
from cpuinfo import get_cpu_info
except ImportError:
get_cpu_info = None
try:
from pycuda.autoinit import device
import cupy
except ImportError:
device = None
cupy = None
try:
from psutil import virtual_memory
......@@ -113,9 +113,11 @@ def memory_sizes_of_current_machine():
if 'l3_cache_size' in cpu_info:
result['L3'] = cpu_info['l3_cache_size']
if device:
size = device.total_memory() / (1024 * 1024)
result['GPU'] = "{0:.0f} MB".format(size)
if cupy:
for i in range(cupy.cuda.runtime.getDeviceCount()):
device = cupy.cuda.Device(i)
size = device.mem_info[1] / (1024 * 1024)
result[f'GPU{i}'] = "{0:.0f} MB".format(size)
if virtual_memory:
mem = virtual_memory()
......@@ -124,7 +126,7 @@ def memory_sizes_of_current_machine():
if not result:
warnings.warn("Couldn't query for any local memory size."
"Install py-cpuinfo to get cache sizes, psutil for RAM size and pycuda for GPU memory size.")
"Install py-cpuinfo to get cache sizes, psutil for RAM size and cupy for GPU memory size.")
return result
......
......@@ -76,38 +76,49 @@ def discrete_maxwellian_equilibrium(stencil, rho=sp.Symbol("rho"), u=sp.symbols(
"""
weights = get_weights(stencil, c_s_sq)
assert stencil.Q == len(weights)
u = u[:stencil.D]
res = [discrete_equilibrium(e_q, u, rho, w_q, order, c_s_sq, compressible) for w_q, e_q in zip(weights, stencil)]
return tuple(res)
def discrete_equilibrium(v=sp.symbols("v_:3"), u=sp.symbols("u_:3"), rho=sp.Symbol("rho"), weight=sp.Symbol("w"),
order=2, c_s_sq=sp.Symbol("c_s") ** 2, compressible=True):
"""
Returns the common discrete LBM equilibrium depending on the mesoscopic velocity and the directional lattice weight
Args:
v: symbols for mesoscopic velocity
u: symbols for macroscopic velocity
rho: sympy symbol for the density
weight: symbol for stencil weights
order: highest order of velocity terms (for hydrodynamics order 2 is sufficient)
c_s_sq: square of speed of sound
compressible: compressibility
"""
rho_outside = rho if compressible else sp.Rational(1, 1)
rho_inside = rho if not compressible else sp.Rational(1, 1)
res = []
for w_q, e_q in zip(weights, stencil):
e_times_u = 0
for c_q_alpha, u_alpha in zip(e_q, u):
e_times_u += c_q_alpha * u_alpha
fq = rho_inside + e_times_u / c_s_sq
e_times_u = 0
for c_q_alpha, u_alpha in zip(v, u):
e_times_u += c_q_alpha * u_alpha
if order <= 1:
res.append(fq * rho_outside * w_q)
continue
fq = rho_inside + e_times_u / c_s_sq
u_times_u = 0
for u_alpha in u:
u_times_u += u_alpha * u_alpha
fq += sp.Rational(1, 2) / c_s_sq**2 * e_times_u ** 2 - sp.Rational(1, 2) / c_s_sq * u_times_u
if order <= 1:
return fq * rho_outside * weight
if order <= 2:
res.append(fq * rho_outside * w_q)
continue
u_times_u = 0
for u_alpha in u:
u_times_u += u_alpha * u_alpha
fq += sp.Rational(1, 2) / c_s_sq**2 * e_times_u ** 2 - sp.Rational(1, 2) / c_s_sq * u_times_u
fq += sp.Rational(1, 6) / c_s_sq**3 * e_times_u**3 - sp.Rational(1, 2) / c_s_sq**2 * u_times_u * e_times_u
if order <= 2:
return fq * rho_outside * weight
res.append(sp.expand(fq * rho_outside * w_q))
fq += sp.Rational(1, 6) / c_s_sq**3 * e_times_u**3 - sp.Rational(1, 2) / c_s_sq**2 * u_times_u * e_times_u
return tuple(res)
return sp.expand(fq * rho_outside * weight)
@disk_cache
......
......@@ -67,8 +67,8 @@ class AbstractLbMethod(abc.ABC):
return d
@property
def subs_dict_relxation_rate(self):
"""returns a dictonary which maps the replaced numerical relaxation rates to its original numerical value"""
def subs_dict_relaxation_rate(self):
"""returns a dictionary which maps the replaced numerical relaxation rates to its original numerical value"""
result = dict()
for i in range(self._stencil.Q):
result[self.symbolic_relaxation_matrix[i, i]] = self.relaxation_matrix[i, i]
......@@ -99,26 +99,33 @@ class AbstractLbMethod(abc.ABC):
# -------------------------------- Helper Functions ----------------------------------------------------------------
def _generate_symbolic_relaxation_matrix(self, relaxation_rates=None):
def _generate_symbolic_relaxation_matrix(self, relaxation_rates=None, relaxation_rates_modifier=None):
"""
This function replaces the numbers in the relaxation matrix with symbols in this case, and returns also
the subexpressions, that assign the number to the newly introduced symbol
"""
rr = relaxation_rates if relaxation_rates is not None else self.relaxation_rates
unique_relaxation_rates = set()
subexpressions = {}
symbolic_relaxation_rates = list()
for relaxation_rate in rr:
relaxation_rate = sp.sympify(relaxation_rate)
if relaxation_rate not in unique_relaxation_rates:
# special treatment for zero, sp.Zero would be an integer ..
if isinstance(relaxation_rate, sp.Symbol):
symbolic_relaxation_rate = relaxation_rate
else:
if isinstance(relaxation_rate, Zero):
relaxation_rate = 0.0
if not isinstance(relaxation_rate, sp.Symbol):
rt_symbol = sp.Symbol(f"rr_{len(subexpressions)}")
subexpressions[relaxation_rate] = rt_symbol
unique_relaxation_rates.add(relaxation_rate)
if relaxation_rate in subexpressions:
symbolic_relaxation_rate = subexpressions[relaxation_rate]
else:
symbolic_relaxation_rate = sp.Symbol(f"rr_{len(subexpressions)}")
subexpressions[relaxation_rate] = symbolic_relaxation_rate
symbolic_relaxation_rates.append(symbolic_relaxation_rate)
new_rr = [subexpressions[sp.sympify(e)] if sp.sympify(e) in subexpressions else e for e in rr]
substitutions = [Assignment(e[1], e[0]) for e in subexpressions.items()]
if relaxation_rates_modifier is not None:
symbolic_relaxation_rates = [r * relaxation_rates_modifier for r in symbolic_relaxation_rates]
else:
for srr in symbolic_relaxation_rates:
assert isinstance(srr, sp.Symbol)
return substitutions, sp.diag(*new_rr)
return substitutions, sp.diag(*symbolic_relaxation_rates)
......@@ -162,7 +162,7 @@ def create_with_continuous_maxwellian_equilibrium(stencil, moment_to_relaxation_
def create_from_equilibrium(stencil, equilibrium, conserved_quantity_computation,
moment_to_relaxation_rate_dict,
collision_space_info=CollisionSpaceInfo(CollisionSpace.POPULATIONS),
zero_centered=False, force_model=None):
zero_centered=False, force_model=None, fraction_field=None):
r"""
Creates a lattice Boltzmann method in either population, moment, or central moment space, from a given
discrete velocity set and equilibrium distribution.
......@@ -202,18 +202,22 @@ def create_from_equilibrium(stencil, equilibrium, conserved_quantity_computation
if cspace.collision_space == CollisionSpace.POPULATIONS:
return MomentBasedLbMethod(stencil, equilibrium, mom_to_rr_dict, conserved_quantity_computation=cqc,
force_model=force_model, zero_centered=zero_centered,
fraction_field=fraction_field,
moment_transform_class=None)
elif cspace.collision_space == CollisionSpace.RAW_MOMENTS:
return MomentBasedLbMethod(stencil, equilibrium, mom_to_rr_dict, conserved_quantity_computation=cqc,
force_model=force_model, zero_centered=zero_centered,
fraction_field=fraction_field,
moment_transform_class=cspace.raw_moment_transform_class)
elif cspace.collision_space == CollisionSpace.CENTRAL_MOMENTS:
return CentralMomentBasedLbMethod(stencil, equilibrium, mom_to_rr_dict, conserved_quantity_computation=cqc,
force_model=force_model, zero_centered=zero_centered,
fraction_field=fraction_field,
central_moment_transform_class=cspace.central_moment_transform_class)
elif cspace.collision_space == CollisionSpace.CUMULANTS:
return CumulantBasedLbMethod(stencil, equilibrium, mom_to_rr_dict, conserved_quantity_computation=cqc,
force_model=force_model, zero_centered=zero_centered,
fraction_field=fraction_field,
central_moment_transform_class=cspace.central_moment_transform_class,
cumulant_transform_class=cspace.cumulant_transform_class)
......@@ -304,7 +308,7 @@ def create_trt_with_magic_number(stencil, relaxation_rate, magic_number=sp.Ratio
relaxation_rate_odd_moments=rr_odd, **kwargs)
def create_mrt_raw(stencil, relaxation_rates, continuous_equilibrium=True, **kwargs):
def create_mrt_raw(stencil, relaxation_rates, continuous_equilibrium=True, conserved_moments=True, **kwargs):
r"""
Creates a MRT method using non-orthogonalized moments.
......@@ -318,6 +322,7 @@ def create_mrt_raw(stencil, relaxation_rates, continuous_equilibrium=True, **kwa
relaxation_rates: relaxation rates (inverse of the relaxation times) for each moment
continuous_equilibrium: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments.
conserved_moments: If lower order moments are conserved or not.
Returns:
:class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
"""
......@@ -325,7 +330,7 @@ def create_mrt_raw(stencil, relaxation_rates, continuous_equilibrium=True, **kwa
check_and_set_mrt_space(CollisionSpace.RAW_MOMENTS)
moments = get_default_moment_set_for_stencil(stencil)
nested_moments = [(c,) for c in moments]
rr_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D)
rr_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D, conserved_moments)
if continuous_equilibrium:
return create_with_continuous_maxwellian_equilibrium(stencil, rr_dict, **kwargs)
else:
......@@ -333,7 +338,7 @@ def create_mrt_raw(stencil, relaxation_rates, continuous_equilibrium=True, **kwa
def create_central_moment(stencil, relaxation_rates, nested_moments=None,
continuous_equilibrium=True, **kwargs):
continuous_equilibrium=True, conserved_moments=True, **kwargs):
r"""
Creates moment based LB method where the collision takes place in the central moment space.
......@@ -346,6 +351,7 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
nested_moments: a list of lists of modes, grouped by common relaxation times.
continuous_equilibrium: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments.
conserved_moments: If lower order moments are conserved or not.
Returns:
:class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` instance
"""
......@@ -367,7 +373,11 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
if not nested_moments:
nested_moments = cascaded_moment_sets_literature(stencil)
rr_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D)
rr_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D, conserved_moments)
if 'fraction_field' in kwargs and kwargs['fraction_field'] is not None:
relaxation_rates_modifier = (1.0 - kwargs['fraction_field'])
rr_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D,
relaxation_rates_modifier=relaxation_rates_modifier)
if continuous_equilibrium:
return create_with_continuous_maxwellian_equilibrium(stencil, rr_dict, **kwargs)
......@@ -455,7 +465,7 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
def create_mrt_orthogonal(stencil, relaxation_rates, continuous_equilibrium=True, weighted=None,
nested_moments=None, **kwargs):
nested_moments=None, conserved_moments=True, **kwargs):
r"""
Returns an orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
These MRT methods are just one specific version - there are many MRT methods possible for all these stencils
......@@ -476,6 +486,7 @@ def create_mrt_orthogonal(stencil, relaxation_rates, continuous_equilibrium=True
nested_moments: a list of lists of modes, grouped by common relaxation times. If this argument is not provided,
Gram-Schmidt orthogonalization of the default modes is performed. The default modes equal the
raw moments except for the separation of the shear and bulk viscosity.
conserved_moments: If lower order moments are conserved or not.
"""
continuous_equilibrium = _deprecate_maxwellian_moments(continuous_equilibrium, kwargs)
check_and_set_mrt_space(CollisionSpace.RAW_MOMENTS)
......@@ -507,7 +518,8 @@ def create_mrt_orthogonal(stencil, relaxation_rates, continuous_equilibrium=True
nested_moments[2] = shear_moments
nested_moments.insert(3, bulk_moment)
moment_to_relaxation_rate_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D)
moment_to_relaxation_rate_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments,
stencil.D, conserved_moments)
if continuous_equilibrium:
return create_with_continuous_maxwellian_equilibrium(stencil,
......@@ -518,8 +530,7 @@ def create_mrt_orthogonal(stencil, relaxation_rates, continuous_equilibrium=True
# ----------------------------------------- Cumulant method creators ---------------------------------------------------
def create_cumulant(stencil, relaxation_rates, cumulant_groups, **kwargs):
def create_cumulant(stencil, relaxation_rates, cumulant_groups, conserved_moments=True, **kwargs):
r"""Creates a cumulant-based lattice Boltzmann method.
Args:
......@@ -531,12 +542,19 @@ def create_cumulant(stencil, relaxation_rates, cumulant_groups, **kwargs):
that the force is applied correctly to the momentum groups
cumulant_groups: Nested sequence of polynomial expressions defining the cumulants to be relaxed. All cumulants
within one group are relaxed with the same relaxation rate.
conserved_moments: If lower order moments are conserved or not.
kwargs: See :func:`create_with_continuous_maxwellian_equilibrium`
Returns:
:class:`lbmpy.methods.cumulantbased.CumulantBasedLbMethod` instance
"""
cumulant_to_rr_dict = _get_relaxation_info_dict(relaxation_rates, cumulant_groups, stencil.D)
cumulant_to_rr_dict = _get_relaxation_info_dict(relaxation_rates, cumulant_groups, stencil.D, conserved_moments)
if 'fraction_field' in kwargs and kwargs['fraction_field'] is not None:
relaxation_rates_modifier = (1.0 - kwargs['fraction_field'])
cumulant_to_rr_dict = _get_relaxation_info_dict(relaxation_rates, cumulant_groups, stencil.D,
relaxation_rates_modifier=relaxation_rates_modifier)
kwargs.setdefault('collision_space_info', CollisionSpaceInfo(CollisionSpace.CUMULANTS))
if kwargs['collision_space_info'].collision_space != CollisionSpace.CUMULANTS:
......@@ -547,7 +565,7 @@ def create_cumulant(stencil, relaxation_rates, cumulant_groups, **kwargs):
**kwargs)
def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
def create_with_monomial_cumulants(stencil, relaxation_rates, conserved_moments=True, **kwargs):
r"""Creates a cumulant lattice Boltzmann model using the given stencil's canonical monomial cumulants.
Args:
......@@ -557,6 +575,7 @@ def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
used for determine the viscosity of the simulation. All other cumulants are relaxed with one.
If a cumulant force model is provided the first order cumulants are relaxed with two to ensure
that the force is applied correctly to the momentum groups
conserved_moments: If lower order moments are conserved or not.
kwargs: See :func:`create_cumulant`
Returns:
......@@ -565,10 +584,10 @@ def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
# Get monomial moments
cumulants = get_default_moment_set_for_stencil(stencil)
cumulant_groups = [(c,) for c in cumulants]
return create_cumulant(stencil, relaxation_rates, cumulant_groups, **kwargs)
return create_cumulant(stencil, relaxation_rates, cumulant_groups, conserved_moments, **kwargs)
def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs):
def create_with_default_polynomial_cumulants(stencil, relaxation_rates, conserved_moments=True, **kwargs):
r"""Creates a cumulant lattice Boltzmann model based on the default polynomial set of :cite:`geier2015`.
Args: See :func:`create_cumulant`.
......@@ -578,10 +597,11 @@ def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs
"""
# Get polynomial groups
cumulant_groups = cascaded_moment_sets_literature(stencil)
return create_cumulant(stencil, relaxation_rates, cumulant_groups, **kwargs)
return create_cumulant(stencil, relaxation_rates, cumulant_groups, conserved_moments, **kwargs)
def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim,
conserved_moments=True, relaxation_rates_modifier=None):
r"""Creates a dictionary where each moment is mapped to a relaxation rate.
Args:
......@@ -591,6 +611,7 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
in the moment group.
nested_moments: list of lists containing the moments.
dim: dimension
conserved_moments: If lower order moments are conserved or not.
"""
result = OrderedDict()
......@@ -599,7 +620,9 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
rr = iter(relaxation_rates(group))
for moment in group:
result[moment] = next(rr)
if relaxation_rates_modifier is not None:
for key in result:
result[key] *= relaxation_rates_modifier
return result
number_of_moments = 0
......@@ -618,12 +641,18 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
if len(relaxation_rates) == 1:
for group in nested_moments:
for moment in group:
if get_order(moment) <= 1:
result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = relaxation_rates[0]
if conserved_moments:
if get_order(moment) <= 1:
result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = relaxation_rates[0]
else:
result[moment] = 1.0
else:
result[moment] = 1.0
if is_shear_moment(moment, dim) or get_order(moment) <= 1:
result[moment] = relaxation_rates[0]
else:
result[moment] = 1.0
# if relaxation rate for each moment is specified they are all used
if len(relaxation_rates) == number_of_moments:
......@@ -647,15 +676,25 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
rr = next(rr_iter)
next_rr = False
for moment in group:
if get_order(moment) <= 1:
result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = shear_rr
elif is_bulk_moment(moment, dim):
result[moment] = bulk_rr
if conserved_moments:
if get_order(moment) <= 1:
result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = shear_rr
elif is_bulk_moment(moment, dim):
result[moment] = bulk_rr
else:
next_rr = True
result[moment] = rr
else:
next_rr = True
result[moment] = rr
if is_shear_moment(moment, dim) or get_order(moment) <= 1:
result[moment] = shear_rr
elif is_bulk_moment(moment, dim):
result[moment] = bulk_rr
else:
next_rr = True
result[moment] = rr
except StopIteration:
raise ValueError("Not enough relaxation rates are specified. You can either specify one relaxation rate, "
"which is used as the shear relaxation rate. In this case, conserved moments are "
......@@ -664,6 +703,9 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
"moment. In this case, conserved moments are also "
"relaxed with 0. The last possibility is to specify a relaxation rate for each moment, "
"including conserved moments")
if relaxation_rates_modifier is not None:
for key in result:
result[key] *= relaxation_rates_modifier
return result
......
from .cumulantbasedmethod import CumulantBasedLbMethod
from .galilean_correction import add_galilean_correction
from .fourth_order_correction import add_fourth_order_correction
__all__ = ['CumulantBasedLbMethod', 'add_galilean_correction']
__all__ = ['CumulantBasedLbMethod', 'add_galilean_correction', 'add_fourth_order_correction']
......@@ -16,11 +16,8 @@ def insert_logs(ac, **kwargs):
def insert_log_products(ac, **kwargs):
def callback(asm):
rhs = asm.rhs
if isinstance(rhs, sp.log):
if rhs.find(sp.log):
return True
if isinstance(rhs, sp.Mul):
if any(isinstance(arg, sp.log) for arg in rhs.args):
return True
return False
return insert_subexpressions(ac, callback, **kwargs)
......
......@@ -23,7 +23,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
This class implements cumulant-based lattice boltzmann methods which relax all the non-conserved quantities
as either monomial or polynomial cumulants. It is mostly inspired by the work presented in :cite:`geier2015`.
This method is implemented modularily as the transformation from populations to central moments to cumulants
This method is implemented modularly as the transformation from populations to central moments to cumulants
is governed by subclasses of :class:`lbmpy.moment_transforms.AbstractMomentTransform`
which can be specified by constructor argument. This allows the selection of the most efficient transformation
for a given setup.
......@@ -47,7 +47,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
def __init__(self, stencil, equilibrium, relaxation_dict,
conserved_quantity_computation=None,
force_model=None, zero_centered=False,
force_model=None, zero_centered=False, fraction_field=None,
central_moment_transform_class=BinomialChimeraTransform,
cumulant_transform_class=CentralMomentsToCumulantsByGeneratingFunc):
assert isinstance(conserved_quantity_computation,
......@@ -63,6 +63,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
self._cqc = conserved_quantity_computation
self._force_model = force_model
self._zero_centered = zero_centered
self._fraction_field = fraction_field
self._weights = None
self._cumulant_transform_class = cumulant_transform_class
self._central_moment_transform_class = central_moment_transform_class
......@@ -72,6 +73,10 @@ class CumulantBasedLbMethod(AbstractLbMethod):
"""Force model employed by this method."""
return self._force_model
@property
def fraction_field(self):
return self._fraction_field
@property
def relaxation_info_dict(self):
"""Dictionary mapping this method's cumulants to their relaxation rates and equilibrium values.
......@@ -124,7 +129,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
@property
def zeroth_order_equilibrium_moment_symbol(self, ):
"""Returns a symbol referring to the zeroth-order moment of this method's equilibrium distribution,
which is the area under it's curve
which is the area under its curve
(see :attr:`lbmpy.equilibrium.AbstractEquilibrium.zeroth_order_moment_symbol`)."""
return self._equilibrium.zeroth_order_moment_symbol
......
import sympy as sp
from lbmpy.moment_transforms import PRE_COLLISION_MONOMIAL_CUMULANT, POST_COLLISION_CUMULANT
from lbmpy.moments import MOMENT_SYMBOLS, statistical_quantity_symbol
from lbmpy.stencils import Stencil, LBStencil
from pystencils import Assignment, AssignmentCollection
from .cumulantbasedmethod import CumulantBasedLbMethod
X, Y, Z = MOMENT_SYMBOLS
CORRECTED_FOURTH_ORDER_POLYNOMIALS = [X ** 2 * Y ** 2 - 2 * X ** 2 * Z ** 2 + Y ** 2 * Z ** 2,
X ** 2 * Y ** 2 + X ** 2 * Z ** 2 - 2 * Y ** 2 * Z ** 2,
X ** 2 * Y ** 2 + X ** 2 * Z ** 2 + Y ** 2 * Z ** 2,
X ** 2 * Y * Z,
X * Y ** 2 * Z,
X * Y * Z ** 2]
FOURTH_ORDER_CORRECTION_SYMBOLS = sp.symbols("corr_fourth_:6")
FOURTH_ORDER_RELAXATION_RATE_SYMBOLS = sp.symbols("corr_rr_:10")
def add_fourth_order_correction(collision_rule, shear_relaxation_rate, bulk_relaxation_rate, limiter):
"""Adds the fourth order correction terms (:cite:`geier2017`, eq. 44-49) to a given polynomial D3Q27
cumulant collision rule."""
method = collision_rule.method
if not isinstance(method, CumulantBasedLbMethod) or method.stencil != LBStencil(Stencil.D3Q27):
raise ValueError("Fourth-order correction is only defined for D3Q27 cumulant methods.")
polynomials = method.cumulants
rho = method.zeroth_order_equilibrium_moment_symbol
if not (set(CORRECTED_FOURTH_ORDER_POLYNOMIALS) < set(polynomials)):
raise ValueError("Fourth order correction requires polynomial cumulants "
f"{', '.join(CORRECTED_FOURTH_ORDER_POLYNOMIALS)} to be present")
# Call
# PC1 = (x^2 * y^2 - 2 * x^2 * z^2 + y^2 * z^2),
# PC2 = (x^2 * y^2 + x^2 * z^2 - 2 * y^2 * z^2),
# PC3 = (x^2 * y^2 + x^2 * z^2 + y^2 * z^2),
# PC4 = (x^2 * y * z),
# PC5 = (x * y^2 * z),
# PC6 = (x * y * z^2)
poly_symbols = [sp.Symbol(f'{POST_COLLISION_CUMULANT}_{polynomials.index(poly)}')
for poly in CORRECTED_FOURTH_ORDER_POLYNOMIALS]
a_symb, b_symb = sp.symbols("a_corr, b_corr")
a, b = get_optimal_additional_parameters(shear_relaxation_rate, bulk_relaxation_rate)
correction_terms = get_fourth_order_correction_terms(method.relaxation_rate_dict, rho, a_symb, b_symb)
optimal_parametrisation = get_optimal_parametrisation_with_limiters(collision_rule, shear_relaxation_rate,
bulk_relaxation_rate, limiter)
subexp_dict = collision_rule.subexpressions_dict
subexp_dict = {**subexp_dict,
a_symb: a,
b_symb: b,
**correction_terms.subexpressions_dict,
**optimal_parametrisation.subexpressions_dict,
**correction_terms.main_assignments_dict,
**optimal_parametrisation.main_assignments_dict}
for sym, cor in zip(poly_symbols, FOURTH_ORDER_CORRECTION_SYMBOLS):
subexp_dict[sym] += cor
collision_rule.set_sub_expressions_from_dict(subexp_dict)
collision_rule.topological_sort()
return collision_rule
def get_optimal_additional_parameters(shear_relaxation_rate, bulk_relaxation_rate):
"""
Calculates the optimal parametrization for the additional parameters provided in :cite:`geier2017`
Equations 115-116.
"""
omega_1 = shear_relaxation_rate
omega_2 = bulk_relaxation_rate
omega_11 = omega_1 * omega_1
omega_12 = omega_1 * omega_2
omega_22 = omega_2 * omega_2
two = sp.Float(2)
three = sp.Float(3)
four = sp.Float(4)
nine = sp.Float(9)
denom_ab = (omega_1 - omega_2) * (omega_2 * (two + three * omega_1) - sp.Float(8) * omega_1)
a = (four * omega_11 + two * omega_12 * (omega_1 - sp.Float(6)) + omega_22 * (
omega_1 * (sp.Float(10) - three * omega_1) - four)) / denom_ab
b = (four * omega_12 * (nine * omega_1 - sp.Float(16)) - four * omega_11 - two * omega_22 * (
two + nine * omega_1 * (omega_1 - two))) / (three * denom_ab)
return a, b
def get_fourth_order_correction_terms(rrate_dict, rho, a, b):
pc1, pc2, pc3, pc4, pc5, pc6 = CORRECTED_FOURTH_ORDER_POLYNOMIALS
omega_1 = rrate_dict[X * Y]
omega_2 = rrate_dict[X ** 2 + Y ** 2 + Z ** 2]
try:
omega_6 = rrate_dict[pc1]
assert omega_6 == rrate_dict[pc2], \
"Cumulants (x^2 - y^2) and (x^2 - z^2) must have the same relaxation rate"
omega_7 = rrate_dict[pc3]
omega_8 = rrate_dict[pc4]
assert omega_8 == rrate_dict[pc5] == rrate_dict[pc6]
except IndexError:
raise ValueError("For the fourth order correction, all six polynomial cumulants"
"(x^2 * y^2 - 2 * x^2 * z^2 + y^2 * z^2),"
"(x^2 * y^2 + x^2 * z^2 - 2 * y^2 * z^2),"
"(x^2 * y^2 + x^2 * z^2 + y^2 * z^2),"
"(x^2 * y * z), (x * y^2 * z) and (x * y * z^2) must be present!")
dxu, dyv, dzw, dxvdyu, dxwdzu, dywdzv = sp.symbols('Dx, Dy, Dz, DxvDyu, DxwDzu, DywDzv')
c_xy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 1, 0))
c_xz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 0, 1))
c_yz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 1, 1))
c_xx = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (2, 0, 0))
c_yy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 2, 0))
c_zz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 0, 2))
cor1, cor2, cor3, cor4, cor5, cor6 = FOURTH_ORDER_CORRECTION_SYMBOLS
one = sp.Float(1)
two = sp.Float(2)
three = sp.Float(3)
# Derivative Approximations
subexpressions = [
Assignment(dxu, - omega_1 / (two * rho) * (two * c_xx - c_yy - c_zz)
- omega_2 / (two * rho) * (c_xx + c_yy + c_zz - rho)),
Assignment(dyv, dxu + (three * omega_1) / (two * rho) * (c_xx - c_yy)),
Assignment(dzw, dxu + (three * omega_1) / (two * rho) * (c_xx - c_zz)),
Assignment(dxvdyu, - three * omega_1 / rho * c_xy),
Assignment(dxwdzu, - three * omega_1 / rho * c_xz),
Assignment(dywdzv, - three * omega_1 / rho * c_yz)]
one_half = sp.Rational(1, 2)
one_over_three = sp.Rational(1, 3)
two_over_three = sp.Rational(2, 3)
four_over_three = sp.Rational(4, 3)
# Correction Terms
main_assignments = [
Assignment(cor1, two_over_three * (one / omega_1 - one_half) * omega_6 * a * rho * (dxu - two * dyv + dzw)),
Assignment(cor2, two_over_three * (one / omega_1 - one_half) * omega_6 * a * rho * (dxu + dyv - two * dzw)),
Assignment(cor3, - four_over_three * (one / omega_1 - one_half) * omega_7 * a * rho * (dxu + dyv + dzw)),
Assignment(cor4, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dywdzv),
Assignment(cor5, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dxwdzu),
Assignment(cor6, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dxvdyu)]
return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)
def get_optimal_parametrisation_with_limiters(collision_rule, shear_relaxation_rate, bulk_relaxation_rate, limiter):
"""
Calculates the optimal parametrization for the relaxation rates provided in :cite:`geier2017`
Equations 112-114.
"""
# if omega numbers
# assert omega_1 != omega_2, "Relaxation rates associated with shear and bulk viscosity must not be identical."
# assert omega_1 > 7/4
omega_1 = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[0]
omega_2 = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[1]
non_limited_omegas = sp.symbols("non_limited_omega_3:6")
limited_omegas = sp.symbols("limited_omega_3:6")
omega_11 = omega_1 * omega_1
omega_12 = omega_1 * omega_2
omega_22 = omega_2 * omega_2
one = sp.Float(1)
two = sp.Float(2)
three = sp.Float(3)
five = sp.Float(5)
six = sp.Float(6)
seven = sp.Float(7)
eight = sp.Float(8)
nine = sp.Float(9)
omega_3 = (eight * (omega_1 - two) * (omega_2 * (three * omega_1 - one) - five * omega_1)) / (
eight * (five - two * omega_1) * omega_1 + omega_2 * (eight + omega_1 * (nine * omega_1 - sp.Float(26))))
omega_4 = (eight * (omega_1 - two) * (omega_1 + omega_2 * (three * omega_1 - seven))) / (
omega_2 * (sp.Float(56) - sp.Float(42) * omega_1 + nine * omega_11) - eight * omega_1)
omega_5 = (sp.Float(24) * (omega_1 - two) * (sp.Float(4) * omega_11 + omega_12 * (
sp.Float(18) - sp.Float(13) * omega_1) + omega_22 * (two + omega_1 * (
six * omega_1 - sp.Float(11))))) / (sp.Float(16) * omega_11 * (omega_1 - six) - two * omega_12 * (
sp.Float(216) + five * omega_1 * (nine * omega_1 - sp.Float(46))) + omega_22 * (omega_1 * (
three * omega_1 - sp.Float(10)) * (sp.Float(15) * omega_1 - sp.Float(28)) - sp.Float(48)))
rho = collision_rule.method.zeroth_order_equilibrium_moment_symbol
# add limiters to improve stability
c_xyy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 2, 0))
c_xzz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 0, 2))
c_xyz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 1, 1))
abs_c_xyy_plus_xzz = sp.Abs(c_xyy + c_xzz)
abs_c_xyy_minus_xzz = sp.Abs(c_xyy - c_xzz)
abs_c_xyz = sp.Abs(c_xyz)
limited_omega_3 = non_limited_omegas[0] + ((one - non_limited_omegas[0]) * abs_c_xyy_plus_xzz) / \
(rho * limiter + abs_c_xyy_plus_xzz)
limited_omega_4 = non_limited_omegas[1] + ((one - non_limited_omegas[1]) * abs_c_xyy_minus_xzz) / \
(rho * limiter + abs_c_xyy_minus_xzz)
limited_omega_5 = non_limited_omegas[2] + ((one - non_limited_omegas[2]) * abs_c_xyz) / (rho * limiter + abs_c_xyz)
subexpressions = [
Assignment(non_limited_omegas[0], omega_3),
Assignment(non_limited_omegas[1], omega_4),
Assignment(non_limited_omegas[2], omega_5),
Assignment(limited_omegas[0], limited_omega_3),
Assignment(limited_omegas[1], limited_omega_4),
Assignment(limited_omegas[2], limited_omega_5)]
# Correction Terms
main_assignments = [
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[0], shear_relaxation_rate),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[1], bulk_relaxation_rate),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[2], limited_omegas[0]),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[3], limited_omegas[1]),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[4], limited_omegas[2]),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[5], one),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[6], one),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[7], one),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[8], one),
Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[9], one),
]
return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)
import sympy as sp
from pystencils.simp.assignment_collection import AssignmentCollection
from pystencils import Assignment
from pystencils import Assignment, AssignmentCollection
from lbmpy.stencils import Stencil, LBStencil
from lbmpy.moments import MOMENT_SYMBOLS, statistical_quantity_symbol
......
......@@ -55,7 +55,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
def __init__(self, stencil, equilibrium, relaxation_dict,
conserved_quantity_computation=None,
force_model=None, zero_centered=False,
force_model=None, zero_centered=False, fraction_field=None,
central_moment_transform_class=BinomialChimeraTransform):
assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
super(CentralMomentBasedLbMethod, self).__init__(stencil)
......@@ -65,6 +65,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
self._cqc = conserved_quantity_computation
self._force_model = force_model
self._zero_centered = zero_centered
self._fraction_field = fraction_field
self._weights = None
self._central_moment_transform_class = central_moment_transform_class
......@@ -73,6 +74,10 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
"""Force model employed by this method."""
return self._force_model
@property
def fraction_field(self):
return self._fraction_field
@property
def relaxation_info_dict(self):
"""Dictionary mapping this method's moments to their relaxation rates and equilibrium values.
......
......@@ -286,7 +286,7 @@ def _get_relaxation_rates(collision_rule):
omega_s = get_shear_relaxation_rate(method)
# if the shear relaxation rate is not specified as a symbol look for its symbolic counter part in the subs dict
for symbolic_rr, rr in method.subs_dict_relxation_rate.items():
for symbolic_rr, rr in method.subs_dict_relaxation_rate.items():
if omega_s == rr:
omega_s = symbolic_rr
......
......@@ -16,8 +16,8 @@ from lbmpy.moment_transforms import PdfsToMomentsByChimeraTransform
class MomentBasedLbMethod(AbstractLbMethod):
"""
Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
These methods work by transforming the pdfs into moment space using a linear transformation. In the moment
space each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
These methods work by transforming the pdfs into moment space using a linear transformation. In the moment-space
each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
Parameters:
......@@ -39,7 +39,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
def __init__(self, stencil, equilibrium, relaxation_dict,
conserved_quantity_computation=None, force_model=None, zero_centered=False,
moment_transform_class=PdfsToMomentsByChimeraTransform):
fraction_field=None, moment_transform_class=PdfsToMomentsByChimeraTransform):
assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
super(MomentBasedLbMethod, self).__init__(stencil)
......@@ -48,6 +48,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
self._cqc = conserved_quantity_computation
self._force_model = force_model
self._zero_centered = zero_centered
self._fraction_field = fraction_field
self._weights = None
self._moment_transform_class = moment_transform_class
......@@ -56,6 +57,10 @@ class MomentBasedLbMethod(AbstractLbMethod):
"""Force model employed by this method."""
return self._force_model
@property
def fraction_field(self):
return self._fraction_field
@property
def relaxation_info_dict(self):
"""Dictionary mapping this method's moments to their relaxation rates and equilibrium values.
......@@ -174,7 +179,14 @@ class MomentBasedLbMethod(AbstractLbMethod):
def get_collision_rule(self, conserved_quantity_equations: AssignmentCollection = None,
pre_simplification: bool = True) -> LbmCollisionRule:
rr_sub_expressions, d = self._generate_symbolic_relaxation_matrix()
if self._fraction_field is not None:
relaxation_rates_modifier = (1.0 - self._fraction_field)
rr_sub_expressions, d = self._generate_symbolic_relaxation_matrix(
relaxation_rates_modifier=relaxation_rates_modifier)
else:
rr_sub_expressions, d = self._generate_symbolic_relaxation_matrix()
ac = self._collision_rule_with_relaxation_matrix(d=d,
additional_subexpressions=rr_sub_expressions,
include_force_terms=True,
......@@ -378,7 +390,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
subexpressions += m_post_to_f_post_eqs.subexpressions
main_assignments = m_post_to_f_post_eqs.main_assignments
else:
# For SRT, TRT by default, and whenever customly required, derive equations entirely in
# For SRT, TRT by default, and whenever customarily required, derive equations entirely in
# population space
if self._zero_centered and not self._equilibrium.deviation_only:
raise Exception("Can only derive population-space equations for zero-centered storage"
......