Commit e40b6b2f authored by Gidon Lucian Bauer's avatar Gidon Lucian Bauer 🏳
Browse files

example 1 funktioniert anscheinend

parent 33f63501
%% 3D Finite Volume method
function [t, U] = FVM_3D(u0, flux, f, xspan, tspan, N)
function [t, U] = FVM_3D(u0, flux, f, xspan, tspan, N, courant_number)
U = zeros(N, 1, 3); % Nx?x3
dx = (xspan(2) - xspan(1)) / N;
......@@ -10,15 +10,24 @@ function [t, U] = FVM_3D(u0, flux, f, xspan, tspan, N)
t(1) = tspan(1);
time = 1;
while t(end) < tspan(2)
dt = dx/max(abs(eigs_euler(U(1, time, :)))); % needs looking into
dt = calc_dt(dx, U(:,time,:), courant_number);
for space = 1:N
% assuming periodic boundary conditions
if space == 1
fluxl = flux(f, U(end, time, :), U(space, time, :));
fluxr = flux(f, U(space, time, :), U(space+1, time, :));
% assuming periodic boundary conditions
% fluxl = flux(f, U(end, time, :), U(space, time, :));
% fluxr = flux(f, U(space, time, :), U(space+1, time, :));
% not assuming periodic boundary conditions
fluxl = flux(f, U(space, time, :), U(space+1, time, :));
fluxr = flux(f, U(space+1, time, :), U(space+2, time, :));
elseif space == N
fluxl = flux(f, U(space-1, time, :), U(space, time, :));
fluxr = flux(f, U(space, time, :), U(1, time, :));
% assuming periodic boundary conditions
% fluxl = flux(f, U(space-1, time, :), U(space, time, :));
% fluxr = flux(f, U(space, time, :), U(1, time, :));
% not assuming periodic boundary conditions
fluxl = flux(f, U(space-2, time, :), U(space-1, time, :));
fluxr = flux(f, U(space-1, time, :), U(space, time, :));
else
fluxl = flux(f, U(space-1, time, :), U(space, time, :));
fluxr = flux(f, U(space, time, :), U(space+1, time, :));
......
......@@ -3,14 +3,13 @@ function [F] = HLL_flux(func, ul, ur)
eigs = zeros(6,1);
eigs(1:3, 1) = eigs_euler(ul);
eigs(4:6, 1) = eigs_euler(ur);
disp(eigs);
al = min(eigs);
ar = max(abs(eigs));
if (al >= 0)
if (0 <= al)
F = func(ul);
elseif ((al <= 0) && (ar >= 0))
elseif ((al <= 0) && (0 <= ar))
F = (ar * func(ul) - al * func(ur) + al * ar * (ur - ul)) / (ar - al);
else
F = func(ur);
......
%% Calculates the dt from dx and the eigenvalues of the jacobian
function [dt] = calc_dt(dx, U, courant_number)
[si, ~, ~] = size(U);
max_ev = 0;
for ii = 1:si
eigs = eigs_euler(U(ii, 1, :));
max_ev = max(max_ev, max(abs(eigs)));
end
dt = (courant_number * dx)/max_ev;
end
\ No newline at end of file
......@@ -14,10 +14,12 @@ f = @f_euler;
flux = @HLL_flux;
%% Setup
N = 100;
N = 50;
xspan = [0;1];
tspan = [0;0.1];
tspan = [0;0.25];
courant_number = 0.8;
%% Set initial conditions
% Note: this ul and ur contains [rho, v, p], while all other u further in
......@@ -46,7 +48,7 @@ end
u0 = init_cond(ul, ur, N);
%% Solving the PDE
[t, U] = FVM_3D(u0, flux, f, xspan, tspan, N);
[t, U] = FVM_3D(u0, flux, f, xspan, tspan, N, courant_number);
%% Plotting
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment