Function differential

phi/physics/advect.py:78–130

Computes the differential advection term using the differentiation Scheme indicated by `order`, ´implicit´ and `upwind`. For a velocity field u, the advection term as it appears on the right-hand-side of a PDE is -u·∇u, including the negative sign. For unstructured meshes, computes -1

(u: Field,
                 velocity: Field,
                 density: float = 1.,
                 order=2,
                 implicit: Solve = None,
                 upwind=True)

def differential(u: Field,
                 velocity: Field,
                 density: float = 1.,
                 order=2,
                 implicit: Solve = None,
                 upwind=True) -> Field:
    """
    Computes the differential advection term using the differentiation Scheme indicated by `order`, ´implicit´ and `upwind`.

    For a velocity field u, the advection term as it appears on the right-hand-side of a PDE is -u·∇u, including the negative sign.

    For unstructured meshes, computes -1/V ∑_f (n·u_prev) u ρ A

    Args:
        u: Scalar or vector-valued `Field` sampled on a `CenteredGrid`, `StaggeredGrid` or `Mesh`.
        velocity: `Field` that can be sampled at the elements of `u`.
            For FVM, the advection term is typically linearized by setting `velocity = previous_velocity`.
            Passing `velocity=u` yields non-linear terms which cannot be traced inside linear functions.
        order: Spatial order of accuracy.
            Higher orders entail larger stencils and more computation time but result in more accurate results assuming a large enough resolution.
            Supported for grids: 2 explicit, 4 explicit, 6 implicit (inherited from `phi.field.spatial_gradient()` and resampling).
            Passing order=4 currently uses 2nd-order resampling. This is work-in-progress.
            For FVM, the order is used when interpolating centroid values to faces if needed.
        implicit: When a `Solve` object is passed, performs an implicit operation with the specified solver and tolerances.
            Otherwise, an explicit stencil is used.
        upwind: Whether to use upwind interpolation. Only supported for FVM at the moment.

    Returns:
        Differential convection term as `Field` on the same geometry.
    """
    if u.is_grid and u.is_staggered:
        grad_list = [spatial_gradient(field_component, stack_dim=channel('grad_dim'), order=order, implicit=implicit) for field_component in u.vector]
        grad_grid = u.with_values(math.stack([component.values for component in grad_list], channel(velocity).as_dual()))
        if order == 4:
            amounts = [grad * vel.at(grad, order=2) for grad, vel in zip(grad_grid.grad_dim, velocity.vector)]  # ToDo resampling does not yet support order=4
        else:
            amounts = [grad * vel.at(grad, order=order, implicit=implicit) for grad, vel in zip(grad_grid.grad_dim, velocity.vector)]
        amount = sum(amounts)
        return u.with_values(- amount)
    elif u.is_grid and u.is_centered:
        grad_tensor = math.stack(
            [spatial_gradient(component, stack_dim=channel('gradient'), order=order, implicit=implicit).values for
             component in u.vector], dim=channel('vector'))
        velocity_tensor = math.stack(math.unstack(velocity.values, dim='vector'), dim=channel('gradient'))
        amounts = velocity_tensor * grad_tensor
        amount = sum(amounts.gradient)
        return velocity.with_values(- amount)
    elif u.is_mesh:
        u = u.at_faces(boundary=NONE, order=order, upwind=velocity if upwind is True else upwind)
        velocity = velocity.at_faces(boundary=NONE, order=order, upwind=velocity if upwind is True else upwind)
        conv = density * u.mesh.integrate_surface(u.values * (velocity.values.vector @ velocity.face_normals.vector)) / u.mesh.volume
        return Field(u.geometry, -conv, 0)
    raise NotImplementedError(u)


finite_difference = differential