[Python] Fluid simulation: From linear to non-linear

Introduction

In addition to studying computational fluid dynamics (CFD), I would like to summarize the knowledge necessary for constructing computational fluid dynamics (CFD) for liquids. I feel that computational fluid dynamics is a rather difficult discipline, so I would like to write it in a way that is easy for beginners to understand. It seems that there are many mistakes, so please contact us if you find it. Also, I would appreciate it if you could comment on what is difficult to understand. We will update it from time to time.

Target audience

series

Rough content of this article

We will explain how to handle non-linear equations as a preliminary step to construct the basic equations of fluid required for water simulation. Finally, we perform a numerical simulation of the Burgers equation, which is a nonlinear equation.

Library used: Numpy, Scipy, Matplotlib

I will make this 2D version with this guy.

burgers.gif

table of contents

chapter title Remarks
1. Linear and non-linear
1-2. Features of linear and nonlinear equations
1-3. How to distinguish between linear and non-linear Trivia
2. Implementation of linear convection-diffusion equation
2-1. Linear convection-diffusion equation
2-2. Implementation
3. Implementation of Burgers equation
3-1. What is Burgers' equation?? The equivalent of a non-linear convection-diffusion equation
3-2. Implementation
3-3. Supplement of discretization formula
4. Extend to 2D
4-1. Two-dimensional advection equation
4-2. Two-dimensional diffusion equation
4-3. 2D Burgers equation

1. Linear and non-linear

1-1. Characteristics of linear and non-linear equations

The following is a brief summary of linear and non-linear equations.

equation Feature
linear 1.Proportional relationship holds
2. Can find a solution
non-linear 1.Proportional relationship does not hold
2. Basically I can't find a solution(Only a few equations such as the kdV equation can be solved)

Fluid equations to be dealt with in the future are classified as non-linear equations, and basically no solution can be obtained. Therefore, when finding the solution of a nonlinear equation, the basic policy is to perform numerical calculation using the discretized equations introduced in the previous articles and find the solution approximately. ..

1-2. How to distinguish between linear and non-linear

Apart from the subject, I will explain how to distinguish between linear and non-linear equations. To distinguish between linear and non-linear equations, determine whether the solution superposition holds.

I would like to explain using the advection equation introduced in Previous article as an example.

The advection equation is $ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $ It was represented by. However, the advection velocity c is a constant (meaning that the substance moves at a constant velocity of c).

When there are two solutions $ u_1 $ and $ u_2 $ in this equation, respectively

\frac{\partial u_1}{\partial t} + c \frac{\partial u_1}{\partial x} = 0
\frac{\partial u_2}{\partial t} + c \frac{\partial u_2}{\partial x} = 0

Is established. Therefore, for $ u_1 + u_2 , which is the sum of the two solutions, $ \frac{\partial (u_1+u_2)}{\partial t} + c \frac{\partial (u_1+u_2)}{\partial x} = 0 $$ Is established, and we can see that this equation is linear.

Now, let's change the advection velocity c, which is a constant, to the variable u.

\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0
\frac{\partial (u_1+u_2)}{\partial t} + (u_1+u_2) \frac{\partial (u_1+u_2)}{\partial x} \neq \frac{\partial (u_1+u_2)}{\partial t} + u_1 \frac{\partial u_1}{\partial x}+ u_2 \frac{\partial u_2}{\partial x} = 0

Then, in this equation, $ u_1 + u_2 $, which is the sum of the two solutions, is not a solution, so we know that it is a non-linear equation.

In this article, I would like to deal with such nonlinear equations.

2. Implementation of linear convection-diffusion equation

2-1. Linear convection-diffusion equation

First of all, I would like to deal with the linear convection-diffusion equation, which also serves as the developmental meaning of the series so far. This equation is like a combination of the convection equation (the equation that means mass transfer) and the diffusion equation (the equation that means the uniformity of physical quantities) as shown in the equation below (c and). $ \ nu $ is a constant). By the way, the non-linear version of this is the Burgers equation that will appear in the next chapter.

2-2. Implementation

For the time being, we will implement it using the CIP method and the non-stationary iterative method. Simply put, the CIP method is a method of predicting future values from the current value and its derivative value. For details, see the article and this pdf I wrote so far. Please see jp / netlab / mhd_summer2001 / cip-intro.pdf). The unsteady iterative method is a method for finding the solution of one-dimensional simultaneous equations while changing the coefficients. If you want to know more, I wrote the Article on Diffusion Equations. See / items / 97daa509991bb9777a4a) and Article on Python's Transsteady Iterative Library.

It's not that difficult to do. In the CIP method, the advection phase and the non-advection phase are considered separately. In short, we just calculate u from the advection equation and use the calculation result u to solve the diffusion equation. In terms of image, the speed of advection and the speed of uniformization by diffusion? (Frequency?) Are different, so I think it's like thinking separately. The convection-diffusion equation can be expressed as follows by dividing it into an advection phase and a non-advection phase. When using the CIP method, information on the differential value $ \ frac {\ partial u} {\ partial x} $ is also required, so it is performed at the same time as the calculation of u.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 \\
    \frac{\partial (\partial u / \partial x)}{\partial t} + c \frac{\partial (\partial u / \partial x)}{\partial x} = 0
 \end{array}
\right.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t}  = \nu \frac{\partial^2 u}{\partial x^2} \\
    \frac{\partial (\partial u / \partial x)}{\partial t} = \frac{\partial \left(\nu \frac{\partial^2 u}{\partial x^2}\right)}{\partial x}
 \end{array}
\right.

linear_adv_diff.gif

import sys
import os
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.sparse
import scipy.sparse.linalg as spla
from scipy.sparse import csr_matrix
from matplotlib.animation import ArtistAnimation
from mpl_toolkits.mplot3d import Axes3D

# Make stencils
# Creat square wave
Num_stencil_x = 101
x_array = np.arange(Num_stencil_x)
u_array = np.where(((x_array >= 30) |(x_array < 10)), 0.0, 1.0)
u_lower_boundary = 0.0
u_upper_boundary = 0.0
Time_step = 300
Delta_x = max(x_array) / (Num_stencil_x-1)
C = 1
Nu = 0.5
Delta_t = 0.2
CFL = C * Delta_t / Delta_x
diffusion_number = Nu * Delta_t / Delta_x**2
total_movement = C * Delta_t * (Time_step+1)
exact_u_array = np.where(((x_array >= 30 + total_movement) |(x_array < 10 + total_movement)), 0.0, 1.0)
plt.plot(x_array, u_array, label="Initial condition")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)
print("Δx:", Delta_x, "Δt:", Delta_t, "CFL:", CFL, "ν:", Nu, "d:", diffusion_number)

def solve_advection_cip(u_cip, partial_u_cip, lower_boundary, velocity=C):
    u_old = u_cip.copy()
    partial_u_old = partial_u_cip.copy()
    u_cip[0] = lower_boundary
    partial_u_cip[0] = 0  #The boundary value was set to 0
    a = (partial_u_old[1:] + partial_u_old[:-1]) / Delta_x**2 - 2.0 * np.diff(u_old, n=1) / Delta_x**3
    b = 3 * (u_old[:-1] - u_cip[1:]) / Delta_x**2 + (2.0*partial_u_old[1:] + partial_u_old[:-1]) / Delta_x
    c = partial_u_old[1:]
    d = u_old[1:]
    xi = - velocity * Delta_t  # C > 0
    u_cip[1:] = a * xi**3 + b * xi**2 + c * xi + d
    partial_u_cip[1:] = 3 * a * xi**2 + 2 * b * xi + c
    return u_cip, partial_u_cip

def solve_diffusion(u_array, upper_boundary, lower_boundary, diffusion_number):
    a_matrix = np.identity(len(u_array)) * 2 *(1/diffusion_number+1) \
                - np.eye(len(u_array), k=1) \
                - np.eye(len(u_array), k=-1)
    temp_u_array = np.append(np.append(
                        lower_boundary, 
                        u_array), upper_boundary)
    b_array = 2 * (1/diffusion_number - 1) * u_array + temp_u_array[2:] + temp_u_array[:-2]
    b_array[0] += lower_boundary
    b_array[-1] += upper_boundary
    a_csr = scipy.sparse.csr_matrix(a_matrix)
    u_array = spla.isolve.bicgstab(a_csr, b_array)[0]
    return u_array

fig = plt.figure()
images = []
u_cip= u_array.copy()
partial_u_cip = ((np.append(u_cip[1:], u_upper_boundary) + u_cip)/2 - (np.append(u_lower_boundary, u_cip[:-1]) + u_cip)/2)/ Delta_x
#Evolve time
for n in range(Time_step):
    u_cip_star, partial_u_cip_star = solve_advection_cip(u_cip, partial_u_cip, u_lower_boundary)
    u_cip = solve_diffusion(u_cip_star, u_upper_boundary, u_lower_boundary, diffusion_number)
    partial_u_cip = solve_diffusion(partial_u_cip_star, 0, 0, diffusion_number)  #The boundary value was set to 0
    if n % 10 == 0:
        line, = plt.plot(x_array, u_cip, "r")
        images.append([line])
# plt.plot(x_array, u_array, label="Initial condition")
# plt.plot(x_array, u_cip, label="CIP method")
animatoin_gif = ArtistAnimation(fig, images)
animatoin_gif.save('linear_adv_diff.gif', writer="pillow")
plt.show()

3. Implementation of Burgers equation

3-1. What is Burgers' equation?

One of the non-linear PDEs, which is a non-linear version of the convection-diffusion equation explained in the previous chapter. If you make assumptions about the fluid equation and simplify it considerably, it becomes Burgers' equation. This equation is a non-linear equation, but it is known that an exact solution can be obtained by using a call hop transformation. Therefore, it is often used to examine the validity of numerical calculation methods.

3-2. Implementation

Implements Burgers' equation. This is also implemented by the CIP method. It is as follows when it is divided into an advection phase and a non-advection phase.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0 \\
    \frac{\partial (\partial u / \partial x)}{\partial t} + u \frac{\partial (\partial u / \partial x)}{\partial x} = 0
 \end{array}
\right.

\left\{
\begin{array}{l}
    \frac{\partial u}{\partial t}  = \nu \frac{\partial^2 u}{\partial x^2} \\
    \frac{\partial (\partial u / \partial x)}{\partial t} = \frac{\partial \left(\nu \frac{\partial^2 u}{\partial x^2}\right)}{\partial x} - \left(\frac{\partial u}{\partial x} \right)^2
 \end{array}
\right.

There are not many changes from the linear convection-diffusion equation. Where to change

I think about it. We've made a few modifications to the discretization, see the next section for more details.

A brief discussion of the calculation results reveals that the higher the height, the faster the velocity, so the upper side of the wave catches up with the lower side and forms a steep discontinuity similar to a shock wave. After that, you can see that the solution gradually becomes dull due to the effect of diffusion.

burgers.gif

def solve_advection_cip(u_cip, partial_u_cip, lower_boundary, dt=Delta_t, velocity=C):
    """
Solving the advection equation by the CIP method
    """
    u_old = u_cip.copy()
    partial_u_old = partial_u_cip.copy()
    u_cip[0] = lower_boundary
    partial_u_cip[0] = (u_cip[0] - lower_boundary) / Delta_x
    a = (partial_u_old[1:] + partial_u_old[:-1]) / Delta_x**2 - 2.0 * np.diff(u_old, n=1) / Delta_x**3
    b = 3 * (u_old[:-1] - u_cip[1:]) / Delta_x**2 + (2.0*partial_u_old[1:] + partial_u_old[:-1]) / Delta_x
    c = partial_u_old[1:]
    d = u_old[1:]
    xi = - velocity * dt  # C > 0
    u_cip[1:] = a * xi**3 + b * xi**2 + c * xi + d
    partial_u_cip[1:] = 3 * a * xi**2 + 2 * b * xi + c
    return u_cip, partial_u_cip

def solve_diffusion(u_array, lower_boundary, upper_boundary, diffusion_number, 
                    update_partial=False, u_array_integral=None, prop_lambda=1/2):
    """
Solve the diffusion equation.
    
    Parameters
    ----------
    u_array : array-like
n vector
    upper_boundary : numpy.float64
        u_array[n]Next to
    lower_boundary : numpy.float64
        u_array[0]Next to
    diffusion_number : numpy.float64
Diffusion number. Positive number.
    update_partial : Bool, default False
Set to True when updating the differential expression.
    u_array_integral : array-like, default None
Used when updating the differential formula.[lower_boundary, u_array, upper_boudary]N in the form of+Vector of 2.
    prop_lambda : numpy.float64, default 1/2
        0: Explicit method. Centered sapce difference.
        1/2:Semi-implicit method. Crank-Nicolson scheme
        1: Fully implicit method.Time retreat difference.

    Returns
    -------
    u_array : array-like
n-row matrix
    """
    a_matrix = np.identity(len(u_array)) * (1/diffusion_number+2 * prop_lambda) \
                - prop_lambda * np.eye(len(u_array), k=1) \
                - prop_lambda * np.eye(len(u_array), k=-1)
    temp_u_array = np.append(np.append(
                             lower_boundary, 
                             u_array), upper_boundary)
    b_array = (1/diffusion_number - 2 * (1-prop_lambda)) * u_array + (1-prop_lambda)*temp_u_array[2:] + (1-prop_lambda)*temp_u_array[:-2]
    b_array[0] += prop_lambda * lower_boundary
    b_array[-1] += prop_lambda * upper_boundary
    if update_partial:
        #Do this when updating the derivative
        b_array -= ((u_array_integral[2:] - u_array_integral[:-2]) / 2 / Delta_x)**2
    a_csr = scipy.sparse.csr_matrix(a_matrix)
    u_array = spla.isolve.bicgstab(a_csr, b_array)[0]
    return u_array

fig = plt.figure(1)
images = []
u_cip= u_array.copy()
partial_u_cip = ((np.append(u_cip[1:], u_upper_boundary) + u_cip)/2 
                 - (np.append(u_lower_boundary, u_cip[:-1]) + u_cip)/2)/ Delta_x
#Setting stability conditions
cfl_condition = 1
diffusion_number_restriction = 1/2
#Evolve time
for n in range(Time_step):
    delta_t = Delta_t
    cfl_max = np.max(np.abs(u_cip)) * delta_t / Delta_x
    diffusion_max = Nu * delta_t / Delta_x**2
    if cfl_max > cfl_condition:
        #Judgment of CFL conditions
        # cfl_If it is larger than the condition, reduce the time step width dt so that the CFL condition is satisfied.
        delta_t = cfl_condition * Delta_x / np.max(np.abs(u_cip))
    if diffusion_number > diffusion_number_restriction:
        #Judgment of diffusion number
        # diffusion_number_If it is larger than the restriction, reduce the time step width dt so that the diffusion number condition is satisfied.
        delta_t = diffusion_max * Delta_x ** 2 / Nu
    #Solve the advection equation
    u_cip_star, partial_u_cip_star = solve_advection_cip(u_cip, partial_u_cip, u_lower_boundary, delta_t, u_cip[1:])
    #Solve the diffusion equation
    u_cip = solve_diffusion(u_cip_star, u_lower_boundary, u_upper_boundary, diffusion_number)
    u_cip_with_boundary = np.append(np.append(
                                     u_lower_boundary, 
                                     u_cip_star), 
                                     u_upper_boundary)
    partial_u_cip = solve_diffusion(partial_u_cip_star, 0, 0, diffusion_number, 
                                    update_partial=True, 
                                    u_array_integral=u_cip_with_boundary)  #The boundary value was set to 0
    if n % 10 == 0:
        line, = plt.plot(x_array, u_cip, "r")
        images.append([line])
# plt.plot(x_array, u_array, label="Initial condition")
# plt.plot(x_array, u_cip, label="Burgers")
animatoin_gif = ArtistAnimation(fig, images)
animatoin_gif.save('burgers.gif', writer="pillow")

3-3. Discretization formula of diffusion phase

Diffusion equation

  \frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2}

Is expressed by a discretized equation

\frac{u_j^{n+1} - u_j^n}{\Delta t} = \nu \left((1 - \lambda) \frac{u_{j+1}^n - 2 u_j^n + u_{j-1}^n }{\Delta x^2} + \lambda \frac{u_{j+1}^{n+1} - 2 u_j^{n+1} + u_{j-1}^{n+1} }{\Delta x^2}\right)

It will be. Considering the future, unlike Article on Diffusion Equation, the discretization equation is expressed using $ \ lambda $. To summarize the value of this $ \ lambda $

\lambda Discretization method Feature
0 Center difference Explicit method. Calculate the time difference using only the current time value.
1/2 Crank Nicholson's implicit method Semi-implicit method. Calculate the time difference using half the current value and half the value at the next time.
1 Time retreat difference Completely implicit method. Calculate the time difference using only the next time value.

It will be. The bottom line is that you can use $ \ lambda $ to implement three types of discretization techniques.

If you transform this and bring the value of time n + 1 to the left side $ -\lambda u_{j+1}^{n+1} +\left(\frac{1}{d} + 2 \lambda \right) u_j^{n+1} - \lambda u_{j-1}^{n+1} = (1-\lambda) u_{j+1}^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_j^{n} + (1 - \lambda) u_{j-1}^{n} \\ d = \nu \frac{\Delta t}{\Delta x^2} $

Assuming that the grid points exist in the range of 1 to M, the boundary value is represented by $ u_0, u_ {M + 1} $, and it is converted to a matrix display.

  \left(
    \begin{array}{cccc}
      \frac{1}{d} + 2 \lambda & -\lambda                      &  0                      & \ldots   & \ldots                   & 0 \\
      -\lambda                & \frac{1}{d} + 2 \lambda       & -\lambda                & 0        & \ldots                   & 0 \\
      0                       &-\lambda                       & \frac{1}{d} + 2 \lambda & -\lambda & 0                         & \ldots  \\
      \vdots                  & \ddots                        & \ddots                  & \ddots   & \ddots                   & \ddots\\
      0                       & 0                             & \ddots                  & -\lambda & \frac{1}{d} + 2 \lambda & -\lambda \\
      0                       & 0                             & \ldots                  & 0        & -\lambda                 & \frac{1}{d} + 2 \lambda
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      u_1^{n+1}  \\
      u_2^{n+1}  \\
      u_3^{n+1}    \\
      \vdots \\
      u_{M-1}^{n+1} \\
      u_M^{n+1} 
    \end{array}
  \right)

   = 
   \left(
    \begin{array}{c}
      (1 - \lambda) u_2^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_1^{n} + \left((1-\lambda) u_0^n + \lambda u_0^{n+1}\right) \\
      (1-\lambda) u_3^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_2^{n} + (1-\lambda) u_1^n  \\
      (1-\lambda) u_4^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_3^{n} + (1-\lambda)u_2^n  \\
      \vdots \\
      (1-\lambda)u_M^{n} + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_{M-1}^{n} + (1-\lambda)u_{M-2}^n  \\
      \left((1-\lambda)u_{M+1}^n + \lambda u_{M+1}^{n+1}\right) + \left(\frac{1}{d} - 2 (1 - \lambda) \right) u_{M}^{n} + (1-\lambda)u_{M-1}^n
    \end{array}
  \right)

If you want to update $ \ frac {\ partial u} {\ partial x} $, change u in this matrix to $ \ frac {\ partial u} {\ partial x} $ and move it to the right-hand side.

\left(
    \begin{array}{c}
      - \left(\frac{u_2^{n} - u_0^{n}}{2 \Delta x} \right)^2  \\
      - \left(\frac{u_3^{n} - u_1^{n}}{2 \Delta x} \right)^2  \\
      - \left(\frac{u_4^{n} - u_2^{n}}{2 \Delta x} \right)^2    \\
      \vdots \\
      - \left(\frac{u_M^{n} - u_{M-1}^{n}}{2 \Delta x} \right)^2 \\
      - \left(\frac{u_{M+1}^{n} - u_M^{n}}{2 \Delta x} \right)^2 
    \end{array}
  \right)

Add.

This is implemented in the previous section. The BiCG STAB method was used as the iterative solution method. The reason is that the name is cool.

4. Extend to 2D

I want to draw a cool picture, so I will make it two-dimensional in this chapter. It's just mathematical expansion.

4-1. Two-dimensional advection equation

As in the example, we will implement it using the CIP method. For detailed formula expansion, refer to 2D CIP method in the reference. Since there are 10 unknowns, we only need to make 10 equations. You don't have to think about the (i + 1, j + 1) coordinates at the time of differentiation. In the references, I've done it well to reduce the amount of calculation, but it's difficult to understand, so I'll do it honestly here.

Just like in one dimension

F(X,Y) = a_3 X^3 + a_2 X^2 + a_1 X + b_3 Y^3 + b_2 Y^2 + b_1 Y + c_3 X^2 Y + c_2 X Y^2 + c_1 XY + d \quad ,where \quad X = x - x_j

Define $ F (X, Y) $ with a third-order completion like this. From the condition that the function value and the derivative value are continuous on the grid point **

F(0, 0) = f_{i,j}, \quad F(\Delta x, 0) = f_{i+1, j}, \quad F(0, \Delta y) = f_{i, j+1}, \quad F(\Delta x, \Delta y) = f_{i+1, j+1}\\ \frac{\partial F(0,0)}{\partial x} = \frac{\partial f_{i,j}}{\partial x}, \quad \frac{\partial F(\Delta x, 0)}{\partial x} = \frac{\partial f_{i+1, j}}{\partial x},\quad \frac{\partial F(0, \Delta y)}{\partial x} = \frac{\partial f_{i, j+1}}{\partial x} \\ \frac{\partial F(0,0)}{\partial y} = \frac{\partial f_{i,j}}{\partial y}, \quad \frac{\partial F(\Delta x, 0)}{\partial y} = \frac{\partial f_{i+1, j}}{\partial y},\quad \frac{\partial F(0, \Delta y)}{\partial y} = \frac{\partial f_{i, j+1}}{\partial y} \\ , where \quad \Delta x = x_{i+1} - x_i, \quad \Delta y = y_{j+1} - y_{j}

The 10 formulas of However, the value at the grid point (i, j) of the function f and the differential value are expressed as $ f_ {i, j}, \ frac {\ partial f_ {i, j}} {\ partial x} $, respectively. I will. Substituting this into the formula of $ F (X, Y) $ to find the coefficient, $ a_3 = \frac{\partial f_{i+1,j} / \partial x + \partial f_{i, j} / \partial x }{\Delta x^2} + \frac{2 (f_{i,j} - f_{i+1, j})}{\Delta x^3} \\ a_2 = \frac{3 (f_{i+1, j} - f_{i,j})}{\Delta x^2} - \frac{(2 \partial f_{i,j} / \partial x + \partial f_{i+1, j} / \partial x )}{\Delta x} \\ a_1 = \frac{\partial f_{i,j}}{\partial x}\\ b_3 = \frac{\partial f_{i,j+1} / \partial y + \partial f_{i, j} / \partial y }{\Delta y^2} + \frac{2 (f_{i,j} - f_{i, j+1})}{\Delta y^3} \\ b_2 = \frac{3 (f_{i, j+1} - f_{i,j})}{\Delta y^2} - \frac{(2 \partial f_{i,j} / \partial y + \partial f_{i, j+1} / \partial y )}{\Delta y} \\ b_1 = \frac{\partial f_{i,j}}{\partial y}\\ c_3 = \frac{\partial f_{i+1,j}/\partial y - \partial f_{i,j}/\partial y }{(\Delta x)^2} - \frac{c_1}{\Delta x}\\ c_2 = \frac{\partial f_{i,j+1}/\partial x - \partial f_{i,j}/\partial x }{(\Delta y)^2} - \frac{c_1}{\Delta y}\\ c_1 = - \frac{f_{i+1,j+1} - f_{i, j+1} - f_{i+1,j} + f_{i, j}}{\Delta x \Delta y} + \frac{\partial f_{i+1,j}/\partial y - \partial f_{i,j}/\partial y }{\Delta x} + \frac{\partial f_{i,j+1}/\partial x - \partial f_{i,j}/\partial x }{\Delta y}\\ d = f_{i,j} $

From the above, when the speed is constant, $ \ frac {\ partial u} {\ partial t} + c_1 \ frac {\ partial u} {\ partial x} + c_2 \ frac {\ partial u} {\ partial y} Since = 0 $ also satisfies the differential value, the value and the differential value after the time step width $ \ Delta t (n + 1 step) $ seconds are

    u_{i,j}^{n+1} = F(x_i - c_1 \Delta t, y_j -c_2 \Delta t)= a_3 \xi^3 + a_2 \xi^2 + a_1 \xi + b_3 \zeta^3 + b_2 \zeta^2 + b_1 \zeta + c_3 \xi^2 \zeta + c_2 \xi \zeta^2 + c_1 \xi \zeta + d \\
    \frac{\partial u_{i,j}^{n+1}}{\partial x} = \frac{\partial F(x_i - c_1 \Delta t, y_j - c_2 \Delta t)}{\partial x} = 3 a_3 \xi^2 + 2 a_2 \xi + a_1 + 2 c_3 \xi \zeta c_2 \zeta^2 + c_1 \zeta \\
    \frac{\partial u_{i,j}^{n+1}}{\partial y} = \frac{\partial F(x_i - c_1 \Delta t, y_j - c_2 \Delta t)}{\partial y} = 3 b_3 \zeta^2 + 2 b_2 \zeta + b_1 + 2 c_2 \xi \zeta c_3 \xi^2 + c_1 \xi \\
    , where \quad \xi = - c_1 \Delta t, \quad \zeta = - c_2 \Delta t \quad (c_1, c_2 < 0)

If the speed is not constant $ \ frac {\ partial u} {\ partial t} + u \ frac {\ partial u} {\ partial x} + v \ frac {\ partial u} {\ partial y} = 0 $ The extra terms that appear after spatial differentiation are calculated as non-advection terms.

2d-advection.gif

4-2. Two-dimensional diffusion equation

Two-dimensional diffusion equation

    \frac{\partial u}{\partial t} = \nu \left\{ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right\}

Is expressed by a discretized equation

  \frac{u_{i,j}^{n+1} - u_{i,j}^n}{\Delta t} = \nu \left\{(1 - \lambda) \left( \frac{u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n }{\Delta x^2} + \frac{u_{i,j+1}^n - 2 u_{i,j}^n + u_{i,j-1}^n }{\Delta y^2} \right) + \lambda \left( \frac{u_{i+1,j}^{n+1} - 2 u_{i,j}^{n+1} + u_{i-1,j}^{n+1} }{\Delta x^2} + \frac{u_{i,j+1}^{n+1} - 2 u_{i,j}^{n+1} + u_{i,j-1}^{n+1} }{\Delta y^2} \right) \right\}

Sort the grids into a row to make this easier to implement in your code. Set the grid point of the entire system to $ M = NX \ times NY $, and set the position of the grid point (i, j) to the mth (start 0 according to the Python notation) counting from the grid point (0,0). I will. Then, the grid points (i + 1, j) are m + 1th, and the grid points (i, j + 1) are m + NXth. Applying this to the above formula, if you bring the value of time n + 1 to the left side

- d_2 \lambda u_{m+NX}^{n+1} - d_1 \lambda u_{m+1}^{n+1}  +s u_{m}^{n+1} - d_1 \lambda u_{m-1}^{n+1} - d_2 \lambda u_{m-NX}^{n+1} = d_2 (1-\lambda) u_{m+NX}^{n} + d_1 (1-\lambda) u_{m+1}^{n} + t u_m^{n} + d_1 (1 - \lambda) u_{m-1}^{n} + d_2 (1 - \lambda) u_{m-NX}^{n} \\
,where \quad s = 1 + 2 \lambda (d_1 + d_2), \quad t = 1 - 2 (1-\lambda) (d_1 + d_2)\\
d_1 = \nu \frac{\Delta t}{\Delta x^2}, \quad d_2 = \nu \frac{\Delta t}{\Delta y^2}

So, if you implement it including the boundary value, you're done. Since the matrix display is long, I will skip it here.

2d-diffusion.gif

4-3. Two-dimensional Burgers equation

Let's make it two-dimensional for the time being. Change the notation and change u to f to make a distinction from speed.

\left\{
\begin{array}{l}
    \frac{\partial f}{\partial t} + u \frac{\partial f}{\partial x} + v \frac{\partial f}{\partial y} = 0 \\
    \frac{\partial (\partial f / \partial x)}{\partial t} + u \frac{\partial (\partial f / \partial x)}{\partial x} + v \frac{\partial (\partial f / \partial x)}{\partial x}= 0  \\
    \frac{\partial (\partial f / \partial y)}{\partial t} + u \frac{\partial (\partial f / \partial y)}{\partial x} + v \frac{\partial (\partial f / \partial y)}{\partial x}= 0
 \end{array}
\right.

\left\{
\begin{array}{l}
    \frac{\partial f}{\partial t}  = \nu \left( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}\right) \\
    \frac{\partial (\partial f / \partial x)}{\partial t} = \frac{\partial \left\{ \nu \left( \frac{\partial^2 f}{\partial x^2}　+ \frac{\partial^2 f}{\partial y^2}\right)\right\}}{\partial x} - \left(\frac{\partial u}{\partial x} \right)\left(\frac{\partial f}{\partial x} \right) - \left(\frac{\partial v}{\partial x} \right) \left(\frac{\partial f}{\partial y} \right)  \\
      \frac{\partial (\partial f / \partial y)}{\partial t} = \frac{\partial \left\{ \nu \left( \frac{\partial^2 f}{\partial x^2}　+ \frac{\partial^2 f}{\partial y^2}\right)\right\}}{\partial y} - \left(\frac{\partial u}{\partial y} \right)\left(\frac{\partial f}{\partial x} \right) - \left(\frac{\partial v}{\partial y} \right) \left(\frac{\partial f}{\partial y} \right) 
\end{array}
\right.

Just implement this. As shown below, I don't know if it matches, but I get a calculation result like that.

2d-burgers_short.gif

# Make stencils
# Creat square wave
Num_stencil_x = 101
Num_stencil_y = 101
x_stencil = np.arange(Num_stencil_x)
y_stencil = np.arange(Num_stencil_y)
x_array, y_array = np.meshgrid(x_stencil, y_stencil)
u_array = np.where(((30 >x_array) & (x_array >= 10) & (30 > y_array) & (y_array >= 10)), 1.0, 0.0)
u_x_lower_boundary = np.zeros(Num_stencil_x)
u_y_lower_boundary = np.zeros(Num_stencil_y)
u_x_upper_boundary = np.zeros(Num_stencil_x)
u_y_upper_boundary = np.zeros(Num_stencil_y)
ux_x_lower_boundary = np.zeros(Num_stencil_x)
ux_y_lower_boundary = np.zeros(Num_stencil_y)
ux_x_upper_boundary = np.zeros(Num_stencil_x)
ux_y_upper_boundary = np.zeros(Num_stencil_y)
uy_x_lower_boundary = np.zeros(Num_stencil_x) 
uy_y_lower_boundary = np.zeros(Num_stencil_y)
uy_x_upper_boundary = np.zeros(Num_stencil_x)
uy_y_upper_boundary = np.zeros(Num_stencil_y)
coner_xlow_ylow = np.zeros(1)
coner_xup_ylow = np.zeros(1)
coner_xlow_yup = np.zeros(1)
coner_xup_yup = np.zeros(1)
Time_step = 300
Delta_x = max(x_stencil) / (Num_stencil_x-1)
Delta_y = max(y_stencil) / (Num_stencil_y-1)
C = 1
Nu = 0.5
Delta_t = 0.2
CFL_x = C * Delta_t / Delta_x
CFL_y = C * Delta_t / Delta_y
CFL = min(CFL_x, CFL_y)
diffusion_number_x = Nu * Delta_t / Delta_x**2
diffusion_number_y = Nu * Delta_t / Delta_y**2
total_movement = C * Delta_t * (Time_step+1)
print("Δx:", Delta_x, "Δt:", Delta_t, "CFL:", CFL_x, "ν:", Nu, "d:", diffusion_number_x)
fig = plt.figure(2)
axe = fig.add_subplot(111, projection='3d')
surface = axe.plot_surface(x_array, y_array, u_array, cmap="bwr", label="Initail condition")
fig.colorbar(surface)
axe.set_xlabel("x")
axe.set_ylabel("y")
axe.set_zlabel("u")
axe.set_xlim(0, 100)
axe.set_ylim(0, 100)
axe.set_zlim(0, 1)
axe.set_title("2d initial condition")
plt.show()

def solve_advection_cip_2d(u_cip, ux_cip, uy_cip, 
                           x_lower_boundary, x_upper_boundary, 
                           y_lower_boundary, y_upper_boundary,
                           ux_x_lower_boundary, ux_x_upper_boundary,
                           ux_y_lower_boundary, ux_y_upper_boundary,
                           uy_x_lower_boundary, uy_x_upper_boundary,
                           uy_y_lower_boundary, uy_y_upper_boundary,
                           coner_ll, coner_lu, coner_ul, coner_uu,
                           dt=Delta_t, 
                           velocity_x=C, velocity_y=0.0):
    """
    Solve the advection equation by using CIP method.
    
       -  x  +
     - 5333337
       1*****2
     y 1*****2
       1*****2
       1*****2
     + 6444448
     
    *: u_cip or ux_cip or uy_cip
    1: x_lower_boundary or ux_x_lower_boundary or uy_x_lower_boundary
    2: x_upper_boundary or ux_x_upper_boundary or uy_x_upper_boundary
    3: y_lower_boundary or ux_y_lower_boundary or uy_y_lowerboundary
    4: y_upper_boundary or ux_y_upper_boundary or uy_y_upper_boundary
    5: coner_ll
    6: coner_lu
    7: coner_ul
    8: coner_uu
    """
    if type(velocity_x) is not np.ndarray:
        velocity_x = np.ones(u_cip.shape) * velocity_x
    if type(velocity_y) is not np.ndarray:
        velocity_y = np.ones(u_cip.shape) * velocity_y
    # Memory the present values
    u_old = u_cip.copy()
    ux_old = ux_cip.copy()
    uy_old = uy_cip.copy()

    # Prepare for calculation
    xi = - np.sign(velocity_x) * velocity_x * dt
    zeta = - np.sign(velocity_y) * velocity_y * dt
    dx = - np.sign(velocity_x) * Delta_x
    dy = - np.sign(velocity_y) * Delta_y
    
    # Set the neighbourhood values for reducing the for loop
    u_with_xlower_bc = np.concatenate([x_lower_boundary[:, np.newaxis], u_old[:, :-1]], axis=1)
    u_with_xupper_bc = np.concatenate([u_old[:, 1:], x_upper_boundary[:, np.newaxis]], axis=1)
    u_with_ylower_bc = np.concatenate([y_lower_boundary[np.newaxis, :], u_old[:-1, :]], axis=0)
    u_with_yupper_bc = np.concatenate([u_old[1:, :], y_upper_boundary[np.newaxis, :]], axis=0)
    temp_u_with_all_bc = np.concatenate([x_lower_boundary[:, np.newaxis], u_old, 
                                         x_upper_boundary[:, np.newaxis]], axis=1)
    u_with_all_bc = np.concatenate([np.concatenate([coner_ll, y_lower_boundary,coner_ul])[np.newaxis, :],
                                    temp_u_with_all_bc,
                                    np.concatenate([coner_lu, y_upper_boundary, coner_uu])[np.newaxis, :]], axis=0)
    
    u_next_x = np.where(velocity_x > 0.0, u_with_xlower_bc, u_with_xupper_bc)
    u_next_y = np.where(velocity_y > 0.0, u_with_ylower_bc, u_with_yupper_bc)
    u_next_xy = np.where(((velocity_x > 0.0) & (velocity_y > 0.0)), u_with_all_bc[:-2, :-2], 
                    np.where(((velocity_x < 0.0) & (velocity_y > 0.0)), u_with_all_bc[:-2, 2:],
                        np.where(((velocity_x > 0.0) & (velocity_y < 0.0)), u_with_all_bc[2:, :-2], u_with_all_bc[2:, 2:])))
    ux_next_x = np.where(velocity_x > 0, 
                        np.concatenate([ux_x_lower_boundary[:, np.newaxis], ux_old[:, :-1]], axis=1), 
                        np.concatenate([ux_old[:, 1:], ux_x_upper_boundary[:, np.newaxis]], axis=1))
    ux_next_y = np.where(velocity_y > 0, 
                        np.concatenate([ux_y_lower_boundary[np.newaxis, :], ux_old[:-1, :]], axis=0), 
                        np.concatenate([ux_old[1:, :], ux_y_upper_boundary[np.newaxis, :]], axis=0))
    uy_next_x = np.where(velocity_x > 0, 
                        np.concatenate([uy_x_lower_boundary[:, np.newaxis], uy_old[:, :-1]], axis=1), 
                        np.concatenate([uy_old[:, 1:], uy_x_upper_boundary[:, np.newaxis]], axis=1))
    uy_next_y = np.where(velocity_y > 0, 
                        np.concatenate([uy_y_lower_boundary[np.newaxis, :], uy_old[:-1, :]], axis=0), 
                        np.concatenate([uy_old[1:, :], uy_y_upper_boundary[np.newaxis, :]], axis=0))
    u_next_x = np.where(velocity_x == 0.0, u_old, u_next_x)
    u_next_y = np.where(velocity_y == 0.0, u_old, u_next_y)
    u_next_xy = np.where(((velocity_x == 0.0) & (velocity_y != 0.0)), u_next_y, 
                    np.where(((velocity_x != 0.0) & (velocity_y == 0.0)), u_next_x,
                        np.where(((velocity_x == 0.0) & (velocity_y == 0.0)), u_old, u_next_xy)))
    ux_next_x = np.where(velocity_x == 0, ux_old, ux_next_x)
    ux_next_y = np.where(velocity_y == 0, ux_old, ux_next_y)
    uy_next_x = np.where(velocity_x == 0, uy_old, uy_next_x)
    uy_next_y = np.where(velocity_y == 0, uy_old, uy_next_y)
    dx = np.where(velocity_x == 0, 1.0, dx)
    dy = np.where(velocity_y == 0, 1.0, dy)
    
    # Calculate the cip coefficient
    c_1 = - (u_next_xy - u_next_x - u_next_y + u_old) / dx / dy \
            + (ux_next_y - ux_old) / dy + (uy_next_x - uy_old) / dx
    c_3 = (uy_next_x - uy_old) / dx ** 2 - c_1 / dx
    c_2 = (ux_next_y - ux_old) / dy ** 2 - c_1 / dy
    a_1 = ux_old
    a_2 = 3 * (u_next_x - u_old) / dx**2 - (ux_next_x + 2 * ux_old) / dx
    a_3 = (ux_next_x - ux_old) / 3 / dx**2 - 2 * a_2 / 3 / dx
    b_1 = uy_old
    b_2 = 3 * (u_next_y - u_old) / dy**2 - (uy_next_y + 2 * uy_old) / dy
    b_3 = (uy_next_y - uy_old) / 3 / dy**2 - 2 * b_2 / 3 / dy
    d = u_old
    
    # Calculate the values
    u_cip = a_3 * xi**3 + a_2 * xi**2 + a_1 * xi + \
            b_3 * zeta**3 + b_2 * zeta**2 + b_1 * zeta + \
            c_3 * xi**2 * zeta + c_2 * xi * zeta**2 + c_1 * xi * zeta + d
    ux_cip = 3 * a_3 * xi**2   + 2 * a_2 * xi + a_1 + 2 * c_3 * xi * zeta + c_2 * zeta**2 + c_1 * zeta
    uy_cip = 3 * b_3 * zeta**2 + 2 * b_2 * zeta + b_1 + 2 * c_2 * xi * zeta + c_2 * xi**2 + c_1 * xi
    return u_cip, ux_cip, uy_cip

def update_boundary_condition(u_cip, x_lower_boundary, x_upper_boundary, 
                             y_lower_boundary, y_upper_boundary):
    x_lower_boundary = min(0.0, np.min(u_cip[:, 0]))
    x_upper_boundary = max(0.0, np.max(u_cip[:, -1]))
    y_lower_boundary = min(0.0, np.min(u_cip[0, :]))
    y_upper_boundary = max(0.0, np.max(u_cip[-1, :]))

def solve_diffusion_2d(u_2d, 
                       x_lower_boundary, x_upper_boundary, 
                       y_lower_boundary, y_upper_boundary,
                       d_number_x, d_number_y,
                       v_x=None, v_y=None,
                       v_x_lower_bc=None, v_x_upper_bc=None,
                       v_y_lower_bc=None, v_y_upper_bc=None,
                       update_partial=0, u_array_integral=None, 
                       prop_lambda=1/2):
    """
Solve the diffusion equation. It also contains elements of non-advection terms.
    
    Parameters
    ----------
    u_2d : array-like
nx × ny row matrix
    upper_boundary : numpy.float64
        u_array[n]Next to
    lower_boundary : numpy.float64
        u_array[0]Next to
    d_number_* : numpy.float64
Diffusion number. Positive number.
    update_partial : Bool, default False
Set to True when updating the differential expression.
    u_array_integral : array-like, default None
Used when updating the differential formula.[lower_boundary, u_array, upper_boudary]N in the form of+Vector of 2.
    prop_lambda : numpy.float64, default 1/2
        0: Explicit method. Centered sapce difference.
        1/2:Semi-implicit method. Crank-Nicolson scheme
        1: Fully implicit method.Time retreat difference.

    Returns
    -------
    u_2d : array-like
nx × ny row matrix
    """
    #Get matrix size
    nx, ny = u_2d.shape
    matrix_size = nx * ny
    # Ax=Create A and b for b
    s_param = 1 + 2 * prop_lambda * (d_number_x + d_number_y)
    t_param = 1 - 2 * (1 - prop_lambda) * (d_number_x + d_number_y)
    east_matrix = np.eye(matrix_size, k=1)
    east_matrix[np.arange(nx-1, matrix_size-nx, nx), np.arange(nx, matrix_size-nx+1, nx)] *= 0
    west_matrix = np.eye(matrix_size, k=-1)
    west_matrix[np.arange(nx, matrix_size-nx+1, nx), np.arange(nx-1, matrix_size-nx, nx)] *= 0
    a_matrix = np.identity(matrix_size) * s_param \
                - prop_lambda * d_number_x * east_matrix \
                - prop_lambda * d_number_x * west_matrix \
                - prop_lambda * d_number_y * np.eye(matrix_size, k=nx) \
                - prop_lambda * d_number_y * np.eye(matrix_size, k=-nx)
    b_array = t_param * u_2d.flatten()
    temp_csr = d_number_x * (1 - prop_lambda) * (csr_matrix(east_matrix) + csr_matrix(west_matrix)) \
               + d_number_y * (1 - prop_lambda) * (csr_matrix(np.eye(matrix_size, k=nx)) + csr_matrix(np.eye(matrix_size, k=-nx)))
    b_array += temp_csr.dot(b_array)
    #Set boundary conditions/ Setting of boundary condition
    b_array[:nx] += ((1 - prop_lambda) * y_lower_boundary + prop_lambda * y_lower_boundary) * d_number_y
    b_array[matrix_size-nx:] += ((1 - prop_lambda) * y_upper_boundary + prop_lambda * y_upper_boundary) * d_number_y
    b_array[np.arange(nx-1, matrix_size, nx)] += ((1 - prop_lambda) * x_upper_boundary + prop_lambda * x_upper_boundary) * d_number_x
    b_array[np.arange(0, matrix_size-nx+1, nx)] += ((1 - prop_lambda) * x_lower_boundary + prop_lambda * x_lower_boundary) * d_number_x
    if update_partial == 1:
        #Executed when updating the differential expression in the x direction. Not used when solving ordinary diffusion equations.
        if type(v_x) is not np.ndarray:
            v_x = np.ones(u_2d.shape) * v_x
        if type(v_y) is not np.ndarray:
            v_y = np.ones(u_2d.shape) * v_y
        v_x_with_bc = np.concatenate([v_x_lower_bc[:, np.newaxis], v_x, v_x_upper_bc[:, np.newaxis]], axis=1)
        v_y_with_bc = np.concatenate([v_y_lower_bc[:, np.newaxis], v_y, v_y_upper_bc[:, np.newaxis]], axis=1)
        u_with_xbc = np.concatenate([x_lower_boundary[:, np.newaxis], u_array_integral, x_upper_boundary[:, np.newaxis]], axis=1)
        u_with_ybc = np.concatenate([y_lower_boundary[np.newaxis, :], u_array_integral, y_upper_boundary[np.newaxis, :]], axis=0)
        temp_b_array = - ((v_x_with_bc[:, 2:] - v_x_with_bc[:, :-2]) / 2 / Delta_x) * \
                             ((u_with_xbc[:, 2:] - u_with_xbc[:, :-2]) / 2 / Delta_x)
        temp_b_array -= ((v_y_with_bc[:, 2:] - v_y_with_bc[:, :-2]) / 2 / Delta_x) * \
                               ((u_with_ybc[2:, :] - u_with_ybc[:-2, :]) / 2 / Delta_y)
        b_array += temp_b_array.flatten()
    elif update_partial == 2:
        #Executed when updating the differential expression in the y direction. Not used when solving ordinary diffusion equations.
        if type(v_x) is not np.ndarray:
            v_x = np.ones(u_2d.shape) * v_x
        if type(v_y) is not np.ndarray:
            v_y = np.ones(u_2d.shape) * v_y
        v_x_with_bc = np.concatenate([v_x_lower_bc[np.newaxis, :], v_x, v_x_upper_bc[np.newaxis, :]], axis=0)
        v_y_with_bc = np.concatenate([v_y_lower_bc[np.newaxis, :], v_y, v_y_upper_bc[np.newaxis, :]], axis=0)
        u_with_xbc = np.concatenate([x_lower_boundary[:, np.newaxis], u_array_integral, x_upper_boundary[:, np.newaxis]], axis=1)
        u_with_ybc = np.concatenate([y_lower_boundary[np.newaxis, :], u_array_integral, y_upper_boundary[np.newaxis, :]], axis=0)
        
        temp_b_array = - ((v_x_with_bc[2:, :] - v_x_with_bc[:-2, :]) / 2 / Delta_y) * \
                             ((u_with_xbc[:, 2:] - u_with_xbc[:, :-2]) / 2 / Delta_x) \
                       - ((v_y_with_bc[2:, :] - v_y_with_bc[:-2, :]) / 2 / Delta_y) * \
                               ((u_with_ybc[2:, :] - u_with_ybc[:-2, :]) / 2 / Delta_y)
        b_array += temp_b_array.flatten()
    a_csr = csr_matrix(a_matrix)
    u_2d = spla.isolve.bicgstab(a_csr, b_array)[0]
    return u_2d.reshape(nx, ny)

#2d Bugers equation
images = []
fig_5 = plt.figure(5, figsize=(5,5))
axe = fig_5.add_subplot(111, projection='3d')
# surface = axe.plot_surface(x_array, y_array, u_array, cmap="bwr")
# fig.colorbar(surface)
axe.set_xlabel("x")
axe.set_ylabel("y")
axe.set_zlabel("u")
axe.set_xlim(0, 100)
axe.set_ylim(0, 100)
axe.set_zlim(0, 1)
axe.set_title("2d burgers")

u_cip = u_array.copy()
partial_u_cip_x \
    = (np.concatenate([u_cip[:, 1:], u_x_upper_boundary[:, np.newaxis]], axis=1) \
       - np.concatenate([u_x_lower_boundary[:, np.newaxis], u_cip[:, :-1]], axis=1)) / 2 / Delta_x
partial_u_cip_y \
    = (np.concatenate([u_cip[1:, :], u_y_upper_boundary[np.newaxis, :]], axis=0) \
       - np.concatenate([u_y_lower_boundary[np.newaxis, :], u_cip[:-1, :]], axis=0)) / 2 / Delta_y
u_boundary = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary)

#Setting stability conditions
cfl_condition = 1
diffusion_number_restriction = 1/2
#Evolve time
for n in range(21):
    update_boundary_condition(u_cip, *u_boundary)
    delta_t = Delta_t
    cfl_max = np.max(np.abs(u_cip)) * delta_t / min(Delta_x, Delta_y)
    if cfl_max > cfl_condition:
        #Judgment of CFL conditions
        # cfl_If it is larger than the condition, reduce the time step width dt so that the CFL condition is satisfied.
        delta_t = CFL * min(Delta_x, Delta_y) / np.max(np.abs(u_cip))
    diffusion_max = Nu * delta_t / (Delta_x**2 + Delta_y**2)
    if diffusion_max > diffusion_number_restriction:
        #Judgment of diffusion number
        # diffusion_number_If it is larger than the restriction, reduce the time step width dt so that the diffusion number condition is satisfied.
        delta_t = diffusion_max * (Delta_x ** 2 + Delta_y**2) / Nu
        diffusion_number_x *= delta_t / Delta_t
        diffusion_number_y *= delta_t / Delta_t
    u_props = (u_cip, partial_u_cip_x, partial_u_cip_y)
    u_boundary = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary)
    ux_boundary = (ux_x_lower_boundary, ux_x_upper_boundary, ux_y_lower_boundary, ux_y_upper_boundary)
    uy_boundary = (uy_x_lower_boundary, uy_x_upper_boundary, uy_y_lower_boundary, uy_y_upper_boundary)
    u_coner = (coner_xlow_ylow, coner_xlow_yup, coner_xup_ylow, coner_xup_yup)
    diffusion_numbers = (diffusion_number_x, diffusion_number_y)
    u_cip_star, partial_u_cip_x_star, partial_u_cip_y_star \
        = solve_advection_cip_2d(*u_props, *u_boundary, *ux_boundary, *uy_boundary, *u_coner, dt=delta_t)
    velocities = (u_cip_star, 0.0)
    velocity_x_boundaries = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary)
    velocity_y_boundaries = (u_x_lower_boundary, u_x_upper_boundary, u_y_lower_boundary, u_y_upper_boundary) #For the time being
    u_cip = solve_diffusion_2d(u_cip_star, *u_boundary, *diffusion_numbers)
    partial_u_cip_x = solve_diffusion_2d(partial_u_cip_x_star, *ux_boundary, *diffusion_numbers,
                                         *velocities, *velocity_x_boundaries,
                                         update_partial=1, u_array_integral=u_cip_star)
    partial_u_cip_y = solve_diffusion_2d(partial_u_cip_y_star, *uy_boundary, *diffusion_numbers,
                                         *velocities, *velocity_y_boundaries,
                                         update_partial=2, u_array_integral=u_cip_star)
    
    if n % 1 == 0 and n > 0:
        img = axe.plot_wireframe(x_array, y_array, u_cip, cmap="bwr")
        images.append([img])
ani = animation.ArtistAnimation(fig_5, images, interval=100)
# ani.save('wf_anim_art.mp4',writer='ffmpeg',dpi=100)
ani.save('2d-burgers.gif', writer="pillow")
HTML(ani.to_html5_video())

Summary

Next time, I would like to deal with basic equations for fluids. Also, the calculation has become quite heavy, so I would like to speed it up if I feel like it.

References

Recommended Posts

[Python] Fluid simulation: From linear to non-linear
Changes from Python 3.0 to Python 3.5
Changes from Python 2 to Python 3.0
Post from Python to Slack
Cheating from PHP to Python
Anaconda updated from 4.2.0 to 4.3.0 (python3.5 updated to python3.6)
Switch from python2.7 to python3.6 (centos7)
Connect to sqlite from python
Call Matlab from Python to optimize
Create folders from '01' to '12' with python
Post from python to facebook timeline
[Lambda] [Python] Post to Twitter from Lambda!
Connect to utf8mb4 database from python
Python (from first time to execution)
Post images from Python to Tumblr
[Python] Fluid simulation: Incompressible Navier-Stokes equation
How to access wikipedia from python
Python to switch from another language
Did not change from Python 2 to 3
Update Python on Mac from 2 to 3
From Python to using MeCab (and CaboCha)
Introduction to Discrete Event Simulation Using Python # 1
[Python] Fluid simulation: Implement the advection equation
How to update Google Sheets from Python
Send a message from Python to Slack
Private Python handbook (updated from time to time)
I want to use jar from python
Convert from katakana to vowel kana [python]
Push notification from Python server to Android
Connecting from python to MySQL on CentOS 6.4
Porting and modifying doublet-solver from python2 to python3.
How to access RDS from Lambda (python)
[Python] Fluid simulation: Implement the diffusion equation
Python> Output numbers from 1 to 100, 501 to 600> For csv
Introduction to Discrete Event Simulation Using Python # 2
Convert from Markdown to HTML in Python
[Amazon Linux] Switching from Python 2 series to Python 3 series
API explanation to touch mastodon from python
Introduction to Vectors: Linear Algebra in Python <1>
Connect to coincheck's Websocket API from Python
Send a message from Slack to a Python server
Edit Excel from Python to create a PivotTable
Updated to Python 2.7.9
How to open a web browser from python
Introduction to OPTIMIZER ~ From Linear Regression to Adam to Eve
Study from Python Hour7: How to use classes
[Python] Convert from DICOM to PNG or CSV
Import Excel file from Python (register to DB)
I want to email from Gmail using Python.
[Python] I want to manage 7DaysToDie from Discord! 1/3
From file to graph drawing in Python. Elementary elementary
[Python] How to read data from CIFAR-10 and CIFAR-100
[python] Create table from pandas DataFrame to postgres
[Bash] Use here-documents to get python power from bash
How to generate a Python object from JSON
Sum from 1 to 10
sql from python
How to handle Linux commands well from Python
ambisonics simulation python
[Python] Flow from web scraping to data analysis
MeCab from Python