2021-11-10

fvm_1d_swe

理论主要是参考：【教学向】从零开始的流体力学（一维浅水方程+有限体积法） - 知乎 (zhihu.com)

而且我话的不是动图，而且耗散项系数为0时，也太震荡了。

function fvm_test_1d
clear all ;
clc;
L = 1;
div = 100;%space number
dx = L/div;
x0 = 0;x1 = 1;
y00 = 0;y1 = 0.25;
t0 = 0.0;
Final_time = 5;

g = 1.0;

x = (1:div)*dx - 0.5*dx;
x = x';
%initial height
ini_h = u0(x);
%initial velocity vector
ini_m = zeros(length(x),1);
ini_U = zeros(length(x),2);
ini_U(:,1) = ini_h;
ini_U(:,2) = ini_m; 
[T,Y] = rk4(0,Final_time,ini_U,0.01,dx);
save Y Y
%initial height profile 
function result = u0(x)
    result = 0.1+0.1*exp(-64*(x-0.25).*(x-0.25));
end
function result = f(H,M,i)
    result(1) = M(i);
    result(2) = M(i)*M(i)/H(i) + (g/2)*H(i)*H(i);
end
function result = flux(H,M,L,R,div)
    if(L==0)
        result = f(H,M,R);%只用到第二个返回值
        result(1) = 0;%将第一个返回值重写 
        return;
    end
    if(R==div+1)
        result = f(H,M,L);
        result(1) = 0;
        return;
    end
    fl = f(H,M,L);
    fr = f(H,M,R);
    d = 0;%0.5;%0.2;%[0.5,0.2,0]
    result(1) = (fl(1) + fr(1))/2 + d*(H(L) - H(R))/2;
    result(2) = (fl(2) + fr(2))/2 + d*(M(L) - M(R))/2;   
end 
function curr_du = dUdt(t,U,dx)%传入的t参数用不到
    curr_dhdt = [];
    curr_dmdt = [];
    curr_h = U(:,1);
    curr_m = U(:,2);
    div = size(U,1);
    for i = 1:div
        flux_l = flux(curr_h,curr_m,i-1,i,div);
        flux_r = flux(curr_h,curr_m,i,i+1,div);
        curr_dhdt = [curr_dhdt;(flux_l(1) - flux_r(1))/dx];
        curr_dmdt = [curr_dmdt;(flux_l(2) - flux_r(2))/dx];        
    end
    curr_du = [curr_dhdt,curr_dmdt];
end
function [T,Y] = rk4(t0,t1,y0,h,dx)
    M = (t1-t0)/h;
    T=t0:h:t1; 
    Y = y0;
    for j=1:M
        k1=h*feval(@dUdt,T(j),Y,dx);
        k2=h*feval(@dUdt,T(j)+h/2,Y+k1/2,dx);
        k3=h*feval(@dUdt,T(j),Y,dx);
        k4=h*feval(@dUdt,T(j)+h,Y+k3,dx);
        Y=Y+(k1+2*k2+2*k3+k4)/6;
        Save_Indicator = 1;
        if (Save_Indicator)
            fprintf('Plot ... \n') ; plot(x, Y(:,1)); %grid on
            title(['t = ', num2str(T(j), '%1.2f')]);
            xlim([x0,x1])
            ylim([y00,y1])
            pause(0.001)
        end
    end
    %R = [T,Y];
%     plot(x,Y)
end
end