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

Programmieraufgabe 4 -> läuft noch nicht

parent f7b07415
%% multidimensional Finite Volume method
function [U] = FVM_3D(u0, flux, f, xspan, tspan, N)
temp = 4;
U = zeros(N, temp*N, 3); % Nx(2N)x3
dx = (xspan(2) - xspan(1)) / N;
dt = dx/temp; % this is currently arbitrarily choosen
U(:, 1, :) = u0;
for time = 1:temp*N-1
for space = 2:N-1
fluxl = flux(f, U(space-1, time, :), U(space, time, :));
fluxr = flux(f, U(space, time, :), U(space+1, time, :));
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
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);
if (0 <= al)
F = func(ul);
elseif ((al <= 0) && (0 <= ar))
F = (ar * func(ul) - al * func(ur) + al * ar * (ur - ul)) / (ar - al);
else
F = func(ur);
end
end
\ No newline at end of file
%% Non linear Euler equation
% u(1) ^= rho : density
% u(2) ^= v : velocity
% 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);
y(1) = u(1) * u(2);
y(2) = u(1) * u(2)^2 + p;
y(3) = u(2) * (u(3) + p);
end
\ No newline at end of file
%% Set 3D inititial conditions with discontinuity at x = 0.5
function [u0] = init_cond(ul, ur, N)
u0 = zeros(N,3);
gamma = 1.4;
tot_energy = @(rho, v, p, gamma) (p/(gamma-1) + 0.5 * rho * v^2);
rhol = ul(1);
vl = ul(2);
pl = ul(3);
rhor = ur(1);
vr = ur(2);
pr = ur(3);
El = tot_energy(rhol, vl, pl, gamma);
Er = tot_energy(rhor, vr, pr, gamma);
if (mod(N,2) == 0)
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), 3) = El;
u0((N/2)+1:end, 3) = Er;
else
u0(1:floor(N/2), 1) = rhol;
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), 3) = El;
u0(ceil(N/2), 3) = (Er + El)/2;
u0(ceil(N/2)+1:end, 3) = Er;
end
end
\ No newline at end of file
%% 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;
fh = func(u + h.*ej);
f = func(u);
jac(i,j) = (fh(i) - f(i))/h;
end
end
end
\ No newline at end of file
%% Mathe V - Progammieraufgabe 04
% Stefan Basermann, Adrien Bordes & Gidon Bauer
%% Clear the workspace
clear; clc; close all;
%% Choose the part you want to test
choice = input('Choose which part you want to test. (1-2) ');
%% Define PDE
f = @f_euler;
%% Define Flux
flux = @HLL_flux;
%% Setup
N = 5;
xspan = [0;1];
tspan = [0;0.25];
%% 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]
ul = zeros(3,1);
ur = zeros(3,1);
if (choice == 1)
ul(1) = 1;
ul(2) = 0;
ul(3) = 1;
ur(1) = 0.125;
ur(2) = 0;
ur(3) = 0.1;
else
ul(1) = 1;
ul(2) = -1;
ul(3) = 0.4;
ur(1) = 1;
ur(2) = 1;
ur(3) = 0.4;
end
u0 = init_cond(ul, ur, N);
U = FVM_3D(u0, flux, f, xspan, tspan, N);
\ No newline at end of file
%% Calculates the maximal eigenvalue of Jacobian at ul and ur
function [ar] = max_eig(func, ul, ur)
h = 1e-8; % 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
disp(ar);
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-8; % 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);
if (al < -1e5)
al = -1e5;
end
disp(al);
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