Module `phi.physics.fluid`

Functions for simulating incompressible fluids, both grid-based and particle-based.

The main function for incompressible fluids (Eulerian as well as FLIP / PIC) is make_incompressible() which removes the divergence of a velocity field.

Functions

def apply_boundary_conditions(velocity: phi.field._field.Field, obstacles: Obstacle)

Enforces velocities boundary conditions on a velocity grid. Cells inside obstacles will get their velocity from the obstacle movement. Cells outside far away will be unaffected.

Args

velocity: Velocity Grid. obstacles: Obstacle or Geometry or tuple/list thereof to specify boundary conditions inside the domain.

Returns

Velocity of same type as velocity

def boundary_push(particles: , obstacles: tuple, separation: float = 0.5) ‑>

Enforces boundary conditions by correcting possible errors of the advection step and shifting particles out of obstacles or back into the domain.

Args

particles: PointCloud holding particle positions as elements
obstacles: List of Obstacle or Geometry objects where any particles inside should get shifted outwards
separation: Minimum distance between particles and domain boundary / obstacle surface after particles have been shifted.

Returns

PointCloud where all particles are inside the domain / outside of obstacles.

def incompressible_rk4(pde: Callable, velocity: phi.field._field.Field, pressure: phi.field._field.Field, dt, pressure_order=4, pressure_solve=CG with tolerance None (rel), None (abs), max_iterations=1000, **pde_aux_kwargs)

Implements the 4th-order Runge-Kutta time advancement scheme for incompressible vector fields. This approach is inspired by Kampanis et. al., 2006 and incorporates the pressure treatment into the time step.

Args

pde: Momentum equation. Function that computes all PDE terms not related to pressure, e.g. diffusion, advection, external forces.
velocity: Velocity grid at time t.
pressure: Pressure at time t.
dt: Time increment to integrate.
pressure_order: spatial order for derivative computations. For Higher-order schemes, the laplace operation is not conducted with a stencil exactly corresponding to the one used in divergence calculations but a smaller one instead. While this disrupts the formal correctness of the method it only induces insignificant errors and yields considerable performance gains. supported: explicit 2/4th order - implicit 6th order (obstacles are only supported with explicit 2nd order)
pressure_solve: Solve object specifying method and tolerances for the implicit pressure solve.
**pde_aux_kwargs: Auxiliary arguments for pde. These are considered constant over time.

Returns

velocity: Velocity at time t+dt, same type as velocity.
pressure: Pressure grid at time t+dt, CenteredGrid.

def make_incompressible(velocity: phi.field._field.Field, obstacles: Obstacle = (), solve: phiml.math._optimize.Solve = auto with tolerance None (rel), None (abs), max_iterations=1000, active: = None, order: int = 2, correct_skew=False, wide_stencil: bool = None) ‑> Tuple[phi.field._field.Field, phi.field._field.Field]

Projects the given velocity field by solving for the pressure and subtracting its spatial_gradient.

This method is similar to :func:field.divergence_free() but differs in how the boundary conditions are specified.

Args

velocity: Vector field sampled on a grid.
obstacles: Obstacle or Geometry or tuple/list thereof to specify boundary conditions inside the domain.
solve: Solve object specifying method and tolerances for the implicit pressure solve.
active: (Optional) Mask for which cells the pressure should be solved. If given, the velocity may take NaN values where it does not contribute to the pressure. Also, the total divergence will never be subtracted if active is given, even if all values are 1.
order: spatial order for derivative computations. For Higher-order schemes, the laplace operation is not conducted with a stencil exactly corresponding to the one used in divergence calculations but a smaller one instead. While this disrupts the formal correctness of the method it only induces insignificant errors and yields considerable performance gains. supported: explicit 2/4th order - implicit 6th order (obstacles are only supported with explicit 2nd order)

Returns

velocity: divergence-free velocity of type type(velocity)
pressure: solved pressure field, CenteredGrid

Classes

class Obstacle (geometry, velocity=0, angular_velocity=0)

An obstacle defines boundary conditions inside a geometry. It can also have a linear and angular velocity.

Args

geometry: Physical shape and size of the obstacle.
velocity: Linear velocity vector of the obstacle.
angular_velocity: Rotation speed of the obstacle. Scalar value in 2D, vector in 3D.

Expand source code

class Obstacle:
    """
    An obstacle defines boundary conditions inside a geometry.
    It can also have a linear and angular velocity.
    """

    def __init__(self, geometry, velocity=0, angular_velocity=0):
        """
        Args:
            geometry: Physical shape and size of the obstacle.
            velocity: Linear velocity vector of the obstacle.
            angular_velocity: Rotation speed of the obstacle. Scalar value in 2D, vector in 3D.
        """
        self.geometry = geometry
        self.velocity = wrap(velocity, channel(geometry)) if isinstance(velocity, (tuple, list)) else velocity
        self.angular_velocity = angular_velocity
        self.shape = shape(geometry) & non_channel(self.velocity) & non_channel(angular_velocity)

    @property
    def is_stationary(self):
        """ Test whether the obstacle is completely still, i.e. not moving or rotating. """
        return not self.is_moving and not self.is_rotating

    @property
    def is_rotating(self):
        """
        Checks whether this obstacle might be rotating.
        This also evaluates to `True` if the angular velocity is unknown at this time.
        """
        return not math.always_close(self.angular_velocity, 0)

    @property
    def is_moving(self):
        """
        Checks whether this obstacle might be moving.
        This also evaluates to `True` if the velocity is unknown at this time.
        """
        return not math.always_close(self.velocity, 0)

    def copied_with(self, **kwargs):
        warnings.warn("Obstacle.copied_with is deprecated. Use math.copy_with instead.", DeprecationWarning, stacklevel=2)
        return math.copy_with(self, **kwargs)

    def __variable_attrs__(self) -> Tuple[str, ...]:
        return 'geometry', 'velocity', 'angular_velocity'

    def with_geometry(self, geometry):
        return Obstacle(geometry, self.velocity, self.angular_velocity)

    def shifted(self, delta: Tensor):
        return self.with_geometry(self.geometry.shifted(delta))

    def at(self, position: Tensor):
        return self.with_geometry(self.geometry.at(position))

    def rotated(self, angle: Union[float, Tensor]):
        return self.with_geometry(self.geometry.rotated(angle))

    def __eq__(self, other):
        if not isinstance(other, Obstacle):
            return False
        return self.geometry == other.geometry and self.velocity == other.velocity and self.angular_velocity == other.angular_velocity

Instance variables

prop is_moving

Checks whether this obstacle might be moving. This also evaluates to True if the velocity is unknown at this time.

Expand source code

@property
def is_moving(self):
    """
    Checks whether this obstacle might be moving.
    This also evaluates to `True` if the velocity is unknown at this time.
    """
    return not math.always_close(self.velocity, 0)

prop is_rotating

Checks whether this obstacle might be rotating. This also evaluates to True if the angular velocity is unknown at this time.

Expand source code

@property
def is_rotating(self):
    """
    Checks whether this obstacle might be rotating.
    This also evaluates to `True` if the angular velocity is unknown at this time.
    """
    return not math.always_close(self.angular_velocity, 0)

prop is_stationary

Test whether the obstacle is completely still, i.e. not moving or rotating.

Expand source code

@property
def is_stationary(self):
    """ Test whether the obstacle is completely still, i.e. not moving or rotating. """
    return not self.is_moving and not self.is_rotating

Methods

def at(self, position: phiml.math._tensors.Tensor)
def copied_with(self, **kwargs)
def rotated(self, angle: Union[phiml.math._tensors.Tensor, float])
def shifted(self, delta: phiml.math._tensors.Tensor)
def with_geometry(self, geometry)