Finite Difference Method - Upwind forward Euler

Krigo大约 4 分钟

1D linear advection equation

\begin{equation} \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0, x_{L} \leq x \leq x_{R}, t\geq 0 \tag{1} \end{equation}

where the independent variables are $t$ (time) and $x$ (space)

$x$ is restricted to the finite interval $[p,q]$ which is called the computational domain.
$a$ is a constant and the denpend variable
where $u$ = $u(x,t)$

Initial Condition

u(0,x) = \sin(x), p \leq x \leq q \tag{2}

Step 1: Spatial Discretisation

First we replace the computational domain by a finite set. Because the computaional domain contains infinte number of $x$ .

This process we called spatial discretisation.

The computational domain is repalced by a grid of $N$ equally spaced grid points.

The first grid point is $x = p$
The last grid point is $x = q$
The constant grid spacing, $\Delta x$

\Delta x = \frac{(q-p)}{N-1}

The values of $x$ in the discretised computational domain are indexed by the subscripts to give,

x_1 = p, \quad, x_2 = p+\Delta x, \quad \cdots \quad, x_N = p + (N-1)\Delta x = q

Since the grid spacing is a constant,

x_{i+1} = x_i + \Delta x

Fixing $t$ at $t = t_n$ , we approximate the, $u_x$ , in (1) for each point $(t_n, x_i)$ using the forward difference formula.

u_{x}(t_n,x_i) \approx \frac{u_{i+1}^{n} - u^{n}}{\Delta x}

Then, replacing $u_x$ in (1) by its approximation,

u_{t} + a\frac{u_{i+1}^n - u_i^n}{\Delta x} = 0

which is said to be in semi-discrete from since only the spatial derivative has been discretised 这里是只有空间被离散了

注

The grid is also called the mesh and the operation of discretising the computational domain is called gridding or meshing

Step 2: Time Discretisation

Then, we fix $x$ at $x = x_i$ and we approximate the temporal partial derivative, $u_t$ , in (1) for each point $(t_n, x_i)$ using the first order forward difference formula to gives,

u_t \approx \frac{u_{i}^{n+1} - u_{i}^{n}}{\Delta t}

Then, substituting the formula and we can get

\frac{u_i^{n+1} - u_i^n}{\Delta t} + a\frac{u^{n}_{i+1} - u^{n}_{i}}{\Delta t} = 0

Which rearranges to gives,

u^{n+1}_{j} = u^{n}_{j} - \frac{a\Delta t}{\Delta x}(u^{n}_{j} - u^{n}_{j-1}) \tag{3}

Code Implementation

The main file is below

program ying_feng_1D

    use globals1d
    use module_data1d
    implicit none

    ! Variables
    call parameters

    ! Body of ying_feng_1D
    print *, 'Start Run'
    do kkkk = 1,2!7
        !call cpu_time(start)
        
        call setup
        call allocate_var
        
        call init
        
        time = 0.
        
        do while( time < time_final )
            call setdt(time_final)
            call get_stream_tn
        enddo
        
        call deallocate_var
        !call cpu_time(finish)
    enddo
    
    print *, 'End Run'
    pause
    
    contains
    include "setup.f90"
    include "init.f90"
    include "allocate_var.f90"
    include "parameters.f90"
    
    include "setdt.f90"
    include "get_upstream_tn.f90"
    include "order_dg.f90"
    
end program ying_feng_1D

Step 1: Initial Parameters

First, we need to initial some variables which I mentioned above.

It should be a subroutine

subroutine init
    implicit none
    integer :: i
    
    dx = (xright-xleft)/nx
    
    do i = 1 - nghost,nx + nghost
        xgrid(i) = xleft + i * dx
    enddo
    
    do i = 1 - nghost,nx + nghost
        vertex(i)%coor = fun_init(xgrid(i))
    enddo
end subroutine init

dx : $\Delta x$
xgrid(i) : grid points
vertex(i)%coor : the value of initial value at $t_n$

The, we need to set up some value for our initial variables like this,

subroutine setup
    implicit none
    
    nx = 10 * 2 ** kkkk
    
    xleft = 0.
    xright = 2. * pi
    
    nghost = 1.
    time_final = 20.
    
    cfl = 0.1 !0.5
    
end subroutine setup
    
    
real function ax(x)
implicit none
real, intent(in) :: x

ax = 1.

end function ax
!***********************
real function exact(x,t)
implicit none
real, intent(in) :: x,t

exact = sin(x-t)

end function exact
!***********************
real function fun_init(x)
implicit none
real, intent(in) :: x

fun_init = sin(x)

end function fun_init

Varibales

nx : is the number of grid/mesh
xleft, right: $q,p$ in (1)
nghost : 这个是因为第一个值的预测值需要周期性信息
time_final : Time distance
cfl: a number for some pythical conservation and make the time move on.

Functions

ax(x): This is $a$ in the (1)
exact(x) : the exact solution of the PDE
fun_init: the $u(x,t_0)$

Step 2: Data structures

I think for this method, we at least need to set two different data structures.

For this file, It should be a module

module module_data1d
    
    !data structure
    !***********************
    type, public :: point1d
        sequence
        real :: coor
    end type
    
    type(point1d),allocatable,target,public :: vertex(:)

    !**********************
    type, public :: element1d
        sequence
        !type(point1d) :: porigin, pend
        real :: coor
        !integer :: id
    end type
    
    type(element1d),allocatable,target,public :: element(:)
    !*********************

end module module_data1d

The first one is point1d which is a public data structure and we have a real variable, coor.

coor: We need to store the initial value.

The second one is element1d which is a public data structure and we have... 编不下去了，当初写着写代码发现这个数据格式没意义，但是留下来了

Step 3: Start to calculate

First, we need have a for-loop in the main file and It looks like this in my code,

do kkkk = 1,2!7
    !call cpu_time(start)
    
    call setup
    call allocate_var
    
    call init
    
    time = 0.
    
    do while( time < time_final )
        call setdt(time_final)
        call get_stream_tn
    enddo
    
    call deallocate_var
    !call cpu_time(finish)
enddo

There have two main file for making the predicted value, setdt(time_final) and get_stream_tn. 这两个名字感觉命名的不是很好

This setdt is used for time forward.

subroutine setdt(time_f)
    implicit none
    real,intent(in) :: time_f
    
    ! ste dt and get tn+1
    dt = cfl*dx
    if(time+dt>time_f) dt = time_f - time

    time = time + dt
    
    print *,"Time��", time
    pause
    
end subroutine setdt

The get_stream_tn is used to execute the equation (3)

subroutine get_stream_tn
    implicit none
    integer i

    element(0)%coor = vertex(0)%coor - dt * (vertex(0)%coor - vertex(nx)%coor)/dx

    do i = 1,nx + nghost

        element(i)%coor = vertex(i)%coor - dt * (vertex(i)%coor - vertex(i-1)%coor)/dx
    enddo
    
    do i = 1 - nghost, nx + nghost
        
        vertex(i)%coor = element(i)%coor
    enddo
end subroutine get_stream_tn

Step 4: Calculate the Error

subroutine order_dg
    implicit none
    
    integer :: kk0
    real :: error1,error2,error3
    real :: rr1,rr2,rr3
    real :: exe
    
    do kk0 = 1 - nghost, nx + nghost
        write(1,*) xgrid(kk0), element(kk0)%coor
    enddo
    close(1)
    
    open(101, file='error_dg.txt')
    error1=0.0
    error2=0.0
    error3=0.0
    do kk0 = 1 - nghost, nx + nghost
        
        exe = exact( xgrid(kk0) , time_final ) - vertex(kk0)%coor
        
        error1 = error1 + abs(exe)
        error2 = error2 + exe**2
        error3 = max(error3,abs(exe))
    enddo
    
    error1=error1/nx
    error2=sqrt(error2/nx)
    
    if(kkkk.eq.1) write(101,103), nx, error1, error2, error3
    
    write(*,*) 'error', error1, error2

    if(kkkk.gt.1) then
        rr1 = log(er1/error1)/log(2.)
        rr2 = log(er2/error2)/log(2.)
        rr3 = log(er3/error3)/log(2.)
        write(101,102) nx,error1,rr1,error2,rr2,error3,rr3
        write(*,*) nx, rr1, rr2, rr3
    endif
    
    er1 = error1
    er2 = error2
    er3 = error3
    
102 format(i6,1x,3('&',1x, es12.2e2,1x,'&' 1x,f8.2 ,1x),'\\',1x,'\hline')
103 format(i6,1x,3('&',1x,es12.2E2,1x,'&',1x),'\\',1x,'\hline')
123 format(4(1x,f16.6))
    
    open(2, file="exact.plt")
    do kk0 = 1 - nghost, nx + nghost
        exe = 0.
        exe = exact(xgrid(kk0), time_final)
        write(2,123) xgrid(kk0), exe
    enddo
    close(2)

end subroutine order_dg