.. sectionauthor::   Yasuyuki SHIMIZU <yasu@i-ric.org>

.. _2d_advection_eq:

######################################
The 2D advection equation
######################################

Basic formulas
================

Here I show the 2-dimensional advection equation.

.. math::
    :label: 2d_advection

    {{\partial f}\over{\partial t}} +
    u{{\partial f}\over{\partial x}}+
    v{{\partial f}\over{\partial y}} = 0

Where, :math:`t` is the time, :math:`x, y` are the axes of Cartesian coordinates, :math:`f` is a variable to be advected, and :math:`u,v` are advection velocities in the directions of :math:`x, y`. 

Similar to the 1D advection equation, while the distribution form is maintained, the theoretical solutions to this equation, :math:`f`, move in the direction of :math:`x` at a velocity of :math:`u` and in the direction of :math:`ｙ` at a velocity of :math:`v`.

.. _2dkazakami:

The upwind difference scheme
===============================

The upwind difference equation for the advection term is as follows.

.. math::
    :label: 2d_ad_1

    u>0 \mbox{When}; \; \; 
    u\cfrac{\partial f}{\partial x}=u\cfrac
    {f(i,j)-f(i-1,j)}{\Delta x}    \\
    u<0 \mbox{When}; \; \; 
    u\cfrac{\partial f}{\partial x}=u\cfrac
    {f(i+1,j)-f(i,j)}{\Delta x}    

.. math::
    :label: 2d_ad_2

    v>0 \mbox{When}; \; \; 
    v\cfrac{\partial f}{\partial x}=v\cfrac
    {f(i,j)-f(i,j-1)}{\Delta y}    \\
    v<0 \mbox{When}; \; \; 
    v\cfrac{\partial f}{\partial x}=v\cfrac
    {f(i,j+1)-f(i,j)}{\Delta y}    

Where, :math:`\Delta x, \Delta y, \Delta t` are the calculation step sizes in the directions of :math:`x, y, t`, respectively, and :math:`(i,j)` are the grid numbers in the directions of :math:`(x,y)`, respectively.
Consequently, for example, when :math:`u>0, \; v>0`, the difference scheme for Equation :eq:`2d_advection` is as follows.

.. math::
    :label: 2d_diff_1

    \cfrac{f_n(i,j)-f(i,j)}{\Delta t}+
    u\cfrac{f(i,j)-f(i-1,j)}{\Delta x}+
    v\cfrac{f(i,j)-f(i,j-1)}{\Delta y}=0

Where, :math:`f_n(i,j)` is the value of :math:`f(i,j)` at :math:`t=t+\Delta t`.
From Equation :eq:`2d_diff_1`, :math:`f_n(i,j)` can be obtained by repeating calculations for the following equation for a step of :math:`\Delta t`.
In this case, the scheme gives a backward difference in the direction of :math:`x` and in the direction of :math:`y`.

.. math::
    :label: 2d_diff_2

    f_n(i,j)=f(i,j)-\Delta t \left\{
    u\cfrac{f(i,j)-f(i-1,j)}{\Delta x}+
    v\cfrac{f(i,j)-f(i,j-1)}{\Delta y}\right\}

An example of calculations by backward difference scheme
----------------------------------------------------------------

This is an example of calculation for when :math:`u>0, \; v>0`. The initial distribution of :math:`f` is conical, as shown in Figure :numref:`3dcone`.


        .. _3dcone:

.. figure:: images/04/3dcone.png
    :width: 70%

    : Initial conditions for the 2D advection calculations

A conical distribution is given by following equation.

.. math::
    :label: cone_profile

    
    r<R \; \mbox{When;} &  f(x,y)=\cfrac{R-r}{R} f_p \\
    r>R \; \mbox{When;} &  f(x,y)=0

Where, 

.. math::
    :label: teimen

    r=\sqrt{(X_p-x)^2+(Y_p-y)^2}
    
Where, :math:`X_p, Y_p` are the locations of peaks, :math:`f_p` is the value of the peak, and :math:`R` is the radius of the bottom.

The Python code and a calculation example 
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The Python code using :eq:`2d_diff_2` is as follows.
https://github.com/YasuShimizu/Python-Code/blob/main/2d_backward_p3d.py
The following figure shows the calculation results calculated by using the code.
In this example, :math:`u=v=0.5` (m/s) and:math:`f_p=0.5, R=20` (m).

        .. _3dcone_backward:

.. figure:: images/04/3dcone_backward.gif
    :width: 95%

    : The calculation results for the 2D advection equation (backward difference scheme)

As with the 1D case, the peak values gradually decrease due to numerical diffusion and the cone gradually flattens while it moves.
The results of this calculation are not entirely accurate.

The upwind difference scheme with varying advection direction
-------------------------------------------------------------------------

The updating formula for :math:`f` in which the advection terms of Equation :eq:`2d_ad_1` and Equation :eq:`2d_ad_2` are expressed with one equation is as follows.

.. Math::

    f_n(i,j)=f(i,j)
    
.. Math::

    -\cfrac{\Delta t}{2 \Delta x}
    \left[(u+|u|)\left\{f(i,j)-f(i-1,j)\right\}+(u-|u|)\left\{f(i+1,j)-f(i,j)\right\} \right. 

.. Math::
    :label: kazakami2d  

    +\left.(v+|v|)\left\{f(i,j)-f(i,j-1)\right\}+(v-|v|)\left\{f(i,j+1)-f(i,j)\right\} \right]

Equation :eq:`kazakami2d` is an upwind difference scheme that can be applied regardless of the sizes or signs of of :math:`u, v`.

Now, let us try to test this 2D upwind difference scheme under the condition where :math:`u, v` vary. We assume that :math:`u,v` is expressed by the following equations.

.. Math::
    :label: circular_uv

    u = & -r_t \omega \sin \theta \\
    v = &  r_t \omega \cos \theta 

The motion expressed by these equations is circular, where :math:`r_t` is the radius of curvature for circular motion, :math:`\omega` is the angular velocity, and :math:`\theta` is the phase angle. 
Assuming that the peak of the cone shown in Figure :numref:`3dcone` moves according to Equation :eq:`circular_uv`, the calculation results are shown by the following figure (:numref:`3d_cone_exact_solution`). 
In this example, :math:`T` is the period of circular motion, and when we assume that :math:`T=` = 100 sec, :math:`r_t` = 30 m and `f_p` = 0.5. 
The relationship between :math:`T` and :math:`\omega` is given by the following equation.

.. Math::
    :label: shuuuki

    \omega=\cfrac{2\pi}{T}

.. _3d_cone_exact_solution:

.. figure:: images/04/CircularMovementCone.gif
    :align: center
    :width: 110%

    : Circular motions of the cone

:numref:`3d_cone_exact_solution` is given merely by moving the distribution of the cone with the velocity of Equation :eq:`circular_uv`.
It is not the difference calculation　result for Equation :eq:`2d_advection`.
However, in a sense, this is the target exact solution for the numerical calculations we will be performing. For reference, the Python code for calculating the animation of :numref:`3d_cone_exact_solution` is below.

https://github.com/YasuShimizu/Python-Code/blob/main/3dcone_circle_m3d.py

I am sorry for the long introduction, but now let us try the upwind difference scheme of the 2D advection equation, which is the goal of this section.

The Python code and a calculation example
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Below is the Python code using Equation :eq:`2d_advection`.

https://github.com/YasuShimizu/Python-Code/blob/main/2d_upwind_circle.py

Using this code, under the initial conditions of Equation :eq:`cone_profile`, :math:`u,v` are varied using Equation :eq:`circular_uv`.
The results are shown in the figure below.

.. _3d_upwind_circular:

.. figure:: images/04/2d_upwind_circle.gif
    :align: center
    :width: 100%

    : Calculating the circular motion for a conical distribution by using the 2D upwind difference scheme

As expected, the upwind difference scheme results in the gradual diffusion of the initial shape and in the attenuation of the peak values.
Although the calculation results are stable, as they were the 1D calculation, the upwind difference scheme is proven not to give accurate calculation results.

Problem (4)
^^^^^^^^^^^^^^

.. list-table:: 

    * - :math:`f` is a variable whose initial condition is the cone-shape distribution. 

Calculate the changes in :math:`f` with time when :math:`f` continuously moves along the sides of the rectangular area.
    
The higher-order difference scheme
===============================================

Here, as I explained for the 1D higher-order difference scheme, I describe the 2D version of the CIP method as a 2D higher-order difference scheme.

Calculating the advection term by using the 2D CIP method
-----------------------------------------------------------------

To simplify the explanation, first, we calculate for when :math:`u<0` and :math:`v<0`.

    .. _2duv:

.. figure:: images/04/2duv.png
   :width: 90%
   :align: center

   : Estimating :math:`f` on the calculation plane

As shown in :numref:`2duv`, set the calculation point numbers in the :math:`x, y` directions as :math:`i,j`.
We assume that points within the range of :math:`i, i+1, j, j+1` (shown by the white dot) at time :math:`t` move during the time of :math:`\Delta t` to the black dot in the figure.
前

As was explained for the 1D calculation, this is a problem for the estimation of the distribution of :math:`f` regarding the black dot in :numref:`2duv`.
How can we estimate the distribution form of :math:`f`?
Therefore, assuming that when... 

.. math::
    :label: 2d_1

    X=-u \Delta t, \; \;  Y=-v \Delta t

When we can estimate the distribution form of :math:`F(X,Y)`, we have... 

.. math::
    :label: 2d_2

    f(i,j,t+\Delta t) = F(X, Y)

With this, we can update the time of :math:`f`.  
Note that, here, :math:`f` is a discrete value for the :math:`(i,j)` grid, and :math:`F` is a continuous function.

Although :math:`F` was assumed to be a curve of the third order for the 1D problem, here we assume it to be a surface of the third order. 
:math:`F(X,Y)` must satisfy the following relationships when :math:`X=0, Y=0`.

.. math::
    :label: 2d_3

    F(X,Y)=f(i,j), \; \; \; F_x(X,Y)=f_x(i,j), \; \; \;F_y(X,Y) = f_y(i,j) 

(Where, :math:`F_x(X,Y)=\displaystyle{{\partial F}\over{\partial x}}`, 
:math:`F_y(X,Y)=\displaystyle{{\partial F}\over{\partial y}} (X,Y)`,
:math:`f_x(i,j)=\displaystyle{{\partial f}\over{\partial x}} (i,j)`,
:math:`f_y(i,j)=\displaystyle{{\partial f}\over{\partial y}} (i,j)`.)

Therefore, we set :math:`F(X,Y)` as follows.

.. math::
    :label: 2d_5

    F(X,Y) =
    \left[ (a_1 X + c_1 Y + e_1) X + g_1 Y + f_x(i,j) \right] X +  \\
    \left[ (b_1 Y + d_1 X + f_1) Y + f_y(i,j) \right] Y + f(i,j)


This equation automatically satisfies Equation :eq:`2d_3`. Other conditions that need to be satisfied are...

.. math::
    :label: 2d_6

    F(0,-\Delta y)=f(i,j+1), \; \; \;
    F(-\Delta x, 0)=f(i+1,j) \\
    F(-\Delta x, -\Delta y) = f(i+1,j+1), \; \; \; 
    F_x(0, -\Delta y)=f_x(i,j+1) \\
    F_x(-\Delta x, 0)=f_x(i+1,j), \; \; \;
    F_x(-\Delta x, -\Delta y) = f_x(i+1,j+1) \\
    F_y(0, -\Delta y)=f_y(i,j+1), \; \; \;
    F_y(-\Delta x, 0)=f_y(i+1,j)

By using these relationships, the factors in Equation :eq:`2d_5` are as follows.

.. math::

    a_1 = {{
    \left[ f_x(i+1,j)+f_x(i,j) \right] \Delta x
    +2 \left[ f(i,j)-f(i+1,j) \right]}\over{\Delta x^3}}

.. math::

    b_1 = {{
    \left[ f_y(i,j+1)+f_y(i,j) \right] \Delta y
    +2 \left[ f(i,j)-f(i,j+1) \right]}\over{\Delta y^3}}

.. math::

    c_1 = {{-\left\{
    f(i,j)-f(i,j+1)-f(i+1,j)+f(i+1,j+1)\right\}-\left[f_x(i,j+1)-f_x(i,j) \right] \Delta x
    }\over{\Delta x^2 \Delta y}}

.. math::

    d_1 = {{-\left\{f(i,j)-f(i,j+1)-f(i+1,j)+f(i+1,j+1)\right\}-
    \left[ f_y(i+1,j)-f_y(i,j) \right] \Delta y}\over{\Delta x \Delta y^2}}

.. math::

    e_1 = {{
    3 \left[ f(i+1,j)-f(i,j) \right] +
    \left[ f_x(i+1,j)+2 f_x(i,j) \right] \Delta x
    }\over{\Delta x^2}}

.. math::

    f_1 = {{
    3 \left[ f(i,j+1)-f(i,j) \right] +
    \left[ f_y(i,j+1)+2 f_y(i,j) \right] \Delta y
    }\over{\Delta y^2}}

.. math::
    :label: coefs

    g_1 =
    {{-\left\{f_y(i+1,j)-f_y(i,j)\right\}+c_1 \Delta x^2
    }\over{\Delta x}}

Problem (5)
^^^^^^^^^^^^^^

.. list-table:: 

    * - Derive Equation :eq:`coefs` to check whether the expression of :math:`a_1 \sim g_1` is right.

Calculating the advection term with the 2D CIP method (generalization)
-----------------------------------------------------------------------------------------------

The explanation thus far is for :math:`u<0` and :math:`v<0`. Actually, there are eight possible combinations for the combinations of signs of :math:`u, v`.
However, it would be very complicated to formulate and program all the combinations.
Therefore, we introduce a sign function (:math:`\mbox{sign}`) as we did for the ID calculation.


.. math::
    :label: advect2d_1

    i_s = \mbox{sign}(u), \; \; \;  & j_s = \mbox{sign}(v) \\
    i_m = i-i_s, \; \; \;  & j_m = j-j_s 


:math:`F(X,Y)`, the distribution form of the variables between grids, is defined as :eq:`advect2d_2` by defining :math:`X, Y` as :numref:`2duv_1`, and by transforming Equation :eq:`fxnew` to the second order. 

    .. _2duv_1:

.. figure:: images/04/2duv_1.png
   :width: 80%
   :align: center

   : Defining :math:`X, Y` on the calculation plane 

.. math::
    :label: advect2d_2

    F(X,Y) = \left[ (a_1 X + c_1 Y + e_1) X + g_1 Y + f_x(i,j) \right] X + \\
    \left[ (b_1 Y + d_1 X + f_1) Y + f_y(i,j) \right] Y + f(i,j)


Where, the factors :math:`a_1 \sim g_1` are expressed as follows by using :math:`i_s, j_s, i_m, j_m`, regardless of the signs of :math:`u, v`. 

.. math::

    a_1 = {{
    i_s \left[ f_x(i_m,j)+f_x(i,j) \right] \Delta x
    -2 \left[ f(i,j)-f(i_m,j) \right]}\over{i_s \Delta x^3}}

.. math::

    b_1 = {{
    j_s \left[ f_y(i,j_m)+f_y(i,j) \right] \Delta y
    -2 \left[ f(i,j)-f(i,j_m) \right]}\over{j_s \Delta y^3}}

.. math::

    c_1 = {{
    -\left\{f(i,j)-f(i,j_m)-f(i_m,j)+f(i_m,j_m)\right\}-i_s \left[f_x(i,j_m)-f_x(i,j) \right] \Delta x
    }\over{j_s \Delta x^2 \Delta y}}

.. math::
 
    d_1 = {{
    -\left\{f(i,j)-f(i,j_m)-f(i_m,j)+f(i_m,j_m)\right\}-
    j_s \left[ f_y(i_m,j)-f_y(i,j) \right] \Delta y}\over{i_s \Delta x \Delta y^2}}

.. math::

    e_1 = {{
    3 \left[ f(i_m,j)-f(i,j) \right] +
    i_s \left[ f_x(i_m,j)+2 f_x(i,j) \right] \Delta x
    }\over{\Delta x^2}}

.. math::

    f_1 = {{
    3 \left[ f(i,j_m)-f(i,j) \right] +
    j_s \left[ f_y(i,j_m)+2 f_y(i,j) \right] \Delta y
    }\over{\Delta y^2}}

.. math::
    :label: coefs_1
    
    g_1 =
    {{-\left\{f_y(i_m,j)-f_y(i,j)\right\}+c_1 \Delta x^2
    }\over{i_s \Delta x}}

Problem (6)
^^^^^^^^^^^^^^

.. list-table:: 

    * - Derive Equation :eq:`coefs_1` to check whether the expression of :math:`a_1 \sim g_1` is right.

Updating the differentiation amount
-----------------------------------------------------------

Under the CIP method, in addition to the variable :math:`f`, its differentiation, :math:`f_x, f_y`, also needs to be updated from time to time.
Now, we perform partial differentiation for both sides of Equation :eq:`2d_advection` with respect to :math:`x`.

.. math::
    :label: bibun2_1

    {{\partial^2 f}\over{\partial t \partial x}}
    + {{\partial u}\over{\partial x}} {{\partial f}\over{\partial x}}
    + u {{\partial^2 f}\over{\partial x^2}}
    + {{\partial v}\over{\partial x}} {{\partial f}\over{\partial y}}
    + v {{\partial^2 f}\over{\partial x \partial x}} =0


Assuming that :math:`\displaystyle{{{\partial f}\over{\partial x}}}` is :math:`f_x`, we have...

.. math::
    :label: bibun2_2

    {{\partial f_x}\over{\partial t}}
    + u{{\partial f_x}\over{\partial x}}
    + v{{\partial f_x}\over{\partial y}}
    = -\left( f_x {{\partial u}\over{\partial x}} +
    f_y {{\partial v}\over{\partial x}} \right)

As we can see, the left side is a 2D advection equation with respect to :`f_x`.
This means that, here, we can directly apply the 2D CIP method.
When :math:`u, v` are constants, the right side is 0.
When they are variables, we can apply [ :ref:`label1` ], which is explained under the 1D CIP method.
The calculation should apply the split method. The central difference scheme can be applied for the difference form on the right side.

Similarly, when we perform partial differentiation for Equation :eq:`2d_advection` with respect to :math:`y`, and we assume that :math:`\displaystyle{{{\partial f}\over{\partial y}}}` is :math:`f_y`, then we have...

.. math::
    :label: bibun2_3

    {{\partial f_y}\over{\partial t}}
    + u{{\partial f_y}\over{\partial x}}
    + v{{\partial f_y}\over{\partial y}}
    = -\left( f_x {{\partial u}\over{\partial y}} +
    f_y {{\partial v}\over{\partial y}} \right)

This is the same form as Equation :eq:`bibun2_2`. 
Therefore, the advection term on the left side can be calculated with the CIP method, and the source term on the right side can be calculated with a central difference scheme.

The Python code and a calculation example
-----------------------------------------------------------------------

The Python code for the 2D CIP method mentioned above is found at...

https://github.com/YasuShimizu/2d_advection_cip

.. _3dcone_cip:

The calculation that results from the code is shown below. The calculation conditions are the same as the aforementioned :ref:`2dkazakami`; that is, a cone distribution with circular motion.

.. figure:: images/04/2d_cip_circle.gif
    :width: 95%
    :align: center

    : The calculation results for the 2D advection equation (the CIP method)

The calculation results show that the circular motion continues and that the cone shape is largely preserved. 
The CIP method is a numerical calculation method for the 2D advection equation and gives highly accurate calculation results.