.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_ch3\ch3_4_6.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_ch3_ch3_4_6.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_ch3_ch3_4_6.py:


3.4.6 Leray-orthogonal operator
===============================

The Leray-orthogonal operator plays a fundamental role in fluid dynamics,
particularly in the mathematical formulation and numerical modeling of
incompressible flows governed by the Euler or Navier–Stokes equations.

Overview
--------

The Leray projection is used to decompose a vector field into divergence-free
(incompressible) and curl-free components. In CodPy, this is achieved using
kernel-based differential operators.

We define the **Leray-orthogonal operator** as:

$$
L_k(X, Y)^{\perp} =
\nabla_k(X, Y) \cdot
\Delta_k(X, Y)^{-1} \cdot
\nabla_k^T(X, Y)
$$

or alternatively:

$$
L_k(\cdot)^{\perp} =
\nabla_k(\cdot) \cdot
\nabla_k^{-1}(\cdot)
$$

This operator projects any vector field onto the space orthogonal to the
gradient field—i.e., the divergence-free subspace.

Action on a Vector Field
------------------------

Given a vector field $f(Z) \in \mathbb{R}^{D \times N_z \times D_f}$,
the Leray-orthogonal operator acts as:

$$
L_k(Z)^{\perp} f(Z)
\in \mathbb{R}^{D \times N_z \times D_f}
$$

That is, it returns the projection of $f(Z)$ onto the orthogonal complement
of the gradient image in RKHS, which corresponds to the divergence-free component.

Use Case in Fluid Dynamics
--------------------------

This operator enables **Helmholtz–Hodge decomposition** of a vector field
into divergence-free and gradient parts:

$$
f(Z) = f^{\perp}(Z) + \nabla_k h(Z)
$$

where $f^{\perp} = L_k^{\perp} f$ is divergence-free and $\nabla_k h$
is the gradient component.

CodPy Implementation
--------------------

To compute the Leray projection operator in CodPy:

.. GENERATED FROM PYTHON SOURCE LINES 68-206

.. code-block:: Python


    # Importing necessary modules
    import os
    import sys

    from matplotlib import pyplot as plt

    curr_f = os.path.join(os.getcwd(), "codpy-book", "utils")
    sys.path.insert(0, curr_f)


    import numpy as np

    # from codpy.plotting import plot1D
    # Lets import multi_plot function from codpy utils
    from codpy.plot_utils import multi_plot


    # Define the sinusoidal function
    def periodic_fun(x):
        """
        A sinusoidal function that generates a sum of sines based on the input ``x``.
        """
        from math import pi

        sinss = np.cos(2 * x * pi)
        if x.ndim == 1:
            sinss = np.prod(sinss, axis=0)
            ress = np.sum(x, axis=0)
        else:
            sinss = np.prod(sinss, axis=1)
            ress = np.sum(x, axis=1)
        return ress + sinss


    def nabla_my_fun(x):
        from math import pi

        import numpy as np

        sinss = np.cos(2 * x * pi)
        if x.ndim == 1:
            sinss = np.prod(sinss, axis=0)
            D = len(x)
            out = np.ones((D))

            def helper(d):
                out[d] += 2.0 * sinss * pi * np.sin(2 * x[d] * pi) / np.cos(2 * x[d] * pi)

            [helper(d) for d in range(0, D)]
        else:
            sinss = np.prod(sinss, axis=1)
            N = x.shape[0]
            D = x.shape[1]
            out = np.ones((N, D))

            def helper(d):
                out[:, d] += (
                    2.0 * sinss * pi * np.sin(2 * x[:, d] * pi) / np.cos(2 * x[:, d] * pi)
                )

            [helper(d) for d in range(0, D)]
        return out


    # Function to generate periodic data
    def generate_periodic_data_cartesian(size_x, size_z, fun=None, nabla_fun=None):
        """
        Generates 2D structured Cartesian grid data for x and z domains,
        and evaluates a given function and optionally its gradient.

        Parameters:
        - size_x: number of points per axis for x (grid will be size_x^2)
        - size_z: number of points per axis for z (grid will be size_z^2)
        - fun: function to evaluate at each point
        - nabla_fun: optional gradient function to evaluate

        Returns:
        - x, z: 2D Cartesian grids of shape (N, 2)
        - fx, fz: function values at x and z
        - nabla_fx, nabla_fz (if nabla_fun is provided)
        """

        def cartesian_grid(size, box):
            lin = [np.linspace(box[0, d], box[1, d], size) for d in range(2)]
            X, Y = np.meshgrid(*lin)
            return np.stack([X.ravel(), Y.ravel()], axis=1)

        # Define domain boxes
        X_box = np.array([[-1, -1], [1, 1]])
        Z_box = np.array([[-1.5, -1.5], [1.5, 1.5]])

        # Generate Cartesian grids
        x = cartesian_grid(size_x, X_box)
        z = cartesian_grid(size_z, Z_box)

        # Function evaluations
        fx = fun(x).reshape(-1, 1) if fun else None
        fz = fun(z).reshape(-1, 1) if fun else None

        if nabla_fun:
            nabla_fx = nabla_fun(x)
            nabla_fz = nabla_fun(z)
            return x, fx, z, fz, nabla_fx, nabla_fz

        return x, fx, z, fz


    # Lets define helper function to plot 3D projection of the function
    def plot_trisurf(xfx, ax, legend="", elev=90, azim=-100, **kwargs):
        from matplotlib import cm

        """
        Helper function to plot a 3D surface using a trisurf plot.

        Parameters:
        - xfx: A tuple containing the x-coordinates (2D points) and their 
          corresponding function values.
        - ax: The matplotlib axis object for plotting.
        - legend: The legend/title for the plot.
        - elev, azim: Elevation and azimuth angles for the 3D view.
        - kwargs: Additional keyword arguments for further customization.
        """

        xp, fxp = xfx[0], xfx[1]
        x, fx = xp, fxp

        X, Y = x[:, 0], x[:, 1]
        Z = fx.flatten()
        ax.plot_trisurf(X, Y, Z, antialiased=False, cmap=cm.jet)
        ax.view_init(azim=azim, elev=elev)
        ax.title.set_text(legend)


    # import CodPy's core module and Kernel class
    from codpy import core
    from codpy.kernel import Kernel








.. GENERATED FROM PYTHON SOURCE LINES 207-210

Leray operator and Helmholtz-Hodge decomposition

########################################################################

.. GENERATED FROM PYTHON SOURCE LINES 210-299

.. code-block:: Python



    def fun_LerayT(size_x=50, size_y=50):
        x, fx, z, fz, _, nabla_fz = generate_periodic_data_cartesian(
            size_x, size_y, periodic_fun, nabla_fun=nabla_my_fun
        )

        nabla_fz = nabla_fz.reshape(-1, 2, 1)

        kernel_ptr = Kernel(
            x=x, fx=fx, 
        ).get_kernel()
        LerayT_fz = core.DiffOps.nabla(
            x=x,
            y=x,
            z=z,
            fx=core.DiffOps.nabla_inv(x=x, y=x, z=x, fz=nabla_fz, kernel_ptr=kernel_ptr),
            kernel_ptr=kernel_ptr,
            order=2,
            regularization=1e-8,
        )
        multi_plot(
            [
                (z, nabla_fz[:, 0, :]),
                (z, LerayT_fz[:, 0, :]),
                (z, nabla_fz[:, 1, :]),
                (z, LerayT_fz[:, 1, :]),
            ],
            plot_trisurf,
            projection="3d",
            mp_nrows=1,
            mp_figsize=(12, 3),
            mp_title=[
                "Comparing f(z) and the transpose of the Leray operator on each direction"
            ],
        )
        plt.show()


    fun_LerayT()


    def fun_LerayT_nabla_nablainv(size_x=50, size_y=50):
        x, fx, z, fz, _, nabla_fz = generate_periodic_data_cartesian(
            size_x, size_y, periodic_fun, nabla_fun=nabla_my_fun
        )

        kernel_ptr = Kernel(
            x=x, fx=fx, set_kernel=core.kernel_setter("tensornorm", "scale_to_unitcube")
        ).get_kernel()
        nabla_fx = core.DiffOps.nabla(
            x=x,
            y=x,
            z=z,
            fx=fx,
            kernel_ptr=kernel_ptr,
            order=2,
            regularization=1e-8,
        )
        nabla_fz = core.DiffOps.leray_t(x=x, y=x, fx=nabla_fx, kernel_ptr=kernel_ptr)

        LerayT_fz = core.DiffOps.nabla(
            x=x,
            y=x,
            z=z,
            fx=core.DiffOps.nabla_inv(x=x, y=x, z=z, fz=nabla_fz, kernel_ptr=kernel_ptr),
            kernel_ptr=kernel_ptr,
            order=2,
            regularization=1e-8,
        )
        multi_plot(
            [
                (z, nabla_fz[:, 0, :]),
                (z, LerayT_fz[:, 0, :]),
                (z, nabla_fz[:, 1, :]),
                (z, LerayT_fz[:, 1, :]),
            ],
            plot_trisurf,
            projection="3d",
            mp_nrows=1,
            mp_figsize=(12, 3),
            mp_title=[
                "Comparing $\nabla \nabla^{-1}f(z)$ and the transpose of the Leray operator $L_k \nabla f(z)$ on each direction"
            ],
        )
        plt.show()


    fun_LerayT_nabla_nablainv()



.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /auto_ch3/images/sphx_glr_ch3_4_6_001.png
         :alt: ['Comparing f(z) and the transpose of the Leray operator on each direction']
         :srcset: /auto_ch3/images/sphx_glr_ch3_4_6_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_ch3/images/sphx_glr_ch3_4_6_002.png
         :alt: ['Comparing $\nabla \nabla^{-1}f(z)$ and the transpose of the Leray operator $L_k \nabla f(z)$ on each direction']
         :srcset: /auto_ch3/images/sphx_glr_ch3_4_6_002.png
         :class: sphx-glr-multi-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 5.676 seconds)


.. _sphx_glr_download_auto_ch3_ch3_4_6.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: ch3_4_6.ipynb <ch3_4_6.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: ch3_4_6.py <ch3_4_6.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: ch3_4_6.zip <ch3_4_6.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_