Азбучный вопрос по методу крупных частиц

[Roy Rogers]

Здравствуйте.

Я пытаюсь запрограммировать одномерный метод крупных частиц первого порядка без модификаций в точности по книжке Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике.
Решил протестировать на задаче Сода (то есть задача о распаде разрыва с начальными условиями слева $p = 1$, $\rho = 1$, $u = 0$ и справа $p = 0.1$, $\rho = 0.125$, $u=0$). Получается не то, чтобы очень хорошо, но терпимо



Хотя провал давления в центре меня несколько смущает. Явно от аналитического решения отличается. Можно ли это как-то улучшить?

И более важно то, что когда ударная волна доходит до правой границы расчетной области, схему разносит. Я там ставлю условие отражения $v(n, Nx+1) = - v(n, Nx)$.

Как быть? Буду очень благодарен за любую помощь.

Прилагаю свой код на фортране.

Используется синтаксис Fortran

program KRUCHA_1D_Sode

    implicit none

    ! Variables
    integer, parameter:: Nt=2500, Nx=300
    real Lx, dx, dt, gamma, Ku, TIME, pr, pl, ur, ul, flowR, flowL, dMR, dML, Xij, vc, vmax
    real u(0:Nt,0:Nx+1), rho(0:Nt,0:Nx+1), p(0:Nt,0:Nx+1), E(0:Nt,0:Nx+1)
    real ui(0:Nx+1), rhoi(0:Nx+1), pi(0:Nx+1), Ei(0:Nx+1)
    integer n, i, DL, DR

    ! Body of KRUCHA_1D_Sode
    Lx=2.0 
    gamma=1.4
    
    ! Initial conditions
    do i=0, Nx+1
        if(i.le.(Nx/2)) then
            u(0,i)=0.0
            rho(0,i)=1.0
            p(0,i) = 1.0
            E(0,i) = p(0,i)/(rho(0,i)*(gamma-1))
            else
            u(0,i)=0.0
            rho(0,i)=0.125
            p(0,i) = 0.1
            E(0,i) = p(0,i)/(rho(0,i)*(gamma-1))
         endif
    end do
    
    ! Boundary conditions at initial step
    u(0,0) = - u(0,1); u(0,Nx+1) = -u(0,Nx)
    
    ! Time & scheme parameters
    TIME = 0.0
    Ku = 0.1
    dt = Ku*dx/(sqrt(gamma*p(0,1)/rho(0,1)))
    dx=Lx/Nx
    
    ! Main cycle
    do n=0, Nt-1
        TIME = TIME + dt
        vmax = 0.0
        ! I. Euler step
        do i=1,Nx
            pr = 0.5*(p(n,i+1)+p(n,i)); pl = 0.5*(p(n,i-1)+p(n,i))
            ur = 0.5*(u(n,i+1)+u(n,i)); ul = 0.5*(u(n,i-1)+u(n,i))
            ui(i) = u(n,i) - (pr - pl)*dt/(dx*rho(n,i))
            Ei(i) = E(n,i) - (pr*ur - pl*ul)*dt/(dx*rho(n,i))
        end do
        
        ! Boundary conditions at Euler step
        ui(0) = -ui(1); Ei(0) = Ei(1)
        ui(Nx+1) = -ui(Nx); Ei(Nx+1) = Ei(Nx)
        
        ! Lagrange & final step cycle
        do i=1,Nx
            ! II. Lagrange step
            flowR = 0.5*(ui(i) + ui(i+1))
            if (flowR.ge.0) then
                dMR = rho(n,i)*flowR*dt; DR = 0
                else
                dMR = rho(n,i+1)*flowR*dt; DR = 1
            end if
            
            flowL = 0.5*(ui(i) + ui(i-1))
            if (flowL.ge.0) then
                dML = rho(n,i-1)*flowL*dt; DL = 1
                else
                dML = rho(n,i)*flowL*dt; DL = 0
            end if
            ! III. Final step
            rho(n + 1,i) = (rho(n,i)*dx + (dML - dMR))/dx
            Xij = rho(n+1,i)*dx - (1-DR)*dMR - (1-DL)*dML
            u(n+1,i) = (DL*dML*ui(i-1) - DR*dMR*ui(i+1) + ui(i)*Xij)/(rho(n+1,i)*dx)
            E(n+1,i) = (DL*dML*Ei(i-1) - DR*dMR*Ei(i+1) + Ei(i)*Xij)/(rho(n+1,i)*dx)
            p(n+1,i) = (gamma - 1)*rho(n+1,i)*(E(n+1,i) - 0.5*u(n+1,i)*u(n+1,i))
            
            ! Maximum velocity
            vc = Abs(u(n+1,i))+Sqrt(gamma*p(n+1,i)/rho(n+1,i))
            if (vc.gt.vmax) then ; vmax = vc; endif
        end do
        
        ! Boundary conditions at final step
        rho(n+1,0) = rho(n+1,1); u(n+1,0) = -u(n+1,1); p(n+1,0) = p(n+1,1); E(n+1,0) = E(n+1,1)
        rho(n+1,Nx+1) = rho(n+1,Nx); u(n+1,Nx+1) = -u(n+1,Nx); p(n+1,Nx+1) = p(n+1,Nx); E(n+1,Nx+1) = E(n+1,Nx);
        
        ! Time step
        dt = Ku*dx/vmax
    end do
    
    end program

[Утундрий]

Roy Rogers в сообщении #687151 писал(а):

Можно ли это как-то улучшить?

Можно, конечно. У вас грубо говоря три дороги: TVD, ENO, SPH.

P.S. А сам доисторический метод Белоцерковского - незамедлительно в печку, разумеется.

[Zai]

Цитата:

метод Белоцерковского - незамедлительно в печку...

Утундрий, то что в методе Харлоу-Белоцерковского использовано, обгорело в сотнях возвращаемых тел на Махе 20 без видимых расхождений с экспериментом не сопряжено с Вашими суждениями.

[Утундрий]

Боевая колесница ездить, конечно, могёт, но приходит время и вместе с ним танки.

[Zai]

Ваши "TVD, ENO, SPH" - убогая реинкарнация по времени, без повышения порядка аппроксимации на разрывах по пространственным координатам( не говоря о излучении и ряду других важных физических процессов), с заметным увеличением вычислительной вязкости и дикой недокументированной дисперсии. К тому же они сейчас уже бордаты.

[Утундрий]

Zai
Думайте как хотите.

[kor_viktor]

я тоже осваиваю метод крупных частиц
Скажите, а в каком компиляторе fortran вы работаете?, потому что у меня при вводе вашего кода кучу ошибок выкидывает

[Утундрий]

kor_viktor
У него стандарт 95 или выше.

[margullina]

Здравствуйте, на основании чего Вами был выбран шаг?

 dt = Ku*dx/(sqrt(gamma*p(0,1)/rho(0,1)))

и что Вы здесь вычисляете:

! Maximum velocity
            vc = Abs(u(n+1,i))+Sqrt(gamma*p(n+1,i)/rho(n+1,i))
            if (vc.gt.vmax) then ; vmax = vc; endif