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

.. only:: html

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

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

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

.. _sphx_glr_auto_examples_plot_advdif.py:


.. currentmodule:: fdg

N-D Advection Diffusion
=======================

This example demonstrates system for N-dimensional advection-diffusion equation can be set
up.

The equation is defined in the weak form as:

.. math::
    :label: examples_nd_advdiff_1

    \left( p^{(n - 1)}, q^{(n - 1)} \right)_\Omega + \left( \mathrm{d} p^{(n - 1)},
    u^{(n)} \right)_\Omega = \int_{\partial \Omega} p^{(n - 1)} \wedge \star u^{(n)}

.. math::
    :label: examples_nd_advdiff_2

    \left( v^{(n)}, \mathrm{d} q^{(n - 1)} \right)_\Omega + \left( i_{\vec{a}} v^{(n)},
    q^{(n - 1)} \right)_\Omega = \left( v^{(n)}, f^{(n)} \right)_\Omega

.. GENERATED FROM PYTHON SOURCE LINES 25-53

.. code-block:: Python


    from functools import partial
    from time import perf_counter
    from typing import Sequence

    import matplotlib.pyplot as plt
    import numpy as np
    import numpy.typing as npt
    from fdg import (
        BasisSpecs,
        BasisType,
        CoordinateMap,
        DegreesOfFreedom,
        FunctionSpace,
        IntegrationMethod,
        IntegrationSpace,
        IntegrationSpecs,
        KForm,
        KFormSpecs,
        SpaceMap,
        compute_kform_interior_product_matrix,
        compute_kform_mass_matrix,
        incidence_kform_operator,
        projection_kform_l2_dual,
        reconstruct,
        transform_kform_to_target,
    )








.. GENERATED FROM PYTHON SOURCE LINES 54-77

The manufactured solution for the general N-dimensional case uses
the following for the manufactured solution:

.. math::
    :label: examples_nd_advdiff_man_sol

    u^{(n)}(x_1, \dots, x_n) = k \left(\prod\limits_{i=1}^n \cos\left( \frac{\pi}{2} x_i
    \right)\right) \mathrm{d} x_1 \wedge \dots \wedge \mathrm{d} x_n


This gives the forcing function:

.. math::
    :label: examples_nd_advdiff_man_for

    f^{(n)}(x_1, \dots, x_n) = k \left( - n \left(\frac{\pi}{2}\right)^2 \left(
    \prod\limits_{i=1}^n \cos\left( \frac{\pi}{2} x_i
    \right)\right)  +
    \sum\limits_{i = 1}^n a_i \frac{\pi}{2} \sin \frac{\pi x_i}{2} \prod\limits_{j=1,
    j \ne i}^n \cos \frac{\pi x_j}{2} \right)
    \mathrm{d} x_1 \wedge \dots \wedge \mathrm{d} x_n



.. GENERATED FROM PYTHON SOURCE LINES 78-135

.. code-block:: Python


    SCALE = 0.1
    ALPHA = 0.05


    def manufactured_solution(*x: npt.NDArray[np.double]) -> npt.NDArray[np.double]:
        """Exact manufactured solution."""
        res = np.cos(x[0] * np.pi / 2)
        for v in x[1:]:
            res *= np.cos(v * np.pi / 2)
        return res * SCALE


    def vector_field(
        coeffs: npt.NDArray[np.double], *x: npt.NDArray[np.double]
    ) -> tuple[npt.NDArray[np.double], ...]:
        """Vector field used for advecion."""
        vals: list[npt.NDArray[np.double]] = list()
        for c_row in coeffs:
            vals.append(np.sum([c * v for c, v in zip(c_row, x, strict=True)], axis=0))
        return tuple(vals)


    def manufactured_source_laplacian(*x: npt.NDArray[np.double]) -> npt.NDArray[np.double]:
        """Exact Laplacian source term."""
        res = np.cos(x[0] * np.pi / 2)
        for v in x[1:]:
            res *= np.cos(v * np.pi / 2)
        res *= -((np.pi / 2) ** 2) * len(x)
        return res * SCALE


    def manufactured_source_advection(
        coeffs: npt.NDArray[np.double], *x: npt.NDArray[np.double]
    ) -> npt.NDArray[np.double]:
        """Exact advection source term."""
        vf = vector_field(coeffs, *x)
        res = np.zeros_like(x[0])
        for idx, v in enumerate(vf):
            grad_i = -np.pi / 2 * np.sin(np.pi / 2 * x[idx])
            for j, s in enumerate(x):
                if j == idx:
                    continue
                grad_i *= np.cos(np.pi / 2 * s)
            res += grad_i * v
        return res * SCALE


    def manufactured_source(
        coeffs: npt.NDArray[np.double], *x: npt.NDArray[np.double]
    ) -> npt.NDArray[np.double]:
        """Exact source term."""
        return ALPHA * manufactured_source_laplacian(*x) + manufactured_source_advection(
            coeffs, *x
        )









.. GENERATED FROM PYTHON SOURCE LINES 136-148

First the geometry of the space this will be solved on will be defined. For this case,
we use the unit square, where the interior is deformed, while boundaries are the same.
The mapping for each coordinate is based on Equation
:eq:`examples_nd_advdiff_deformation`, with :math:`c` being the parameter that
determines the scale of deformation.


.. math::
    :label: examples_nd_advdiff_deformation

    x_i = \xi_i + c \prod\limits_{j=1}^n \left( 1 - {x_j}^2 \right) \sin \pi x_j


.. GENERATED FROM PYTHON SOURCE LINES 149-179

.. code-block:: Python



    def disturbed_mapping(
        c: float, idx: int, *x: npt.NDArray[np.double]
    ) -> npt.NDArray[np.double]:
        """Return a perturbed map, where the boundaries are not affected, but the interior is.

        Parameters
        ----------
        c : float
            Strenght of the disturbance.

        idx : int
            Index of the input to base the mapping on.

        *x : array
            Coordinates where the mapping should be computed.

        Returns
        -------
        array
            Input coordinate ``idx``, but somewhat.
        """
        base = x[idx]
        d = np.full_like(base, c)
        for v in x:
            d *= (1 - v**2) * np.sin(np.pi * v)
        return base + d









.. GENERATED FROM PYTHON SOURCE LINES 180-183

Mappings for every coordinate are collected into a joined :class:`SpaceMap`,
which is then used to map :math:`k`-form components between the reference domain
and the physical domain.

.. GENERATED FROM PYTHON SOURCE LINES 184-208

.. code-block:: Python



    def create_space_map(
        c: float, orders: Sequence[int], space: IntegrationSpace
    ) -> SpaceMap:
        """Create space map that are is disturbed."""
        func_space = FunctionSpace(
            *(BasisSpecs(BasisType.LAGRANGE_UNIFORM, order) for order in orders)
        )
        ndim = len(orders)
        points = np.meshgrid(
            *[np.linspace(-1, +1, order + 1) for order in orders], indexing="ij"
        )
        return SpaceMap(
            *[
                CoordinateMap(
                    DegreesOfFreedom(func_space, disturbed_mapping(c, idim, *points)),
                    space,
                )
                for idim in range(ndim)
            ]
        )









.. GENERATED FROM PYTHON SOURCE LINES 209-218

With these two utilities, we can start working on computing convergence
of the FEM discretization in the :math:`L^2`-norm

First we must define how integration will be done. Here two
:class:`IntegrationSpace` objects are defined - one for computing our
results and another, finer, for computing the error.

With these discretizations defined, we can use previously written functions to
create the space mappings.

.. GENERATED FROM PYTHON SOURCE LINES 219-234

.. code-block:: Python



    def create_space_maps(order_integration, type_integration, ndim, dp):
        """Create integration spaces."""
        int_space = IntegrationSpace(
            *((IntegrationSpecs(order_integration, type_integration),) * ndim)
        )
        int_space_higher = IntegrationSpace(
            *((IntegrationSpecs(order_integration + dp, type_integration),) * ndim)
        )
        space_map = create_space_map(0.1, [5] * ndim, int_space)
        space_map_high = create_space_map(0.1, [5] * ndim, int_space_higher)
        return space_map, space_map_high









.. GENERATED FROM PYTHON SOURCE LINES 235-242

Along with the :class:`IntegrationSpace` objects to define integration we
must define the discretization of the :math:`k`-forms using a :class:`FunctionSpace`
object.

With base function space defined, we can define some :math:`k`-form
specifications using :class:`KFormSpecs`. These do not contain any degrees of
freedom themselves, but provide information about the order and function spaces.

.. GENERATED FROM PYTHON SOURCE LINES 243-254

.. code-block:: Python



    def create_kform_specs(type_basis, order_basis, ndim):
        """Create k-form specifications."""
        base_space = FunctionSpace(*((BasisSpecs(type_basis, order_basis),) * ndim))

        specs_u = KFormSpecs(ndim, base_space)
        specs_q = KFormSpecs(ndim - 1, base_space)
        return specs_u, specs_q









.. GENERATED FROM PYTHON SOURCE LINES 255-259

Since we aim to solve a (linear) advection-diffusion problem, we need
the vector field. In this example it is a linear function of the coordinates.
While we could use completely coefficients, by making them depend only on
the number of dimensions we can fairly compare subsequent solves.

.. GENERATED FROM PYTHON SOURCE LINES 260-272

.. code-block:: Python



    def creat_vector_filed_coefficients(ndim):
        """Create vector field coefficients."""
        return np.array(
            [
                [((-1) ** i * (2 * i - j**2) / (1 + i**2 + j**2)) for i in range(ndim)]
                for j in range(ndim)
            ]
        )









.. GENERATED FROM PYTHON SOURCE LINES 273-275

To be able to use the vector field, we must evaluate the components in the
physical domain at integration points.

.. GENERATED FROM PYTHON SOURCE LINES 275-285

.. code-block:: Python



    def evaluate_vector_field(coeffs, sm: SpaceMap):
        """Evaluate vector field at integration points."""
        vf = vector_field(
            coeffs, *[sm.coordinate_map(idx).values for idx in range(sm.input_dimensions)]
        )
        return vf









.. GENERATED FROM PYTHON SOURCE LINES 286-294

While for left side symmetry could be exploited, there's no harm to compute
it in full for this example. To that end, we first compute the two mass
matrices for the two :math:`k`-forms, then apply the incidence operators
as needed.

Now that we have the required blocks, we can assemble the system matrix.
Alternatively we could have used Schur's complement to compute the solution,
but for the sake of simplicity here we use the full dense solve.

.. GENERATED FROM PYTHON SOURCE LINES 294-324

.. code-block:: Python



    def assemble_lhs(sm, specs_q, specs_u, vf):
        """Crate the system matrix."""
        mq = compute_kform_mass_matrix(
            sm, specs_q.order, specs_q.base_space, specs_q.base_space
        )
        mu = compute_kform_mass_matrix(
            sm, specs_u.order, specs_u.base_space, specs_u.base_space
        )
        mqip = compute_kform_interior_product_matrix(
            sm,
            specs_u.order,
            specs_u.base_space,
            specs_q.base_space,
            np.asarray(vf, np.double),
        )

        mu_e = incidence_kform_operator(specs_q, mu, right=True)
        et_mu = incidence_kform_operator(specs_q, mu, transpose=True)

        system_matrix = np.block(
            [
                [mq, et_mu],
                [ALPHA * mu_e + mqip.T, np.zeros_like(mu)],
            ]
        )
        return system_matrix









.. GENERATED FROM PYTHON SOURCE LINES 325-327

The right side of the Poisson equation can be computed from the "dual projection"
of the manufactured source term on the function space.

.. GENERATED FROM PYTHON SOURCE LINES 327-341

.. code-block:: Python



    def assemble_rhs(specs_u, specs_q, sm_h, coeffs):
        """Assemble the system's RHS."""
        source_vals = projection_kform_l2_dual(
            [partial(manufactured_source, coeffs)], specs_u, sm_h
        )[0]

        rhs = np.concatenate(
            (np.zeros(sum(specs_q.component_dof_counts)), source_vals.flatten())
        )
        return rhs









.. GENERATED FROM PYTHON SOURCE LINES 342-346

From here we can split solution vector into degrees of freedom of individual
:math:`k`-forms represented by :class:`KForm` objects.
To compute the :math:`L^2` error, we need to reconstruct the computed solution,
subtract the manufactured solution, then integrate the square of the error.

.. GENERATED FROM PYTHON SOURCE LINES 346-378

.. code-block:: Python

    def reconstruct_error_l2(specs_q, specs_u, solution_dofs, ndim, sm_h):
        """Reconstruct the solution and compute the L2 error."""
        sol_q = KForm(specs_q)
        sol_u = KForm(specs_u)

        nq = sum(specs_q.component_dof_counts)

        sol_q.values[:] = solution_dofs[:nq]
        sol_u.values[:] = solution_dofs[nq:]

        u_dofs = DegreesOfFreedom(
            specs_u.get_component_function_space(0), sol_u.get_component_dofs(0)
        )
        # K-form computed values at integration nodes
        computed_values = transform_kform_to_target(
            ndim,
            sm_h,
            [reconstruct(u_dofs, *sm_h.integration_space.nodes())],
        )[0]
        # K-form exact values at integration nodes
        real_values = manufactured_solution(
            *[sm_h.coordinate_map(idx).values for idx in range(ndim)]
        )
        # Error
        err_l2 = np.sum(
            (computed_values - real_values) ** 2
            * sm_h.determinant
            * sm_h.integration_space.weights()
        )
        return err_l2









.. GENERATED FROM PYTHON SOURCE LINES 379-381

All the small building blocks discussed before can now be put together
to form the error calculation function.

.. GENERATED FROM PYTHON SOURCE LINES 381-412

.. code-block:: Python



    def compute_l2_error(
        order_integration: int,
        type_integration: IntegrationMethod,
        order_basis: int,
        type_basis: BasisType,
        ndim: int,
        dp: int,
    ) -> float:
        """Solve the N-dimensional Poisson equation and compute the L^2 error."""
        # Space maps
        sm, sm_h = create_space_maps(order_integration, type_integration, ndim, dp)
        # Vector coefficients
        coeffs = creat_vector_filed_coefficients(ndim)
        # Vector field values
        vf = evaluate_vector_field(coeffs, sm)
        # K-form specs
        specs_u, specs_q = create_kform_specs(type_basis, order_basis, ndim)
        # LHS of the system
        lhs = assemble_lhs(sm, specs_q, specs_u, vf)
        # RHS
        rhs = assemble_rhs(specs_u, specs_q, sm_h, coeffs)
        # Solve
        solution_dofs = np.linalg.solve(lhs, rhs)
        # Compute error squared
        err_l2 = reconstruct_error_l2(specs_q, specs_u, solution_dofs, ndim, sm_h)
        # Retur the error
        return float(np.sqrt(err_l2))









.. GENERATED FROM PYTHON SOURCE LINES 413-415

For this test, we will use Bernstein basis, Gauss integration rule, and the order
difference of 1 between the lower and higher order integration rules.

.. GENERATED FROM PYTHON SOURCE LINES 415-455

.. code-block:: Python

    BTYPE = BasisType.BERNSTEIN
    ITYPE = IntegrationMethod.GAUSS
    DP = 1

    for ndim, plimit in zip(range(1, 4), [15, 11, 6], strict=True):
        pvals = np.arange(1, plimit)
        evals = np.zeros(pvals.size)
        tvals = np.zeros(pvals.size)
        for ip, p in enumerate(pvals):
            pv = int(p)
            t0 = perf_counter()
            l2 = compute_l2_error(pv + 0, ITYPE, pv, BTYPE, ndim, DP)
            t1 = perf_counter()
            evals[ip] = l2
            tvals[ip] = t1 - t0

        k1, k0 = np.polyfit(pvals, np.log(evals), deg=1)
        c = np.exp(k0)
        b = np.exp(k1)

        fig, ax = plt.subplots()

        ax.scatter(pvals, evals)
        ax.plot(
            pvals,
            c * b**pvals,
            linestyle="dashed",
            label=f"$\\varepsilon = {c:.2g} \\cdot {b:.2g}^{{p}}$",
        )
        ax.set(
            yscale="log",
            xlabel="$p$",
            ylabel="$\\left|\\left| \\varepsilon \\right|\\right|_{ L^2 }$",
        )
        ax.grid()
        ax.legend()
        ax.set_title(f"{ndim}-dimensional Poisson equation convergence")
        fig.tight_layout()

    plt.show()



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


    *

      .. image-sg:: /auto_examples/images/sphx_glr_plot_advdif_001.png
         :alt: 1-dimensional Poisson equation convergence
         :srcset: /auto_examples/images/sphx_glr_plot_advdif_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_examples/images/sphx_glr_plot_advdif_002.png
         :alt: 2-dimensional Poisson equation convergence
         :srcset: /auto_examples/images/sphx_glr_plot_advdif_002.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_examples/images/sphx_glr_plot_advdif_003.png
         :alt: 3-dimensional Poisson equation convergence
         :srcset: /auto_examples/images/sphx_glr_plot_advdif_003.png
         :class: sphx-glr-multi-img






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

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


.. _sphx_glr_download_auto_examples_plot_advdif.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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