Implemented zero center logic

lbmpy/creationfunctions.py

@@ -412,6 +412,7 @@ def create_lb_method(lbm_config=None, **params):
 
     common_params = {
         'compressible': lbm_config.compressible,
+        'zero_centered_pdfs': lbm_config.zero_centered_pdfs,
         'equilibrium_order': lbm_config.equilibrium_order,
         'force_model': lbm_config.force_model,
         'maxwellian_moments': lbm_config.maxwellian_moments,
@@ -421,6 +422,7 @@ def create_lb_method(lbm_config=None, **params):
     }
 
     cumulant_params = {
+        'zero_centered_pdfs': lbm_config.zero_centered_pdfs,
         'equilibrium_order': lbm_config.equilibrium_order,
         'force_model': lbm_config.force_model,
         'c_s_sq': lbm_config.c_s_sq,
@@ -579,6 +581,8 @@ class LBMConfig:
         relaxation_rates: List of relexation rates for the method. Can contain symbols, ints and floats
         relaxation_rate: For TRT and KBC methods the remaining relaxation rate is determined via magic parameter
         compressible: If Helmholtz decomposition is used to better approximate the incompressible Navier Stokes eq
+        zero_centered_pdfs: If True the density is subtracted from the PDFs to center them around zero.
+                            This is important for simulations in single precision
         equilibrium_order: order of the highest term in the equilibrium formulation
         c_s_sq: speed of sound. Usually 1 / 3
         weighted: For MRT methods. If True the moments are weigted orthogonal. Usually desired.
@@ -618,6 +622,7 @@ class LBMConfig:
     relaxation_rates: Any = None
     relaxation_rate: Union[int, float, Type[sp.Symbol]] = None
     compressible: bool = False
+    zero_centered_pdfs: bool = True
     equilibrium_order: int = 2
     c_s_sq: sympy.core.numbers.Rational = sp.Rational(1, 3)
     weighted: bool = True

lbmpy/macroscopic_value_kernels.py

@@ -67,7 +67,9 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
         iteration_slice: if not None, iteration is done only over this slice of the field
         field_layout: layout for output field, also used for pdf field if pdf_arr is not given
         target: `Target.CPU` or `Target.GPU`
-        previous_step_accessor: The accessor used by the streaming pattern of the previous timestep
+        streaming_pattern: define the streaming pattern for the getter kernel. Default is pull
+        previous_timestep: timestep which is used. Timestep.BOTH for two population methods.
+                           Timestep.EVEN to use the even order timestep, Timestep.ODD to use the odd order timestep
 
     Returns:
         a function to compute macroscopic values:
@@ -80,8 +82,8 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
     cqc = lb_method.conserved_quantity_computation
     unknown_quantities = [oq for oq in output_quantities if oq not in cqc.conserved_quantities]
     if unknown_quantities:
-        raise ValueError("No such conserved quantity: %s, conserved quantities are %s" %
-                         (str(unknown_quantities), str(cqc.conserved_quantities.keys())))
+        raise ValueError(f"No such conserved quantity: {str(unknown_quantities)}, "
+                         f"conserved quantities are {str(cqc.conserved_quantities.keys())}")
 
     if pdf_arr is None:
         pdf_field = Field.create_generic('pdfs', lb_method.dim, index_dimensions=1, layout=field_layout)
@@ -126,15 +128,15 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
         kernel = gpu.make_python_function(gpu.create_cuda_kernel(
             eqs, ghost_layers=ghost_layers, iteration_slice=iteration_slice))
     else:
-        raise ValueError("Unknown target '%s'. Possible targets are `Target.CPU` and `Target.GPU`" % (target,))
+        raise ValueError(f"Unknown target '{target}'. Possible targets are `Target.CPU` and `Target.GPU`")
 
     def getter(pdfs, **kwargs):
         if pdf_arr is not None:
             assert pdfs.shape == pdf_arr.shape and pdfs.strides == pdf_arr.strides, \
                 "Pdf array not matching blueprint which was used to compile" + str(pdfs.shape) + str(pdf_arr.shape)
         if not set(output_quantities).issubset(kwargs.keys()):
-            raise ValueError("You have to specify the output field for each of the following quantities: %s"
-                             % (str(output_quantities),))
+            raise ValueError(f"You have to specify the output field for "
+                             f"each of the following quantities: {str(output_quantities)}")
         kernel(pdfs=pdfs, **kwargs)
 
     return getter
@@ -158,7 +160,9 @@ def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None
         iteration_slice: if not None, iteration is done only over this slice of the field
         field_layout: layout of the pdf field if pdf_arr was not given
         target: `Target.CPU` or `Target.GPU`
-        previous_step_accessor: The accessor used by the streaming pattern of the previous timestep
+        streaming_pattern: define the streaming pattern for the getter kernel. Default is pull
+        previous_timestep: timestep which is used. Timestep.BOTH for two population methods.
+                           Timestep.EVEN to use the even order timestep, Timestep.ODD to use the odd order timestep
 
     Returns:
         function taking pdf array as single argument and which sets the field to the given values
@@ -211,12 +215,12 @@ def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None
         kernel = gpu.make_python_function(gpu.create_cuda_kernel(eq))
         kernel = functools.partial(kernel, **fixed_kernel_parameters)
     else:
-        raise ValueError("Unknown target '%s'. Possible targets are `Target.CPU` and `Target.GPU`" % (target,))
+        raise ValueError(f"Unknown target '{target}'. Possible targets are `Target.CPU` and `Target.GPU`")
 
     def setter(pdfs, **kwargs):
         if pdf_arr is not None:
             assert pdfs.shape == pdf_arr.shape and pdfs.strides == pdf_arr.strides, \
-                "Pdf array not matching blueprint which was used to compile" + str(pdfs.shape) + str(pdf_arr.shape)
+                f"Pdf array not matching blueprint which was used to compile {str(pdfs.shape)} {str(pdfs.strides)}"
         kernel(pdfs=pdfs, **kwargs)
 
     return setter

lbmpy/methods/centeredcumulant/centeredcumulantmethod.py

@@ -177,7 +177,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
 
     @property
     def central_moment_transform_class(self):
-        self._central_moment_transform_class
+        return self._central_moment_transform_class
 
     @property
     def cumulants(self):
@@ -189,7 +189,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
 
     @property
     def cumulant_transform_class(self):
-        self._cumulant_transform_class
+        return self._cumulant_transform_class
 
     @property
     def first_order_equilibrium_moment_symbols(self, ):
@@ -211,6 +211,14 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
     def relaxation_rates(self):
         return tuple([e.relaxation_rate for e in self._cumulant_to_relaxation_info_dict.values()])
 
+    @property
+    def compressible(self):
+        return self._conserved_quantity_computation.compressible
+
+    @property
+    def zero_centered_pdfs(self):
+        return self._conserved_quantity_computation.zero_centered_pdfs
+
     @property
     def zeroth_order_equilibrium_moment_symbol(self, ):
         return self._conserved_quantity_computation.zeroth_order_moment_symbol

lbmpy/methods/conservedquantitycomputation.py

@@ -82,13 +82,14 @@ class AbstractConservedQuantityComputation(abc.ABC):
 
 
 class DensityVelocityComputation(AbstractConservedQuantityComputation):
-    def __init__(self, stencil, compressible, force_model=None,
+    def __init__(self, stencil, compressible, zero_centered_pdfs,  force_model=None,
                  zeroth_order_moment_symbol=sp.Symbol("rho"),
                  first_order_moment_symbols=sp.symbols("u_:3"),
                  second_order_moment_symbols=sp.symbols("p_:9")):
         dim = stencil.D
         self._stencil = stencil
         self._compressible = compressible
+        self._zero_centered_pdfs = zero_centered_pdfs
         self._forceModel = force_model
         self._symbolOrder0 = zeroth_order_moment_symbol
         self._symbolsOrder1 = first_order_moment_symbols[:dim]
@@ -103,6 +104,10 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
     def compressible(self):
         return self._compressible
 
+    @property
+    def zero_centered_pdfs(self):
+        return self._zero_centered_pdfs
+
     def defined_symbols(self, order='all'):
         if order == 'all':
             return {'density': self._symbolOrder0,
@@ -114,10 +119,6 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
         else:
             return None
 
-    @property
-    def zero_centered_pdfs(self):
-        return not self._compressible
-
     @property
     def zeroth_order_moment_symbol(self):
         return self._symbolOrder0
@@ -149,7 +150,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
 
     def equilibrium_input_equations_from_init_values(self, density=1, velocity=(0, 0, 0), force_substitution=True):
         dim = self._stencil.D
-        zeroth_order_moment = density if self._compressible else density - sp.Rational(1, 1)
+        zeroth_order_moment = density - sp.Rational(1, 1) if self._zero_centered_pdfs else density
         first_order_moments = velocity[:dim]
         vel_offset = [0] * dim
 
@@ -180,7 +181,8 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
 
         if self._compressible:
             ac = divide_first_order_moments_by_rho(ac, dim)
-        else:
+
+        if self._zero_centered_pdfs:
             ac = add_density_offset(ac)
 
         if self._forceModel is not None:

lbmpy/methods/creationfunctions.py

@@ -27,12 +27,13 @@ from lbmpy.moments import (
     is_bulk_moment, is_shear_moment, get_order, set_up_shift_matrix)
 
 from lbmpy.relaxationrates import relaxation_rate_from_magic_number
-from lbmpy.stencils import get_stencil
+from lbmpy.stencils import LBStencil
 from pystencils.stencil import have_same_entries
 from pystencils.sympyextensions import common_denominator
 
 
-def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, compressible=False,
+def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict,
+                                               compressible=False, zero_centered_pdfs=True,
                                                force_model=None, equilibrium_order=2, c_s_sq=sp.Rational(1, 3),
                                                central_moment_space=False,
                                                moment_transform_class=None,
@@ -42,7 +43,7 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
     These moments are relaxed against the moments of the discrete Maxwellian distribution.
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See `LBStencil`
         moment_to_relaxation_rate_dict: dict that has as many entries as the stencil. Each moment, which can be
                                         represented by an exponent tuple or in polynomial form
                                         (see `lbmpy.moments`), is mapped to a relaxation rate.
@@ -50,6 +51,8 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
                       where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
                       This approximates the incompressible Navier-Stokes equations better than the standard
                       compressible model.
+        zero_centered_pdfs: If True the density is subtracted from the PDFs to center them around zero.
+                            This is important for simulations in single precision
         force_model: force model instance, or None if no external forces
         equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
         c_s_sq: Speed of sound squared
@@ -64,12 +67,13 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
         :class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` 
     """
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
     mom_to_rr_dict = OrderedDict(moment_to_relaxation_rate_dict)
     assert len(mom_to_rr_dict) == len(stencil), \
         "The number of moments has to be the same as the number of stencil entries"
 
-    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
+    density_velocity_computation = DensityVelocityComputation(stencil, compressible,
+                                                              zero_centered_pdfs, force_model)
 
     moments = tuple(mom_to_rr_dict.keys())
     eq_values = get_moments_of_discrete_maxwellian_equilibrium(stencil, moments,
@@ -89,7 +93,8 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
         return MomentBasedLbMethod(stencil, rr_dict, density_velocity_computation, force_model, moment_transform_class)
 
 
-def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, compressible=False,
+def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict,
+                                                 compressible=False, zero_centered_pdfs=True,
                                                  force_model=None, equilibrium_order=2, c_s_sq=sp.Rational(1, 3),
                                                  central_moment_space=False,
                                                  moment_transform_class=None,
@@ -101,7 +106,7 @@ def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_r
     By using the continuous Maxwellian we automatically get a compressible model.
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See `LBStencil`
         moment_to_relaxation_rate_dict: dict that has as many entries as the stencil. Each moment, which can be
                                         represented by an exponent tuple or in polynomial form
                                         (see `lbmpy.moments`), is mapped to a relaxation rate.
@@ -109,6 +114,8 @@ def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_r
                       where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
                       This approximates the incompressible Navier-Stokes equations better than the standard
                       compressible model.
+        zero_centered_pdfs: If True the density is subtracted from the PDFs to center them around zero.
+                            This is important for simulations in single precision
         force_model: force model instance, or None if no external forces
         equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
         c_s_sq: Speed of sound squared
@@ -123,11 +130,12 @@ def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_r
         :class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` 
     """
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
     mom_to_rr_dict = OrderedDict(moment_to_relaxation_rate_dict)
     assert len(mom_to_rr_dict) == stencil.Q, "The number of moments has to be equal to the number of stencil entries"
     dim = stencil.D
-    density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
+    density_velocity_computation = DensityVelocityComputation(stencil, compressible,
+                                                              zero_centered_pdfs, force_model)
     moments = tuple(mom_to_rr_dict.keys())
 
     if compressible and central_moment_space:
@@ -191,7 +199,7 @@ def create_from_equilibrium(stencil, equilibrium, moment_to_relaxation_rate_dict
         force_model: see create_with_discrete_maxwellian_eq_moments
     """
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
     if any(isinstance(moment_to_relaxation_rate_dict, t) for t in (sp.Symbol, float, int)):
         moment_to_relaxation_rate_dict = {m: moment_to_relaxation_rate_dict
                                           for m in get_default_moment_set_for_stencil(stencil)}
@@ -212,7 +220,7 @@ def create_srt(stencil, relaxation_rate, maxwellian_moments=False, **kwargs):
     r"""Creates a single relaxation time (SRT) lattice Boltzmann model also known as BGK model.
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
         relaxation_rate: relaxation rate (inverse of the relaxation time)
                         usually called :math:`\omega` in LBM literature
         maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
@@ -222,7 +230,7 @@ def create_srt(stencil, relaxation_rate, maxwellian_moments=False, **kwargs):
         :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
     """
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
     moments = get_default_moment_set_for_stencil(stencil)
     rr_dict = OrderedDict([(m, relaxation_rate) for m in moments])
     if maxwellian_moments:
@@ -243,7 +251,7 @@ def create_trt(stencil, relaxation_rate_even_moments, relaxation_rate_odd_moment
     If unsure how to choose the odd relaxation rate, use the function :func:`lbmpy.methods.create_trt_with_magic_number`
     """
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
     moments = get_default_moment_set_for_stencil(stencil)
     rr_dict = OrderedDict([(m, relaxation_rate_even_moments if is_even(m) else relaxation_rate_odd_moments)
                            for m in moments])
@@ -260,7 +268,7 @@ def create_trt_with_magic_number(stencil, relaxation_rate, magic_number=sp.Ratio
     For possible parameters see :func:`lbmpy.methods.create_trt`
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
         relaxation_rate: relaxation rate (inverse of the relaxation time)
                         usually called :math:`\omega` in LBM literature
         magic_number: magic number which is used to calculate the relxation rate for the odd moments.
@@ -278,7 +286,7 @@ def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs
     Creates a MRT method using non-orthogonalized moments.
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
         relaxation_rates: relaxation rates (inverse of the relaxation times) for each moment
         maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
                         used to compute the equilibrium moments.
@@ -286,7 +294,7 @@ def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs
         :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
     """
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
     moments = get_default_moment_set_for_stencil(stencil)
     rr_dict = OrderedDict(zip(moments, relaxation_rates))
     if maxwellian_moments:
@@ -301,7 +309,7 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
     Creates moment based LB method where the collision takes place in the central moment space.
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
         relaxation_rates: relaxation rates (inverse of the relaxation times) for each moment
         maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
                         used to compute the equilibrium moments.
@@ -309,7 +317,7 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
         :class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` instance
     """
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
 
     dim = len(stencil)
 
@@ -415,7 +423,7 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
                        [omega_s] * len(shear_part) + \
                        [omega_h] * len(rest_part)
 
-    stencil = get_stencil("D2Q9") if dim == 2 else get_stencil("D3Q27")
+    stencil = LBStencil("D2Q9") if dim == 2 else LBStencil("D3Q27")
     all_moments = rho + velocity + shear_part + rest_part
     moment_to_rr = OrderedDict((m, rr) for m, rr in zip(all_moments, relaxation_rates))
 
@@ -434,7 +442,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
     To create a generic MRT method use `create_with_discrete_maxwellian_eq_moments`
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See `LBStencil`
         relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
         maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
                                                used to compute the equilibrium moments
@@ -445,7 +453,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
     """
     dim = len(stencil[0])
     if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
+        stencil = LBStencil(stencil)
 
     if weighted is None and not nested_moments:
         raise ValueError("Please specify whether you want weighted or unweighted orthogonality using 'weighted='")
@@ -493,7 +501,7 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted):
     This is for commonly used MRT models found in literature.
 
     Args:
-        stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+        stencil: nested tuple defining the discrete velocity space. See `LBStencil`
         is_weighted: whether to use weighted or unweighted orthogonality
 
     MRT schemes as described in the following references are used
@@ -501,7 +509,7 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted):
     x, y, z = MOMENT_SYMBOLS
     one = sp.Rational(1, 1)
 
-    if have_same_entries(stencil, get_stencil("D2Q9")) and is_weighted:
+    if have_same_entries(stencil, LBStencil("D2Q9")) and is_weighted:
         # Reference:
         # Duenweg, B., Schiller, U. D., & Ladd, A. J. (2007). Statistical mechanics of the fluctuating