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

Neue Struktur -> läuft noch nicht

parent 0f789684
%% 3D Finite Volume method
function [U] = FVM_3D(u0, flux, f, xspan, tspan, N)
U = zeros(N, N, 3); % NxNx3
function [t, U] = FVM_3D(u0, flux, f, xspan, tspan, N)
U = zeros(N, 1, 3); % Nx?x3
dx = (xspan(2) - xspan(1)) / N;
dt = (tspan(2) - tspan(1)) / N;
U(:, 1, :) = u0;
for time = 1:N-1
t = [];
t(1) = tspan(1);
time = 1;
while t(end) < tspan(2)
dt = dx/max(abs(eigs_euler(U(1, time, :)))); % needs looking into
for space = 1:N
% assuming periodic boundary conditions
if space == 1
......@@ -19,10 +22,12 @@ function [U] = FVM_3D(u0, flux, f, xspan, tspan, N)
else
fluxl = flux(f, U(space-1, time, :), U(space, time, :));
fluxr = flux(f, U(space, time, :), U(space+1, time, :));
end
end
U(space, time+1, 1) = U(space, time, 1) + (dt/dx) * (fluxl(1) - fluxr(1));
U(space, time+1, 2) = U(space, time, 2) + (dt/dx) * (fluxl(2) - fluxr(2));
U(space, time+1, 3) = U(space, time, 3) + (dt/dx) * (fluxl(3) - fluxr(3));
end
t(time+1) = t(time) + dt;
time = time + 1;
end
end
\ No newline at end of file
%% HLL Flux function
function [F] = HLL_flux(func, ul, ur)
al = min_eig(func, ul, ur);
ar = max_eig(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)
F = func(ul);
......
%% Calculates the eigenvalues for the given euler function
% u(1) ^= rho : density
% u(2) ^= rho*v : momentum
% u(3) ^= E : total energy
% pressure p = (gamma - 1) * (E - 0.5 * rho * v^2)
% speed of sound c = sqrt((gamma * p) / rho)
function [eigs] = eigs_euler(u)
gamma = 1.4;
p = (gamma - 1) * (u(3) - 0.5 * (u(2)^2/u(1)));
c = sqrt((gamma * p) / u(1));
eigs = (u(2) / u(1)) * ones(3,1);
eigs(1) = eigs(1) - c;
eigs(3) = eigs(3) + c;
end
\ No newline at end of file
%% Non linear Euler equation
% u(1) ^= rho : density
% u(2) ^= v : velocity
% u(2) ^= rho*v : momentum
% u(3) ^= E : total energy
% pressure p = (gamma - 1) * (E - 0.5 * rho * v^2)
function [y] = f_euler(u)
gamma = 1.4;
p = (gamma - 1) * (u(3) - 0.5 * u(1) * u(2)^2);
p = (gamma - 1) * (u(3) - 0.5 * (u(2)^2/u(1)));
y(1) = u(1) * u(2);
y(2) = u(1) * u(2)^2 + p;
y(3) = u(2) * (u(3) + p);
y(1) = u(2);
y(2) = u(2)^2/u(1) + p;
y(3) = (u(2)/u(1)) * (u(3) + p);
end
\ No newline at end of file
%% Set 3D inititial conditions with discontinuity at x = 0.5
% Note: this ul and ur contains [rho, v, p], while all other u further in
% this programm will contain [rho, v, E]
% this programm will contain [rho, rho*v, E]
function [u0] = init_cond(ul, ur, N)
u0 = zeros(N,3);
......@@ -22,8 +22,8 @@ function [u0] = init_cond(ul, ur, N)
u0(1:(N/2), 1) = rhol;
u0((N/2)+1:end, 1) = rhor;
u0(1:(N/2), 2) = vl;
u0((N/2)+1:end, 2) = vr;
u0(1:(N/2), 2) = rhol * vl;
u0((N/2)+1:end, 2) = rhor * vr;
u0(1:(N/2), 3) = El;
u0((N/2)+1:end, 3) = Er;
......@@ -32,9 +32,9 @@ function [u0] = init_cond(ul, ur, N)
u0(ceil(N/2), 1) = (rhor + rhol)/2;
u0(ceil(N/2)+1:end, 1) = rhor;
u0(1:floor(N/2), 2) = vl;
u0(ceil(N/2), 2) = (vr + vl)/2;
u0(ceil(N/2)+1:end, 2) = vr;
u0(1:floor(N/2), 2) = rhol * vl;
u0(ceil(N/2), 2) = (rhor * vr + rhol * vl)/2;
u0(ceil(N/2)+1:end, 2) = rhor * vr;
u0(1:floor(N/2), 3) = El;
u0(ceil(N/2), 3) = (Er + El)/2;
......
%% Evaluates the jacobian of a function at point u; assumes square
function [jac] = jacobi(func, u, h)
si = length(u);
jac = zeros(si);
for i = 1:si
for j = 1:si
ej = zeros(si,1);
ej(j) = 1;
% First derivative using central finite differences: O(h^2)
f_plus_h = func(u + h.*ej);
f_minus_h = func(u - h.*ej);
jac(i,j) = (f_plus_h(i) - f_minus_h(i)) / (2*h);
end
end
% disp(jac);
end
\ No newline at end of file
......@@ -14,15 +14,14 @@ f = @f_euler;
flux = @HLL_flux;
%% Setup
N = 10;
N = 100;
xspan = [0;1];
tspan = [0;0.15];
% tspan = [0;0.25];
tspan = [0;0.1];
%% Set initial conditions
% Note: this ul and ur contains [rho, v, p], while all other u further in
% this programm will contain [rho, v, E]
% this programm will contain [rho, rho*v, E]
ul = zeros(3,1);
ur = zeros(3,1);
......@@ -46,12 +45,13 @@ end
u0 = init_cond(ul, ur, N);
U = FVM_3D(u0, flux, f, xspan, tspan, N);
%% Solving the PDE
[t, U] = FVM_3D(u0, flux, f, xspan, tspan, N);
%% Plotting
x = linspace(xspan(1), xspan(2), N);
t = linspace(tspan(1), tspan(2), N);
% t = linspace(tspan(1), tspan(2), N);
figure(1);
tiledlayout(2,2);
......@@ -66,7 +66,7 @@ ylabel('x');
nexttile;
pcolor(t, x, U(:,:,2));
% surf(t, x, U(:,:,2));
title('Velocity');
title('Momentum');
xlabel('t');
ylabel('x');
......
%% Calculates the maximal eigenvalue of Jacobian at ul and ur
function [ar] = max_eig(func, ul, ur)
h = 1e-7; % for numerical approximation of the jacobian
si = length(ul);
Jacl = jacobi(func, ul, h);
Jacr = jacobi(func, ur, h);
eigs = zeros(2*si, 1);
eigs(1:si) = eig(Jacl);
eigs(si+1:end) = eig(Jacr);
ar = max(abs(eigs));
if ar > 1e5
ar = 1e5;
end
end
\ No newline at end of file
%% Calculates the minimal eigenvalue of Jacobian at ul and ur
function [al] = min_eig(func, ul, ur)
h = 1e-7; % for numerical approximation of the jacobian
si = length(ul);
Jacl = jacobi(func, ul, h);
Jacr = jacobi(func, ur, h);
eigs = zeros(2*si, 1);
eigs(1:si) = eig(Jacl);
eigs(si+1:end) = eig(Jacr);
% al = min(eigs);
al = min(real(eigs)); % only using the real part of the eigenvalue
if (al < -1e5)
al = -1e5;
end
end
\ No newline at end of file
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