Blame view

output/advection.m 858 Bytes
d380ec6c   Jérémy Baudry   second commit
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function E = advection(E,c,dx,dt)
%ADVECTION  is a 1d advection code Lax Wendroff scheme with Superbee flux limiter
%
%   E = advection(E,c,dx,dt)
%
% INPUTS:
%   E  = vector of the thing to be advected;
%   c  = scalar speed;
%   dx = the spatial resolution;
%   dt = the temporal resolution;
% 
% OUTPUT:
%   E  = vector of the thing to be advected - after advection;

if ~isequal(size(E),size(c))
    c = c';
end

f  = flux(E,c,dx,dt);
E  = E-dt*diffl(f)/dx;

function f = flux(E,c,h,dt)
   % Lax-Wendroff with Superbee flux limiting;
   theta = diffl(E)./(diffr(E)+3e-14);
   phi   = limiter(theta);
   f     = c.*E+c/2.*(1-c*dt/h).*diffr(E).*phi;
end


function y = diffl(x)
y = [x(1);diff(x)];
end

function y = diffr(x)
y = [diff(x);-x(end)];
end

 % Superbee
function phi = limiter(r)
    phi = max(0,max(min(1,2*r),min(r,2)));              
end
end