lbmpy/methods/momentbased/momentbasedmethod.py

-from collections import OrderedDict
-
-import sympy as sp
-
-from lbmpy.maxwellian_equilibrium import get_weights
-from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
-from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
-from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix
-from pystencils import Assignment
-from pystencils.sympyextensions import subs_additive
-
-
-class MomentBasedLbMethod(AbstractLbMethod):
-    def __init__(self, stencil, moment_to_relaxation_info_dict, conserved_quantity_computation=None, force_model=None):
-        """
-        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
-        equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
-
-        Args:
-            stencil: see :func:`lbmpy.stencils.get_stencil`
-            moment_to_relaxation_info_dict: a dictionary mapping moments in either tuple or polynomial formulation
-                                            to a RelaxationInfo, which consists of the corresponding equilibrium moment
-                                            and a relaxation rate
-            conserved_quantity_computation: instance of :class:`lbmpy.methods.AbstractConservedQuantityComputation`.
-                                            This determines how conserved quantities are computed, and defines
-                                            the symbols used in the equilibrium moments like e.g. density and velocity
-            force_model: force model instance, or None if no forcing terms are required
-        """
-        assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
-        super(MomentBasedLbMethod, self).__init__(stencil)
-
-        self._forceModel = force_model
-        self._momentToRelaxationInfoDict = OrderedDict(moment_to_relaxation_info_dict.items())
-        self._conservedQuantityComputation = conserved_quantity_computation
-        self._weights = None
-
-    @property
-    def force_model(self):
-        return self._forceModel
-
-    @property
-    def relaxation_info_dict(self):
-        return self._momentToRelaxationInfoDict
-
-    @property
-    def conserved_quantity_computation(self):
-        return self._conservedQuantityComputation
-
-    @property
-    def moments(self):
-        return tuple(self._momentToRelaxationInfoDict.keys())
-
-    @property
-    def moment_equilibrium_values(self):
-        return tuple([e.equilibrium_value for e in self._momentToRelaxationInfoDict.values()])
-
-    @property
-    def relaxation_rates(self):
-        return tuple([e.relaxation_rate for e in self._momentToRelaxationInfoDict.values()])
-
-    @property
-    def zeroth_order_equilibrium_moment_symbol(self, ):
-        return self._conservedQuantityComputation.zeroth_order_moment_symbol
-
-    @property
-    def first_order_equilibrium_moment_symbols(self, ):
-        return self._conservedQuantityComputation.first_order_moment_symbols
-
-    @property
-    def weights(self):
-        if self._weights is None:
-            self._weights = self._compute_weights()
-        return self._weights
-
-    def override_weights(self, weights):
-        assert len(weights) == len(self.stencil)
-        self._weights = weights
-
-    def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False):
-        relaxation_matrix = sp.eye(len(self.relaxation_rates))
-        return self._collision_rule_with_relaxation_matrix(relaxation_matrix,
-                                                           conserved_quantity_equations=conserved_quantity_equations,
-                                                           include_force_terms=include_force_terms)
-
-    def get_equilibrium_terms(self):
-        equilibrium = self.get_equilibrium()
-        return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
-
-    def get_collision_rule(self, conserved_quantity_equations=None, pre_simplification=True):
-        d = self.relaxation_matrix
-        relaxation_rate_sub_expressions, d = self._generate_relaxation_matrix(d, pre_simplification)
-        ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
-                                                         True, conserved_quantity_equations)
-        return ac
-
-    def set_zeroth_moment_relaxation_rate(self, relaxation_rate):
-        e = sp.Rational(1, 1)
-        prev_entry = self._momentToRelaxationInfoDict[e]
-        new_entry = RelaxationInfo(prev_entry[0], relaxation_rate)
-        self._momentToRelaxationInfoDict[e] = new_entry
-
-    def set_first_moment_relaxation_rate(self, relaxation_rate):
-        for e in MOMENT_SYMBOLS[:self.dim]:
-            assert e in self._momentToRelaxationInfoDict, "First moments are not relaxed separately by this method"
-        for e in MOMENT_SYMBOLS[:self.dim]:
-            prev_entry = self._momentToRelaxationInfoDict[e]
-            new_entry = RelaxationInfo(prev_entry[0], relaxation_rate)
-            self._momentToRelaxationInfoDict[e] = new_entry
-
-    def set_conserved_moments_relaxation_rate(self, relaxation_rate):
-        self.set_zeroth_moment_relaxation_rate(relaxation_rate)
-        self.set_first_moment_relaxation_rate(relaxation_rate)
-
-    def set_force_model(self, force_model):
-        self._forceModel = force_model
-
-    @property
-    def collision_matrix(self):
-        pdfs_to_moments = self.moment_matrix
-        d = self.relaxation_matrix
-        return pdfs_to_moments.inv() * d * pdfs_to_moments
-
-    @property
-    def inverse_collision_matrix(self):
-        pdfs_to_moments = self.moment_matrix
-        inverse_relaxation_matrix = self.relaxation_matrix.inv()
-        return pdfs_to_moments.inv() * inverse_relaxation_matrix * pdfs_to_moments
-
-    @property
-    def moment_matrix(self):
-        return moment_matrix(self.moments, self.stencil)
-
-    @property
-    def relaxation_matrix(self):
-        d = sp.zeros(len(self.relaxation_rates))
-        for i in range(0, len(self.relaxation_rates)):
-            d[i, i] = self.relaxation_rates[i]
-        return d
-
-    @property
-    def is_orthogonal(self):
-        return (self.moment_matrix * self.moment_matrix.T).is_diagonal()
-
-    @property
-    def is_weighted_orthogonal(self):
-        w = get_weights(self.stencil, sp.Rational(1, 3))
-        return (sp.matrix_multiply_elementwise(self.moment_matrix, sp.Matrix([w] * len(w))) * self.moment_matrix.T
-                ).is_diagonal()
-
-    def __getstate__(self):
-        # Workaround for a bug in joblib
-        self._momentToRelaxationInfoDictToPickle = [i for i in self._momentToRelaxationInfoDict.items()]
-        return self.__dict__
-
-    def _repr_html_(self):
-        table = """
-        <table style="border:none; width: 100%">
-            <tr {nb}>
-                <th {nb} >Moment</th>
-                <th {nb} >Eq. Value </th>
-                <th {nb} >Relaxation Rate</th>
-            </tr>
-            {content}
-        </table>
-        """
-        content = ""
-        for moment, (eq_value, rr) in self._momentToRelaxationInfoDict.items():
-            vals = {
-                'rr': sp.latex(rr),
-                'moment': sp.latex(moment),
-                'eq_value': sp.latex(eq_value),
-                'nb': 'style="border:none"',
-            }
-            content += """<tr {nb}>
-                            <td {nb}>${moment}$</td>
-                            <td {nb}>${eq_value}$</td>
-                            <td {nb}>${rr}$</td>
-                         </tr>\n""".format(**vals)
-        return table.format(content=content, nb='style="border:none"')
-
-    def _compute_weights(self):
-        replacements = self._conservedQuantityComputation.default_values
-        ac = self.get_equilibrium(include_force_terms=False)
-        ac = ac.new_with_substitutions(replacements, substitute_on_lhs=False).new_without_subexpressions()
-
-        new_assignments = [Assignment(e.lhs,
-                                      subs_additive(e.rhs, sp.sympify(1), sum(self.pre_collision_pdf_symbols),
-                                                    required_match_replacement=1.0))
-                           for e in ac.main_assignments]
-        ac = ac.copy(new_assignments)
-
-        weights = []
-        for eq in ac.main_assignments:
-            value = eq.rhs.expand()
-            assert len(value.atoms(sp.Symbol)) == 0, "Failed to compute weights " + str(value)
-            weights.append(value)
-        return weights
-
-    def _collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
-                                               conserved_quantity_equations=None):
-        f = sp.Matrix(self.pre_collision_pdf_symbols)
-        pdf_to_moment = self.moment_matrix
-        m_eq = sp.Matrix(self.moment_equilibrium_values)
-
-        collision_rule = f + pdf_to_moment.inv() * d * (m_eq - pdf_to_moment * f)
-        collision_eqs = [Assignment(lhs, rhs) for lhs, rhs in zip(self.post_collision_pdf_symbols, collision_rule)]
-
-        if conserved_quantity_equations is None:
-            conserved_quantity_equations = self._conservedQuantityComputation.equilibrium_input_equations_from_pdfs(f)
-
-        simplification_hints = conserved_quantity_equations.simplification_hints.copy()
-        simplification_hints.update(self._conservedQuantityComputation.defined_symbols())
-        simplification_hints['relaxation_rates'] = [d[i, i] for i in range(d.rows)]
-
-        all_subexpressions = list(additional_subexpressions) + conserved_quantity_equations.all_assignments
-
-        if self._forceModel is not None and include_force_terms:
-            force_model_terms = self._forceModel(self)
-            force_term_symbols = sp.symbols("forceTerm_:%d" % (len(force_model_terms,)))
-            force_subexpressions = [Assignment(sym, force_model_term)
-                                    for sym, force_model_term in zip(force_term_symbols, force_model_terms)]
-            all_subexpressions += force_subexpressions
-            collision_eqs = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
-                             for eq, force_term_symbol in zip(collision_eqs, force_term_symbols)]
-            simplification_hints['force_terms'] = force_term_symbols
-
-        return LbmCollisionRule(self, collision_eqs, all_subexpressions,
-                                simplification_hints)
-
-    @staticmethod
-    def _generate_relaxation_matrix(relaxation_matrix, keep_rr_symbolic):
-        """
-        For SRT and TRT the equations can be easier simplified if the relaxation times are symbols, not numbers.
-        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_matrix[i, i] for i in range(relaxation_matrix.rows)]
-        if keep_rr_symbolic <= 2:
-            unique_relaxation_rates = set(rr)
-            subexpressions = {}
-            for rt in unique_relaxation_rates:
-                rt = sp.sympify(rt)
-                if not isinstance(rt, sp.Symbol):
-                    rt_symbol = sp.Symbol("rr_%d" % (len(subexpressions),))
-                    subexpressions[rt] = rt_symbol
-
-            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()]
-            return substitutions, sp.diag(*new_rr)
-        else:
-            return [], relaxation_matrix

lbmpy/phasefield_allen_cahn/analytical.py

-
-def analytic_rising_speed(gravitational_acceleration, bubble_diameter, viscosity_gas):
-    r"""
-    Calculated the analytical rising speed of a bubble. This is the expected end rising speed.
-    Args:
-        gravitational_acceleration: the gravitational acceleration acting in the simulation scenario. Usually it gets
-                                    calculated based on dimensionless parameters which describe the scenario
-        bubble_diameter: the diameter of the bubble at the beginning of the simulation
-        viscosity_gas: the viscosity of the fluid inside the bubble
-    """
-    result = -(gravitational_acceleration * bubble_diameter * bubble_diameter) / (12.0 * viscosity_gas)
-    return result

lbmpy/phasefield_allen_cahn/force_model.py

-import sympy as sp
-import numpy as np
-
-from lbmpy.forcemodels import Simple
-
-
-class MultiphaseForceModel:
-    r"""
-    A force model based on PhysRevE.96.053301. This model realises the modified equilibrium distributions meaning the
-    force gets shifted by minus one half multiplied with the collision operator
-    """
-    def __init__(self, force, rho=1):
-        self._force = force
-        self._rho = rho
-
-    def __call__(self, lb_method):
-        stencil = lb_method.stencil
-
-        force_symp = sp.symbols("force_:{}".format(lb_method.dim))
-        simple = Simple(force_symp)
-        force = [f / self._rho for f in simple(lb_method)]
-
-        moment_matrix = lb_method.moment_matrix
-        relaxation_rates = sp.Matrix(np.diag(lb_method.relaxation_rates))
-        mrt_collision_op = moment_matrix.inv() * relaxation_rates * moment_matrix
-
-        result = -0.5 * mrt_collision_op * sp.Matrix(force) + sp.Matrix(force)
-
-        for i in range(0, len(stencil)):
-            result[i] = result[i].simplify()
-
-        subs_dict = dict(zip(force_symp, self._force))
-
-        for i in range(0, len(stencil)):
-            result[i] = result[i].subs(subs_dict)
-
-        return result

lbmpy/phasefield_allen_cahn/kernel_equations.py

-from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
-from pystencils import Assignment, AssignmentCollection
-
-from lbmpy.maxwellian_equilibrium import get_weights
-from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAccessor
-
-import sympy as sp
-import numpy as np
-
-
-def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
-    r"""
-    Get a symbolic expression for the chemical potential according to equation (5) in PhysRevE.96.053301.
-    Args:
-        phi_field: the phase-field on which the chemical potential is applied
-        stencil: stencil to derive the finite difference for the laplacian (2nd order isotropic)
-        beta: coefficient related to surface tension and interface thickness
-        kappa: coefficient related to surface tension and interface thickness
-    """
-    dimensions = len(stencil[0])
-    q = len(stencil)
-    lap = sp.simplify(0)
-    for i in range(dimensions):
-        deriv = FiniteDifferenceStencilDerivation((i, i), stencil)
-        for j in range(dimensions):
-            # assume the stencil is symmetric
-            deriv.assume_symmetric(dim=j, anti_symmetric=False)
-
-        # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic
-        if q == 9:
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-        elif q == 15:
-            deriv.set_weight((0, 0, 0), sp.Rational(-32, 27))
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-        elif q == 19:
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-        else:
-            deriv.set_weight((0, 0, 0), sp.Rational(-38, 27))
-            res = deriv.get_stencil(isotropic=True)
-            lap += res.apply(phi_field.center)
-
-    # get the chemical potential
-    mu = 4.0 * beta * phi_field.center * (phi_field.center - 1.0) * (phi_field.center - 0.5) - \
-        kappa * lap
-    return mu
-
-
-def isotropic_gradient_symbolic(phi_field, stencil):
-    r"""
-    Get a symbolic expression for the isotropic gradient of the phase-field
-    Args:
-        phi_field: the phase-field on which the isotropic gradient is applied
-        stencil: stencil to derive the finite difference for the gradient (2nd order isotropic)
-    """
-    dimensions = len(stencil[0])
-    q = len(stencil)
-    deriv = FiniteDifferenceStencilDerivation((0,), stencil)
-
-    deriv.assume_symmetric(0, anti_symmetric=True)
-    deriv.assume_symmetric(1, anti_symmetric=False)
-    if dimensions == 3:
-        deriv.assume_symmetric(2, anti_symmetric=False)
-
-    # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic
-    # furthermore the stencils gets rotated to get the y and z components
-    if q == 9:
-        res = deriv.get_stencil(isotropic=True)
-        grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
-    elif q == 15:
-        res = deriv.get_stencil(isotropic=True)
-        grad = [res.apply(phi_field.center),
-                res.rotate_weights_and_apply(phi_field.center, (0, 1)),
-                res.rotate_weights_and_apply(phi_field.center, (1, 2))]
-    elif q == 19:
-        deriv.set_weight((0, 0, 0), sp.sympify(0))
-        deriv.set_weight((1, 0, 0), sp.Rational(1, 6))
-
-        res = deriv.get_stencil(isotropic=True)
-        grad = [res.apply(phi_field.center),
-                res.rotate_weights_and_apply(phi_field.center, (0, 1)),
-                res.rotate_weights_and_apply(phi_field.center, (1, 2))]
-    else:
-        deriv.set_weight((0, 0, 0), sp.sympify(0))
-        deriv.set_weight((1, 0, 0), sp.Rational(2, 9))
-
-        res = deriv.get_stencil(isotropic=True)
-        grad = [res.apply(phi_field.center),
-                res.rotate_weights_and_apply(phi_field.center, (0, 1)),
-                res.rotate_weights_and_apply(phi_field.center, (1, 2))]
-
-    return grad
-
-
-def normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil=None):
-    r"""
-    Get a symbolic expression for the normalized isotropic gradient of the phase-field
-    Args:
-        phi_field: the phase-field on which the normalized isotropic gradient is applied
-        stencil: stencil of the lattice Boltzmann step
-        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
-        field. If it is not given the stencil of the LB method will be applied.
-    """
-    if fd_stencil is None:
-        fd_stencil = stencil
-
-    tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, fd_stencil))) + 1.e-32) ** 0.5
-
-    result = [x / tmp for x in isotropic_gradient_symbolic(phi_field, fd_stencil)]
-    return result
-
-
-def pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil=None):
-    r"""
-    Get a symbolic expression for the pressure force
-    Args:
-        phi_field: phase-field
-        stencil: stencil of the lattice Boltzmann step
-        density_heavy: density of the heavier fluid
-        density_light: density of the lighter fluid
-        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
-        field. If it is not given the stencil of the LB method will be applied.
-    """
-    if fd_stencil is None:
-        fd_stencil = stencil
-
-    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
-    result = list(map(lambda x: sp.Rational(-1, 3) * sp.symbols("rho") * (density_heavy - density_light) * x, iso_grad))
-    return result
-
-
-def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light, fd_stencil=None):
-    r"""
-    Get a symbolic expression for the viscous force
-    Args:
-        lb_velocity_field: hydrodynamic distribution function
-        phi_field: phase-field
-        mrt_method: mrt lattice boltzmann method used for hydrodynamics
-        tau: relaxation time of the hydrodynamic lattice boltzmann step
-        density_heavy: density of the heavier fluid
-        density_light: density of the lighter fluid
-        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
-        field. If it is not given the stencil of the LB method will be applied.
-    """
-    stencil = mrt_method.stencil
-    dimensions = len(stencil[0])
-
-    if fd_stencil is None:
-        fd_stencil = stencil
-
-    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
-
-    non_equilibrium = lb_velocity_field.center_vector - mrt_method.get_equilibrium_terms()
-
-    stress_tensor = [0] * 6
-    # Calculate Stress Tensor MRT
-    for i, d in enumerate(stencil):
-        stress_tensor[0] = sp.Add(stress_tensor[0], non_equilibrium[i] * (d[0] * d[0]))
-        stress_tensor[1] = sp.Add(stress_tensor[1], non_equilibrium[i] * (d[1] * d[1]))
-
-        if dimensions == 3:
-            stress_tensor[2] = sp.Add(stress_tensor[2], non_equilibrium[i] * (d[2] * d[2]))
-            stress_tensor[3] = sp.Add(stress_tensor[3], non_equilibrium[i] * (d[1] * d[2]))
-            stress_tensor[4] = sp.Add(stress_tensor[4], non_equilibrium[i] * (d[0] * d[2]))
-
-        stress_tensor[5] = sp.Add(stress_tensor[5], non_equilibrium[i] * (d[0] * d[1]))
-
-    density_difference = density_heavy - density_light
-
-    # Calculate Viscous Force MRT
-    fmx = (0.5 - tau) * (stress_tensor[0] * iso_grad[0]
-                         + stress_tensor[5] * iso_grad[1]
-                         + stress_tensor[4] * iso_grad[2]) * density_difference
-
-    fmy = (0.5 - tau) * (stress_tensor[5] * iso_grad[0]
-                         + stress_tensor[1] * iso_grad[1]
-                         + stress_tensor[3] * iso_grad[2]) * density_difference
-
-    fmz = (0.5 - tau) * (stress_tensor[4] * iso_grad[0]
-                         + stress_tensor[3] * iso_grad[1]
-                         + stress_tensor[2] * iso_grad[2]) * density_difference
-
-    return [fmx, fmy, fmz]
-
-
-def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None):
-    r"""
-    Get a symbolic expression for the surface tension force
-    Args:
-        phi_field: the phase-field on which the chemical potential is applied
-        stencil: stencil of the lattice Boltzmann step
-        beta: coefficient related to surface tension and interface thickness
-        kappa: coefficient related to surface tension and interface thickness
-        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
-        field. If it is not given the stencil of the LB method will be applied.
-    """
-    if fd_stencil is None:
-        fd_stencil = stencil
-
-    chemical_potential = chemical_potential_symbolic(phi_field, fd_stencil, beta, kappa)
-    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
-    return [chemical_potential * x for x in iso_grad]
-
-
-def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
-                       density_heavy, density_light, kappa, beta, body_force, fd_stencil=None):
-    r"""
-    Get a symbolic expression for the hydrodynamic force
-    Args:
-        lb_velocity_field: hydrodynamic distribution function
-        phi_field: phase-field
-        lb_method: Lattice boltzmann method used for hydrodynamics
-        tau: relaxation time of the hydrodynamic lattice boltzmann step
-        density_heavy: density of the heavier fluid
-        density_light: density of the lighter fluid
-        beta: coefficient related to surface tension and interface thickness
-        kappa: coefficient related to surface tension and interface thickness
-        body_force: force acting on the fluids. Usually the gravity
-        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
-        field. If it is not given the stencil of the LB method will be applied.
-    """
-    stencil = lb_method.stencil
-    dimensions = len(stencil[0])
-
-    if fd_stencil is None:
-        fd_stencil = stencil
-
-    fp = pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil)
-    fm = viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light, fd_stencil)
-    fs = surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil)
-
-    result = []
-    for i in range(dimensions):
-        result.append(fs[i] + fp[i] + fm[i] + body_force[i])
-
-    return result
-
-
-def interface_tracking_force(phi_field, stencil, interface_thickness, fd_stencil=None):
-    r"""
-    Get a symbolic expression for the hydrodynamic force
-    Args:
-        phi_field: phase-field
-        stencil: stencil of the phase-field distribution lattice Boltzmann step
-        interface_thickness: interface thickness
-        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
-        field. If it is not given the stencil of the LB method will be applied.
-    """
-    if fd_stencil is None:
-        fd_stencil = stencil
-
-    dimensions = len(stencil[0])
-    normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
-    result = []
-    for i in range(dimensions):
-        result.append(((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness) * normal_fd[i])
-
-    return result
-
-
-def get_update_rules_velocity(src_field, u_in, lb_method, force, density, sub_iterations=2):
-    r"""
-     Get assignments to update the velocity with a force shift
-     Args:
-         src_field: the source field of the hydrodynamic distribution function
-         u_in: velocity field
-         lb_method: mrt lattice boltzmann method used for hydrodynamics
-         force: force acting on the hydrodynamic lb step
-         density: the interpolated density of the simulation
-         sub_iterations: number of updates of the velocity field
-     """
-    stencil = lb_method.stencil
-    dimensions = len(stencil[0])
-
-    moment_matrix = lb_method.moment_matrix
-    eq = lb_method.moment_equilibrium_values
-
-    first_eqs = lb_method.first_order_equilibrium_moment_symbols
-    indices = list()
-    for i in range(dimensions):
-        indices.append(eq.index(first_eqs[i]))
-
-    src = [src_field.center(i) for i, _ in enumerate(stencil)]
-    m0 = np.dot(moment_matrix.tolist(), src)
-
-    update_u = list()
-    update_u.append(Assignment(sp.symbols("rho"), m0[0]))
-
-    index = 0
-    u_symp = sp.symbols("u_:{}".format(dimensions))
-    aleph = sp.symbols("aleph_:{}".format(dimensions * sub_iterations))
-
-    for i in range(dimensions):
-        update_u.append(Assignment(aleph[i], u_in.center_vector[i]))
-        index += 1
-
-    for k in range(sub_iterations - 1):
-        subs_dict = dict(zip(u_symp, aleph[k * dimensions:index]))
-        for i in range(dimensions):
-            update_u.append(Assignment(aleph[index], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
-            index += 1
-
-    subs_dict = dict(zip(u_symp, aleph[index - dimensions:index]))
-    for i in range(dimensions):
-        update_u.append(Assignment(u_symp[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
-        # update_u.append(Assignment(u_in.center_vector[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
-
-    return update_u
-
-
-def get_collision_assignments_hydro(lb_method, density, velocity_input, force, sub_iterations, symbolic_fields,
-                                    kernel_type):
-    r"""
-     Get collision assignments for the hydrodynamic lattice Boltzmann step. Here the force gets applied in the moment
-     space. Afterwards the transformation back to the pdf space happens.
-     Args:
-         lb_method: moment based lattice Boltzmann method
-         density: the interpolated density of the simulation
-         velocity_input: velocity field for the hydrodynamic and Allen-Chan LB step
-         force: force vector containing a summation of the surface tension-, pressure-, viscous- and bodyforce vector
-         sub_iterations: number of updates of the velocity field
-         symbolic_fields: PDF fields for source and destination
-         kernel_type: collide_stream_push or collide_only
-     """
-
-    stencil = lb_method.stencil
-    dimensions = len(stencil[0])
-
-    src_field = symbolic_fields['symbolic_field']
-    dst_field = symbolic_fields['symbolic_temporary_field']
-
-    moment_matrix = lb_method.moment_matrix
-    rel = lb_method.relaxation_rates
-    eq = lb_method.moment_equilibrium_values
-
-    first_eqs = lb_method.first_order_equilibrium_moment_symbols
-    indices = list()
-    for i in range(dimensions):
-        indices.append(eq.index(first_eqs[i]))
-
-    eq = np.array(eq)
-
-    g_vals = [src_field.center(i) for i, _ in enumerate(stencil)]
-    m0 = np.dot(moment_matrix.tolist(), g_vals)
-
-    mf = np.zeros(len(stencil), dtype=object)
-    for i in range(dimensions):
-        mf[indices[i]] = force[i] / density
-
-    m = sp.symbols("m_:{}".format(len(stencil)))
-
-    update_m = get_update_rules_velocity(src_field, velocity_input, lb_method, force,
-                                         density, sub_iterations=sub_iterations)
-    u_symp = sp.symbols("u_:{}".format(dimensions))
-
-    for i in range(0, len(stencil)):
-        update_m.append(Assignment(m[i], m0[i] - (m0[i] - eq[i] + mf[i] / 2) * rel[i] + mf[i]))
-
-    update_g = list()
-    var = np.dot(moment_matrix.inv().tolist(), m)
-    if kernel_type == 'collide_stream_push':
-        push_accessor = StreamPushTwoFieldsAccessor()
-        post_collision_accesses = push_accessor.write(dst_field, stencil)
-    else:
-        collide_accessor = CollideOnlyInplaceAccessor()
-        post_collision_accesses = collide_accessor.write(src_field, stencil)
-
-    for i in range(0, len(stencil)):
-        update_g.append(Assignment(post_collision_accesses[i], var[i]))
-
-    for i in range(dimensions):
-        update_g.append(Assignment(velocity_input.center_vector[i], u_symp[i]))
-
-    hydro_lb_update_rule = AssignmentCollection(main_assignments=update_g,
-                                                subexpressions=update_m)
-
-    return hydro_lb_update_rule
-
-
-def initializer_kernel_phase_field_lb(lb_phase_field, phi_field, velocity_field, mrt_method, interface_thickness,
-                                      fd_stencil=None):
-    r"""
-    Returns an assignment list for initializing the phase-field distribution functions
-    Args:
-        lb_phase_field: source field of phase-field distribution function
-        phi_field: phase-field
-        velocity_field: velocity field
-        mrt_method: lattice Boltzmann method of the phase-field lattice Boltzmann step
-        interface_thickness: interface thickness
-        fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
-        field. If it is not given the stencil of the LB method will be applied.
-    """
-    stencil = mrt_method.stencil
-    dimensions = len(stencil[0])
-
-    if fd_stencil is None:
-        fd_stencil = stencil
-
-    weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
-    u_symp = sp.symbols("u_:{}".format(dimensions))
-
-    normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
-
-    gamma = mrt_method.get_equilibrium_terms()
-    gamma = gamma.subs({sp.symbols("rho"): 1})
-    gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)})
-    # create the kernels for the initialization of the h field
-    h_updates = list()
-
-    def scalar_product(a, b):
-        return sum(a_i * b_i for a_i, b_i in zip(a, b))
-
-    f = []
-    for i, d in enumerate(stencil):
-        f.append(weights[i] * ((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness)
-                 * scalar_product(d, normal_fd[0:dimensions]))
-
-    for i, _ in enumerate(stencil):
-        h_updates.append(Assignment(lb_phase_field.center(i), phi_field.center * gamma_init[i] - 0.5 * f[i]))
-
-    return h_updates
-
-
-def initializer_kernel_hydro_lb(lb_velocity_field, velocity_field, mrt_method):
-    r"""
-    Returns an assignment list for initializing the velocity distribution functions
-    Args:
-        lb_velocity_field: source field of velocity distribution function
-        velocity_field: velocity field
-        mrt_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step
-    """
-    stencil = mrt_method.stencil
-    dimensions = len(stencil[0])
-    weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
-    u_symp = sp.symbols("u_:{}".format(dimensions))
-
-    gamma = mrt_method.get_equilibrium_terms()
-    gamma = gamma.subs({sp.symbols("rho"): 1})
-    gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)})
-
-    g_updates = list()
-    for i, _ in enumerate(stencil):
-        g_updates.append(Assignment(lb_velocity_field.center(i), gamma_init[i] - weights[i]))
-
-    return g_updates

lbmpy/phasefield_allen_cahn/parameter_calculation.py

-import math
-
-
-def calculate_parameters_rti(reference_length=256,
-                             reference_time=16000,
-                             density_heavy=1.0,
-                             capillary_number=0.26,
-                             reynolds_number=3000,
-                             atwood_number=0.5,
-                             peclet_number=1000,
-                             density_ratio=3,
-                             viscosity_ratio=1):
-    r"""
-    Calculate the simulation parameters for the Rayleigh-Taylor instability.
-    The calculation is done according to the description in part B of PhysRevE.96.053301.
-
-    Args:
-        reference_length: reference length of the RTI
-        reference_time: chosen reference time
-        density_heavy: density of the heavier fluid
-        capillary_number: capillary number of the simulation represents the relative effect of viscous drag
-                          forces versus surface tension forces
-        reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
-                         and the inertial forces
-        atwood_number: atwood number of the simulation is a dimensionless density ratio
-        peclet_number: peclet number of the simulation is the ratio of the rate of advection
-                       of a physical quantity by the flow to the rate of diffusion of the same quantity
-                       driven by an appropriate gradient
-        density_ratio: density ratio of the heavier and the lighter fluid
-        viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
-    """
-
-    # calculate the gravitational acceleration and the reference velocity
-    g = reference_length / (reference_time ** 2 * atwood_number)
-    reference_velocity = math.sqrt(g * reference_length)
-
-    dynamic_viscosity_heavy = (density_heavy * reference_velocity * reference_length) / reynolds_number
-    dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
-
-    density_light = density_heavy / density_ratio
-
-    kinematic_viscosity_heavy = dynamic_viscosity_heavy / density_heavy
-    kinematic_viscosity_light = dynamic_viscosity_light / density_light
-
-    relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy
-    relaxation_time_light = 3.0 * kinematic_viscosity_light
-
-    surface_tension = (dynamic_viscosity_heavy * reference_velocity) / capillary_number
-    mobility = (reference_velocity * reference_length) / peclet_number
-
-    parameters = {
-        "density_light": density_light,
-        "dynamic_viscosity_heavy": dynamic_viscosity_heavy,
-        "dynamic_viscosity_light": dynamic_viscosity_light,
-        "relaxation_time_heavy": relaxation_time_heavy,
-        "relaxation_time_light": relaxation_time_light,
-        "gravitational_acceleration": -g,
-        "mobility": mobility,
-        "surface_tension": surface_tension
-    }
-    return parameters
-
-
-def calculate_dimensionless_rising_bubble(reference_time=18000,
-                                          density_heavy=1.0,
-                                          bubble_radius=16,
-                                          bond_number=1,
-                                          reynolds_number=40,
-                                          density_ratio=1000,
-                                          viscosity_ratio=100):
-    r"""
-    Calculate the simulation parameters for a rising bubble. The parameter calculation follows the ideas of Weber and
-    is based on the Reynolds number. This means the rising velocity of the bubble is implicitly stated with the
-    Reynolds number
-
-    Args:
-        reference_time: chosen reference time
-        density_heavy: density of the heavier fluid
-        bubble_radius: initial radius of the rising bubble
-        bond_number: also called eötvös number is measuring the importance of gravitational forces compared to
-                     surface tension forces
-        reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
-                 and the inertial forces
-        density_ratio: density ratio of the heavier and the lighter fluid
-        viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
-    """
-
-    bubble_diameter = bubble_radius * 2
-    g = bubble_diameter / (reference_time ** 2)
-
-    density_light = density_heavy / density_ratio
-
-    dynamic_viscosity_heavy = (density_heavy * math.sqrt(g * bubble_diameter ** 3)) / reynolds_number
-    dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
-
-    kinematic_viscosity_heavy = dynamic_viscosity_heavy / density_heavy
-    kinematic_viscosity_light = dynamic_viscosity_light / density_light
-
-    relaxation_time_heavy = 3 * kinematic_viscosity_heavy
-    relaxation_time_light = 3 * kinematic_viscosity_light
-
-    surface_tension = (density_heavy - density_light) * g * bubble_diameter ** 2 / bond_number
-    # calculation of the Morton number
-    # Mo = gravitational_acceleration * dynamic_viscosity_heavy / (density_heavy * surface_tension ** 3)
-
-    parameters = {
-        "density_light": density_light,
-        "dynamic_viscosity_heavy": dynamic_viscosity_heavy,
-        "dynamic_viscosity_light": dynamic_viscosity_light,
-        "relaxation_time_heavy": relaxation_time_heavy,
-        "relaxation_time_light": relaxation_time_light,
-        "gravitational_acceleration": -g,
-        "surface_tension": surface_tension
-    }
-    return parameters

lbmpy/simplificationfactory.py

-import sympy as sp
-
-from lbmpy.innerloopsplit import create_lbm_split_groups
-from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
-from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
-from lbmpy.methods.centeredcumulant.simplification import insert_aliases, insert_zeros, insert_constants
-from lbmpy.methods.momentbased.momentbasedsimplifications import (
-    factor_density_after_factoring_relaxation_times, factor_relaxation_rates,
-    replace_common_quadratic_and_constant_term, replace_density_and_velocity, replace_second_order_velocity_products)
-from pystencils.simp import (
-    SimplificationStrategy, add_subexpressions_for_divisions, apply_to_all_assignments,
-    subexpression_substitution_in_main_assignments)
-
-
-def create_simplification_strategy(lb_method, split_inner_loop=False):
-    s = SimplificationStrategy()
-    expand = apply_to_all_assignments(sp.expand)
-
-    if isinstance(lb_method, MomentBasedLbMethod):
-        if len(set(lb_method.relaxation_rates)) <= 2:
-            s.add(expand)
-            s.add(replace_second_order_velocity_products)
-            s.add(expand)
-            s.add(factor_relaxation_rates)
-            s.add(replace_density_and_velocity)
-            s.add(replace_common_quadratic_and_constant_term)
-            s.add(factor_density_after_factoring_relaxation_times)
-            s.add(subexpression_substitution_in_main_assignments)
-            if split_inner_loop:
-                s.add(create_lbm_split_groups)
-            s.add(add_subexpressions_for_divisions)
-        else:
-            s.add(subexpression_substitution_in_main_assignments)
-            if split_inner_loop:
-                s.add(create_lbm_split_groups)
-    elif isinstance(lb_method, CenteredCumulantBasedLbMethod):
-        s.add(insert_zeros)
-        s.add(insert_aliases)
-        s.add(insert_constants)
-        s.add(lambda ac: ac.new_without_unused_subexpressions())
-    return s

lbmpy/stencils.py

-def get_stencil(name, ordering='walberla'):
-    """
-    Stencils are tuples of discrete velocities. They are commonly labeled in the 'DxQy' notation, where d is the
-    dimension (length of the velocity tuples) and y is number of discrete velocities.
-
-    Args:
-        name: DxQy notation
-        ordering: the LBM literature does not use a common order of the discrete velocities, therefore here
-                  different common orderings are available. All orderings lead to the same method, it just has
-                  to be used consistently. Here more orderings are available to compare intermediate results with
-                  the literature.
-    """
-    try:
-        return get_stencil.data[name.upper()][ordering.lower()]
-    except KeyError:
-        err_msg = ""
-        for stencil, ordering_names in get_stencil.data.items():
-            err_msg += "  %s: %s\n" % (stencil, ", ".join(ordering_names.keys()))
-
-        raise ValueError("No such stencil available. "
-                         "Available stencils: <stencil_name>( <ordering_names> )\n" + err_msg)
-
-
-get_stencil.data = {
-    'D2Q9': {
-        'walberla': ((0, 0),
-                     (0, 1), (0, -1), (-1, 0), (1, 0),
-                     (-1, 1), (1, 1), (-1, -1), (1, -1),),
-        'counterclockwise': ((0, 0),
-                             (1, 0), (0, 1), (-1, 0), (0, -1),
-                             (1, 1), (-1, 1), (-1, -1), (1, -1)),
-        'braunschweig': ((0, 0),
-                         (-1, 1), (-1, 0), (-1, -1), (0, -1),
-                         (1, -1), (1, 0), (1, 1), (0, 1)),
-        'uk': ((0, 0),
-               (1, 0), (-1, 0), (0, 1), (0, -1),
-               (1, 1), (-1, -1), (-1, 1), (1, -1),
-               )
-    },
-    'D2V17': {
-        'walberla': (
-            (0, 0), (0, -1), (-1, 0), (1, 0), (0, 1), (-1, -1), (1, -1), (-1, 1), (1, 1), (-2, -2), (2, -2), (-2, 2),
-            (2, 2), (0, -3), (-3, 0), (3, 0), (0, 3)),
-    },
-    'D2V37': {
-        'walberla': (
-            (0, 0), (0, -1), (-1, 0), (1, 0), (0, 1), (-1, -1), (1, -1), (-1, 1), (1, 1), (0, -2), (-2, 0), (2, 0),
-            (0, 2), (-1, -2), (1, -2), (-2, -1), (2, -1), (-2, 1), (2, 1), (-1, 2), (1, 2), (-2, -2), (2, -2), (-2, 2),
-            (2, 2), (0, -3), (-3, 0), (3, 0), (0, 3), (-1, -3), (1, -3), (-3, -1), (3, -1), (-3, 1), (3, 1), (-1, 3),
-            (1, 3))
-    },
-    'D3Q15': {
-        'walberla':
-            ((0, 0, 0),
-             (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
-             (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
-        'premnath': ((0, 0, 0),
-                     (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
-                     (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
-                     (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
-    },
-    'D3Q19': {
-        'walberla': ((0, 0, 0),
-                     (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
-                     (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
-                     (0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
-                     (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1)),
-        'braunschweig': ((0, 0, 0),
-                         (1, 0, 0), (-1, 0, 0),
-                         (0, 1, 0), (0, -1, 0),
-                         (0, 0, 1), (0, 0, -1),
-                         (1, 1, 0), (-1, -1, 0),
-                         (1, -1, 0), (-1, 1, 0),
-                         (1, 0, 1), (-1, 0, -1),
-                         (1, 0, -1), (-1, 0, 1),
-                         (0, 1, 1), (0, -1, -1),
-                         (0, 1, -1), (0, -1, 1)),
-        'premnath': ((0, 0, 0),
-                     (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
-                     (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
-                     (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
-                     (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)),
-    },
-    'D3Q27': {
-        'walberla': ((0, 0, 0),
-                     (0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
-                     (-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
-                     (0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
-                     (0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1),
-                     (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1),
-                     (-1, -1, -1)),
-        'premnath': ((0, 0, 0),
-                     (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
-                     (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
-                     (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
-                     (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1),
-                     (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
-                     (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
-        'fakhari': ((0, 0, 0),
-                    (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
-                    (1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
-                    (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1),
-                    (1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
-                    (1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1), (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)),
-    }
-}

lbmpy/turbulence_models.py

-import sympy as sp
-
-from lbmpy.relaxationrates import get_shear_relaxation_rate
-from pystencils import Assignment
-
-
-def second_order_moment_tensor(function_values, stencil):
-    """Returns (D x D) Matrix of second order moments of the given function where D is the dimension"""
-    assert len(function_values) == len(stencil)
-    dim = len(stencil[0])
-    return sp.Matrix(dim, dim, lambda i, j: sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))
-
-
-def frobenius_norm(matrix, factor=1):
-    """Computes the Frobenius norm of a matrix defined as the square root of the sum of squared matrix elements
-    The optional factor is added inside the square root"""
-    return sp.sqrt(sum(i * i for i in matrix) * factor)
-
-
-def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_field=None):
-    r""" Adds a smagorinsky model to a lattice Boltzmann collision rule. To add the Smagorinsky model to a LB scheme
-        one has to first compute the strain rate tensor $S_{ij}$ in each cell, and compute the turbulent
-        viscosity :math:`nu_t` from it. Then the local relaxation rate has to be adapted to match the total viscosity
-        :math `\nu_{total}` instead of the standard viscosity :math `\nu_0`.
-
-        A fortunate property of LB methods is, that the strain rate tensor can be computed locally from the
-        non-equilibrium part of the distribution function. This is somewhat surprising, since the strain rate tensor
-        contains first order derivatives. The strain rate tensor can be obtained by
-
-        .. math ::
-            S_{ij} = - \frac{3 \omega_s}{2 \rho_{(0)}} \Pi_{ij}^{(neq)}
-
-        where :math `\omega_s` is the relaxation rate that determines the viscosity, :math `\rho_{(0)}` is :math `\rho`
-        in compressible models and :math `1` for incompressible schemes.
-        :math `\Pi_{ij}^{(neq)}` is the second order moment tensor of the non-equilibrium part of
-        the distribution functions
-        :math `f^{(neq)} = f - f^{(eq)}` and can be computed as
-
-        .. math ::
-            \Pi_{ij}^{(neq)} = \sum_q c_{qi} c_{qj} \; f_q^{(neq)}
-
-
-    """
-    method = collision_rule.method
-    omega_s = get_shear_relaxation_rate(method)
-    if isinstance(omega_s, float) or isinstance(omega_s, int):
-        raise ValueError("For the smagorinsky model the shear relaxation rate has to be a symbol")
-    f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
-
-    tau_0 = sp.Symbol("tau_0_")
-    second_order_neq_moments = sp.Symbol("Pi")
-    rho = method.zeroth_order_equilibrium_moment_symbol if method.conserved_quantity_computation.compressible else 1
-    adapted_omega = sp.Symbol("smagorinsky_omega")
-
-    collision_rule = collision_rule.new_with_substitutions({omega_s: adapted_omega})
-    # for derivation see notebook demo_custom_LES_model.ipynb
-    eqs = [Assignment(tau_0, 1 / omega_s),
-           Assignment(second_order_neq_moments,
-                      frobenius_norm(second_order_moment_tensor(f_neq, method.stencil), factor=2) / rho),
-           Assignment(adapted_omega,
-                      2 / (tau_0 + sp.sqrt(18 * smagorinsky_constant ** 2 * second_order_neq_moments + tau_0 ** 2)))]
-    collision_rule.subexpressions += eqs
-    collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
-
-    if omega_output_field:
-        collision_rule.main_assignments.append(Assignment(omega_output_field.center, adapted_omega))
-
-    return collision_rule

lbmpy_tests/advanced_streaming/test_communication.py

-import pystencils as ps
-
-import numpy as np
-
-from lbmpy.stencils import get_stencil
-from pystencils.slicing import get_slice_before_ghost_layer, get_ghost_region_slice
-from lbmpy.creationfunctions import create_lb_update_rule
-from lbmpy.advanced_streaming.communication import get_communication_slices, _fix_length_one_slices, \
-    LBMPeriodicityHandling
-from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep
-
-import pytest
-
-
-@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
-@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
-@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
-def test_slices_not_empty(stencil, streaming_pattern, timestep):
-    stencil = get_stencil(stencil)
-    dim = len(stencil[0])
-    q = len(stencil)
-    arr = np.zeros((4,) * dim + (q,))
-    slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep,
-                                      ghost_layers=1)
-    for _, slices_list in slices.items():
-        for src, dst in slices_list:
-            assert all(s != 0 for s in arr[src].shape)
-            assert all(s != 0 for s in arr[dst].shape)
-
-
-@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
-@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
-@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
-def test_src_dst_same_shape(stencil, streaming_pattern, timestep):
-    stencil = get_stencil(stencil)
-    dim = len(stencil[0])
-    q = len(stencil)
-    arr = np.zeros((4,) * dim + (q,))
-    slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep,
-                                      ghost_layers=1)
-    for _, slices_list in slices.items():
-        for src, dst in slices_list:
-            src_shape = arr[src].shape
-            dst_shape = arr[dst].shape
-            assert src_shape == dst_shape
-
-
-@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
-def test_pull_communication_slices(stencil):
-    stencil = get_stencil(stencil)
-
-    slices = get_communication_slices(
-        stencil, streaming_pattern='pull', prev_timestep=Timestep.BOTH, ghost_layers=1)
-    for i, d in enumerate(stencil):
-        if i == 0:
-            continue
-
-        for s in slices[d]:
-            if s[0][-1] == i:
-                src = s[0][:-1]
-                dst = s[1][:-1]
-                break
-
-        inner_slice = _fix_length_one_slices(get_slice_before_ghost_layer(d, ghost_layers=1))
-        inv_dir = (-e for e in d)
-        gl_slice = _fix_length_one_slices(get_ghost_region_slice(inv_dir, ghost_layers=1))
-        assert src == inner_slice
-        assert dst == gl_slice
-
-
-@pytest.mark.parametrize('stencil_name', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
-def test_optimised_and_full_communication_equivalence(stencil_name):
-    target = 'cpu'
-    stencil = get_stencil(stencil_name)
-    dim = len(stencil[0])
-    domain_size = (4, ) * dim
-
-    dh = ps.create_data_handling(domain_size, periodicity=(True, ) * dim,
-                                 parallel=False, default_target=target)
-
-    pdf = dh.add_array("pdf", values_per_cell=len(stencil), dtype=np.int64)
-    dh.fill("pdf", 0, ghost_layers=True)
-    pdf_tmp = dh.add_array("pdf_tmp", values_per_cell=len(stencil), dtype=np.int64)
-    dh.fill("pdf_tmp", 0, ghost_layers=True)
-
-    gl = dh.ghost_layers_of_field("pdf")
-
-    num = 0
-    for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']):
-        dh.cpu_arrays['pdf'][idx] = num
-        dh.cpu_arrays['pdf_tmp'][idx] = num
-        num += 1
-
-    ac = create_lb_update_rule(stencil=stencil,
-                               optimization={"symbolic_field": pdf,
-                                             "symbolic_temporary_field": pdf_tmp},
-                               kernel_type='stream_pull_only')
-    ast = ps.create_kernel(ac, target=dh.default_target, cpu_openmp=True)
-    stream = ast.compile()
-
-    full_communication = dh.synchronization_function(pdf.name, target=dh.default_target, optimization={"openmp": True})
-    full_communication()
-
-    dh.run_kernel(stream)
-    dh.swap("pdf", "pdf_tmp")
-    pdf_full_communication = np.copy(dh.cpu_arrays['pdf'])
-
-    num = 0
-    for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']):
-        dh.cpu_arrays['pdf'][idx] = num
-        dh.cpu_arrays['pdf_tmp'][idx] = num
-        num += 1
-
-    optimised_communication = LBMPeriodicityHandling(stencil=stencil, data_handling=dh, pdf_field_name=pdf.name,
-                                                     streaming_pattern='pull')
-    optimised_communication()
-    dh.run_kernel(stream)
-    dh.swap("pdf", "pdf_tmp")
-
-    if dim == 3:
-        for i in range(gl, domain_size[0]):
-            for j in range(gl, domain_size[1]):
-                for k in range(gl, domain_size[2]):
-                    for f in range(len(stencil)):
-                        assert dh.cpu_arrays['pdf'][i, j, k, f] == pdf_full_communication[i, j, k, f], print(f)
-    else:
-        for i in range(gl, domain_size[0]):
-            for j in range(gl, domain_size[1]):
-                for f in range(len(stencil)):
-                    assert dh.cpu_arrays['pdf'][i, j, f] == pdf_full_communication[i, j, f]

lbmpy_tests/centeredcumulant/test_central_moment_transform.py

-import sympy as sp
-from lbmpy.moments import get_default_moment_set_for_stencil, extract_monomials
-import pytest
-
-from lbmpy.stencils import get_stencil
-from lbmpy.methods.centeredcumulant.centered_cumulants import statistical_quantity_symbol
-from lbmpy.methods.momentbased.moment_transforms import (
-    PdfsToCentralMomentsByMatrix, FastCentralMomentTransform, PdfsToCentralMomentsByShiftMatrix,
-    PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT
-)
-
-
-@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
-def test_forward_transform(stencil):
-    stencil = get_stencil(stencil)
-    dim = len(stencil[0])
-    q = len(stencil)
-    moment_exponents = list(extract_monomials(get_default_moment_set_for_stencil(stencil)))
-    moment_exponents = sorted(moment_exponents, key=sum)
-    pdfs = sp.symbols(f"f_:{q}")
-    rho = sp.Symbol('rho')
-    u = sp.symbols(f"u_:{dim}")
-
-    matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_exponents, rho, u)
-    fast_transform = FastCentralMomentTransform(stencil, moment_exponents, rho, u)
-    shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_exponents, rho, u)
-
-    f_to_k_matrix = matrix_transform.forward_transform(pdfs, moment_symbol_base=PRE_COLLISION_CENTRAL_MOMENT)
-    f_to_k_matrix = f_to_k_matrix.new_without_subexpressions().main_assignments_dict
-
-    f_to_k_fast = fast_transform.forward_transform(pdfs, moment_symbol_base=PRE_COLLISION_CENTRAL_MOMENT)
-    f_to_k_fast = f_to_k_fast.new_without_subexpressions().main_assignments_dict
-
-    f_to_k_shift = shift_transform.forward_transform(pdfs, moment_symbol_base=PRE_COLLISION_CENTRAL_MOMENT,
-                                                     simplification=False)
-    f_to_k_shift = f_to_k_shift.new_without_subexpressions().main_assignments_dict
-
-    for e in moment_exponents:
-        moment_symbol = statistical_quantity_symbol(PRE_COLLISION_CENTRAL_MOMENT, e)
-        rhs_matrix = f_to_k_matrix[moment_symbol].expand()
-        rhs_fast = f_to_k_fast[moment_symbol].expand()
-        rhs_shift = f_to_k_shift[moment_symbol].expand()
-        assert (rhs_matrix - rhs_fast) == 0, f"Mismatch between matrix and fast transform at {moment_symbol}."
-        assert (rhs_matrix - rhs_shift) == 0, f"Mismatch between matrix and shift-matrix transform at {moment_symbol}."
-
-
-@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
-def test_backward_transform(stencil):
-    stencil = get_stencil(stencil)
-    dim = len(stencil[0])
-    q = len(stencil)
-    moment_exponents = list(extract_monomials(get_default_moment_set_for_stencil(stencil)))
-    moment_exponents = sorted(moment_exponents, key=sum)
-    pdfs = sp.symbols(f"f_:{q}")
-    rho = sp.Symbol('rho')
-    u = sp.symbols(f"u_:{dim}")
-
-    matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_exponents, rho, u)
-    fast_transform = FastCentralMomentTransform(stencil, moment_exponents, rho, u)
-    shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_exponents, rho, u)
-
-    k_to_f_matrix = matrix_transform.backward_transform(pdfs, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT)
-    k_to_f_matrix = k_to_f_matrix.new_without_subexpressions().main_assignments_dict
-
-    k_to_f_fast = fast_transform.backward_transform(pdfs, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT)
-    k_to_f_fast = k_to_f_fast.new_without_subexpressions().main_assignments_dict
-
-    k_to_f_shift = shift_transform.backward_transform(pdfs, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT)
-    k_to_f_shift = k_to_f_shift.new_without_subexpressions().main_assignments_dict
-
-    for f in pdfs:
-        rhs_matrix = k_to_f_matrix[f].expand()
-        rhs_fast = k_to_f_fast[f].expand()
-        rhs_shift = k_to_f_shift[f].expand()
-        assert (rhs_matrix - rhs_fast) == 0, f"Mismatch between matrix and fast transform at {f}."
-        assert (rhs_matrix - rhs_shift) == 0, f"Mismatch between matrix and shift-matrix transform at {f}."

lbmpy_tests/centeredcumulant/test_periodic_pipe_flow.ipynb

-%% Cell type:code id: tags:
-
-``` python
-import pytest
-pytest.importorskip('pycuda')
-```
-
-%% Output
-
-    <module 'pycuda' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/pycuda/__init__.py'>
-
-%% Cell type:code id: tags:
-
-``` python
-%load_ext autoreload
-%autoreload 2
-```
-
-%% Cell type:code id: tags:
-
-``` python
-import sympy as sp
-import numpy as np
-from numpy.testing import assert_allclose
-
-from pystencils.session import *
-from lbmpy.session import *
-
-from lbmpy.stencils import get_stencil
-
-from lbmpy.methods.centeredcumulant import CenteredCumulantForceModel
-
-from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
-```
-
-%% Cell type:markdown id: tags:
-
-## Pipe Flow Scenario
-
-%% Cell type:code id: tags:
-
-``` python
-class PeriodicPipeFlow:
-    def __init__(self, method_params, optimization=dict()):
-
-        wall_boundary = NoSlip()
-
-        self.stencil = method_params.get('stencil', get_stencil("D2Q9"))
-        self.streaming_pattern = method_params.get('streaming_pattern', 'pull')
-        self.target = optimization.get('target', 'cpu')
-
-        self.gpu = self.target in ['gpu', 'opencl']
-
-        #   Stencil
-        self.q = len(self.stencil)
-        self.dim = len(self.stencil[0])
-
-        #   Streaming
-        self.inplace = is_inplace(self.streaming_pattern)
-        self.timesteps = get_timesteps(self.streaming_pattern)
-        self.zeroth_timestep = self.timesteps[0]
-
-        #   Domain, Data Handling and PDF fields
-        self.pipe_length = 60
-        self.pipe_radius = 15
-        self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
-        self.periodicity = (True, ) + (False, ) * (self.dim - 1)
-        force = (0.0001, ) + (0.0,) * (self.dim - 1)
-        self.force = method_params.get('force', force)
-
-        self.dh = create_data_handling(domain_size=self.domain_size,
-                                       periodicity=self.periodicity, default_target=self.target)
-
-        self.pdfs = self.dh.add_array('pdfs', self.q)
-        if not self.inplace:
-            self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
-
-        method_params['force'] = self.force
-        optimization['symbolic_field'] = self.pdfs
-
-        if not self.inplace:
-            optimization['symbolic_temporary_field'] = self.pdfs_tmp
-
-        self.lb_collision = create_lb_collision_rule(optimization=optimization, **method_params)
-        self.lb_method = self.lb_collision.method
-
-        self.lb_kernels = []
-        for t in self.timesteps:
-            self.lb_kernels.append(create_lb_function(collision_rule=self.lb_collision,
-                                                      optimization=optimization,
-                                                      timestep=t,
-                                                      **method_params))
-
-        #   Macroscopic Values
-        self.density = 1.0
-        self.density_field = self.dh.add_array('rho', 1)
-        u_x = 0.0
-        self.velocity = (u_x,) * self.dim
-        self.velocity_field = self.dh.add_array('u', self.dim)
-
-        setter = macroscopic_values_setter(
-            self.lb_method, self.density, self.velocity, self.pdfs,
-            streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
-        self.init_kernel = create_kernel(setter, ghost_layers=1, target=self.target).compile()
-
-        self.getter_kernels = []
-        for t in self.timesteps:
-            getter = macroscopic_values_getter(
-                self.lb_method, self.density_field, self.velocity_field, self.pdfs,
-                streaming_pattern=self.streaming_pattern, previous_timestep=t)
-            self.getter_kernels.append(create_kernel(getter, ghost_layers=1, target=self.target).compile())
-
-        #   Periodicity
-        self.periodicity_handler = LBMPeriodicityHandling(
-            self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
-
-        #   Boundary Handling
-        self.wall = wall_boundary
-        self.bh = LatticeBoltzmannBoundaryHandling(
-            self.lb_method, self.dh, self.pdfs.name,
-            streaming_pattern=self.streaming_pattern, target=self.target)
-
-        self.bh.set_boundary(boundary_obj=self.wall, mask_callback=self.mask_callback)
-
-        self.current_timestep = self.zeroth_timestep
-
-    def mask_callback(self, x, y, z=None):
-        y = y - self.pipe_radius
-        z = z - self.pipe_radius if z is not None else 0
-        return np.sqrt(y**2 + z**2) >= self.pipe_radius
-
-    def init(self):
-        self.current_timestep = self.zeroth_timestep
-        self.dh.run_kernel(self.init_kernel)
-
-    def step(self):
-        #   Order matters! First communicate, then boundaries, otherwise
-        #   periodicity handling overwrites reflected populations
-        # Periodicty
-        self.periodicity_handler(self.current_timestep)
-
-        # Boundaries
-        self.bh(prev_timestep=self.current_timestep)
-
-        # Here, the next time step begins
-        self.current_timestep = self.current_timestep.next()
-
-        # LBM Step
-        self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
-
-        # Field Swaps
-        if not self.inplace:
-            self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
-
-        # Macroscopic Values
-        self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
-
-    def run(self, iterations):
-        for _ in range(iterations):
-            self.step()
-
-    @property
-    def velocity_array(self):
-        if self.gpu:
-            self.dh.to_cpu(self.velocity_field.name)
-        return self.dh.gather_array(self.velocity_field.name)
-
-    def get_trimmed_velocity_array(self):
-        if self.gpu:
-            self.dh.to_cpu(self.velocity_field.name)
-        u = np.copy(self.dh.gather_array(self.velocity_field.name))
-        mask = self.bh.get_mask(None, self.wall)
-        for idx in np.ndindex(u.shape[:-1]):
-            if mask[idx] != 0:
-                u[idx] = np.full((self.dim, ), np.nan)
-        return u
-```
-
-%% Cell type:markdown id: tags:
-
-## General Setup
-
-%% Cell type:code id: tags:
-
-``` python
-stencil = get_stencil('D3Q19')
-dim = len(stencil[0])
-target = 'gpu'
-streaming_pattern = 'aa'
-
-optimization = { 'target': target }
-
-viscous_rr = 1.1
-force = (0.0001, ) + (0.0,) * (dim - 1)
-```
-
-%% Cell type:markdown id: tags:
-
-## 1. Reference: SRT Method
-
-%% Cell type:code id: tags:
-
-``` python
-srt_params = {
-    'stencil': stencil,
-    'method': 'srt',
-    'relaxation_rate': viscous_rr,
-    'force_model': 'guo',
-    'force' : force,
-    'streaming_pattern': streaming_pattern
-}
-
-srt_flow = PeriodicPipeFlow(srt_params, optimization)
-srt_flow.init()
-srt_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-srt_u = srt_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(srt_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x7f20b2f75bb0>
-
-
-
-%% Cell type:markdown id: tags:
-
-## 2. Centered Cumulant Method with Implicit Forcing
-
-The default setup of the centered cumulant method does not use a force model to compute the force contributions on all populations, but instead applies the force in cumulant space simply by setting the momentum relaxation rate $\omega_1 = 2$. Due to the half-force shift of the equilibrium input velocity, the first central moments (which correspond to momentum relative to the moving frame of reference) are not zero, but equal to $-F/2$. The relaxation process causes the first central moments to simply change sign:
-$$
-    \kappa_{100}^* = - \kappa_{100}.
-$$
-Thus, $\kappa_{100}^* = F_x /2$. In total, the entire acceleration given by the force is added onto the momentum.
-
-%% Cell type:code id: tags:
-
-``` python
-cm_method_params = {
-    'method' : 'monomial_cumulant',
-    'stencil': stencil,
-    'relaxation_rate': viscous_rr,  #   Specify viscous relaxation rate only
-    'force_model': CenteredCumulantForceModel(force),
-    'force' : force,
-    'streaming_pattern' : streaming_pattern
-}
-
-optimization = {
-    'target': target,
-    'pre_simplification' : True
-}
-
-cm_impl_f_flow = PeriodicPipeFlow(cm_method_params, optimization)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-lb_method = cm_impl_f_flow.lb_method
-assert all(rr == 2 for rr in lb_method.relaxation_rates[1:4])
-assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_impl_f_flow.init()
-cm_impl_f_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_impl_f_u = cm_impl_f_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(cm_impl_f_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x7f20aa846670>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert_allclose(cm_impl_f_u, srt_u, rtol=1, atol=1e-4)
-```
-
-%% Cell type:markdown id: tags:
-
-## 3. Centered Cumulant Method with Explicit Forcing
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.forcemodels import Schiller
-
-cm_method_params_expl_force = {
-    'method' : 'cumulant',
-    'stencil': stencil,
-    'relaxation_rates': [viscous_rr],  #   Specify viscous relaxation rate only
-    'force_model': Schiller(force),
-    'force' : force,
-    'streaming_pattern' : streaming_pattern
-}
-
-optimization = {
-    'target': target,
-    'pre_simplification' : True
-}
-
-cm_expl_f_flow = PeriodicPipeFlow(cm_method_params_expl_force, optimization)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-lb_method = cm_expl_f_flow.lb_method
-assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
-assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_expl_f_flow.init()
-cm_expl_f_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_expl_f_u = cm_expl_f_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(cm_expl_f_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x7f20aa641070>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert_allclose(cm_expl_f_u, srt_u, rtol=1, atol=1e-4)
-assert_allclose(cm_expl_f_u, cm_impl_f_u, rtol=1, atol=1e-4)
-```

lbmpy_tests/phasefield/test_analytical_3phase_model.ipynb

-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from pystencils.fd import evaluate_diffs
-```
-
-%% Cell type:markdown id: tags:
-
-# Analytical checks for 3-Phase model
-
-Here you can inspect the components of the free energy. View only bulk or interface contributions, before and after transformation from $U \rightarrow (\rho, \phi,\psi)$:
-
-%% Cell type:code id: tags:
-
-``` python
-order_parameters = sp.symbols("rho phi psi")
-rho, phi, psi = order_parameters
-F, _ = free_energy_functional_3_phases(include_bulk=True,
-                                       include_interface=True,
-                                       transformed=True,
-                                       expand_derivatives=False)
-F
-```
-
-%% Output
-
-    $$\frac{\alpha^{2} \kappa_{0}}{2} {\partial (\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{1}}{2} {\partial (- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{2}}{2} {\partial \psi}^{2} + \frac{\kappa_{0}}{2} \left(\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(- \frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{1}}{2} \left(- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(\frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{2} \psi^{2}}{2} \left(- \psi + 1\right)^{2}$$
-    
-    
-     2                            2    2                             2    2
-    α ⋅κ₀⋅D(phi/2 - psi/2 + rho/2)    α ⋅κ₁⋅D(-phi/2 - psi/2 + rho/2)    α ⋅κ₂⋅D(p
-    ─────────────────────────────── + ──────────────────────────────── + ─────────
-                   2                                 2                         2
-    
-                         2                  2                   2                2
-              ⎛φ   ψ   ρ⎞  ⎛  φ   ψ   ρ    ⎞       ⎛  φ   ψ   ρ⎞  ⎛φ   ψ   ρ    ⎞
-       2   κ₀⋅⎜─ - ─ + ─⎟ ⋅⎜- ─ + ─ - ─ + 1⎟    κ₁⋅⎜- ─ - ─ + ─⎟ ⋅⎜─ + ─ - ─ + 1⎟
-    si)       ⎝2   2   2⎠  ⎝  2   2   2    ⎠       ⎝  2   2   2⎠  ⎝2   2   2    ⎠
-    ──── + ────────────────────────────────── + ──────────────────────────────────
-                           2                                    2
-    
-    
-    
-           2         2
-       κ₂⋅ψ ⋅(-ψ + 1)
-     + ───────────────
-              2
-
-%% Cell type:markdown id: tags:
-
-### Analytically checking the phase transition profile
-
-%% Cell type:markdown id: tags:
-
-Automatically deriving chemical potential as functional derivative of free energy
-
-%% Cell type:code id: tags:
-
-``` python
-F, _ = free_energy_functional_3_phases(order_parameters)
-
-mu_diff_eq = chemical_potentials_from_free_energy(F, order_parameters)
-mu_diff_eq[0]
-```
-
-%% Output
-
-    $$- \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \rho}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \rho}} + \frac{\kappa_{0} \phi^{3}}{8} - \frac{3 \kappa_{0}}{8} \phi^{2} \psi + \frac{3 \kappa_{0}}{8} \phi^{2} \rho - \frac{3 \kappa_{0}}{8} \phi^{2} + \frac{3 \kappa_{0}}{8} \phi \psi^{2} - \frac{3 \kappa_{0}}{4} \phi \psi \rho + \frac{3 \kappa_{0}}{4} \phi \psi + \frac{3 \kappa_{0}}{8} \phi \rho^{2} - \frac{3 \kappa_{0}}{4} \phi \rho + \frac{\kappa_{0} \phi}{4} - \frac{\kappa_{0} \psi^{3}}{8} + \frac{3 \kappa_{0}}{8} \psi^{2} \rho - \frac{3 \kappa_{0}}{8} \psi^{2} - \frac{3 \kappa_{0}}{8} \psi \rho^{2} + \frac{3 \kappa_{0}}{4} \psi \rho - \frac{\kappa_{0} \psi}{4} + \frac{\kappa_{0} \rho^{3}}{8} - \frac{3 \kappa_{0}}{8} \rho^{2} + \frac{\kappa_{0} \rho}{4} - \frac{\kappa_{1} \phi^{3}}{8} - \frac{3 \kappa_{1}}{8} \phi^{2} \psi + \frac{3 \kappa_{1}}{8} \phi^{2} \rho - \frac{3 \kappa_{1}}{8} \phi^{2} - \frac{3 \kappa_{1}}{8} \phi \psi^{2} + \frac{3 \kappa_{1}}{4} \phi \psi \rho - \frac{3 \kappa_{1}}{4} \phi \psi - \frac{3 \kappa_{1}}{8} \phi \rho^{2} + \frac{3 \kappa_{1}}{4} \phi \rho - \frac{\kappa_{1} \phi}{4} - \frac{\kappa_{1} \psi^{3}}{8} + \frac{3 \kappa_{1}}{8} \psi^{2} \rho - \frac{3 \kappa_{1}}{8} \psi^{2} - \frac{3 \kappa_{1}}{8} \psi \rho^{2} + \frac{3 \kappa_{1}}{4} \psi \rho - \frac{\kappa_{1} \psi}{4} + \frac{\kappa_{1} \rho^{3}}{8} - \frac{3 \kappa_{1}}{8} \rho^{2} + \frac{\kappa_{1} \rho}{4}$$
-       2                 2                 2                 2                 2
-      α ⋅κ₀⋅D(D(phi))   α ⋅κ₀⋅D(D(psi))   α ⋅κ₀⋅D(D(rho))   α ⋅κ₁⋅D(D(phi))   α ⋅κ
-    - ─────────────── + ─────────────── - ─────────────── + ─────────────── + ────
-             4                 4                 4                 4
-    
-                   2                    3         2           2           2
-    ₁⋅D(D(psi))   α ⋅κ₁⋅D(D(rho))   κ₀⋅φ    3⋅κ₀⋅φ ⋅ψ   3⋅κ₀⋅φ ⋅ρ   3⋅κ₀⋅φ    3⋅κ₀
-    ─────────── - ─────────────── + ───── - ───────── + ───────── - ─────── + ────
-       4                 4            8         8           8          8
-    
-        2                                   2                         3         2
-    ⋅φ⋅ψ    3⋅κ₀⋅φ⋅ψ⋅ρ   3⋅κ₀⋅φ⋅ψ   3⋅κ₀⋅φ⋅ρ    3⋅κ₀⋅φ⋅ρ   κ₀⋅φ   κ₀⋅ψ    3⋅κ₀⋅ψ ⋅
-    ───── - ────────── + ──────── + ───────── - ──────── + ──── - ───── + ────────
-    8           4           4           8          4        4       8         8
-    
-              2           2                         3         2              3
-    ρ   3⋅κ₀⋅ψ    3⋅κ₀⋅ψ⋅ρ    3⋅κ₀⋅ψ⋅ρ   κ₀⋅ψ   κ₀⋅ρ    3⋅κ₀⋅ρ    κ₀⋅ρ   κ₁⋅φ    3
-    ─ - ─────── - ───────── + ──────── - ──── + ───── - ─────── + ──── - ───── - ─
-           8          8          4        4       8        8       4       8
-    
-         2           2           2           2                                   2
-    ⋅κ₁⋅φ ⋅ψ   3⋅κ₁⋅φ ⋅ρ   3⋅κ₁⋅φ    3⋅κ₁⋅φ⋅ψ    3⋅κ₁⋅φ⋅ψ⋅ρ   3⋅κ₁⋅φ⋅ψ   3⋅κ₁⋅φ⋅ρ
-    ──────── + ───────── - ─────── - ───────── + ────────── - ──────── - ─────────
-       8           8          8          8           4           4           8
-    
-                             3         2           2           2
-       3⋅κ₁⋅φ⋅ρ   κ₁⋅φ   κ₁⋅ψ    3⋅κ₁⋅ψ ⋅ρ   3⋅κ₁⋅ψ    3⋅κ₁⋅ψ⋅ρ    3⋅κ₁⋅ψ⋅ρ   κ₁⋅ψ
-     + ──────── - ──── - ───── + ───────── - ─────── - ───────── + ──────── - ────
-          4        4       8         8          8          8          4        4
-    
-           3         2
-       κ₁⋅ρ    3⋅κ₁⋅ρ    κ₁⋅ρ
-     + ───── - ─────── + ────
-         8        8       4
-
-%% Cell type:markdown id: tags:
-
-Checking if expected profile is a solution of the differential equation:
-
-%% Cell type:code id: tags:
-
-``` python
-x = sp.symbols("x")
-expectedProfile = analytic_interface_profile(x)
-expectedProfile
-```
-
-%% Output
-
-
-    $$\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}$$
-        ⎛ x ⎞
-    tanh⎜───⎟
-        ⎝2⋅α⎠   1
-    ───────── + ─
-        2       2
-
-%% Cell type:markdown id: tags:
-
-Checking a phase transition from $C_0$ to $C_2$. This means that $\rho=1$ while $phi$ and $psi$ are the analytical profile or 1-analytical profile
-
-%% Cell type:code id: tags:
-
-``` python
-for eq in mu_diff_eq:
-    eq = eq.subs({rho: 1,
-                  phi: 1 - expectedProfile,
-                  psi: expectedProfile})
-    eq = evaluate_diffs(eq, x).expand()
-    assert eq == 0
-```
-
-%% Cell type:markdown id: tags:
-
-### Checking the surface tensions parameters
-
-%% Cell type:code id: tags:
-
-``` python
-F, _ = free_energy_functional_3_phases(order_parameters)
-F = expand_diff_linear(F, functions=order_parameters)  # expand derivatives using product rule
-two_phase_free_energy = F.subs({rho: 1,
-                                phi: 1 - expectedProfile,
-                                psi: expectedProfile})
-
-two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
-two_phase_free_energy
-```
-
-%% Output
-
-
-    $$\frac{\kappa_{0} + \kappa_{2}}{16 \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}$$
-       κ₀ + κ₂
-    ─────────────
-           4⎛ x ⎞
-    16⋅cosh ⎜───⎟
-            ⎝2⋅α⎠
-
-%% Cell type:code id: tags:
-
-``` python
-gamma = cosh_integral(two_phase_free_energy, x)
-gamma
-```
-
-%% Output
-
-
-    $$\frac{\alpha}{6} \left(\kappa_{0} + \kappa_{2}\right)$$
-    α⋅(κ₀ + κ₂)
-    ───────────
-         6
-
-%% Cell type:code id: tags:
-
-``` python
-alpha, k0, k2 = sp.symbols("alpha, kappa_0, kappa_2")
-assert gamma == alpha/6 * (k0 + k2)
-```

lbmpy_tests/phasefield/test_analytical_n_phase.ipynb

-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from functools import partial
-```
-
-%% Cell type:markdown id: tags:
-
-# Analytical checks for N-Phase model
-
-%% Cell type:markdown id: tags:
-
-### Formulation of free energy
-
-%% Cell type:markdown id: tags:
-
-In the next cell you can inspect the formulation of the free energy functional.
-Bulk and interface terms can be viewed separately.
-
-%% Cell type:code id: tags:
-
-``` python
-num_phases = 3
-order_params = symbolic_order_parameters(num_symbols=num_phases-1)
-f2 = partial(n_phases_correction_function, beta=1)
-F = free_energy_functional_n_phases(order_parameters=order_params,
-                                    include_interface=False, f2=f2,
-                                    include_bulk=True )
-F
-```
-
-%% Output
-
-
-    $$\frac{3}{2 \alpha} \left(\tau_{0 1} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} - \begin{cases} - \left(\phi_{0} + \phi_{1}\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} < 0 \\- \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} > 1 \\\left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{0 2} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{1} + 1\right)^{2} & \text{for}\: - \phi_{1} + 1 < 0 \\- \phi_{1}^{2} & \text{for}\: - \phi_{1} + 1 > 1 \\\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{1 2} \left(\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{0} + 1\right)^{2} & \text{for}\: - \phi_{0} + 1 < 0 \\- \phi_{0}^{2} & \text{for}\: - \phi_{0} + 1 > 1 \\\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} & \text{otherwise} \end{cases}\right)\right)$$
-      ⎛     ⎛                                  ⎛⎧                 2
-      ⎜     ⎜                                  ⎜⎪       -(φ₀ + φ₁)           for φ
-      ⎜     ⎜                                  ⎜⎪
-      ⎜     ⎜  2          2     2          2   ⎜⎪                    2
-    3⋅⎜τ₀ ₁⋅⎜φ₀ ⋅(-φ₀ + 1)  + φ₁ ⋅(-φ₁ + 1)  - ⎜⎨     -(-φ₀ - φ₁ + 1)        for φ
-      ⎜     ⎜                                  ⎜⎪
-      ⎜     ⎜                                  ⎜⎪         2               2
-      ⎜     ⎜                                  ⎜⎪(φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)      ot
-      ⎝     ⎝                                  ⎝⎩
-    ──────────────────────────────────────────────────────────────────────────────
-    
-    
-              ⎞⎞        ⎛                                              ⎛⎧
-    ₀ + φ₁ < 0⎟⎟        ⎜                                              ⎜⎪ -(-φ₁ +
-              ⎟⎟        ⎜                                              ⎜⎪
-              ⎟⎟        ⎜  2          2            2               2   ⎜⎪        2
-    ₀ + φ₁ > 1⎟⎟ + τ₀ ₂⋅⎜φ₀ ⋅(-φ₀ + 1)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  - ⎜⎨     -φ₁
-              ⎟⎟        ⎜                                              ⎜⎪
-              ⎟⎟        ⎜                                              ⎜⎪  2
-    herwise   ⎟⎟        ⎜                                              ⎜⎪φ₁ ⋅(-φ₁
-              ⎠⎠        ⎝                                              ⎝⎩
-    ──────────────────────────────────────────────────────────────────────────────
-                                                            2⋅α
-    
-      2                   ⎞⎞        ⎛
-    1)     for -φ₁ + 1 < 0⎟⎟        ⎜
-                          ⎟⎟        ⎜
-                          ⎟⎟        ⎜  2          2            2               2
-           for -φ₁ + 1 > 1⎟⎟ + τ₁ ₂⋅⎜φ₁ ⋅(-φ₁ + 1)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  -
-                          ⎟⎟        ⎜
-        2                 ⎟⎟        ⎜
-    + 1)      otherwise   ⎟⎟        ⎜
-                          ⎠⎠        ⎝
-    ──────────────────────────────────────────────────────────────────────────────
-    
-    
-     ⎛⎧           2                   ⎞⎞⎞
-     ⎜⎪ -(-φ₀ + 1)     for -φ₀ + 1 < 0⎟⎟⎟
-     ⎜⎪                               ⎟⎟⎟
-     ⎜⎪        2                      ⎟⎟⎟
-     ⎜⎨     -φ₀        for -φ₀ + 1 > 1⎟⎟⎟
-     ⎜⎪                               ⎟⎟⎟
-     ⎜⎪  2          2                 ⎟⎟⎟
-     ⎜⎪φ₀ ⋅(-φ₀ + 1)      otherwise   ⎟⎟⎟
-     ⎝⎩                               ⎠⎠⎠
-    ─────────────────────────────────────
-    
-
-%% Cell type:markdown id: tags:
-
-### Analytically checking the phase transition profile
-
-First we define the order parameters and free energy, including bulk and interface terms:
-
-%% Cell type:code id: tags:
-
-``` python
-F = free_energy_functional_n_phases(order_parameters=symbolic_order_parameters(num_symbols=num_phases-1))
-```
-
-%% Cell type:markdown id: tags:
-
-Then we automatically derive the differential equations for the chemial potential $\mu$
-
-%% Cell type:code id: tags:
-
-``` python
-mu_diff_eq = chemical_potentials_from_free_energy(F, order_params)
-
-# there is one equation less than phases
-assert len(mu_diff_eq) == num_phases - 1
-
-# show the first one
-mu_diff_eq[0]
-```
-
-%% Output
-
-    $$3 \alpha \tau_{0 1} {\partial {\partial \phi_{1}}} - 6 \alpha \tau_{0 2} {\partial {\partial \phi_{0}}} - 3 \alpha \tau_{0 2} {\partial {\partial \phi_{1}}} - 3 \alpha \tau_{1 2} {\partial {\partial \phi_{1}}} + \frac{12 \phi_{0}^{3}}{\alpha} \tau_{0 2} - \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 1} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 2} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{1 2} - \frac{18 \phi_{0}^{2}}{\alpha} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 1} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{1 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 1} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{1 2} + \frac{6 \phi_{0}}{\alpha} \tau_{0 2} - \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 1} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 2} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{1 2} + \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 1} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 2} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{1 2} - \frac{3 \phi_{1}}{\alpha} \tau_{0 1} + \frac{3 \phi_{1}}{\alpha} \tau_{0 2} + \frac{3 \phi_{1}}{\alpha} \tau_{1 2}$$
-    
-    
-    3⋅α⋅τ₀ ₁⋅D(D(phi_1)) - 6⋅α⋅τ₀ ₂⋅D(D(phi_0)) - 3⋅α⋅τ₀ ₂⋅D(D(phi_1)) - 3⋅α⋅τ₁ ₂⋅
-    
-    
-                       3             2                2                2
-                  12⋅φ₀ ⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₀ ₁   18⋅φ₀ ⋅φ₁⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₁ ₂
-    D(D(phi_1)) + ─────────── - ────────────── + ────────────── + ────────────── -
-                       α              α                α                α
-    
-          2                2                2                2
-     18⋅φ₀ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₀ ₁   18⋅φ₀⋅φ₁ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₁ ₂   18⋅φ₀⋅φ₁⋅τ₀
-     ─────────── - ────────────── + ────────────── + ────────────── + ────────────
-          α              α                α                α                α
-    
-                                                        3            3
-    ₁   18⋅φ₀⋅φ₁⋅τ₀ ₂   18⋅φ₀⋅φ₁⋅τ₁ ₂   6⋅φ₀⋅τ₀ ₂   6⋅φ₁ ⋅τ₀ ₁   6⋅φ₁ ⋅τ₀ ₂   6⋅φ₁
-    ─ - ───────────── - ───────────── + ───────── - ────────── + ────────── + ────
-              α               α             α           α            α
-    
-    3            2            2            2
-     ⋅τ₁ ₂   9⋅φ₁ ⋅τ₀ ₁   9⋅φ₁ ⋅τ₀ ₂   9⋅φ₁ ⋅τ₁ ₂   3⋅φ₁⋅τ₀ ₁   3⋅φ₁⋅τ₀ ₂   3⋅φ₁⋅τ
-    ────── + ────────── - ────────── - ────────── - ───────── + ───────── + ──────
-    α            α            α            α            α           α           α
-    
-    
-    ₁ ₂
-    ───
-    
-
-%% Cell type:markdown id: tags:
-
-Next we check, that the $tanh$ is indeed a solution of this equation...
-
-%% Cell type:code id: tags:
-
-``` python
-x = sp.symbols("x")
-expected_profile = analytic_interface_profile(x)
-expected_profile
-```
-
-%% Output
-
-
-    $$\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}$$
-        ⎛ x ⎞
-    tanh⎜───⎟
-        ⎝2⋅α⎠   1
-    ───────── + ─
-        2       2
-
-%% Cell type:markdown id: tags:
-
-... by inserting it for $\phi_0$, and setting all other order parameters to zero
-
-%% Cell type:code id: tags:
-
-``` python
-# zero other parameters
-diff_eq = mu_diff_eq[0].subs({p: 0 for p in order_params[1:]})
-
-# insert analytical solution
-diff_eq = diff_eq.subs(order_params[0], expected_profile)
-diff_eq.factor()
-```
-
-%% Output
-
-    $$- \frac{3 \tau_{0 2}}{2 \alpha} \left(4 \alpha^{2} {\partial {\partial (\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}) }} - \tanh^{3}{\left (\frac{x}{2 \alpha} \right )} + \tanh{\left (\frac{x}{2 \alpha} \right )}\right)$$
-            ⎛   2                                       3⎛ x ⎞       ⎛ x ⎞⎞
-    -3⋅τ₀ ₂⋅⎜4⋅α ⋅D(D(tanh(x/(2*alpha))/2 + 1/2)) - tanh ⎜───⎟ + tanh⎜───⎟⎟
-            ⎝                                            ⎝2⋅α⎠       ⎝2⋅α⎠⎠
-    ────────────────────────────────────────────────────────────────────────
-                                      2⋅α
-
-%% Cell type:markdown id: tags:
-
-finally the differentials have to be evaluated...
-
-%% Cell type:code id: tags:
-
-``` python
-from pystencils.fd import evaluate_diffs
-diff_eq = evaluate_diffs(diff_eq, x)
-diff_eq
-```
-
-%% Output
-
-
-    $$\frac{12}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{3} - \frac{18}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{2} + \frac{6}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right) + \frac{3 \tau_{0 2}}{2 \alpha} \left(- \tanh^{2}{\left (\frac{x}{2 \alpha} \right )} + 1\right) \tanh{\left (\frac{x}{2 \alpha} \right )}$$
-                           3                          2
-            ⎛    ⎛ x ⎞    ⎞            ⎛    ⎛ x ⎞    ⎞           ⎛    ⎛ x ⎞    ⎞
-            ⎜tanh⎜───⎟    ⎟            ⎜tanh⎜───⎟    ⎟           ⎜tanh⎜───⎟    ⎟
-            ⎜    ⎝2⋅α⎠   1⎟            ⎜    ⎝2⋅α⎠   1⎟           ⎜    ⎝2⋅α⎠   1⎟
-    12⋅τ₀ ₂⋅⎜───────── + ─⎟    18⋅τ₀ ₂⋅⎜───────── + ─⎟    6⋅τ₀ ₂⋅⎜───────── + ─⎟
-            ⎝    2       2⎠            ⎝    2       2⎠           ⎝    2       2⎠
-    ──────────────────────── - ──────────────────────── + ────────────────────── +
-               α                          α                         α
-    
-    
-    
-    
-            ⎛      2⎛ x ⎞    ⎞     ⎛ x ⎞
-     3⋅τ₀ ₂⋅⎜- tanh ⎜───⎟ + 1⎟⋅tanh⎜───⎟
-            ⎝       ⎝2⋅α⎠    ⎠     ⎝2⋅α⎠
-     ───────────────────────────────────
-                     2⋅α
-
-%% Cell type:markdown id: tags:
-
-.. and the result simplified...
-
-%% Cell type:code id: tags:
-
-``` python
-assert diff_eq.expand() == 0
-diff_eq.expand()
-```
-
-%% Output
-
-
-    $$0$$
-    0
-
-%% Cell type:markdown id: tags:
-
-...and indeed the expected tanh profile satisfies this differential equation.
-
-Next lets check for the interface between phase 0 and phase 1:
-
-%% Cell type:code id: tags:
-
-``` python
-for diff_eq in mu_diff_eq:
-    eq = diff_eq.subs({order_params[0]: expected_profile,
-                       order_params[1]: 1 - expected_profile})
-    assert evaluate_diffs(eq, x).expand() == 0
-```
-
-%% Cell type:markdown id: tags:
-
-### Checking the surface tensions parameters
-
-Computing the excess free energy per unit area of an interface between two phases.
-This should be exactly the surface tension parameter.
-
-%% Cell type:code id: tags:
-
-``` python
-order_params = symbolic_order_parameters(num_symbols=num_phases-1)
-F = free_energy_functional_n_phases(order_parameters=order_params)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-two_phase_free_energy = F.subs({order_params[0]: expected_profile,
-                                order_params[1]: 1 - expected_profile})
-
-# Evaluate differentials and simplification
-two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-excess_free_energy = sp.Integral(two_phase_free_energy, x)
-excess_free_energy
-```
-
-%% Output
-
-
-    $$\int \frac{3 \tau_{0 1}}{8 \alpha \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}\, dx$$
-    ⌠
-    ⎮     3⋅τ₀ ₁
-    ⎮ ────────────── dx
-    ⎮         4⎛ x ⎞
-    ⎮ 8⋅α⋅cosh ⎜───⎟
-    ⎮          ⎝2⋅α⎠
-    ⌡
-
-%% Cell type:markdown id: tags:
-
-Sympy cannot integrate this automatically - help with a manual substitution $\frac{x}{2\alpha} \rightarrow u$
-
-%% Cell type:code id: tags:
-
-``` python
-coshTerm = list(excess_free_energy.atoms(sp.cosh))[0]
-transformed_int = excess_free_energy.transform(coshTerm.args[0], sp.Symbol("u", real=True))
-transformed_int
-```
-
-%% Output
-
-
-    $$\int \frac{3 \tau_{0 1}}{4 \cosh^{4}{\left (u \right )}}\, du$$
-    ⌠
-    ⎮   3⋅τ₀ ₁
-    ⎮ ────────── du
-    ⎮       4
-    ⎮ 4⋅cosh (u)
-    ⌡
-
-%% Cell type:markdown id: tags:
-
-Now the integral can be done:
-
-%% Cell type:code id: tags:
-
-``` python
-result = sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
-assert result == symmetric_symbolic_surface_tension(0,1)
-result
-```
-
-%% Output
-
-
-    $$\tau_{0 1}$$
-    τ₀ ₁

lbmpy_tests/phasefield/test_n_phase_boyer_noncoupled.ipynb

-%% Cell type:code id: tags:
-
-``` python
-import pytest
-pytest.importorskip('pycuda')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.n_phase_boyer import *
-from lbmpy.phasefield.kerneleqs import *
-from lbmpy.phasefield.contact_angle_circle_fitting import *
-from scipy.ndimage.filters import gaussian_filter
-from pystencils.simp import sympy_cse_on_assignment_list
-one = sp.sympify(1)
-
-import pyximport
-pyximport.install(language_level=3)
-from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
-```
-
-%% Cell type:markdown id: tags:
-
-# Simulation arbitrary surface tension case
-
-%% Cell type:code id: tags:
-
-``` python
-n = 4
-dx, dt = 1, 1
-mobility = 2e-3
-domain_size = (150, 150)
-ε = one * 4
-penalty_factor = 0
-stabilization_factor = 10
-
-κ = (one,  one/2, one/3, one/4)
-sigma_factor = one / 15
-σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-#σ
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_target='gpu')
-c = dh.add_array('c', values_per_cell=n)
-c_tmp = dh.add_array_like('c_tmp', 'c')
-
-μ = dh.add_array('mu', values_per_cell=n)
-
-cvec = c.center_vector
-μvec = μ.center_vector
-```
-
-%% Cell type:code id: tags:
-
-``` python
-α, _ = diffusion_coefficients(σ)
-
-f = lambda c: c**2 * ( 1 - c ) **2
-a, b = compute_ab(f)
-
-capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
-
-f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
-f_if = free_energy_interfacial(cvec, σ, a, ε)
-f = f_bulk + f_if
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#f_bulk
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ_assignments = mu_kernel(f, cvec, c, μ)
-μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
-μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-def lapl(e):
-    return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-rhs = α * μvec
-discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
-                   for i, rhs_i in enumerate(rhs)]
-c_assignments = [Assignment(lhs, rhs)
-                 for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#c_assignments
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ_sync = dh.synchronization_function(μ.name)
-c_sync = dh.synchronization_function(c.name)
-optimization = {'cpu_openmp': 4, 'cpu_vectorize_info': None}
-μ_kernel = create_kernel(μ_assignments, target=dh.default_target, **optimization).compile()
-c_kernel = create_kernel(c_assignments, target=dh.default_target, **optimization).compile()
-
-def set_c(slice_obj, values):
-    for block in dh.iterate(slice_obj):
-        arr = block[c.name]
-        arr[..., : ] = values
-
-def smooth():
-    for block in dh.iterate(ghost_layers=True):
-        c_arr = block[c.name]
-        for i in range(n):
-            gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
-
-def time_loop(steps):
-    dh.all_to_gpu()
-    for t in range(steps):
-        c_sync()
-        dh.run_kernel(μ_kernel)
-        μ_sync()
-        dh.run_kernel(c_kernel)
-        dh.swap(c.name, c_tmp.name)
-        #simplex_projection_2d(dh.cpu_arrays[c.name])
-    dh.all_to_cpu()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-set_c(make_slice[:, :], [0, 0, 0, 0])
-set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
-set_c(make_slice[:, :0.5], [0, 1, 0, 0])
-set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
-smooth()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#dh.load_all('n_phases_state_size200_stab10.npz')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-plt.phase_plot(dh.gather_array(c.name))
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 146.44269023807928, \  117.81813928465394, \  95.73917047726678\right]$
-    [146.44269023807928, 117.81813928465394, 95.73917047726678]
-
-%% Cell type:code id: tags:
-
-``` python
-import time
-for i in range(10):
-    start = time.perf_counter()
-    time_loop(1_000)
-    end = time.perf_counter()
-
-    try:
-        print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
-    except Exception:
-        print(i, end - start, "none found")
-```
-
-%% Output
-
-    0 0.30607624700132874 [83.97888710273061, 100.48794346625529, 175.5331694310141]
-    1 0.2600655169990205 [83.73094685534376, 100.65854574856168, 175.6105073960945]
-    2 0.2601136189987301 [83.49914818603683, 100.82173327601079, 175.67911853795226]
-    3 0.25987518599868054 [83.31519592224448, 100.94468140501989, 175.74012267273554]
-    4 0.2651959220020217 [83.14239972296966, 101.06100094405181, 175.79659933297853]
-    5 0.25910847799968906 [82.984481834461, 101.16731750637399, 175.8482006591651]
-    6 0.259863024999504 [82.84781128433397, 101.2570276449976, 175.89516107066834]
-    7 0.2606479199966998 [82.7456965110742, 101.31687551766585, 175.93742797125986]
-    8 0.25991897900166805 [82.67010885583116, 101.35099855297112, 175.97889259119805]
-    9 0.2590353729974595 [75.9000280154447, 108.9652166787719, 175.1347553057833]
-
-%% Cell type:code id: tags:
-
-``` python
-plt.subplot(1,3,1)
-t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
-plt.plot(t);
-
-plt.subplot(1,3,2)
-plt.phase_plot(dh.gather_array(c.name), linewidth=1)
-
-plt.subplot(1,3,3)
-plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
-plt.colorbar();
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert not np.isnan(dh.max(c.name))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
-plt.hlines(0.5, 0, 30)
-plt.plot(t);
-```
-
-%% Output
-
-

lbmpy_tests/phasefield/test_n_phase_paper_boyer_lbm.ipynb

-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from lbmpy.phasefield.n_phase_boyer import *
-from lbmpy.phasefield.kerneleqs import *
-
-from pystencils.fd.spatial import discretize_spatial
-from pystencils.simp import sympy_cse_on_assignment_list
-from lbmpy.phasefield.cahn_hilliard_lbm import *
-from pystencils.fd.derivation import *
-
-one = sp.sympify(1)
-```
-
-%% Cell type:markdown id: tags:
-
-# Chemical potential
-
-Current state:
-
-- not stable (yet)
-- without LB coupling the model is stable
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (100, 100)
-
-n = 4
-c = sp.symbols(f"c_:{n}")
-simple_potential = False
-omega_h = 1.3
-
-if simple_potential:
-    f = free_energy_functional_n_phases_penalty_term(c, interface_width=1, kappa=(0.05, 0.05/2, 0.05/4))
-else:
-    ε = one * 4
-    mobility = one * 2 / 1000
-    κ = (one,  one/2, one/3, one/4)
-    sigma_factor = one / 15
-    σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-
-    α, _ = diffusion_coefficients(σ)
-    f_b, f_if, mu_b, mu_if = chemical_potential_n_phase_boyer(c, ε, σ, one,
-                                                              assume_nonnegative=True, zero_threshold=1e-10)
-```
-
-%% Cell type:markdown id: tags:
-
-# Data Setup
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_ghost_layers=2)
-p_linearization = symmetric_tensor_linearization(dh.dim)
-
-c_field = dh.add_array("c", values_per_cell=n)
-rho_field = dh.add_array("rho")
-mu_field = dh.add_array("mu", values_per_cell=n, latex_name="\\mu")
-pressure_field = dh.add_array("P", values_per_cell=len(p_linearization))
-force_field = dh.add_array("F", values_per_cell=dh.dim)
-u_field = dh.add_array("u", values_per_cell=dh.dim)
-
-# Distribution functions for each order parameter
-pdf_field = []
-pdf_dst_field = []
-for i in range(n):
-    pdf_field_local = dh.add_array("f%d" % i, values_per_cell=9) # 9 for D2Q9
-    pdf_dst_field_local = dh.add_array("f%d_dst"%i, values_per_cell=9, latex_name="f%d_{dst}" % i)
-    pdf_field.append(pdf_field_local)
-    pdf_dst_field.append(pdf_dst_field_local)
-
-# Distribution functions for the hydrodynamics
-pdf_hydro_field = dh.add_array("fh", values_per_cell=9)
-pdf_hydro_dst_field = dh.add_array("fh_dst", values_per_cell=9, latex_name="fh_{dst}")
-```
-
-%% Cell type:markdown id: tags:
-
-### μ-Kernel
-
-%% Cell type:code id: tags:
-
-``` python
-if simple_potential:
-    mu_assignments = mu_kernel(f, c, c_field, mu_field, discretization='isotropic')
-else:
-    mu_subs = {a: b for a, b in zip(c, c_field.center_vector)}
-    mu_if_discrete = [discretize_spatial(e.subs(mu_subs), dx=1, stencil='isotropic') for e in mu_if]
-    mu_b_discrete = [e.subs(mu_subs) for e in mu_b]
-
-    mu_assignments = [Assignment(mu_field(i),
-                                 mu_if_discrete[i] + mu_b_discrete[i]) for i in range(n)]
-
-    mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
-
-μ_kernel = create_kernel(mu_assignments).compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-second_neighbor_stencil = [(i, j)
-                           for i in (-2, -1, 0, 1, 2)
-                           for j in (-2, -1, 0, 1, 2)
-                          ]
-x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
-x_diff.set_weight((2, 0), sp.Rational(1, 10))
-x_diff.assume_symmetric(0, anti_symmetric=True)
-x_diff.assume_symmetric(1)
-x_diff_stencil = x_diff.get_stencil(isotropic=True)
-
-y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
-y_diff.set_weight((0, 2), sp.Rational(1, 10))
-y_diff.assume_symmetric(1, anti_symmetric=True)
-y_diff.assume_symmetric(0)
-y_diff_stencil = y_diff.get_stencil(isotropic=True)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ = mu_field
-μ_discretization = {}
-
-for i in range(n):
-    μ_discretization.update({Diff(μ(i), 0): x_diff_stencil.apply(μ(i)),
-                             Diff(μ(i), 1): y_diff_stencil.apply(μ(i))})
-force_rhs = force_from_phi_and_mu(order_parameters=c_field.center_vector, dim=dh.dim, mu=mu_field.center_vector)
-force_rhs = force_rhs.subs(μ_discretization)
-force_assignments = [Assignment(force_field(i), force_rhs[i]) for i in range(dh.dim)]
-
-force_kernel = create_kernel(force_assignments).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-## Lattice Boltzmann kernels
-
-For each order parameter a Cahn-Hilliard LB method computes the time evolution.
-
-%% Cell type:code id: tags:
-
-``` python
-if simple_potential:
-    mu_alpha = mu_field.center_vector
-else:
-    mu_alpha = mobility * α * mu_field.center_vector
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# Defining the Cahn-Hilliard Collision assignments
-ch_kernels = []
-
-for i in range(n):
-    ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu_alpha[i],
-                                        relaxation_rate=1.0, gamma=1.0)
-    kernel = create_lb_function(lb_method=ch_method,
-                                velocity_input=u_field.center_vector,
-                                compressible=True,
-                                output={'density': c_field(i)},
-                                optimization={"symbolic_field":pdf_field[i],
-                                              "symbolic_temporary_field": pdf_dst_field[i]})
-
-    ch_kernels.append(kernel)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-hydro_lbm = create_lb_function(relaxation_rate=omega_h, force=force_field,
-                               compressible=True,
-                               optimization={"symbolic_field": pdf_hydro_field,
-                                             "symbolic_temporary_field": pdf_hydro_dst_field},
-                               output={'velocity': u_field}
-                               )
-#hydro_lbm.update_rule
-```
-
-%% Cell type:markdown id: tags:
-
-# Initialization
-
-%% Cell type:code id: tags:
-
-``` python
-def set_c(slice_obj, values):
-    for block in dh.iterate(slice_obj):
-        arr = block[c_field.name]
-        arr[..., : ] = values
-
-def init():
-    dh.fill(u_field.name, 0)
-    dh.fill(mu_field.name, 0)
-    dh.fill(force_field.name, 0)
-
-    set_c(make_slice[:, :], [0, 0, 0, 0])
-    set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
-    set_c(make_slice[:, :0.5], [0, 1, 0, 0])
-    set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
-
-pdf_sync_fns = dh.synchronization_function([f.name for f in pdf_field])
-hydro_sync_fn=dh.synchronization_function([pdf_hydro_field.name])
-c_sync_fn = dh.synchronization_function([c_field.name])
-mu_sync_fn = dh.synchronization_function([mu_field.name])
-
-def time_loop(steps):
-    for i in range(steps):
-        c_sync_fn()
-        dh.run_kernel(μ_kernel)
-
-        mu_sync_fn()
-        dh.run_kernel(force_kernel)
-
-        hydro_sync_fn()
-        dh.run_kernel(hydro_lbm)
-        dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
-
-        pdf_sync_fns()
-        for i in range(n):
-            dh.run_kernel(ch_kernels[i])
-            dh.swap(pdf_field[i].name,pdf_dst_field[i].name)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-plt.phase_plot(dh.gather_array(c_field.name))
-```
-
-%% Output
-
-    /local/bauer/miniconda3/envs/pystencils_dev/lib/python3.7/site-packages/matplotlib/contour.py:1243: UserWarning: No contour levels were found within the data range.
-      warnings.warn("No contour levels were found"
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-time_loop(10)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#plt.scalar_field(dh.cpu_arrays[mu_field.name][..., 0])
-plt.scalar_field(dh.cpu_arrays[u_field.name][..., 0])
-plt.colorbar();
-```
-
-%% Output
-
-