%% Die (nichtviskose) Burgers-Gleichung
% Wir loesen hier eine Erhaltungsgleichung der Form
%
% $$ u_t + f(u)_x = 0 $$
%
% mit der speziellen Wahl
%
% $$ f(u) = \frac12 u^2 $$ (Nichtviskose Burgers-Gleichung)
%

%% Zusammenhang mit der linearen Advektionsgleichung
% Die Code-Struktur ist fast ident mit dem Beispiel zur linearen
% Advektionsgleichung $$ u_t +  v u_x = 0 $$.
% Diese Gleichung ergibt sich als Spezialfall fuer die Wahl $$f(u) = v u $$.
%
% Bei der linearen Advektionsgleichung ist die
% Courant-Friedrichs-Lewy-Zahl $$C = \frac{v \Delta t}{h} $$ (Ausbreitung
% pro Zeitschritt im Verhaeltnis zur Gitterlaenge) die wichtige Kennzahl fuer
% die Stabilitaet expliziter Verfahren. Bei nichtlinearen Gleichungen, die
% in ausdifferenzierter Schreibung $$ u_t +  f'(u) u_x = 0 $$ lauten,
% uebernimmt $$ f' $$ die Rolle von $$ v $$ als einer "lokalen"
% Ausbreitungsgeschwindigkeit.
%
%% Rechengebiet und raeumliches Gitter
clear all
clc

xanf = -1.5; xend = 1.5;        % Intervallgrenzen
n=2100;                           % n innere Gitterpunkte
x=linspace(xanf,xend,n+2);
h=(xend-xanf)/(n+1);            % Gitterweite

%% Anfangsbedingung und Randbedingungen
% Verschiedene Anfangsbedingungen sind moeglich...
init='welle';
switch(init)
    case('rampe')
        u0 = -x;
        u0(x>0) = 0;
        u0(1)=0; %Randbedingung
    case('stufe')
        u0 = 0*x;
        u0((x>-1) & (x<0)) = 0.5;
        u0(1)=0; %Randbedingung
    case('welle')
        u0 = 0*x;
        range=((x>-1) & (x<0));
        u0(range) = 1-cos(2*pi*x(range));
        u0(1)=0; %Randbedingung
end

u = u0; % initialisiere zeitabhaengige Loesung u
unew = u; % zum jeweils alten und neuen Zeitschritt

%% Zeitschritt- und Ausgabe-Steuerung
% Die lokale Courant-Zahl waere  C = f'(u)*dt/h. Wir
% verwenden hier C allerdings als Abkuerzung fuer dt/h. Damit laesst sich der Zeitschritt
% automatisch der Gitterweite h anpassen.
C = 0.5;
dt=h*C;

t=0; tend = 4; % Start- und End-Zeitpunkt
delta_print=tend/20; t_next=delta_print; % Steuerung der Ausgabe-Intervalle

figure(1);
plot(x,u); hold on
pause;


%% Zeitschritt-Schleife
% Verschiedene Formulierungen des Zeitschritts sind moeglich...
tstep='MacCormack';


while (t< tend-0.5*dt)
    
    switch tstep
        case('konservativ')
            % Formulierung in konservativer Form mit
            % einer numerischen Fluss-Funktion fnum:  Es sollte hier die Differenz
            % fnum(i-1/2) - fnum(i+1/2) stehen.
            % Programmtechnisch realisiert als i+1/2 -> i !
            for i=2:length(u)-1
                unew(i)= u(i) + C*( fnum(u,i-1,C) - fnum(u,i,C) );
            end
            unew(end) = unew(end-1); %numerische Randbedingung
            u =unew;
            
        case('twostepLW')
            % Two-Step Lax-Wendroff:
            % Zuerst (erste Schleife) wird unew als Hilfsgroesse berechnet.
            % Dann (zweite Schleife) wird f dort ausgewertet und das tatsaechliche neue u
            % mit einer konservative Formulierung berechnet.
            for i=1:length(u)-1
                unew(i)= 0.5*(u(i) + u(i+1)) -C/2*( f(u(i+1)) - f(u(i)) );
            end
            unew(1)=0;
            for i=2:length(u)
                u(i)= u(i) - C* ( f(unew(i)) - f(unew(i-1)) );
            end
        case('MacCormack')
            % Auch ein Zweischritt-Verfahren.
            % Zuerst (erste Schleife) wird unew als Hilfsgroesse berechnet.
            % Dann (zweite Schleife) wird f dort ausgewertet und das tatsaechliche neue u
            % mit einer konservative Formulierung berechnet.
            for i=1:length(u)-1
                unew(i)= u(i)  +C*( f(u(i)) - f(u(i+1)) );
            end
            unew(1)=0;
            for i=2:length(u)
                u(i)= (u(i)+unew(i))/2 + C/2*( f(unew(i-1)) - f(unew(i)) );
            end
    end
    
    
    %u=unew;
    t = t+dt;
    
    if t>t_next-0.5*dt
        plot(x,unew)
        pause(0.1)
        t_next=t_next+delta_print;
    end
    
end
hold off


