diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 479504565cf07220883e559038680db553169d44..79050dbbd39685496c60e6d82436ab61e0c0f6a4 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -52,6 +52,8 @@ For example, to modify the AST one can run::
func = create_lb_function(ast=ast, ...)
"""
+import copy
+
from dataclasses import dataclass, field, replace
from typing import Union, List, Tuple, Any, Type, Iterable
from warnings import warn, filterwarnings
@@ -65,6 +67,7 @@ from lbmpy.enums import Stencil, Method, ForceModel, CollisionSpace
import lbmpy.forcemodels as forcemodels
from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor
from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule
+from lbmpy.partially_saturated_cells import add_psm_to_collision_rule, PSMConfig
from lbmpy.non_newtonian_models import add_cassons_model, CassonsParameters
from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_central_moment,
create_srt, create_trt, create_trt_kbc)
@@ -79,7 +82,9 @@ from lbmpy.turbulence_models import add_smagorinsky_model
from lbmpy.updatekernels import create_lbm_kernel, create_stream_pull_with_output_kernel
from lbmpy.advanced_streaming.utility import Timestep, get_accessor
from pystencils import CreateKernelConfig, create_kernel
+from pystencils.astnodes import Conditional, Block
from pystencils.cache import disk_cache_no_fallback
+from pystencils.node_collection import NodeCollection
from pystencils.typing import collate_types
from pystencils.field import Field
from pystencils.simp import sympy_cse, SimplificationStrategy
@@ -257,6 +262,12 @@ class LBMConfig:
Temperature for fluctuating lattice Boltzmann methods.
"""
+ psm_config: PSMConfig = None
+ """
+ If a PSM config is specified, (1 - fractionField) is added to the relaxation rates of the collision
+ and to the potential force term, and a solid collision is build and added to the main assignments.
+ """
+
output: dict = field(default_factory=dict)
"""
A dictionary mapping macroscopic quantites e.g. the strings 'density' and 'velocity' to pystencils
@@ -446,6 +457,9 @@ class LBMConfig:
'centralmoment': forcemodels.CentralMoment
}
+ if self.psm_config is not None and self.psm_config.fraction_field is not None:
+ self.force = [(1.0 - self.psm_config.fraction_field.center) * f for f in self.force]
+
if isinstance(self.force_model, str):
new_force_model = ForceModel[self.force_model.upper()]
warn(f'ForceModel "{self.force_model}" as str is deprecated. Use {new_force_model} instead or '
@@ -657,6 +671,11 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
else:
collision_rule = lb_method.get_collision_rule(pre_simplification=pre_simplification)
+ if lbm_config.psm_config is not None:
+ if lbm_config.psm_config.fraction_field is None or lbm_config.psm_config.object_velocity_field is None:
+ raise ValueError("Specify a fraction and object velocity field in the PSM Config")
+ collision_rule = add_psm_to_collision_rule(collision_rule, lbm_config.psm_config)
+
if lbm_config.galilean_correction:
from lbmpy.methods.cumulantbased import add_galilean_correction
collision_rule = add_galilean_correction(collision_rule)
@@ -686,6 +705,7 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
omega_output_field=lbm_config.omega_output_field)
else:
collision_rule = add_entropy_condition(collision_rule, omega_output_field=lbm_config.omega_output_field)
+
elif lbm_config.smagorinsky:
if lbm_config.cassons:
raise ValueError("Cassons model can not be combined with Smagorinsky model")
@@ -738,6 +758,11 @@ def create_lb_method(lbm_config=None, **params):
if isinstance(lbm_config.force, Field):
lbm_config.force = tuple(lbm_config.force(i) for i in range(dim))
+ if lbm_config.psm_config is None:
+ fraction_field = None
+ else:
+ fraction_field = lbm_config.psm_config.fraction_field
+
common_params = {
'compressible': lbm_config.compressible,
'zero_centered': lbm_config.zero_centered,
@@ -747,6 +772,7 @@ def create_lb_method(lbm_config=None, **params):
'continuous_equilibrium': lbm_config.continuous_equilibrium,
'c_s_sq': lbm_config.c_s_sq,
'collision_space_info': lbm_config.collision_space_info,
+ 'fraction_field': fraction_field,
}
cumulant_params = {
@@ -754,6 +780,7 @@ def create_lb_method(lbm_config=None, **params):
'force_model': lbm_config.force_model,
'c_s_sq': lbm_config.c_s_sq,
'collision_space_info': lbm_config.collision_space_info,
+ 'fraction_field': fraction_field,
}
if lbm_config.method == Method.SRT:
@@ -815,6 +842,54 @@ def create_lb_method(lbm_config=None, **params):
return method
+def create_psm_update_rule(lbm_config, lbm_optimisation):
+ node_collection = []
+
+ # Use regular lb update rule for no overlapping particles
+ config_without_psm = copy.deepcopy(lbm_config)
+ config_without_psm.psm_config = None
+ # TODO: the force is still multiplied by (1.0 - self.psm_config.fraction_field.center)
+ # (should not harm if memory bound since self.psm_config.fraction_field.center should always be 0.0)
+ lb_update_rule = create_lb_update_rule(
+ lbm_config=config_without_psm, lbm_optimisation=lbm_optimisation
+ )
+ node_collection.append(
+ Conditional(
+ lbm_config.psm_config.fraction_field.center(0) <= 0.0,
+ Block(lb_update_rule.all_assignments),
+ )
+ )
+
+ # Only one particle, i.e., no individual_fraction_field is provided
+ if lbm_config.psm_config.individual_fraction_field is None:
+ assert lbm_config.psm_config.MaxParticlesPerCell == 1
+ psm_update_rule = create_lb_update_rule(
+ lbm_config=lbm_config, lbm_optimisation=lbm_optimisation
+ )
+ node_collection.append(
+ Conditional(
+ lbm_config.psm_config.fraction_field.center(0) > 0.0,
+ Block(psm_update_rule.all_assignments),
+ )
+ )
+ else:
+ for p in range(lbm_config.psm_config.MaxParticlesPerCell):
+ # Add psm update rule for p overlapping particles
+ config_with_p_particles = copy.deepcopy(lbm_config)
+ config_with_p_particles.psm_config.MaxParticlesPerCell = p + 1
+ psm_update_rule = create_lb_update_rule(
+ lbm_config=config_with_p_particles, lbm_optimisation=lbm_optimisation
+ )
+ node_collection.append(
+ Conditional(
+ lbm_config.psm_config.individual_fraction_field.center(p) > 0.0,
+ Block(psm_update_rule.all_assignments),
+ )
+ )
+
+ return NodeCollection(node_collection)
+
+
# ----------------------------------------------------------------------------------------------------------------------
def update_with_default_parameters(params, opt_params=None, lbm_config=None, lbm_optimisation=None, config=None):
# Fix CreateKernelConfig params
diff --git a/lbmpy/methods/abstractlbmethod.py b/lbmpy/methods/abstractlbmethod.py
index 096a3dd6f264a0a15b12b813a24a251d2abd1514..7bc4f5c920e0c8e3474a052a8842c1a70a040042 100644
--- a/lbmpy/methods/abstractlbmethod.py
+++ b/lbmpy/methods/abstractlbmethod.py
@@ -99,7 +99,7 @@ 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
@@ -120,5 +120,7 @@ class AbstractLbMethod(abc.ABC):
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:
+ new_rr = [r * relaxation_rates_modifier for r in new_rr]
return substitutions, sp.diag(*new_rr)
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index 78f35192b95d1fe718b02d75b830f40be6c82f65..e7440bd82a81b68c457db4b298ed404e9728c891 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -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.
@@ -201,11 +201,11 @@ 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,
+ 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,
+ 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,
@@ -333,7 +333,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, fraction_field=None, **kwargs):
r"""
Creates moment based LB method where the collision takes place in the central moment space.
@@ -368,6 +368,10 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
nested_moments = cascaded_moment_sets_literature(stencil)
rr_dict = _get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D)
+ if fraction_field is not None:
+ relaxation_rates_modifier = (1.0 - fraction_field.center)
+ 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)
@@ -519,7 +523,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, fraction_field=None, **kwargs):
r"""Creates a cumulant-based lattice Boltzmann method.
Args:
@@ -537,6 +541,12 @@ def create_cumulant(stencil, relaxation_rates, cumulant_groups, **kwargs):
:class:`lbmpy.methods.cumulantbased.CumulantBasedLbMethod` instance
"""
cumulant_to_rr_dict = _get_relaxation_info_dict(relaxation_rates, cumulant_groups, stencil.D)
+
+ if fraction_field is not None:
+ relaxation_rates_modifier = (1.0 - fraction_field.center)
+ 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:
@@ -581,7 +591,7 @@ def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs
return create_cumulant(stencil, relaxation_rates, cumulant_groups, **kwargs)
-def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
+def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim, relaxation_rates_modifier=None):
r"""Creates a dictionary where each moment is mapped to a relaxation rate.
Args:
@@ -599,7 +609,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
@@ -664,6 +676,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
diff --git a/lbmpy/methods/momentbased/momentbasedmethod.py b/lbmpy/methods/momentbased/momentbasedmethod.py
index 79aee5c813b71c358e93964b52b84b6b1532269f..740834072a0210772d407a7e441b9e1ef46d0f7a 100644
--- a/lbmpy/methods/momentbased/momentbasedmethod.py
+++ b/lbmpy/methods/momentbased/momentbasedmethod.py
@@ -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
@@ -174,7 +175,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.center)
+ 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,
diff --git a/lbmpy/partially_saturated_cells.py b/lbmpy/partially_saturated_cells.py
new file mode 100644
index 0000000000000000000000000000000000000000..0798ac69501b46ffc22251ff13e7f5ecc036d3b9
--- /dev/null
+++ b/lbmpy/partially_saturated_cells.py
@@ -0,0 +1,115 @@
+import sympy as sp
+from dataclasses import dataclass
+
+from lbmpy.methods.abstractlbmethod import LbmCollisionRule
+from pystencils import Assignment, AssignmentCollection
+from pystencils.field import Field
+
+
+@dataclass
+class PSMConfig:
+ fraction_field: Field = None
+ """
+ Fraction field for PSM
+ """
+
+ object_velocity_field: Field = None
+ """
+ Object velocity field for PSM
+ """
+
+ SC: int = 1
+ """
+ Solid collision option for PSM
+ """
+
+ MaxParticlesPerCell: int = 1
+ """
+ Maximum number of particles overlapping with a cell
+ """
+
+ individual_fraction_field: Field = None
+ """
+ Fraction field for each overlapping particle in PSM
+ """
+
+ particle_force_field: Field = None
+ """
+ Force field for each overlapping particle in PSM
+ """
+
+
+def add_psm_to_collision_rule(collision_rule, psm_config):
+ if psm_config.individual_fraction_field is None:
+ psm_config.individual_fraction_field = psm_config.fraction_field
+
+ method = collision_rule.method
+ pre_collision_pdf_symbols = method.pre_collision_pdf_symbols
+ stencil = method.stencil
+
+ # Get equilibrium from object velocity for solid collision
+ forces_rhs = [0] * psm_config.MaxParticlesPerCell * stencil.D
+ solid_collisions = [0] * stencil.Q
+ for p in range(psm_config.MaxParticlesPerCell):
+ equilibrium_fluid = method.get_equilibrium_terms()
+ equilibrium_solid = []
+ for eq in equilibrium_fluid:
+ eq_sol = eq
+ for i in range(stencil.D):
+ eq_sol = eq_sol.subs(sp.Symbol("u_" + str(i)),
+ psm_config.object_velocity_field.center(p * stencil.D + i), )
+ equilibrium_solid.append(eq_sol)
+
+ # Build solid collision
+ for i, (eqFluid, eqSolid, f, offset) in enumerate(
+ zip(equilibrium_fluid, equilibrium_solid, pre_collision_pdf_symbols, stencil)):
+ inverse_direction_index = stencil.stencil_entries.index(stencil.inverse_stencil_entries[i])
+ if psm_config.SC == 1:
+ solid_collision = psm_config.individual_fraction_field.center(p) * (
+ (
+ pre_collision_pdf_symbols[inverse_direction_index]
+ - equilibrium_fluid[inverse_direction_index]
+ )
+ - (f - eqSolid)
+ )
+ elif psm_config.SC == 2:
+ # TODO get relaxation rate vector from method and use the right relaxation rate [i]
+ solid_collision = psm_config.individual_fraction_field.center(p) * (
+ (eqSolid - f) + (1.0 - method.relaxation_rates[0]) * (f - eqFluid)
+ )
+ elif psm_config.SC == 3:
+ solid_collision = psm_config.individual_fraction_field.center(p) * (
+ (
+ pre_collision_pdf_symbols[inverse_direction_index]
+ - equilibrium_solid[inverse_direction_index]
+ )
+ - (f - eqSolid)
+ )
+ else:
+ raise ValueError("Only SC=1, SC=2 and SC=3 are supported.")
+ solid_collisions[i] += solid_collision
+ for j in range(stencil.D):
+ forces_rhs[p * stencil.D + j] -= solid_collision * int(offset[j])
+
+ # Add solid collision to main assignments of collision rule
+ collision_assignments = []
+ for main, sc in zip(collision_rule.main_assignments, solid_collisions):
+ collision_assignments.append(Assignment(main.lhs, main.rhs + sc))
+
+ # Add hydrodynamic force calculations to collision assignments if two-way coupling is used
+ # (i.e., force field is not None)
+ if psm_config.particle_force_field is not None:
+ for p in range(psm_config.MaxParticlesPerCell):
+ for i in range(stencil.D):
+ collision_assignments.append(
+ Assignment(
+ psm_config.particle_force_field.center(p * stencil.D + i),
+ forces_rhs[p * stencil.D + i],
+ )
+ )
+
+ collision_assignments = AssignmentCollection(collision_assignments)
+ ac = LbmCollisionRule(method, collision_assignments, collision_rule.subexpressions,
+ collision_rule.simplification_hints)
+ ac.topological_sort()
+ return ac