<commit_before>"""
Advection-Diffusion Equation
============================

Here we use the discretize package to model the advection-diffusion
equation. The goal of this tutorial is to demonstrate:

    - How to solve time-dependent PDEs
    - How to apply Neumann boundary conditions


Derivation
----------

If we assume the fluid is incompressible and that diffusivity is spatially,
the advection-diffusion equation with Neumann boundary conditions is given by:

.. math::
    p_t = \\alpha \\nabla \\cdot \\nabla p
    - \\mathbf{u} \\cdot \\nabla p + s \n
    \\textrm{s.t.} \;\;\; \\frac{\\partial p}{\\partial n} \Bigg |_{\partial \Omega} = 0

where :math:`p` is an unknown variable, :math:`\\alpha` is the diffusivity
constant, :math:`\\mathbf{u}` the velocity field and :math:`s` is the source
term. We will consider the case where there are two point sources within our
domain. Thus:

.. math::
    s = s_0 \\big [ \\delta ( \\mathbf{r_1}) + \\delta (\\mathbf{r_2} ) \\big ]

To solve this problem numerically, we use the weak formulation; that is, we
take the inner product of the equation with an appropriate test function
(:math:`\\psi`):

.. math::
    \int_\\Omega \\psi \\, p_t \\, dv =
    \\alpha \int_\\Omega \\psi (\\nabla \\cdot \\nabla p) \\, dv
    - \int_\\Omega \\psi \\, \mathbf{u} \\cdot \\nabla p \\, dv
    + \int_\\Omega \\psi \\, s \\, dv

The source term is a volume integral containing the Dirac delta function, thus:

.. math::
    \int_\\Omega \\psi \\, s \\, dv = \\psi \\, q

where :math:`q=s_0` at the location of either point source and zero everywhere
else. In the weak formulation, we generally remove divergence terms. If we use
the identity
:math:`\\phi \\nabla \\cdot \\mathbf{a} = \\nabla \\cdot (\\phi \\mathbf{a}) - \\mathbf{a} \\cdot (\\nabla \\phi )`
on the diffusion term and apply the divergence theorem we obtain:

.. math::
    \int_\\Omega \\psi \\, p_t \\, dv =
    \\alpha \int_{\\partial \\Omega} \\mathbf{n} \\cdot ( \\psi \\nabla p ) \\, da
    - \\alpha \int_\\Omega (\\nabla \\psi ) \\cdot (\\nabla p ) \\, dv
    - \int_\\Omega (\\psi \\, \mathbf{u} ) \\cdot \\nabla p dv
    + \\psi \\, q

Since the flux on the faces is zero at the boundaries, we can eliminate
the first term on the right-hand side. By evaluating the inner products
according to the finite volume approach we obtain:

.. math::
    \\mathbf{\\psi^T M_c p_t} = - \\alpha \\mathbf{\\psi^T G^T M_f G \\, p}
    - \\mathbf{\\psi^T M_c A_{fc}} \\textrm{diag} ( \\mathbf{u} ) \\mathbf{G \\, p}
    + \\mathbf{\\psi^T M_c s}

where :math:`\\mathbf{\\psi}`, :math:`\\mathbf{p}` and :math:`\\mathbf{p_t}`
live at cell centers and :math:`\\mathbf{u}` lives on faces. :math:`\\mathbf{G}`
is the discrete gradient operators. And :math:`\\mathbf{A_{cf}}`
averages the component of :math:`\\mathbf{u} \\cdot \\nabla p` to cell centers
and sums them. :math:`\\mathbf{M_c}` and :math:`\\mathbf{M_f}` are the cell
center and face inner product matricies, respectively.

By eliminating :math:`\\psi^T` and multiplying both sides by
:math:`\\mathbf{M_c^{-1}}` we obtain:

.. math::
    \\mathbf{p_t} = - \\mathbf{M} \\mathbf{p} + \\mathbf{s}

where

.. math::
    \\mathbf{M} = \\alpha \\mathbf{M_c^{-1} G^T M_f G}
    + \\mathbf{A_{fc}} \\textrm{diag}(\\mathbf{u}) \\mathbf{G}

For the example, we will discretize in time using backward Euler. This results
in the following system which must be solve at every time step :math:`k`.
Where :math:`\\Delta t` is the step size:

.. math::
    \\big [ \\mathbf{I} + \\Delta t \\, \\mathbf{M} \\big ] \\mathbf{p}^{k+1} =
    \\mathbf{p}^k + \\Delta t \\, \\mathbf{s}



"""

###################################################
#
# Import Packages
# ---------------
#
# Here we import the packages required for this tutorial.
#

from discretize import TensorMesh
from pymatsolver import SolverLU
import matplotlib.pyplot as plt
import numpy as np
from discretize.utils.matutils import sdiag


# Create a tensor mesh
h = np.ones(75)
mesh = TensorMesh([h, h], 'CC')

# Define diffusivity constant
a = 25

# Define velocity vector u = [ux, uy] = [1, -1] on faces
ux = 10*np.ones((mesh.nFx))
uy = - 10*np.ones((mesh.nFy))
u = np.r_[ux, uy]

# Define source term diag(v)*s<commit_msg>Add subheader<commit_after>"""
Advection-Diffusion Equation
============================

Here we use the discretize package to model the advection-diffusion
equation. The goal of this tutorial is to demonstrate:

    - How to solve time-dependent PDEs
    - How to apply Neumann boundary conditions


Derivation
----------

If we assume the fluid is incompressible and that diffusivity is spatially,
the advection-diffusion equation with Neumann boundary conditions is given by:

.. math::
    p_t = \\alpha \\nabla \\cdot \\nabla p
    - \\mathbf{u} \\cdot \\nabla p + s \n
    \\textrm{s.t.} \;\;\; \\frac{\\partial p}{\\partial n} \Bigg |_{\partial \Omega} = 0

where :math:`p` is an unknown variable, :math:`\\alpha` is the diffusivity
constant, :math:`\\mathbf{u}` the velocity field and :math:`s` is the source
term. We will consider the case where there are two point sources within our
domain. Thus:

.. math::
    s = s_0 \\big [ \\delta ( \\mathbf{r_1}) + \\delta (\\mathbf{r_2} ) \\big ]

To solve this problem numerically, we use the weak formulation; that is, we
take the inner product of the equation with an appropriate test function
(:math:`\\psi`):

.. math::
    \int_\\Omega \\psi \\, p_t \\, dv =
    \\alpha \int_\\Omega \\psi (\\nabla \\cdot \\nabla p) \\, dv
    - \int_\\Omega \\psi \\, \mathbf{u} \\cdot \\nabla p \\, dv
    + \int_\\Omega \\psi \\, s \\, dv

The source term is a volume integral containing the Dirac delta function, thus:

.. math::
    \int_\\Omega \\psi \\, s \\, dv = \\psi \\, q

where :math:`q=s_0` at the location of either point source and zero everywhere
else. In the weak formulation, we generally remove divergence terms. If we use
the identity
:math:`\\phi \\nabla \\cdot \\mathbf{a} = \\nabla \\cdot (\\phi \\mathbf{a}) - \\mathbf{a} \\cdot (\\nabla \\phi )`
on the diffusion term and apply the divergence theorem we obtain:

.. math::
    \int_\\Omega \\psi \\, p_t \\, dv =
    \\alpha \int_{\\partial \\Omega} \\mathbf{n} \\cdot ( \\psi \\nabla p ) \\, da
    - \\alpha \int_\\Omega (\\nabla \\psi ) \\cdot (\\nabla p ) \\, dv
    - \int_\\Omega (\\psi \\, \mathbf{u} ) \\cdot \\nabla p dv
    + \\psi \\, q

Since the flux on the faces is zero at the boundaries, we can eliminate
the first term on the right-hand side. By evaluating the inner products
according to the finite volume approach we obtain:

.. math::
    \\mathbf{\\psi^T M_c p_t} = - \\alpha \\mathbf{\\psi^T G^T M_f G \\, p}
    - \\mathbf{\\psi^T M_c A_{fc}} \\textrm{diag} ( \\mathbf{u} ) \\mathbf{G \\, p}
    + \\mathbf{\\psi^T M_c s}

where :math:`\\mathbf{\\psi}`, :math:`\\mathbf{p}` and :math:`\\mathbf{p_t}`
live at cell centers and :math:`\\mathbf{u}` lives on faces. :math:`\\mathbf{G}`
is the discrete gradient operators. And :math:`\\mathbf{A_{cf}}`
averages the component of :math:`\\mathbf{u} \\cdot \\nabla p` to cell centers
and sums them. :math:`\\mathbf{M_c}` and :math:`\\mathbf{M_f}` are the cell
center and face inner product matricies, respectively.

By eliminating :math:`\\psi^T` and multiplying both sides by
:math:`\\mathbf{M_c^{-1}}` we obtain:

.. math::
    \\mathbf{p_t} = - \\mathbf{M} \\mathbf{p} + \\mathbf{s}

where

.. math::
    \\mathbf{M} = \\alpha \\mathbf{M_c^{-1} G^T M_f G}
    + \\mathbf{A_{fc}} \\textrm{diag}(\\mathbf{u}) \\mathbf{G}

For the example, we will discretize in time using backward Euler. This results
in the following system which must be solve at every time step :math:`k`.
Where :math:`\\Delta t` is the step size:

.. math::
    \\big [ \\mathbf{I} + \\Delta t \\, \\mathbf{M} \\big ] \\mathbf{p}^{k+1} =
    \\mathbf{p}^k + \\Delta t \\, \\mathbf{s}
    


"""

###################################################
#
# Import Packages
# ---------------
#
# Here we import the packages required for this tutorial.
#

from discretize import TensorMesh
from pymatsolver import SolverLU
import matplotlib.pyplot as plt
import numpy as np
from discretize.utils.matutils import sdiag

###############################################
#
# Solving the Problem
# -------------------
#

# Create a tensor mesh
h = np.ones(75)
mesh = TensorMesh([h, h], 'CC')

# Define diffusivity constant
a = 25

# Define velocity vector u = [ux, uy] = [1, -1] on faces
ux = 10*np.ones((mesh.nFx))
uy = - 10*np.ones((mesh.nFy))
u = np.r_[ux, uy]

# Define source term diag(v)*s
{'fraction_non_alphanumeric': 0.20037322136692326, 'fraction_numerical': 0.004198740377886634, 'mean_word_length': 2.7808641975308643, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6629679', 'n_tokens_mistral': 2855, 'n_tokens_neox': 2789, 'n_words': 1033}